IDL 所会的重投影方法整合

最新推荐文章于 2024-05-30 10:13:55 发布

include_dx

最新推荐文章于 2024-05-30 10:13:55 发布

阅读量760

点赞数 2

文章标签：经验分享图像处理

本文链接：https://blog.csdn.net/m0_59609950/article/details/123700164

版权

前言

本人遥感小白，仅仅记录所会的重投影方法，若在某些地方说得有错误，请各位指正。

正文

在使用遥感影像的时候，往往需要对各类下载的文件读取遥感影像进行重投影。在我接触遥感的三年中也陆续接触了几种重投影方法（此处仅记录IDL方法），envi接口glt重投影，IDL中warp方法，包括自己也写过底层glt(后来因为总觉得有误差就放弃使用)。

前几天在刚学会warp方法后便萌生了重投影方法整合的念头，下面的IDL代码也因此诞生了。

代码

下列代码由几个功能函数组成，主函数为re_projection_by_lon_lat，因此各位保存的时候也要主要文件名称要与主函数一致。

该程序本质也为一个重投影函数，需要从原始资料提取到待重投影影像数组后，调用本函数进行重投影。以modis数据为例，所需重投影类型大致为两类：swath和grid，swath数据在hdf中会提供经度和纬度数据集，而grid则需要通过读取hdf源信息进行计算获得经纬度。详细处理方法在主程序帮助中也有说明。

;;coding=utf-8
pro warp_re_projection,input_img,lon_data,lat_data,out_path,pixel_size=pixel_size,percent=percent,method=method ;;两类warp方法，通过method指定，建议使用poly_warp方法
  data_size=size(input_img)
  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_i
  for line_i=0,data_size[2]-1 do line_ori[*,line_i]=line_i

  out_col_num=ceil((max(lon_data)-min(lon_data))/pixel_size)
  out_line_num=ceil((max(lat_data)-min(lat_data))/pixel_size)
  col_out=floor((lon_data-min(lon_data))/pixel_size)
  line_out=floor((max(lat_data)-lat_data)/pixel_size)
  
  if KEYWORD_SET(percent) eq 1 then begin
    control_point_num=floor(n_elements(input_img)*percent)
    if control_point_num lt 36 then begin
      print,'当前影像被选择影像像元个数不满足重投影条件,当前影像跳过'
      method='error 1'
    endif else begin
    selected_pos=findgen(control_point_num)/float(control_point_num)
    pixel_pos=floor(selected_pos*n_elements(input_img))
    col_ori=col_ori[pixel_pos]
    line_ori=line_ori[pixel_pos]
    col_out=col_out[pixel_pos]
    line_out=line_out[pixel_pos]
    endelse
  endif
  
  
  if method eq 'warp_tri' then begin
    out_img=warp_tri(col_out,line_out,col_ori,line_ori,input_img,OUTPUT_SIZE=[out_col_num,out_line_num])
  endif
  
  if method eq 'poly_warp' then begin
    polywarp,col_ori,line_ori,col_out,line_out,5,coe_x,coe_y
    out_img=poly_2d(input_img,coe_x,coe_y,0,out_col_num,out_line_num,missing=0.0)
  endif
  
  if out_img ne !null then begin
    geo_info={$
      MODELPIXELSCALETAG:[pixel_size,pixel_size,0.0],$
      MODELTIEPOINTTAG:[0.0,0.0,0.0,min(lon_data),max(lat_data),0.0],$
      GTMODELTYPEGEOKEY:2,$
      GTRASTERTYPEGEOKEY:1,$
      GEOGRAPHICTYPEGEOKEY:4326,$
      GEOGCITATIONGEOKEY:'GCS_WGS_1984',$
      GEOGANGULARUNITSGEOKEY:9102}
    WRITE_TIFF,out_path,out_img,/float,GEOTIFF=geo_info
  endif
  
end

