IDL遥感应用入门(14):IDL读取netcdf数据(.nc)

netcdf数据也是遥感应用中经常见到的一种图像格式,关于这种影像格式的一些介绍,你可以看这篇文章:关于netcdf文件(.nc)的一点概述,在这篇文章中,有一个pdf文件,专门介绍IDL读取netcdf的,本文还是以实际实例为基础,对IDL读取netcdf做一个简单的介绍。

一些说明

本文所用到的netcdf数据,你可以从这里下载的到。链接:http://pan.baidu.com/s/1kTgLjIf 密码:x53w

详细代码

;Author:Sailor
;2013-12-25
PRO Course_13

;定义文件路径
MyRootDir='D:\10\'

;遍历文件夹
filearr = file_search(MyRootDir,'*.nc',count=num);
  FOR fileindex=0,num-1,1 DO BEGIN
    nid = ncdf_open(filearr[fileindex], /nowrite )

    ; inquire about this file; returns structure
    file_info = ncdf_inquire( nid )

    ; print out the dimensions of this file
    FOR dimid=0, file_info.ndims -1 DO BEGIN
      ncdf_diminq, nid, dimid, name, size
      print, ' ---> dimension ' + name  + ' is: ', size 
    ENDFOR

    FOR varid=0, file_info.nvars-1 DO BEGIN
      ; inquire about the variable; returns structure
      var = ncdf_varinq( nid, varid )
      print,var
      print,'========================'
      ;read all attributes
      FOR var_att_id=0,var.natts -1 DO BEGIN
        att_name = ncdf_attname( nid, varid, var_att_id )
        print,att_name
        ncdf_attget, nid, varid, att_name, tematt
        print,string(tematt)
      ENDFOR
    ENDFOR

    ;read sst
    sstid = ncdf_varid(nid, 'analysed_sst')
    ncdf_varget, nid, sstid, sst
    sst=sst*0.01+273.150
    landindex=where(sst eq min(sst))
    sst[landindex]=max(sst)+1
    Img=image(sst,rgb_table=39,title=filearr(fileindex),grid_units=1,POSITION=[0.1,0.2,0.9,0.9])
    xaxis=axis('X',LOCATION=[0,0],AXIS_RANGE=[-180,180],MINOR=0, MAJOR=19,COORD_TRANSFORM=[-180,360.0/4096.0],title='Longitude(°)')
    yaxis=axis('Y',LOCATION=[0,0],AXIS_RANGE=[-90,90],MINOR=0, MAJOR=7,COORD_TRANSFORM=[-90,180.0/2048.0],title='Latitude(°)')
    c1 = COLORBAR(TARGET=Img, ORIENTATION=0,TITLE='(K) ',POSITION=[0.1,0.1,0.9,0.15])

  ENDFOR
END

效果图

代码解释

;定义文件路径
MyRootDir='D:\10\'

这里定义了文件路径,当然你要改成你自己的。

filearr = file_search(MyRootDir,'*.nc',count=num);
FOR fileindex=0,num-1,1 DO BEGIN

这句代码段是对上面的路径做一个遍历,找出所有的.nc文件,之所以这么写是为了方便批处理,因为不可能只处理一个文件。不过如果仅仅是做测试,这里可以改为num-num只处理一个文件测试效果,效果满意再做批处理。关于file_search函数的详细使用,你可以看这篇文章:IDL使用file_search函数遍历某个文件夹下的所有文件

; inquire about this file; returns structure
file_info = ncdf_inquire( nid )

这里是读取了整个nc文件的文件信息,返回的是一个结构体。详细参数有:

Ndims:文件的维数
Nvars:netcdf文件定义的变量
Ngatts:netcdf文件的全局属性
RecDim:未限制维数的ID,如果没有则为-1

; print out the dimensions of this file
FOR dimid=0, file_info.ndims -1 DO BEGIN
  ncdf_diminq, nid, dimid, name, size
  print, ' ---> dimension ' + name  + ' is: ', size 
ENDFOR

这段代码是读取文件的的所有维的信息,包括名字和大小。

FOR varid=0, file_info.nvars-1 DO BEGIN
  ; inquire about the variable; returns structure
  var = ncdf_varinq( nid, varid )
  print,var
  print,'========================'
  ;read all attributes
  FOR var_att_id=0,var.natts -1 DO BEGIN
    att_name = ncdf_attname( nid, varid, var_att_id )
    print,att_name
    ncdf_attget, nid, varid, att_name, tematt
    print,string(tematt)
  ENDFOR
ENDFOR

读取所有变量和变量的属性信息。

;read sst
sstid = ncdf_varid(nid, 'analysed_sst')
ncdf_varget, nid, sstid, sst
sst=sst*0.01+273.150
landindex=where(sst eq min(sst))
sst[landindex]=max(sst)+1
Img=image(sst,rgb_table=39,title=filearr(fileindex),grid_units=1,POSITION=[0.1,0.2,0.9,0.9])
xaxis=axis('X',LOCATION=[0,0],AXIS_RANGE=[-180,180],MINOR=0, MAJOR=19,COORD_TRANSFORM=[-180,360.0/4096.0],title='Longitude(°)')
yaxis=axis('Y',LOCATION=[0,0],AXIS_RANGE=[-90,90],MINOR=0, MAJOR=7,COORD_TRANSFORM=[-90,180.0/2048.0],title='Latitude(°)')
c1 = COLORBAR(TARGET=Img, ORIENTATION=0,TITLE='(K) ',POSITION=[0.1,0.1,0.9,0.15])

最后出图,其中这里数据处理所做的处理来自文件本身,也就是在上面读取所有变量和变量的属性信息所读取出的数据。

一些说明

这里只是对IDL读取netcdf数据的读取做一个示例,如有问题,欢迎交流。

关于netcdf文件(.nc)的一点概述 利用IDL中的SHIFT函数实现经度翻转

作者:,GIS爱好者。
分享本文,请您带上本文链接
分享到:

已有 4 条评论

  1. 亖木子
    1#
    亖木子  · 2015-01-13 17:44

    很有帮助。我想知道,如果我想把每一层单独存储成tif,后面怎么写。我写出来的没有坐标信息,而且纬度是颠倒的。

  2. hellosyndy
    2#
    hellosyndy  · 2016-02-21 14:22

    请问,在给nid赋值时nid = ncdf_open(filearr[fileindex], /nowrite ),运行会报Attempt to subscript FILEARR with FILEINDEX is out of range错误怎么解决?

  3. hikari
    3#
    hikari  · 2016-11-21 00:36

    请问 sst[landindex]=max(sst)+1 这部分有什么意义

      1. Sailor
        Sailor  · 2016-11-21 09:00

        landindex表示陆地,这里表示将陆地的值统一成所以数据的最大值+1,可以将陆地对比显示出来。

发表评论