前言

关于MODIS正弦投影产品的重投影，其实在我的B站讲IDL的视频里有详细介绍，写这篇文档的目的主要是分享另一种不常见的方法。
IDL里其实是有自带的重投影函数的，函数名叫map_proj_image。只是用这个函数容易踩坑，从而导致做出来的结果不正确。我个人也是尝试了很久才意识到问题的关键。

map_proj_image函数使用关键

IDL的帮助里，map_proj_image函数的语法如下：
Result = MAP_PROJ_IMAGE(Image [, Range] [, /BILINEAR] [, DIMENSIONS=vector] [, IMAGE_STRUCTURE=structure] [, MAP_STRUCTURE=structure] [, MASK=variable] [, MAX_VALUE=value] [, MIN_VALUE=value] [, MISSING=value] [, UVRANGE=variable] [, XINDEX=variable] [, YINDEX=variable] )
其中比较关键的输入：
Image：待重投影影像数组
Range：影像数组的角点投影坐标，4元素数组，x方向最小坐标、y方向最小坐标、x方向最大坐标、y方向最大坐标
DIMENSIONS：重投影后结果数组的大小（与输出分辨率有关）
IMAGE_STRUCTURE：输入数据投影信息，由map_proj_init函数创建
MAP_STRUCTURE：输出数据投影信息，由map_proj_init函数创建
其实这个函数的用法，就是把输入数据Image的角点信息获取到，再通过map_proj_init创建输入输出投影参数，就能实现对Image的重投影。

比较坑的点在于：通过hdf函数读出来的modis数据，一般都是上北下南的方向，而MAP_PROJ_IMAGE这个函数，要求上南下北，也就是需要对读出来的数据做rotate上下翻转。

单个文件的重投影示例

以下代码，data_reprojection即基于map_proj_image实现的重投影函数，modis_grid_corner_get是读取modis grid产品中正弦角点信息的函数，hdf4_data_get是读取hdf4数据集的函数。modis_grid_to_geographic是主程序，功能是把modis 6公里地表温度数据重投影为经纬度格网，可根据自己的实际情况修改主程序。其余关键点在程序里已有相关说明，请自行研究。