pro swath_by_GLT,output_name,target_data,modis_lon,modis_lat,pixel_size=pixel_size  ;;通过envi接口实现GLT重投影，缺点速度较慢，优点则是内存优化很好
  compile_opt idl2
  envi,/restore_base_save_files
  envi_batch_init,/no_status_window

  success_value=1

  name=file_basename(output_name)
  output_directory=strmid(output_name,0,strlen(output_name)-strlen(name))

  out_lon=output_directory+'lon_out.tiff'
  out_lat=output_directory+'lat_out.tiff'
  out_target=output_directory+'out_target.tiff'

  write_tiff,out_lon,modis_lon,/float
  write_tiff,out_lat,modis_lat,/float
  write_tiff,out_target,target_data,/float

  envi_open_file,out_lon,r_fid=x_fid
  envi_open_file,out_lat,r_fid=y_fid
  envi_open_file,out_target,r_fid=target_fid

  out_name_glt=output_directory+'_glt.tiff'
  out_name_glt_hdr=output_directory+'_glt.hdr'

  int_proj=envi_proj_create(/geographic)
  out_proj=envi_proj_create(/geographic)

  envi_glt_doit,i_proj=int_proj,x_fid=x_fid,y_fid=y_fid,x_pos=0,y_pos=0,$
    o_proj=out_proj,pixel_size=pixel_size,rotation=0.0,out_name=out_name_glt,r_fid=glt_fid

  out_name_geo=output_directory+'_georef.img'
  out_name_geo_hdr=output_directory+'_georef.hdr'

  envi_georef_from_glt_doit,glt_fid=glt_fid,fid=target_fid,pos=0,out_name=out_name_geo,r_fid=geo_fid

  reslut_name=output_directory+'_geotiff.tiff'
  map_info=envi_get_map_info(fid=geo_fid)
  geo_zs=map_info.(1)
  pix_size=map_info.(2)
  envi_file_query,geo_fid,dims=dims
  target_data=envi_get_data(fid=geo_fid,pos=0,dims=dims)
  geo_info={$
    MODELPIXELSCALETAG:[pix_size[0],pix_size[1],0.0],$
    MODELTIEPOINTTAG:[0.0,0.0,0.0,geo_zs[2],geo_zs[3],0.0],$
    GTMODELTYPEGEOKEY:2,$
    GTRASTERTYPEGEOKEY:1,$
    GEOGRAPHICTYPEGEOKEY:4326,$
    GEOGCITATIONGEOKEY:'GCS_WGS_1984',$
    GEOGANGULARUNITSGEOKEY:9102,$
    GEOGSEMIMAJORAXISGEOKEY:6378137.0,$
    GEOGINVFLATTENINGGEOKEY:298.25722}

  write_tiff,output_name,target_data,/float,geotiff=geo_info

  envi_file_mng,id=x_fid,/remove
  envi_file_mng,id=y_fid,/remove
  envi_file_mng,id=target_fid,/remove
  envi_file_mng,id=glt_fid,/remove
  envi_file_mng,id=geo_fid,/remove

  file_delete,[out_lon,out_lat,out_target,out_name_glt,out_name_glt_hdr,out_name_geo, out_name_geo_hdr]
end

