IDL学习——处理自带经纬度文件的遥感影像，以哨兵5P数据为例

自带经纬度文件处理流程
- 1.寻找经纬度文件
- 2.构建GLT文件
- 3.使用GLT文件进行几何校正
- 4.批处理过程中遇到的问题

最近一直在做大气方面相关遥感数据的影像处理，由于大气遥感数据的幅宽一般较大，原始数据一般会自带经纬度文件，手动一景一景处理效率极低，故编写IDL代码进行处理。那么这篇文章就带各位在面对不同大气遥感数据时候，如何快速对该数据的处理流程有个大致了解。

自带经纬度文件处理流程

其实在我这篇博客里面https://mp.csdn.net/mp_blog/creation/editor/107604176 已经对自带经纬度文件的遥感数据处理在ENVI里面的流程大致叙述了一遍。如果有还不知道的小伙伴，可以翻看一下。主要就是：

找到数据集中的经纬度文件。
将经纬度文件制作成GLT文件。
根据GLT对目标数据集进行校正。
以哨兵-5P数据为例逐步进行处理

1.寻找经纬度文件

首先你得知道需要处理数据的格式，一般有两种格式–NC文件和HDF文件，知道了数据格式，就知道如何在ENVI里面打开数据，如下图：

因为哨兵5P数据格式为NetCDF-4格式，所以使用 Open as ---- Generic Formats ---- HDF5/NetCDF4方式打开nc数据。接下来就是寻找经纬度文件，一般就找 latitude or longitude or lat or lon 这几个关键字就能找到，像S5P里面就是 latitude和longitude。

那么在IDL代码里面是如何打开经纬度数据呢？
先上代码：

  file = 'G:\CH4\S5P_OFFL_L2__CH4____20200529T043518_20200529T061648_13603_01_010302_20200530T210438.nc'latitude= e.openraster(file,DATASET_NAME='/PRODUCT/latitude')longitude= e.openraster(file,DATASET_NAME='/PRODUCT/longitude')

这里主要用到了 openraster这个函数里面的一个参数---- DATASET_NAME,先看一下官方解释：

这里可以去看一下NC格式文件和HDF文件时如何存储数据的。可以简单这么理解：把NC数据看成一个大文件夹，这个大文件夹下面可能会有若干级子文件夹，DATASET_NAME就是去寻找底层文件夹里面放置的数据，相当于一个相对路径。回过头看我代码：

  latitude= e.openraster(file,DATASET_NAME='/PRODUCT/latitude')

已经知道DATASET_NAME里面写的是一个路径，路径由两部分组成： Path 和 Name。这里纬度文件的Path 就是下图里面 Path后面那一串，Name 就是 Path上面有个Name后面那一串。经度数据和目标数据都是同理打开。我这里处理的是哨兵5P——甲烷数据，目标数据集就为methane_mixing_ratio

FILE = 'G:\CH4\S5P_OFFL_L2__CH4____20200529T043518_20200529T061648_13603_01_010302_20200530T210438.nc'
raster=e.openraster(file,DATASET_NAME='/PRODUCT/methane_mixing_ratio')   ;打开目标数据集

2.构建GLT文件

现在已经找到了经纬度文件，那么构建GLT文件就很简单啦，直接上代码：

  latitude= e.openraster(file,DATASET_NAME='/PRODUCT/latitude')longitude= e.openraster(file,DATASET_NAME='/PRODUCT/longitude')out_name=e.gettemporaryfilename() ;选择创建临时文件输出路径lat_id=ENVIRASTERTOFID(latitude)lon_id=ENVIRASTERTOFID(longitude)envi_file_query,lat_id,nb=nb_laty_pos=lindgen(nb_lat)envi_file_query,lon_id,nb=nb_lonx_pos=lindgen(nb_lon)             ;提取需要用到的参数input_prj=ENVI_PROJ_CREATE(/geographic,datum='WGS-84')  output_prj=ENVI_PROJ_CREATE(/geographic,datum='WGS-84')  ;设置输入输出的坐标系为 WGS-84坐标系ENVI_DOIT,'ENVI_GLT_DOIT',dims=dims,I_PROJ=input_prj,O_PROJ=output_prj,pixel_size=0.08113,$OUT_NAME=out_name,ROTATION=0,X_FID=lon_id,$X_POS=x_pos[0], Y_FID=lat_id, Y_POS=y_pos[0]

这里需要注意一点就是在ENVI_DOIT里面有个参数是pixel_size,这个参数是需要自己去指定的，那么怎么指定呢？我一般是手动在ENVI里面构建GLT文件时，在设置GLT参数时，会自己有个默认值。也可以自己算：公式如下：这个有点不准，还是建议手动
1°111.11KM=卫星分辨率(十进制度)卫星分辨率(KM)\frac{1°}{111.11KM}=\frac{卫星分辨率(十进制度)}{卫星分辨率(KM)} 111.11KM1°=卫星分辨率(KM)卫星分辨率(十进制度)

3.使用GLT文件进行几何校正

这里主要就是找到目标raster，从数据集中找目标raster和找经纬度文件是一样的方法，这里就不赘述了。校正部分代码如下：

      glt_raster = e.openraster(out_name)  ;打开GLT文件glt_fid=envirastertofid(glt_raster)     envi_file_query,glt_fid,dims=dimsraster=e.openraster(file[i],DATASET_NAME='/PRODUCT/methane_mixing_ratio')    ;打开目标数据集raster_id=envirastertofid(raster) envi_file_query,raster_id,ns=ns,nb=nbpos=indgen(nb)                    envi_doit,'ENVI_GEOREF_FROM_GLT_DOIT',FID=raster_id,GLT_DIMS=dims,GLT_FID=glt_fid,out_name=out_path,pos=pos

