卫星数据对应站点的匹配

对于卫星数据与对应站点的匹配过去我一直使用分别检索最近的经度纬度的办法找点然而在某现情况下出现缺失值于是询问老师采用更加稳健的3*3网格匹配的方法1.3*3网格法阈值解释请查看下面的HTML因为图片超过5MB上传不了查看的话可以先复制下面的链接写入记事本接着修改后缀为html就可以用浏览器打开了img srchttps://img.remit.ee/api/file/BQACAgUAAyEGAASHRsPbAAEUF_hqBoOeka9lEGBZlwABV2RH8cHmYpwAAo8eAALdWjBUwExnXQ9P-9s7BA.jpg altBQACAgUAAyEGAASHRsPbAAEUF_hqBoOeka9lEGBZlwABV2RH8cHmYpwAAo8eAALdWjBUwExnXQ9P-9s7BA.jpg /1° 经纬度 ≈ 111 km 4km 像素 ≈ 4/111 ≈ 0.04°%使用panoply打开nc文件查看里面变量名称 lon ncread(fullfile(AQUA_MODIS.20231101_20231108.L3m.8D.RRS.Rrs_412.9km.nc), lon); lat ncread(fullfile(AQUA_MODIS.20231101_20231108.L3m.8D.RRS.Rrs_412.9km.nc), lat); data ncread(fullfile(AQUA_MODIS.20231101_20231108.L3m.8D.RRS.Rrs_412.9km.nc), Rrs_412); %输入目标站点的经纬度lat_target与lon_target %由于我们直接输入的是nc文件自带的lon和lat他们的维度分别是4320*1和2160*1 %这时候我们需要将其转换成网格才能计算 [LON, LAT] meshgrid(lon, lat); dist (LAT-lat_target).^2(LON-lon_target).^2; %获取最短距离 min_dist min(dist(:)); %设置距离阈值 if min_dist0.001%这是针对于4km分辨率来讲 val NAN; else [row, col] find(dist min(min(dist))); row row(1); col col(1);%因为和最小距离相同的可能有多个行列所以我们选取第一个最小值的行列 %提取3*3窗口 %因为data为4320*2160row是纬度最多2160行是data的第二维 %col是经度最多4320是data的第一维 %经过下面的运算就可以获得row-1,row1,col-1,col_1并且不超过边界 r1 max(1, row-1); r2 min(size(data,2), row1); c1 max(1, col-1); c2 min(size(data,1), col1); %计算均值 data_patch data(c1:c2, r1:r2); valid_vals data_patch(~isnan(data_patch)); val ~isempty(valid_vals) ? mean(valid_vals) : NaN; %在运算时如果上面一行代码被提示运算符无效可以尝试使用 if ~isempty(valid) pp_hot mean(valid); else pp_hot NaN; end end这样就可以减少NAN的数量2.最近经纬度的办法%使用panoply打开nc文件查看里面变量名称 lon ncread(fullfile(AQUA_MODIS.20231101_20231108.L3m.8D.RRS.Rrs_412.9km.nc), lon); lat ncread(fullfile(AQUA_MODIS.20231101_20231108.L3m.8D.RRS.Rrs_412.9km.nc), lat); Rrs_412 ncread(fullfile(AQUA_MODIS.20231101_20231108.L3m.8D.RRS.Rrs_412.9km.nc), Rrs_412); %输入目标站点的经纬度lat_target与lon_target station_idxzeros(1,2)%经度索引和纬度索引 [~, ix] min(abs(lon - lon_target)); [~, iy] min(abs(lat - lat_target)); station_idx(s,:) [ix, iy];%保存经纬度索引 data Rrs_412(ix,iy);%得到目标站点Rrs_412结果