pro data_reprojection,input_projection,input_data,input_corner,output_projection,output_resolution,reprojected_data,reprojected_geo_info,helpme=helpmeif (keyword_set(helpme)) then beginprint,'函数功能为实现输入数据的重投影（设计目标不局限为经纬度，但geotiff结构体部分还未完工，目前仅能完整实现重投影为经纬度格网）'print,'输入参数：'print,'input_projection：输入数据自身的投影参数，需由map_proj_init创建'print,'input_data：待重投影数组'print,'input_corner：输入数据的角点坐标数组，4个元素分别为最小x坐标、最小y坐标、最大x坐标、最大y坐标'print,'output_projection：重投影结果的目标投影参数，需由map_proj_init创建'print,'output_resolution：重投影结果的像元分辨率目标大小（最终实际分辨率会与设定值有少量偏差，因此不推荐将该函数用于数据重投影后的分块拼接，而是整体重投影。比如，直接按modis正弦的正弦分块将数据进行拼接，之后使用该函数将拼接数组转换为经纬度格网或其他投影，好处是不会产生缝隙，也不会因分辨率设置产生空缺值）'print,'输出参数:'print,'reprojected_data：投影转换后的结果数组'print,'reprojected_geo_info：投影转换后结果数组的geotiff结构体'print,'使用代码示例:'print,"input_directory='O:/coarse_data/chapter_3/modis_grid/'"+string(10b)+$"output_directory='O:/coarse_data/chapter_3/modis_grid/rpj_out/'"+string(10b)+$"if ~file_test(output_directory,/directory) then file_mkdir,output_directory"+string(10b)+$"file_list=file_search(input_directory+'MOD11B3*.hdf',count=file_n)"+string(10b)+$"in_proj=map_proj_init('Sinusoidal',/gctp,sphere_radius=6371007.181,center_longitude=0.0,false_easting=0.0,false_northing=0.0)"+string(10b)+$"out_proj=map_proj_init('Geographic',/gctp)"+string(10b)+$"out_res=[0.06,0.06]"+string(10b)+$"for file_i=0,file_n-1 do begin"+string(10b)+$"print,file_list[file_i]"+string(10b)+$"result_name=output_directory+file_basename(file_list[file_i],'.hdf')+'_rpj.tiff'"+string(10b)+$"in_cn=modis_grid_corner_get(file_list[file_i])"+string(10b)+$"in_data=float(hdf4_data_get(file_list[file_i],'LST_Day_6km'))"+string(10b)+$"in_data=(in_data gt 0.0)*in_data*0.02"+string(10b)+$"data_reprojection,in_proj,in_data,in_cn,out_proj,out_res,out_data,out_geo_info"+string(10b)+$"write_tiff,result_name,rotate(out_data,7),/float,geotiff=out_geo_info"+string(10b)+$"endfor"returnendifreprojected_data_test=map_proj_image(rotate(input_data,7),input_corner,image_structure=input_projection,map_structure=output_projection,uvrange=reprojected_range)output_dimension=[ceil((reprojected_range[2]-reprojected_range[0])/output_resolution[0]),ceil((reprojected_range[3]-reprojected_range[1])/output_resolution[1])]reprojected_data=map_proj_image(rotate(input_data,7),input_corner,dimensions=output_dimension,image_structure=input_projection,map_structure=output_projection,uvrange=reprojected_range)reprojected_data_size=size(reprojected_data)x_res=(reprojected_range[2]-reprojected_range[0])/reprojected_data_size[1]y_res=(reprojected_range[3]-reprojected_range[1])/reprojected_data_size[2];待完善部分reprojected_geo_info={$MODELPIXELSCALETAG:[x_res,y_res,0.0],$MODELTIEPOINTTAG:[0.0,0.0,0.0,reprojected_range[0],reprojected_range[3],0.0],$GTMODELTYPEGEOKEY:2,$GTRASTERTYPEGEOKEY:1,$GEOGRAPHICTYPEGEOKEY:4326,$GEOGANGULARUNITSGEOKEY:9102}
end

function modis_grid_corner_get,modis_file,helpme=helpmeif keyword_set(helpme) then beginprint,'函数功能为返回MODIS grid数据中包含的正弦投影角点'print,'输入参数：'print,'modis_file：待生成经纬度坐标的modis正弦投影文件名'print,'返回值：'print,'文件包含的正弦投影角点数组，第0到第3个元素分别为x方向最小坐标、y方向最小坐标、x方向最大坐标、y方向最大坐标'return,0endifout_loc=dblarr(4)sd_id=hdf_sd_start(modis_file,/read)gindex=hdf_sd_attrfind(sd_id,'StructMetadata.0')hdf_sd_attrinfo,sd_id,gindex,data=metadatahdf_sd_end,sd_idul_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)out_loc[0]=double(ul_info_spl[1])out_loc[3]=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_spl=strsplit(lr_info,'=(,)',/extract)out_loc[2]=double(lr_info_spl[1])out_loc[1]=double(lr_info_spl[2])return,out_loc
end

function hdf4_data_get,file_name,sds_namesd_id=hdf_sd_start(file_name,/read)sds_index=hdf_sd_nametoindex(sd_id,sds_name)sds_id=hdf_sd_select(sd_id,sds_index)hdf_sd_getdata,sds_id,datahdf_sd_endaccess,sds_idhdf_sd_end,sd_idreturn,data
end

pro modis_grid_to_geographicinput_directory='O:/coarse_data/chapter_3/modis_grid/'output_directory='O:/coarse_data/chapter_3/modis_grid/rpj_out/'if ~file_test(output_directory,/directory) then file_mkdir,output_directoryfile_list=file_search(input_directory+'MOD11B3*.hdf',count=file_n)in_proj=map_proj_init('Sinusoidal',/gctp,sphere_radius=6371007.181,center_longitude=0.0,false_easting=0.0,false_northing=0.0)out_proj=map_proj_init('Geographic',/gctp)out_res=[0.06,0.06];输出分辨率（与公里分辨率接近的°即可）for file_i=0,file_n-1 do beginprint,file_list[file_i]result_name=output_directory+file_basename(file_list[file_i],'.hdf')+'_rpj.tiff'in_cn=modis_grid_corner_get(file_list[file_i])in_data=float(hdf4_data_get(file_list[file_i],'LST_Day_6km'));6 km地表温度数据集名in_data=(in_data gt 0.0)*in_data*0.02;数据预处理（乘scale factor）data_reprojection,in_proj,in_data,in_cn,out_proj,out_res,out_data,out_geo_infowrite_tiff,result_name,rotate(out_data,7),/float,geotiff=out_geo_infoendfor
end

