使用IDL进行netcdf文件读取，出图保存，转geotiff格式

Author:HanDi
上海某高校遥感专业 工科男
相信积微者速成，相信分享的力量
CSDN@这可真是难为我了

这是一个老师布置的小任务。中间包括nc文件的读取，出图样式的调整，色带的自定义花了很长时间去学习，查资料和实验。最终，效果还算满意。希望大家能少走一些弯路，这只是个简单的小示例，更详细的函数参数以及细微调整需要大家自己去动手实践。

初次写技术博客没有什么经验，如有不足，还请大家多多包涵，欢迎大家私信交流 ^ _ ^ !!!一起进步~
以后会分享一些地理数据处理和数据挖掘的示例

下面是完整代码，我的IDL版本是 IDL8.5

这里是原始示例数据
链接：https://pan.baidu.com/s/1iWObxXAVOJwGSqva7Ug9BA
提取码：bb9z

;This demo provides an illustration on how to read and visualize AOD in
; SCHAP dataset and write data into GeoTIFF in IDL.

; The demo was prepared with IDL8.5.

; Last updated: 2021-11-01
PRO read_SCHAP_AOD

  ;Determine compilation rules
  compile_opt idl2

  ;Set output path
  envi, /restore_base_save_files
  envi_batch_init
  CD, 'D:\'
  outpath = 'D:\test\'
  ; read AOD data from NetCDF file


  filearr = dialog_pickfile(/multiple_files,  title = 'Open the original nc file')

  num = n_elements(filearr)


  FOR i=0,num-1,1 DO BEGIN

    file_ID = ncdf_open(filearr, /nowrite )


    latid = NCDF_VARID(file_ID,'lat')       ;read latitude from NetCDF file
    NCDF_VARGET, file_ID, latid, lat
    nlat = N_ELEMENTS(lat)

    lonid = NCDF_VARID(file_ID,'lon')       ;read  longitude from NetCDF file
    NCDF_VARGET, file_ID, lonid, lon
    nlon = N_ELEMENTS(lon)

    AODid = NCDF_VARID(file_ID,'AOD')        ; read AOD data
    NCDF_VARGET, file_ID, AODid, AOD
    nAOD = N_ELEMENTS(AOD)
    ;Obtain the gain, offset, and fill values of AOD data, and perform calibration
    ncdf_attget,file_ID,AODid,'scale_factor',a
    ncdf_attget,file_ID,AODid,'add_offset',b
    ncdf_attget,file_ID,AODid,'_FillValue',fv
    fv_index = where(AOD eq fv)
    output_AOD = (AOD ne fv)*AOD*a+b
    ; define fillvalue
    output_AOD[fv_index] =-999

    ;visualization

    loadct, 33
    TVLCT, r, g, b, /get
    color_table = BYTARR(3, 256)
    color_table[0, *] = r
    color_table[1, *] = g
    color_table[2, *] = b
    color_table[*, 0] = [255, 255, 255]  ;Custom colorbar

    img=image(transpose(output_AOD),rgb_table=color_table,title=file,grid_units=2,POSITION=[0.1,0.15,0.9,0.95],map_projection='geographic', image_dimensions=[max(lon)-min(lon),max(lat)-min(lat)] ,$
      IMAGE_LOCATION=[min(lon),min(lat)],DIMENSIONS=[691,512])

    ;set colormap limits
    img.MAX_VALUE=2
    img.MIN_VALUE=0
    ;Change the figure title, grid type, axis text direction, and color bar form
    img.title = 'Date:'+strmid(file_basename(filearr),16,8)

    img.mapgrid.label_position=0
    img.mapgrid.font_name='Palatino'

    img.mapgrid.linestyle=6
    img.mapgrid.horizon_thick=1


    lons=img.mapgrid.longitudes
    lats=img.mapgrid.latitudes
    for lons_i=0,n_elements(lons)-1 do begin
      lons[lons_i].label_angle=0
      lons[lons_i].label_align=0
    endfor

    ; add colorbar
    c = COLORBAR(TARGET=img, ORIENTATION=1,TITLE='AOD',POSITION=[0.905,0.25,0.925,0.85])
    c.RANGE=[0,2]
    c.BORDER=0
    c.TICKDIR= 1
    c.TEXTPOS = 1

    ;save plot in jpg format with resolution
    Img.save,outpath+'SCHAP.AOD.D001.A'+strmid(file_basename(filearr),16,8)+'.jpg'
    ; save as GeoTIFF
    ;Write geographic information structure
    geo_info={$
      MODELPIXELSCALETAG:[0.01,0.01,0.0],$
      MODELTIEPOINTTAG:[0.0,0.0,0.0,min(lon),max(lat),0.0],$
      GTMODELTYPEGEOKEY:2,$
      GTRASTERTYPEGEOKEY:1,$
      GEOGRAPHICTYPEGEOKEY:4326,$
      GEOGCITATIONGEOKEY:'GCS_WGS_1984'}
      
    ;write GeoTIFF
    WRITE_TIFF,outpath+'SCHAP.AOD.D001.A'+strmid(file_basename(filearr),16,8)+'.tif',reverse(transpose(output_AOD),2),/float, geotiff=geo_info
    ;Close nc file ID
    NCDF_CLOSE, file_ID
  endfor

END

tip

快捷键：Ctrl + Shift+ F
可以根据格式器的设置自动调整编码大小写同时自动对齐

代码运行效果

这里也要感谢
@环境遥感L.F点这里
@ 傻灰点这里，两位博主博客的帮助

如有错误，欢迎评论，私信，会及时改正