饮料网站模板,莞城网页设计,大型网站建设规范,wdcp wordpress搬家IDL学习笔记#xff08;四#xff09; MODIS Grid数据的重投影 正弦投影 是以 米 为单位的 经纬度网格 是以 度 为单位的 但是转换之后#xff0c;不会一一对应#xff0c;所以需要对中间空缺位置需要进行一个填补。
核心问题:
把一个点从一个空间参考系放到另一个空间参… IDL学习笔记四 MODIS Grid数据的重投影 正弦投影 是以 米 为单位的 经纬度网格 是以 度 为单位的 但是转换之后不会一一对应所以需要对中间空缺位置需要进行一个填补。
核心问题:
把一个点从一个空间参考系放到另一个空间参考系
对应的编程关键:
1.数组A的元素应该放到数组B的哪个位置?
2.数组B的行列数量?
3.数组B的空值如何填补
MODIS Grid数据的重投影
1.对应的专业问题:
数组A的元素应该放到数组B的哪个位置投影
2.坐标与经纬度坐标的转换
数组B的行列数量:经纬度格网的范围
3.数组B的空值如何填补:
移动窗口均值计算 首先确定经纬度最值。数组行列数:原图像包含的经纬度范围最值 重投影结果分辨率设置
重投影结果分辨率设置不能比原来更细因为这样的话空缺会更多将带来更大误差。
每个像元坐标从哪里来 知道第一列坐标那么第二列坐标 就是 第一列坐标 像元分辨率。所以其余坐标就可以利用像元分辨率进行 加减乘除 运算 得知。
像元分辨率没有提到需要自己计算。利用左上角、右下角坐标以及行数、列数进行计算得到。 再找到我们需要的“一行”信息再定位到取出里面的坐标信息。为什么不用strmid因为内的信息可能是不固定的长度不唯一。
根据前面的关键字通过strpos定位到起始位置和结束位置获取总共子串长度然后就可以用strsplit分割函数获取所需数值。
strsplit(ul_info, (,), /extract)1.(,)只要遇到引号里的任何一个内容就把他分开 2./extract 如果需要提取内容则输出的是分割的下标。所以需要提取信息的话需要添加关键字。这时候提取结果是一个 字符型 数组里面包含坐标信息。 sd_id hdf_sd_start(file_list[file_i],/read)gindex hdf_sd_attrfind(sd_id,StructMetadata.0)hdf_sd_attrinfo, sd_id, gindex, datametadata;左上角坐标获取ul_start_pos strpos(metadata,UpperLeftPointMtrs)ul_end_pos strpos(metadata,LowerRightMtrs)ul_info strmid(metadata,ul_start_pos, ul_end_pos-ul_start_pos)ul_info_spl strsplit(ul_info,(,),/extract) ; 加 /extract 才返回的是内容而不是切割点的索引ul_prj_x double(ul_info_spl[1])ul_prj_y double(ul_info_spl[2]);右下角坐标获取lr_start_pos strpos(metadata,LowerRightMtrs)lr_end_pos strpos(metadata,Projection)lr_info strmid(metadata,lr_start_pos,le_end_pos - le_start_pos)lr_split strsplit(lr_info,(,),/extract) ;这里一开始忘添加关键字了导致后面查报错查了很久lr_prj_x double(lr_split[1])lr_prj_y double(lr_split[2]);map_proj_inverse 用于将 投影坐标 转为 经纬度坐标sin_prj map_proj_init(sinusoidal, /gctp, sphere_radius 6371007.181,center_longitude0.0,false_easting0.0,false_northing0.0);投影坐标参数从数据里得知geo_loc map_proj_inverse(proj_x, proj_y, map_strcture sin_prj); 进行坐标转换需要知道原来的投影的坐标 lon_min min(geo_x)lon_max max(geo_x)lat_min min(geo_y)lat_max max(geo_y)geo_resolution 0.01 ;1公里的空间分辨率对应的大概是0.01度data_box_geo_col ceil((lon_max - lon_min)/geo_resolution) ; ceil 防止无法整除带来的误差小数部分需要向上取整保证数据全部存放没有丢失部分data_box_geo_line ceil((lat_max - lat_min)/geo_resolution)data_box_geo fltarr(data_box_geo_col, data_box_geo_line) ; 初始化data_box_geo[*,*] -9999.0data_box_geo_col_pos floor((geo_x - lon_min)/geo_resolution) ; 归到哪一列需要向下取整因为data_box_geo_line_pos floor((geo_y - lat_min)/geo_resolution)
对应位置存放
data_box_geo[data_box_geo_col_pos, data_box_geo_line_pos] data ; 对应行列号上应放置的像素是哪一个
原理参考如下 填充条件 1.同时遇到_FillValue这里是0和-9999初始值要忽略不参与运算。 2.有N个以上有效值例如3的时候才进行填充否则就是不填充.不填充结果也设置为0 3.不需要填充的数组保持不动
对空缺值进行填充平滑
data_box_geo_out fltarr(data_box_geo_col, data_box_geo_line)for data_box_geo_col_i1,data_box_geo_col-2 do beginfor data_box_geo_line_i1,data_box_geo_line-2 do beginif data_box_geo[data_box_geo_col_i,data_box_geo_line_i] eq -9999.0 then begintemp_windowdata_box_geo[data_box_geo_col_i-1:data_box_geo_col_i1,data_box_geo_line_i-1:data_box_geo_line_i1]temp_window(temp_window gt 0.0)*temp_windowtemp_window_sumtotal(temp_window)temp_window_numtotal(temp_window gt 0.0)if (temp_window_num gt 3) then begindata_box_geo_out[data_box_geo_col_i,data_box_geo_line_i]temp_window_sum/temp_window_numendifendif else begindata_box_geo_out[data_box_geo_col_i,data_box_geo_line_i]data_box_geo[data_box_geo_col_i,data_box_geo_line_i]endelseendforendfor对geo结构体进行定义 geo_info {$MODELPIXELSCALETAG:[geo_resolution,geo_resolution,0.0],$;x、y、z方向的像元分辨率MODELTIEPOINTTAG:[0.0,0.0,0.0,lon_min,lat_max,0.0],$;坐标转换信息前三个0.0代表栅格图像上的第0,0,0个像元位置z方向一般不存在后面-180.0代表x方向第0个位置对应的经度是-180.0度90.0代表y方向第0个位置对应的纬度是90.0度GTMODELTYPEGEOKEY:2,$GTRASTERTYPEGEOKEY:1,$GEOGRAPHICTYPEGEOKEY:4326,$GEOGCITATIONGEOKEY:GCS_WGS_1984,$GEOGANGULARUNITSGEOKEY:9102,$GEOGSEMIMAJORAXISGEOKEY:6378137.0,$GEOGINVFLATTENINGGEOKEY:298.25722}输出一个文件后打印输出的提示信息可以都设置为以下格式 print,Reprojection time consuming of file file_basename(file_list[file_i]):strcompress(string(end_time-start_time)) s.
完整代码
pro modis_sinusoidal_to_geographicinput_directory E:\Data\IDL\chapter_3\modis_grid\output_directory E:\Data\IDL\chapter_3\modis_grid\geo_out\directory_exist file_test(output_directory,/directory)if directory_exist eq 0 then beginfile_mkdir,output_directoryendiffile_list file_search(input_directory,*.hdf)file_n n_elements(file_list);print, file_nfor file_i 0, file_n - 1 do begin ;file_n - 1 而不是file_nstart_time systime(1)result_name output_directory file_basename(file_list[file_i],.hdf) _geo.tiffsd_id hdf_sd_start(file_list[file_i],/read)gindex hdf_sd_attrfind(sd_id,StructMetadata.0)hdf_sd_attrinfo, sd_id, gindex, datametadata;左上角坐标获取ul_start_pos strpos(metadata,UpperLeftPointMtrs)ul_end_pos strpos(metadata,LowerRightMtrs)ul_info strmid(metadata,ul_start_pos, ul_end_pos-ul_start_pos)ul_info_spl strsplit(ul_info,(,),/extract) ; 加 /extract 才返回的是内容而不是切割点的索引ul_prj_x double(ul_info_spl[1])ul_prj_y double(ul_info_spl[2]);右下角坐标获取lr_start_pos strpos(metadata,LowerRightMtrs)lr_end_pos strpos(metadata,Projection)lr_info strmid(metadata,lr_start_pos,lr_end_pos - lr_start_pos)lr_info_splstrsplit(lr_info,(,),/extract)lr_prj_x double(lr_split[1])lr_prj_y double(lr_split[2])sds_index hdf_sd_nametoindex(sd_id,LST_Day_1km)sds_id hdf_sd_select(sd_id,sds_index)hdf_sd_getdata, sds_id, dataindex hdf_sd_attrfind(sds_id, scale_factor) ; 为什么这次没有计算——fillvalue? hdf_sd_attrinfo,sds_id, index, COUNTcali_num,DATAcali_scaledata data * cali_scale[0]hdf_sd_endaccess,sds_idhdf_sd_end,sd_iddata_size size(data)sin_resolution (lr_prj_x - ul_prj_x)/(data_size[1]) ; 计算原始正弦投影的 像元分辨率 proj_x fltarr(data_size[1], data_size[2]) ; 初始化数组存放坐标proj_y fltarr(data_size[1], data_size[2])for col_i 0, data_size[1] -1 do beginproj_x[col_i,*] ul_prj_x (sin_resolution * col_i) ; 利用 左上坐标、像元分辨率计算每个点坐标endforfor line_i 0, data_size[2] - 1 do beginproj_y[*, line_i] ul_prj_y - (sin_resolution * line_i)endfor;map_proj_inverse 用于将 投影坐标 转为 经纬度坐标sin_prj map_proj_init(sinusoidal, /gctp, sphere_radius 6371007.181,center_longitude0.0,false_easting0.0,false_northing0.0);投影坐标参数从数据里得知geo_loc map_proj_inverse(proj_x, proj_y, map_structure sin_prj); 进行坐标转换需要知道原来的投影的坐标;geo是2列1440000行的数组第一列是精度第二列是纬度geo_x geo_loc[0, *] ; x方向 对应应放入 data_box_geo_col_posgeo_y geo_loc[1, *] ; y方向 对应应放入 data_box_geo_line_pos;到此为止得知了原图上经纬度坐标。lon_min min(geo_x)lon_max max(geo_x)lat_min min(geo_y)lat_max max(geo_y)geo_resolution 0.01 ;1公里的空间分辨率对应的大概是0.01度data_box_geo_col ceil((lon_max - lon_min)/geo_resolution) ; ceil 防止无法整除带来的误差小数部分需要向上取整保证数据全部存放没有丢失部分data_box_geo_line ceil((lat_max - lat_min)/geo_resolution)data_box_geo fltarr(data_box_geo_col, data_box_geo_line) ; 初始化data_box_geo[*,*] -9999.0 ; 为什么赋值-9999.0 是为了区分是原本数据缺失值 还是 本身就有的0值_FillValue,判断是否要做填充data_box_geo_col_pos floor((geo_x - lon_min)/geo_resolution) ; 归到哪一列需要向下取整因为data_box_geo_line_pos floor((geo_y - lat_min)/geo_resolution)data_box_geo[data_box_geo_col_pos, data_box_geo_line_pos] data ; 对应行列号上应放置的像素是哪一个;对空缺值做一个填补操作----平滑data_box_geo_out fltarr(data_box_geo_col, data_box_geo_line)for data_box_geo_col_i1,data_box_geo_col-2 do beginfor data_box_geo_line_i1,data_box_geo_line-2 do beginif data_box_geo[data_box_geo_col_i,data_box_geo_line_i] eq -9999.0 then begintemp_windowdata_box_geo[data_box_geo_col_i-1:data_box_geo_col_i1,data_box_geo_line_i-1:data_box_geo_line_i1]temp_window(temp_window gt 0.0)*temp_windowtemp_window_sumtotal(temp_window)temp_window_numtotal(temp_window gt 0.0)if (temp_window_num gt 3) then begindata_box_geo_out[data_box_geo_col_i,data_box_geo_line_i]temp_window_sum/temp_window_numendifendif else begindata_box_geo_out[data_box_geo_col_i,data_box_geo_line_i]data_box_geo[data_box_geo_col_i,data_box_geo_line_i]endelseendforendforgeo_info {$MODELPIXELSCALETAG:[geo_resolution,geo_resolution,0.0],$;x、y、z方向的像元分辨率MODELTIEPOINTTAG:[0.0,0.0,0.0,lon_min,lat_max,0.0],$;坐标转换信息前三个0.0代表栅格图像上的第0,0,0个像元位置z方向一般不存在后面-180.0代表x方向第0个位置对应的经度是-180.0度90.0代表y方向第0个位置对应的纬度是90.0度GTMODELTYPEGEOKEY:2,$GTRASTERTYPEGEOKEY:1,$GEOGRAPHICTYPEGEOKEY:4326,$GEOGCITATIONGEOKEY:GCS_WGS_1984,$GEOGANGULARUNITSGEOKEY:9102,$GEOGSEMIMAJORAXISGEOKEY:6378137.0,$GEOGINVFLATTENINGGEOKEY:298.25722}write_tiff, result_name, data_box_geo_out,/float, geo_tiff geo_info end_time systime(1)print,Reprojection time consuming of file file_basename(file_list[file_i]):strcompress(string(end_time-start_time)) s.endfor
end
