GLT几何校正

基于GLT方法的风云三号C（FY-3C）卫星影像几何校正_ENVI-IDL技术殿堂_新浪博客科学网-ENVI二次开发时ENVI_GET_DATA等函数的编译报错问题-董彦卿的博文基于GLT方法的风云三号C（FY-3C）卫星影像几何校正_ENVI-IDL技术殿堂_新浪博客

COMPILE_OPT IDL2

玄学第一次可以运行然后一直报错

 

 

pro glt
  
  COMPILE_OPT IDL2

  e=envi(/headless)
  path = 'E:\data\fy'
  files = file_search(path,'*.HDF')
  ;help,files
  ;print,files
  file = files[1]
  ;print,file
  lat = e.OpenRaster(file,DATASET_NAME='/Latitude')
  lon = e.OpenRaster(file,DATASET_NAME='/Longitude')
  ;help,lon
  bt = e.OpenRaster(file,DATASET_NAME='/EARTH_OBSERVE_BT_10_to_89GHz')
  
  
  out_name=e.gettemporaryfilename() ;选择创建临时文件输出路径
  lat_id=ENVIRASTERTOFID(lat)
  lon_id=ENVIRASTERTOFID(lon)
  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.100722,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文件
  glt_fid=envirastertofid(glt_raster)
  envi_file_query,glt_fid,dims=dims
  raster=e.openraster(file,DATASET_NAME='/EARTH_OBSERVE_BT_10_to_89GHz')    ;打开目标数据集

  raster_id=envirastertofid(raster)
  envi_file_query,raster_id,ns=ns,nb=nb
  pos=indgen(nb)
  
  out_path = 'E:\data\ttt\result.tif'
  
  envi_doit,'ENVI_GEOREF_FROM_GLT_DOIT',FID=raster_id,GLT_DIMS=dims,GLT_FID=glt_fid,out_name=out_path,pos=pos


  
  
  
end

pro tiqu
  e=envi(/headless)
  fn='C:\Users\Administrator\Desktop\test.csv';
  data=READ_CSV(fn,count=nsta,header=header);
  point_fid=data.(0);
  Lat=data.(1);站点纬度
  Lon=data.(2);站点经度
  image_fname='‪E:\data\gltfy.tif';文件路径
  ;打开文件
  ENVI_OPEN_FILE,image_fname,r_fid=fid
  ;print ,fid
  map_info=ENVI_GET_MAP_INFO(fid = fid)
  ;print,map_info
  i_proj=ENVI_PROJ_CREATE(/geographic,datum='WGS-84')
  o_proj=map_info.proj
  envi_convert_projection_coordinates,lon,lat,i_proj,$
    xmap,ymap,o_proj
  print,lon,lat,i_proj,xmap,ymap,o_proj
  ;将站点的地图坐标转换为文件坐标
  envi_convert_file_coordinates,fid,xf,yf,xmap,ymap
  xf=floor(xf)
  yf=floor(yf)
  ;循环读取各个站点dn
  value=fltarr(nsta)
  for i=0, nsta-1 do begin
    dims1=[-1, xf[i], xf[i], yf[i], yf[i]]  ;构建空间范围数组DIMS
    value[i]=envi_get_data(fid=fid,dims=dims1,pos=0)
  endfor

  print,value
  
  e.close
end

 

 

 提取这个代码也是突然就一直报错

pro extract_point
  
  COMPILE_OPT IDL2,hidden
  
  e=envi(/headless)
  fn='C:\Users\Administrator\Desktop\tibtation.csv';
  data=READ_CSV(fn,count=nsta,header=header);
  point_fid=data.(0);
  Lat=data.(1);站点纬度
  Lon=data.(2);站点经度
  image_fname='C:\Users\Administrator\Desktop\arc\wendut.tif';文件路径
  ;打开文件
  ENVI_OPEN_FILE,image_fname,r_fid=fid
  ;print ,fid
  map_info=ENVI_GET_MAP_INFO(fid=fid)
  ;print,map_info
  i_proj=ENVI_PROJ_CREATE(/geographic,datum='WGS-84')
  o_proj=map_info.proj
  envi_convert_projection_coordinates,lon,lat,i_proj,$
    xmap,ymap,o_proj
  print,lon,lat,i_proj,xmap,ymap,o_proj
  ;将站点的地图坐标转换为文件坐标
  envi_convert_file_coordinates,fid,xf,yf,xmap,ymap
  xf=floor(xf)
  yf=floor(yf)
  ;循环读取各个站点dn
  value=fltarr(nsta)
  for i=0, nsta-1 do begin
    dims1=[-1, xf[i], xf[i], yf[i], yf[i]]  ;构建空间范围数组DIMS
    value[i]=envi_get_data(fid=fid,dims=dims1,pos=0)
  endfor

  ;***********保存结果*********
  o_fn='C:\Users\Administrator\Desktop\point_IDL_extract_data1.csv';指定存储结果路径
  ;header=[header,strmid(image_fname,28,4)]
  ;print,value
  header=[header,'wendudu']
  data=create_struct(data,'wendu',value);在data中增加一个value域
  write_csv,o_fn,data,header=header;写入csv
  e.close
end

 

pro tiqu

  COMPILE_OPT IDL2,hidden

  e=envi(/headless)
  fn='C:\Users\Administrator\Desktop\test.csv';
  data=READ_CSV(fn,count=nsta,header=header);
  print,data
  point_fid=data.field1;
  Lat=data.(1);站点纬度
  Lon=data.field3;站点经度
  image_fname='‪E:\data\gltfy.tif';文件路径
  ;打开文件
  ENVI_OPEN_FILE,image_fname,r_fid=fid
  ;print ,fid
  map_info=ENVI_GET_MAP_INFO(fid = fid)
  ;print,map_info
  i_proj=ENVI_PROJ_CREATE(/geographic,datum='WGS-84')
  o_proj=map_info.proj
  envi_convert_projection_coordinates,lon,lat,i_proj,xmap,ymap,o_proj
  ;print,lon,lat,i_proj,xmap,ymap,o_proj
  ;将站点的地图坐标转换为文件坐标
  envi_convert_file_coordinates,fid,xf,yf,xmap,ymap
  xf=floor(xf)
  yf=floor(yf)
  ;循环读取各个站点dn
  value=fltarr(nsta)
  for i=0, nsta-1 do begin
    dims1=[-1, xf[i], xf[i], yf[i], yf[i]]  ;构建空间范围数组DIMS
    value[i]=envi_get_data(fid=fid,dims=dims1,pos=0)
  endfor

  print,value
  
  e.close
end

批处理，会报错

批处理正确

pro fy3mwri_deal
  compile_opt idl2
  e=envi(/headless)
  time = systime(1)
  path = 'E:\data\fy'
  file = file_search(path,'*.HDF',count=n)
  for i=0,n-1 do begin
    name = file_basename(file[i],'.HDF')
    print,name
    out_path='E:\data\ttt\'+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='/EARTH_OBSERVE_BT_10_to_89GHz')

      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
      print,file[i]+" is ok!"
      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='/Latitude')
  longitude= e.openraster(file,DATASET_NAME='/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.100722,$
    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