信息技术改变生活

Xiaobo Wu'shared space

NetCDF转成ENVI标准格式

NetCDF(http://www.unidata.ucar.edu/software/netcdf/)格式是比较常见的地学数据存储格式,很多数据都是NetCDF格式存储,读写NetCDF有很多库,如C++,java,Fortran,但是,这些库都是底层,编译比较麻烦;其次,用底层的库将NetCDF格式和其他数据(如,shp矢量数据)结合处理是比较麻烦的,如Fortran,要处理结合shp和NetCDF数据还需要FortranGIS(http://fortrangis.sourceforge.net/)支持,同样需要编译;同时,像matlab和IDL是高级语言,能直接读写NetCDF,并且擅长于栅格矩阵运算,由于和ENVI结合紧密,IDL对shp的支持更优,所以,难免需要将NetCDF格式转成ENVI标准格式,用于显示和后期处理。

1.原始NetCDF【示例数据SA.nc下载地址】

Panoply(http://www.giss.nasa.gov/tools/panoply/)显示的信息如下

属性信息: 属性信息

空间信息: 空间信息

2.转换成ENVI标准格式信息,ENVI显示如下如图:

属性信息: 属性信息

空间信息: 空间信息

3.代码(需要ENVI+IDL运行 支持批处理)

;+
; ;**********************************************************
; Author: wuxb
; Email: wavelet2008@163.com
; NAME:
;    D:\workspace\Tech\Code\IDL\ENVI_IDL4.7\SoilTextureProduce\DataPreprocessing\NCData2ENVI.pro
; PARAMETERS:
;  1、NcPath
;  2、EnviPath
; Write by :
;    2015-4-23 9:29:02
;;MODIFICATION HISTORY:
;;   Modified  :
;***********************************************************
;-
PRO  NCDATA2ENVI, NcDirectory = NcDirectory, ENVIDirectory = ENVIDirectory

  ;**************************************************************************
  ;
  ; 1) Initialize the Common variables
  ;
  ;**************************************************************************
  ;just for Test
  NcDirectory     = 'D:\Temp\Data\Input\'
  ENVIDirectory   = 'D:\Temp\Data\Output\'
  
  ;**************************************************************************
  ;
  ; 2) Get the files to be convert
  ;
  ;**************************************************************************
  
  NCFilePaths = FILE_SEARCH(NcDirectory,'*.NC', COUNT = nCount, /TEST_READ, /FULLY_QUALIFY_PATH)
  IF(nCount LE 0) THEN RETURN
  
  PRINT,'Procedure begins at ' + STRING(SYSTIME(/UTC))
  
  FOR i=0, nCount-1 DO BEGIN
    ;Open Nc File
    NID         = NCDF_OPEN(NCFilePaths[i], /nowrite )
    ;Inquire about this file; returns structure
    File_Info   = NCDF_INQUIRE(NID)
    ;define the dimensions of this file
    NS = 0L
    NL = 0L
    NB = 0L
    ;**************************************************************************
    ; 1) read dim info
    ;**************************************************************************
    FOR dimid=0, FILE_INFO.NDIMS -1 DO BEGIN
      NCDF_DIMINQ, NID, Dimid, Name, Size
      SWITCH Name OF
        "lon": BEGIN
          NS = Size
          BREAK
        END
        "lat": BEGIN
          NL = Size
          BREAK
        END
        "depth": BEGIN
          NB = Size
          BREAK
        END
        ELSE: BEGIN
          PRINT, 'There is no Dim'
        END
      ENDSWITCH
    ENDFOR
    ;**************************************************************************
    ; 2) read data  info
    ;**************************************************************************
    ;define proj request
    StartX             = 0.0;
    StartY             = 0.0;
    EndX               = 0.0;
    EndY               = 0.0;
    GridSize           = 0.0;
    ndepth             = FLTARR(NB)
    FOR Varid=0, FILE_INFO.NVARS-1 DO BEGIN
      IsWRITE  = 1;
      ; inquire about the variable; returns structure
      VAR  =  NCDF_VARINQ( NID, Varid )
      ;read data
      dataID = NCDF_VARID(NID, VAR.NAME)
      NCDF_VARGET, NID, dataID, Data
      SWITCH VAR.NAME OF
        "lon": BEGIN
          StartX   = Data[0];
          nData    = N_ELEMENTS(Data)
          EndX     = Data(nData-1)
          IsWRITE  = 0
          BREAK
        END
        "lat": BEGIN
          StartY   = Data[0];
          nData    = N_ELEMENTS(Data)
          EndY     = Data(nData-1)
          IsWRITE  = 0
          BREAK
        END
        "depth": BEGIN
          ndepth   = Data
          IsWRITE  = 0
          BREAK
        END
      ENDSWITCH
      ;**************************************************************************
      ;1) IF No project then Creat project
      ;**************************************************************************
      ; Projection
      ; GridSize = 0  illustrate  iMap is null
      IF (GridSize EQ 0) AND IsWRITE THEN BEGIN
        GridSize  = (ROUND(EndX)   ROUND(StartX)) / FLOAT(NS)
        ; start point
        iProj     = ENVI_PROJ_CREATE(/GEOGRAPHIC)
        ; Create the map information
        ps        = [GridSize, GridSize]
        mc        = [0D, 0D, StartX, EndY]
        Datum     = 'WGS-84'
        Units     = ENVI_TRANSLATE_PROJECTION_UNITS ('Degrees')
        iMap      = ENVI_MAP_INFO_CREATE(/GEOGRAPHIC,MC=MC,PS=PS,PROJ=iProj,UNITS=Units)
      ENDIF
          ;**************************************************************************
          ; 2) Creat Envi File
          ;**************************************************************************
          IF IsWRITE THEN BEGIN
            OutFilePath = ENVIDirectory + VAR.NAME
            OPENW, HData, OutFilePath, /GET_LUN
            WRITEU, HData,REVERSE(Data,2)
            FREE_LUN, HData
            DATA_TYPE = 4
            DataType  = VAR.DATATYPE
            SWITCH VAR.DATATYPE OF
          "BYTE": BEGIN
            DATA_TYPE    = 1
            BREAK
          END
          "CHAR": BEGIN
            DATA_TYPE    = 1
            BREAK
          END
          "INT": BEGIN
            DATA_TYPE    = 2
            BREAK
          END
          "LONG": BEGIN
            DATA_TYPE    = 3
            BREAK
          END
          "FLOAT": BEGIN
            DATA_TYPE    = 4
            BREAK
          END
          "DOUBLE": BEGIN
           DATA_TYPE    = 5
            BREAK
          END
        ENDSWITCH
        ; Edit the envi header file
        ENVI_SETUP_HEAD,FNAME=OutFilePath,NS=NS,NL=NL,NB=NB,INTERLEAVE=0,$
          DATA_TYPE=DATA_TYPE,OFFSET=0,MAP_INFO=iMap,BNAMES=[VAR.NAME],/WRITE,$
          /OPEN,R_FID=Data_FID
        ;**************************************************************************
        ; 3) read all attributes update Herder file
        ;**************************************************************************
        ;Set the header value
        ENVI_ASSIGN_HEADER_VALUE, FID=Data_FID, KEYWORD ='Author',  VALUE = "Wuxb"
        ENVI_ASSIGN_HEADER_VALUE, FID=Data_FID, KEYWORD ='Soil Depth',  VALUE = ndepth, PRECISION=5
        FOR VarAttID =0, VAR.NATTS -1 DO BEGIN
          AttName = NCDF_ATTNAME( NID, Varid, VarAttID )
          NCDF_ATTGET, NID, varid, AttName, AttValue
          ENVI_ASSIGN_HEADER_VALUE, FID=Data_FID, KEYWORD = AttName,  VALUE = STRING(AttValue)
        ENDFOR
        ;Write the values to the file header
        ENVI_WRITE_FILE_HEADER, Data_FID
        ;remove  file
        ENVI_FILE_MNG, ID=Data_FID, /REMOVE
      ENDIF
    ENDFOR
  ENDFOR
  PRINT, 'Procedure ends at ' + SYSTIME(/UTC)
END