pro grid_data,input_data,file_work,lon_data=lon_data,lat_data=lat_data  ;;当重投影文件是正弦投影时，则可使用grid_file关键字，读取hdf文件内原始坐标信息，计算得到经纬度查找表
  sd_id=hdf_sd_start(file_work,/read)
  gendex=hdf_sd_attrfind(sd_id,'StructMetadata.0')
  hdf_sd_attrinfo,sd_id,gendex,data=yuan_info
  hdf_sd_end,sd_id

  ul_start_pos=strpos(yuan_info,'UpperLeftPointMtrs')
  ul_end_pos=strpos(yuan_info,'LowerRightMtrs')
  ul_info=strmid(yuan_info,ul_start_pos,ul_end_pos-ul_start_pos)
  ul_info_sp=strsplit(ul_info,'=(,)',/extract)
  ul_geo_x=double(ul_info_sp[1])
  ul_geo_y=double(ul_info_sp[2])

  lr_start_pos=strpos(yuan_info,'LowerRightMtrs')
  lr_end_pos=strpos(yuan_info,'Projection')
  lr_info=strmid(yuan_info,lr_start_pos,lr_end_pos-lr_start_pos)
  lr_info_sp=strsplit(lr_info,'=(,)',/extract)
  lr_geo_x=double(lr_info_sp[1])
  lr_geo_y=double(lr_info_sp[2])

  sin_prj=map_proj_init('sinusoidal',/gctp,sphere_radius=6371007.181,center_longitude=0.0,false_easting=0.0,false_northing=0.0)
  data_size=size(input_data)
  sin_pix=(lr_geo_x-ul_geo_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 begin
    proj_x[col_i,*]=ul_geo_x+sin_pix*col_i
  endfor
  for line_i=0,data_size[2]-1 do begin
    proj_y[*,line_i]=ul_geo_y-sin_pix*line_i
  endfor

  geo_loc=map_proj_inverse(proj_x,proj_y,map_structure=sin_prj)

  geo_x=geo_loc[0,*]
  geo_y=geo_loc[1,*]
  
  data_size=size(input_data)
  lon_data=reform(geo_x,data_size[1],data_size[2])
  lat_data=reform(geo_y,data_size[1],data_size[2])
end

pro limit_lon_lat_function,input_img=input_img,lon_data=lon_data,lat_data=lat_data,img_limited=img_limited,success_value=success_value 
    ;;当使用img_limit限制时使用，进行图像经纬度范围的限制，既可方便图像大小存储，也提升重投影速度
  pos=where((lon_data ge img_limited['lon_min']) and (lon_data le img_limited['lon_max']) and $
    (lat_data ge img_limited['lat_min']) and (lat_data le img_limited['lat_max']),pos_n)
  if pos_n ne 0 then begin
    pos_col_line=array_indices(lon_data,pos)
    col_min=min(pos_col_line[0,*])
    col_max=max(pos_col_line[0,*])
    line_min=min(pos_col_line[1,*])
    line_max=max(pos_col_line[1,*])
    if (col_min eq col_max) or (line_min eq line_max) then begin
      success_value=0
    endif else begin
      lon_data=lon_data[col_min:col_max,line_min:line_max]
      lat_data=lat_data[col_min:col_max,line_min:line_max]
      input_img=input_img[col_min:col_max,line_min:line_max]
    endelse
    
  endif
  if pos_n eq 0 then begin
    success_value=0
  endif
end

pro re_projection_by_lon_lat,input_img,lon_data,lat_data,out_img_path=out_img_path,method=method,pixel_size=pixel_size,$
  img_limited=img_limited,help=help,pixel_unit=pixel_unit,percent=percent,grid_file=grid_file  ;;重投影函数的主函数，各类关键字使用方法可查看下列帮助关键字中内容
  
  if KEYWORD_SET(help) eq 1 then begin
    print,'本函数功能在于进行数据重投影,函数输入元素如下'
    print,'input_img为待投影的影像数据'
    print,'lon_data,lat_data分别为重投影的经度数据和纬度数据'
    print,'out_img_path关键字为输出影像的保存路径,此关键字必须使用'
    print,'method关键字为指定重投影方法,默认使用envi接口glt重投影,另还可以使用warp_tri或ploy_warp,建议使用glt或poly_warp方法'
    print,'pixel_size关键字为设置影像分辨率,默认单位为degree,重投影方法中的两种warp均必须规定像元分辨率,若没有设置将使用glt方法重投影'
    print,'img_limited关键字为影像重投影范围字典,字典关键字为lon_max,lon_min,lat_max,lat_min,不使用该关键字时将默认进行全幅重投影'
    print,'pixel_unit关键字为指定空间分辨率单位,degree或km,分别为度、千米,默认为度'
    print,'percent关键字为使用两类wrap时选择参与地形校正公式像元数量的百分比,范围为0.0至1.0'
    print,'grid_file关键字为当重投影类型为正弦投影时使用，输入grid文件路径,当使用该关键字时,则不用输入lon_data,lat_data数据'
  endif
  
  success_value=1
  if KEYWORD_SET(out_img_path) eq 0 then begin
    success_value=-1
    goto,tip2
  endif
  
  method_list=['glt','warp_tri','poly_warp']
  
  if KEYWORD_SET(method) eq 0 then begin
    method='glt'
  endif else if total(method_list eq method) eq 0 then begin
    method='glt'
  endif else begin
    ;; 此行空格保留
  endelse
  
  if KEYWORD_SET(pixel_size) eq 0 then method='glt'
  
  if KEYWORD_SET(pixel_unit) eq 0 then pixel_unit='degree'
  
  if pixel_unit eq 'km' then pixel_size=360.0/40076.0*pixel_size  ;;利用投影分辨率确定经纬度分辨率
  
  if KEYWORD_SET(grid_file) eq 1 then begin
    grid_data,input_img,grid_file,LON_DATA=lon_data,LAT_DATA=lat_data
  endif
  
  
  if (KEYWORD_SET(img_limited) eq 1) then begin
    limit_lon_lat_function,INPUT_IMG=input_img,LON_DATA=lon_data,LAT_DATA=lat_data,IMG_LIMITED=img_limited,success_value=success_value
  endif
  if success_value eq 0 then goto,tip1
  
  
  if (method eq 'poly_warp') or (method eq 'warp_tri') then begin
    if KEYWORD_SET(percent) eq 0 then percent=0.5
    if  (percent gt 1) or (percent le 0) then percent=0.5
    warp_re_projection,input_img,lon_data,lat_data,out_img_path,PIXEL_SIZE=pixel_size,method=method,percent=percent
  endif
  
  if (method eq 'glt') and (input_img ne !null) then begin
    swath_by_GLT,out_img_path,input_img,lon_data,lat_data,pixel_size=PIXEL_SIZE
  endif
  
  tip1:
  if success_value eq 0 then begin
    print,'该限制经纬度范围不满足HDF文件要求，无法输出'+file_basename(out_img_path)
  endif
  tip2:
  if success_value eq -1 then begin
    print,'当前未设置输出保存路径,请设置后运行'
  endif
end