多个文件的重投影+拼接

以下代码解决的是多个正弦投影产品的重投影并拼接的需求。要点已在程序中注释，也用到了上一节中的三个自定义函数，请自行补全。注意中间用到了imsl函数实现年积日到实际日期的转换，确认自己的IDL版本是否支持这个函数（8.5肯定支持）。

pro modis_grid_to_geographic_mosaicinput_directory='O:/coarse_data/chapter_3/modis_grid/';输入路径output_directory='O:/coarse_data/chapter_3/modis_grid/rpj_out/';输出路径product_prefix='MOD11B3*';grid产品文件名前缀dataset_name='LST_Day_6km';grid产品文件中待处理数据集名称out_res=[0.06d,0.06d];输出分辨率（与公里分辨率接近的°即可）threshold_value=0.0;grid产品中的无效值scale_factor=0.02;grid产品数据需要乘以的缩放因子tile_size_x=200L;grid产品数据集本身的列数，常见的为200 1200 2400 4800tile_size_y=200L;grid产品数据集本身的行数，常见的为200 1200 2400 4800in_proj=map_proj_init('Sinusoidal',/gctp,sphere_radius=6371007.181d,center_longitude=0.0d,false_easting=0.0d,false_northing=0.0d)out_proj=map_proj_init('Geographic',/gctp)if ~file_test(output_directory,/directory) then file_mkdir,output_directoryfile_list_all=file_search(input_directory+product_prefix+'.hdf',count=file_n_all)if file_n_all eq 0 then returnfile_prefix=(strmid(file_basename(file_list_all),8,8)).uniq() ;注意文件中的Axxxxxxx的年积日是否是从第8个位置开始取8个字符for prefix_i=0,n_elements(file_prefix)-1 do beginprint,file_prefix[prefix_i]year=fix(strmid(file_prefix[prefix_i],1,4))doy=fix(strmid(file_prefix[prefix_i],5,3))date_julian=imsl_datetodays(31,12,year-1)+doyimsl_daystodate,date_julian,day,month,yearresult_name=output_directory+string(year,format='(I04)')+string(month,format='(I02)')+string(day,format='(I02)')+'_rpj_mosaic.tiff'in_cn_ex=[999999999999999d,999999999999999d,-999999999999999d,-999999999999999d]in_cn_max=in_cn_exin_cn_min=in_cn_exfile_list=file_search(input_directory+product_prefix+file_prefix[prefix_i]+'*.hdf',count=file_n)h_value=fix(strmid(file_basename(file_list),18,2))v_value=fix(strmid(file_basename(file_list),21,2))h_min=fix(min(h_value),type=3)h_max=fix(max(h_value),type=3)v_min=fix(min(v_value),type=3)v_max=fix(max(v_value),type=3)h_num=h_max-h_min+1Lv_num=v_max-v_min+1Lgrid_box=fltarr(h_num*tile_size_x,v_num*tile_size_y)for file_i=0,file_n-1 do beginprint,file_list[file_i]in_cn=modis_grid_corner_get(file_list[file_i])in_cn_max=in_cn>in_cn_maxin_cn_min=in_cn<in_cn_minh_value=fix(strmid(file_basename(file_list[file_i]),18,2))v_value=fix(strmid(file_basename(file_list[file_i]),21,2))data_temp=float(hdf4_data_get(file_list[file_i],dataset_name))data_temp=(data_temp gt threshold_value)*data_temp*scale_factorcol_start_pos=(h_value-h_min)*tile_size_xcol_end_pos=col_start_pos+tile_size_x-1line_start_pos=(v_value-v_min)*tile_size_yline_end_pos=line_start_pos+tile_size_y-1grid_box[col_start_pos:col_end_pos,line_start_pos:line_end_pos]=data_tempendforin_cn_ex[0:1]=in_cn_min[0:1]in_cn_ex[2:3]=in_cn_max[2:3]data_reprojection,in_proj,grid_box,in_cn_ex,out_proj,out_res,out_data,out_geo_infowrite_tiff,result_name,rotate(out_data,7),/float,geotiff=out_geo_infoendfor
end

