IDL实现遥感数据的快速重投影（几何校正）- 以MODIS Swath产品为例

网上关于遥感数据的重投影资料其实已经很多了，比如基于GDAL的实现、基于ENVI GLT的实现、基于一些专用处理工具的实现（MCTK、VCTK、HEG、MRT）等。但实际上，IDL中本身提供了用于数据几何校正的一些函数，如poly_2d、warp_tri。这里主要介绍一种基于poly_2d函数的几何校正，特点是快速（类似于gdalwarp），校正精度相比GLT要低一些，但个人认为对于遥感数据的可视化而言，没有什么差别，毕竟最终目标就只是看一个空间趋势。如果真的要严格定量，可能GLT方法更好。

文章目录

IDL实现遥感数据的快速重投影（几何校正）- 以MODIS Swath产品为例
一、关于poly_2d
二、Xo和Yo的求法
三、代码实现
四、存在的疑问

一、关于poly_2d

poly_2d实际上是基于多项式几何校正的一种函数实现，具体原理在很多讲遥感原理的书里都有，网上也有很多资料，这里就不再提了。poly_2d在IDL中的语法如下：

Result = POLY_2D( Array, P, Q [, Interp [, Dimx, Dimy]] [, CUBIC={-1 to 0}] [, MISSING=value] [, PIXEL_CENTER=value])

其中：
Array是存储待校正数据的数组；
P和Q则是用于几何校正的多项式系数，一般是配合polywarp函数使用来获取；
interp是校正过程中使用的插值方法，0代表最近邻插值，1代表双线性插值，2代表三次卷积插值；
Dimx和Dimy分别代表输出图像的列数和行数；
MISSING代表在图像边界外的填充值，可自定义；
更多参数说明可直接查看IDL中关于poly_2d的详细说明。

上面提到了多项式系数需要用到polywarp函数：

POLYWARP, Xi, Yi, Xo, Yo, Degree, Kx, Ky [, /DOUBLE] [, STATUS=variable]

其中：
Xi和Yi，分别代表待校正图像中像元的列号和行号；
Xo和Yo，分别代表校正后图像像元的列号和行号；
Degree是拟合的自由度（理解不准确，存疑）；
Kx和Ky，分别代表要求的拟合多项式系数；
更多参数说明可直接查看IDL中关于polywarp的详细说明。

从上面可看出，关键的点在于Xo和Yo怎么求。

二、Xo和Yo的求法

这里以MODIS的气溶胶产品为例，假设已经通过HDF读取函数读出了其中的经度数据（modis_lon_data）、纬度数据（modis_lat_data）、AOD数据（modis_aod_data），最终目标想实现经纬度格网的重投影。
所谓的重投影，或几何校正，就是一个像元应该放到对应空间参考格网下哪行哪列的问题，polywarp这个函数实际上就是要用户自己计算出一些点在目标空间参考格网下的行列号，并一一与待校正影像中的原始行列号进行对应，求解出Kx和Ky，也就是原始行列号和目标空间参考格网下行列号的拟合系数，得到一个转换关系式，给poly_2d函数使用。我们以最简单的经纬度格网这种空间参考为例，来简单讲一讲目标空间参考格网下的行列号的计算，实际上是一个很简单的公式：

col_out=floor((modis_lon_data-min(modis_lon_data))/out_res)
line_out=floor((max(modis_lat_data)-modis_lat_data)/out_res)

其中floor代表向下取整，max是求数组中最大值的函数，min是求数组中最小值的函数，out_res是自己设定的空间分辨率（°为单位，我一般会根据数据的实际分辨率大概设置）。这个公式不难理解，假设说最小经度是110°，我现在要计算的像元经度是110.4°，输出分辨率我设置为0.5°，代进去算，得到col_out=0，也就是说这个像元要放到校正后格网中的第0列（只要没超出格网，那都算成属于110.0°到110.5°的第0个格网），纬度的计算也是同理。
有了这个计算公式之后，结合读取到的modis经度和纬度数组，就可以算出每一个像元在目标空间参考格网下的行列号，最终配合poly_2d实现重投影（几何校正）。
其他的一些细节，比如重投影后目标数组大小的计算、geotiff结构体的写法等，本人在B站讲MODIS数据重投影的时候有提及，这里不再多讲。

三、代码实现

用这个代码重投影MODIS的3 km分辨率气溶胶产品，每个文件耗时1.5 s（i9 10980 XE @ 4.8 GHz），Golden Cove架构上会更快（个人笔记本i7 12700H上测试1.2 s），在速度上明显优于MCTK或GLT。

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
endpro modis_swath_warpinput_directory='D:/mod04/'output_directory='D:/mod04/warp_out/'if ~file_test(output_directory,/directory) then file_mkdir,output_directoryfile_list=file_search(input_directory+'*.hdf',count=file_n)for file_i=0,file_n-1 do beginstart_time=systime(1)print,file_list[file_i]modis_lon_data=hdf4_data_get(file_list[file_i],'Longitude')modis_lat_data=hdf4_data_get(file_list[file_i],'Latitude')modis_aod_data=hdf4_data_get(file_list[file_i],'Image_Optical_Depth_Land_And_Ocean')modis_aod_data=(modis_aod_data gt 0)*modis_aod_data*0.001out_res=0.03data_size=size(modis_aod_data)col_ori=intarr(data_size[1],data_size[2])line_ori=intarr(data_size[1],data_size[2])for col_i=0,data_size[1]-1 do col_ori[col_i,*]=col_ifor line_i=0,data_size[2]-1 do line_ori[*,line_i]=line_iout_col_num=ceil((max(modis_lon_data)-min(modis_lon_data))/out_res)out_line_num=ceil((max(modis_lat_data)-min(modis_lat_data))/out_res)col_out=floor((modis_lon_data-min(modis_lon_data))/out_res)line_out=floor((max(modis_lat_data)-modis_lat_data)/out_res)polywarp,col_ori,line_ori,col_out,line_out,5,coe_x,coe_yresult_image=poly_2d(modis_aod_data,coe_x,coe_y,0,out_col_num,out_line_num,missing=0.0)geo_info={$MODELPIXELSCALETAG:[out_res,out_res,0.0],$MODELTIEPOINTTAG:[0.0,0.0,0.0,min(modis_lon_data),max(modis_lat_data),0.0],$GTMODELTYPEGEOKEY:2,$GTRASTERTYPEGEOKEY:1,$GEOGRAPHICTYPEGEOKEY:4326,$GEOGCITATIONGEOKEY:'GCS_WGS_1984',$GEOGANGULARUNITSGEOKEY:9102}write_tiff,output_directory+file_basename(file_list[file_i],'.hdf')+'.tiff',result_image,/float,geotiff=geo_infoend_time=systime(1)print,'Time consumption: '+string(end_time-start_time,format='(F0.4," s.")')endfor
end

四、存在的疑问

polywarp函数中我的degree设置为5，具体取多少合适不知道，我只是对照着MCTK、GLT等方法得到的结果取的“目视最优值”
warp_tri函数也可实现，且校正效果与MCTK、GLT等更接近，但有一些图像边界上的异常值
理论上说，并不需要计算所有像元的位置，均匀分布选择一些样本点即可（创建值域在0.0-1.0的索引数组，和总的样本数相乘即可得出均匀分布的样本点下标）。没有开展尝试，但这对于处理一些行列数量较大的数据比较重要，因为polywarp这个函数内部使用了大量的double型变量，很吃内存。另外，样本点数量少，对于加快polywarp的运算也有很大的作用。