4.批处理过程中遇到的问题

一开始我在主pro里面（就是从头执行到底不插入任何函数）进行上述所有操作，但是频繁给我报错，如下图：

总结下来就是，在for循环里面直接去构建GLT文件这种方法在处理第二景的时候就会报错，有时候提示GLT文件不存在，有时候提示不能连接到指定文件。
然后我便开始找原因，我觉得，可能是第一景文件的GLT文件还存在，导致后续创建GLT文件时，IDL认为GLT文件已存在，所以在处理第二景的时候， ENVI_DOIT,‘ENVI_GLT_DOIT’, 这一步直接给你跳过了，故在后续处理第二景的时候，读取GLT文件的时候读取不到。发现了问题，那就解决问题。按照常理想，你既然以为这个GLT文件存在，那我就在你用完这个GLT文件之后,也就是第一景处理完之后，把这个GLT文件关闭应该就可以了，但是还是不能解决问题。。。附上会报错的代码

pro CH4_dealcompile_opt idl2e=envi(/headless)time = systime(1)path = 'G:\CH4\Oridata'file = file_search(path,'*.nc',count=n)for i=0,n-1 do beginname = file_basename(file[i],'.nc4')print,nameout_path='G:\CH4\result\'+name+'_result.dat'; 重复处理数据不进行处理if file_test(out_path) eq 1 then beginprint,'已处理'endif else begin    ;构建GLT文件 latitude= e.openraster(file[i],DATASET_NAME='/PRODUCT/latitude')longitude= e.openraster(file[i],DATASET_NAME='/PRODUCT/longitude')    out_name=e.gettemporaryfilename()lat_id=ENVIRASTERTOFID(latitude)lon_id=ENVIRASTERTOFID(longitude)envi_file_query,lat_id,nb=nb_laty_pos=lindgen(nb_lat)envi_file_query,lon_id,nb=nb_lonx_pos=lindgen(nb_lon)input_prj=ENVI_PROJ_CREATE(/geographic,datum='WGS-84')output_prj=ENVI_PROJ_CREATE(/geographic,datum='WGS-84')ENVI_DOIT,'ENVI_GLT_DOIT',dims=dims,I_PROJ=input_prj,O_PROJ=output_prj,pixel_size=0.08113,$OUT_NAME=out_name,ROTATION=0,X_FID=lon_id,$X_POS=x_pos[0], Y_FID=lat_id, Y_POS=y_pos[0]; 开始进行几何校正glt_raster = e.openraster(out_name)glt_fid=envirastertofid(glt_raster)envi_file_query,glt_fid,dims=dimsraster=e.openraster(file[i],DATASET_NAME='/PRODUCT/methane_mixing_ratio')   raster_id=envirastertofid(raster) envi_file_query,raster_id,ns=ns,nb=nbpos=indgen(nb)                    envi_doit,'ENVI_GEOREF_FROM_GLT_DOIT',FID=raster_id,GLT_DIMS=dims,GLT_FID=glt_fid,out_name=out_path,pos=pos  all_time = systime(2)print,'共花费时间：',(all_time-time)/60endelseendforend

然后我就尝试把构建GLT文件的代码封装成了一个函数，结果可以执行，哦了，问题解决~，上代码：

pro CH4_dealcompile_opt idl2e=envi(/headless)time = systime(1)path = 'G:\CH4\Oridata'file = file_search(path,'*.nc',count=n)for i=0,n-1 do beginname = file_basename(file[i],'.nc4')print,nameout_path='G:\CH4\result\'+name+'_result.dat'if file_test(out_path) eq 1 then beginprint,'已处理'endif else beginout_name =  build_glt(file[i],name)glt_raster = e.openraster(out_name)glt_fid=envirastertofid(glt_raster)envi_file_query,glt_fid,dims=dimsraster=e.openraster(file[i],DATASET_NAME='/PRODUCT/methane_mixing_ratio')   raster_id=envirastertofid(raster) envi_file_query,raster_id,ns=ns,nb=nbpos=indgen(nb)                    envi_doit,'ENVI_GEOREF_FROM_GLT_DOIT',FID=raster_id,GLT_DIMS=dims,GLT_FID=glt_fid,out_name=out_path,pos=pos  all_time = systime(2)print,'共花费时间：',(all_time-time)/60endelseendforendfunction build_glt,file,namecompile_opt idl2e=envi(/headless)latitude= e.openraster(file,DATASET_NAME='/PRODUCT/latitude')longitude= e.openraster(file,DATASET_NAME='/PRODUCT/longitude')out_name=e.gettemporaryfilename()lat_id=ENVIRASTERTOFID(latitude)lon_id=ENVIRASTERTOFID(longitude)envi_file_query,lat_id,nb=nb_laty_pos=lindgen(nb_lat)envi_file_query,lon_id,nb=nb_lonx_pos=lindgen(nb_lon)input_prj=ENVI_PROJ_CREATE(/geographic,datum='WGS-84')output_prj=ENVI_PROJ_CREATE(/geographic,datum='WGS-84')ENVI_DOIT,'ENVI_GLT_DOIT',dims=dims,I_PROJ=input_prj,O_PROJ=output_prj,pixel_size=0.081130,$OUT_NAME=out_name,ROTATION=0,X_FID=lon_id,$X_POS=x_pos[0], Y_FID=lat_id, Y_POS=y_pos[0]return,out_nameend

如有错误，欢迎指正！