后记

其实已经在data_reprojection函数中写了，最开始的目标并不只是转成经纬度格网。奈何没想出比较好的快速生成其他投影geotiff结构体的方法，所以搁了。主要是平时自己也用不到其他投影，正所谓没需求就没动力，祝各位好运（手动狗头）。

IDL实现MODIS Grid（正弦投影）产品的重投影及拼接处理相关推荐

ENVI_IDL:批量重投影ModisSwath产品(调用二次开发接口)+解析
目录 1. 课堂内容 1. 获取Modis Swath数据(这里只获取Lat.Lon.Aod(气溶胶厚度)三个数据集以及aod数据集的两个属性),并对aod数据进行简单的处理 2. 调用二次开发接口以 ...
地理坐标系、大地坐标系与地图投影与重投影详解
地理坐标系.大地坐标系与地图投影与重投影详解基本概念首先简单介绍一下地理坐标系.大地坐标系以及地图投影的概念: 地理坐标系:为球面坐标. 参考平面地是椭球面,坐标单位:经纬度: 投影坐标系:为平面 ...
地理坐标系、大地坐标系、地图投影与重投影
转自:https://www.cnblogs.com/lishanyang/p/6008256.html 基本概念首先简单介绍一下地理坐标系.大地坐标系以及地图投影的概念: 地理坐标系:为球面坐标. ...
栅格重投影（投影变换）
OpenLayers能够在不同的坐标系统中显示来自WMS.WMTS.静态图像和许多其他源的栅格数据.图像的地图重投影直接发生在web浏览器中.在任何Proj4js支持的坐标参考系统中都是可视的,并且以 ...
IDL实现遥感数据的快速重投影（几何校正）- 以MODIS Swath产品为例
IDL实现遥感数据的快速重投影(几何校正)- 以MODIS Swath产品为例网上关于遥感数据的重投影资料其实已经很多了,比如基于GDAL的实现.基于ENVI GLT的实现.基于一些专用处理工具的实 ...
ENVI_IDL:批量重投影Modis Swath产品并指定范围输出为Geotiff格式+解析
目录 1. 实验内容 2. 知识储备 3. ENVI实操对对应DL代码部分 4. 编程 5. 题外话 5.1 n_element()与一些函数自带的count参数返回的区别 5.2 发现一个难以发现 ...
MODIS产品MCD12Q1数据介绍、下载与拼接重投影格式转换处理
1.MODIS数据介绍 2.MCD12Q1数据介绍 MODIS三级数据土地覆盖类型产品(Land Cover data)是根据一年的Terra和Aqua观测所得的数据经过处理,描述土地覆盖的类型.该土 ...
ENVI_IDL:批量处理Modis Swath数据的重投影并输出为Geotiff格式+详细解析
目录 1. 课堂内容 2. 知识储备 3. 编程 4. 题外话 1. 课堂内容批量处理Modis Swath数据的重投影并输出为Geotiff格式总体思路 1. 先获取Modis Grid产品的数 ...
Python实现VIIRS气溶胶产品重投影-类GLT实现
Python实现VIIRS气溶胶产品重投影-类GLT实现前言代码说明用到的外部python包代码实现注意事项后记前言最近用IDL写了个VIIRS气溶胶产品重投影的代码,之后心想着用py ...

IDL实现MODIS Grid（正弦投影）产品的重投影及拼接处理

IDL实现MODIS Grid（正弦投影）产品的重投影及拼接处理

前言

map_proj_image函数使用关键

单个文件的重投影示例

多个文件的重投影+拼接

后记

IDL实现MODIS Grid（正弦投影）产品的重投影及拼接处理相关推荐

最新文章

热门文章