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

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

自带经纬度文件处理流程

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

1. 找到数据集中的经纬度文件。
2. 将经纬度文件制作成GLT文件。
3. 根据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_lat
  y_pos=lindgen(nb_lat)
  envi_file_query,lon_id,nb=nb_lon
  x_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参数时，会自己有个默认值。也可以自己算：公式如下：这个有点不准，还是建议手动
$$\frac{1°}{111.11KM}=\frac{卫星分辨率(十进制度)}{卫星分辨率(KM)}$$

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

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

      
      glt_raster = e.openraster(out_name)  ;打开GLT文件
      glt_fid=envirastertofid(glt_raster)		
      envi_file_query,glt_fid,dims=dims
      raster=e.openraster(file[i],DATASET_NAME='/PRODUCT/methane_mixing_ratio')    ;打开目标数据集
  
      raster_id=envirastertofid(raster) 
      envi_file_query,raster_id,ns=ns,nb=nb
      pos=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_deal
  compile_opt idl2
  e=envi(/headless)
  time = systime(1)
  path = 'G:\CH4\Oridata'
  file = file_search(path,'*.nc',count=n)
  for i=0,n-1 do begin
    name = file_basename(file[i],'.nc4')
    print,name
    out_path='G:\CH4\result\'+name+'_result.dat'
    ; 重复处理数据不进行处理
    if file_test(out_path) eq 1 then begin
      print,'已处理'
    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_lat
      y_pos=lindgen(nb_lat)
      envi_file_query,lon_id,nb=nb_lon
      x_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=dims
      raster=e.openraster(file[i],DATASET_NAME='/PRODUCT/methane_mixing_ratio')   
      raster_id=envirastertofid(raster) 
      envi_file_query,raster_id,ns=ns,nb=nb
      pos=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)/60
    endelse
  endfor
  
end

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

pro CH4_deal
  compile_opt idl2
  e=envi(/headless)
  time = systime(1)
  path = 'G:\CH4\Oridata'
  file = file_search(path,'*.nc',count=n)
  for i=0,n-1 do begin
    name = file_basename(file[i],'.nc4')
    print,name
    out_path='G:\CH4\result\'+name+'_result.dat'
    if file_test(out_path) eq 1 then begin
      print,'已处理'
    endif else begin
      
      out_name =  build_glt(file[i],name)
  
      glt_raster = e.openraster(out_name)
      glt_fid=envirastertofid(glt_raster)
      envi_file_query,glt_fid,dims=dims
  
      raster=e.openraster(file[i],DATASET_NAME='/PRODUCT/methane_mixing_ratio')   
  
      raster_id=envirastertofid(raster) 
      envi_file_query,raster_id,ns=ns,nb=nb
      pos=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)/60
    endelse
  endfor
  
end

function build_glt,file,name
  compile_opt idl2
  e=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_lat
  y_pos=lindgen(nb_lat)
  envi_file_query,lon_id,nb=nb_lon
  x_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_name
 end

如有错误，欢迎指正！