了解知识
Netcdf是一种数据存储格式,其特点自描述、便携式、易提取、可追加、能共享。其实最大的特点就是易提取数据,打个比方,比如在一座图书馆中,并排放了100个书架,但书架没有编号,让你从中找第8888本书,你只能从第一个书架的第一本书一个一个向下数,直到你找到那本书为止;但如果每个书架、每层、每本都有编号,那就容易找了。前一种找书的方式就是我们通常在二进制或ASCII码文件中用sequent方式读数据,后一种方式就是读Netcdf文件了。当然,前一种方式写程序比较简单,后一种方式因为要知道其是如何编码的,如每个书架有几层,每层放几本,程序写起来要麻烦点。所以,要想找的快,总要付出代价的。
1.读取变量的全部内容
在linux下我通常是用ncdump和ncview先看查看netcdf文件,确定要读取的变量及其维数后,再用fortran读,只要三条语句就可以搞定了(下面程序中的红色部分)。
下面是一个简单的Fortran程序读取wrfout文件,程序名为read_netcdf.f90
! This program is for reading wrfout 
      implicit none
      integer, parameter :: LLm=144, MMm=96, Ndas=25   !定义维数
      integer ierr,ncid,varid,len_file
      real  h(LLm,MMm,Ndas)                              !变量名:地形高度
      character*299 file_wrfout                            !文件名路径
      include 'netcdf.inc'
      file_wrfout='/data/zlzang/fly_wrfout/wrfout_d03_2010-06-03_00:00:00'
 
      len_file=len_trim(file_wrfout)
      ierr=nf_open(trim(file_wrfout),nf_nowrite,ncid)    !打开netcdf文件,获取文件的ID号(ncid)
      ierr=nf_inq_varid (ncid, 'HGT', varid)                  !打开指定变量名'HGT',获取该变量的ID号
      ierr=nf_get_var_real (ncid,varid,h)                     !读取该变量的内容,这里已知h的维数和变量类型
      print*,'ierr'
      print*, h(1,1,1),h(144,1,1)                                   !检查读取的内容
      end
!-----------------------------------------------------------------------------------------------------------------------
对上述read_netcdf.f90程序的编译也只用一行就就行了:
ifort   -I/home/zlzang/local/include read_netcdf.f90 -o a.exe  -L/home/zlzang/local/lib -lnetcdf
运行上面的脚本会生成一个a.exe文件,运行后就可以看到print的结果了,之所以把ierr输出来,是要看看其是否为零,只有ierr=0才表明程序正常执行了,但结果不一定是正确的啊。
 
但如果对变量的维数和变量类型不是很清楚(我就经常对维数的排序和大小犯迷糊),则上面的程序运行时不一定会报错(只有当维数的大小比实际的维数偏小时才会报错),但读出来的数据会有位移,所以无论如何都要对读出来的数据进行检查核对。
如果要知道变量的维数,则要借助于nf_inq_var和nf_inq_dim两个函数,先用第一个函数,将上述程序改写为:
      implicit none
      integer, parameter :: LLm=145, MMm=96, Ndas=25
      integer ierr,ncid,varid,len_file,xtype,ndims,dimids(3),natts
      real  h(LLm,MMm,Ndas)
      character*299 file_wrfout
      character*10 vname
      include 'netcdf.inc'
      file_wrfout='/data/zlzang/fly_wrfout/wrfout_d03_2010-06-03_00:00:00'
 
      len_file=len_trim(file_wrfout)
      ierr=nf_open(trim(file_wrfout),nf_nowrite,ncid)
      ierr=nf_inq_varid (ncid, 'HGT', varid)
      ierr=nf_inq_var(ncid,varid,vname,xtype,ndims,dimids,natts)   !获取变量信息
      print*,'ierr=',ierr                   
      print*,'vname=',vname     !变量名
      print*,'xtype=',xtype          !变量类型,4表示整型,5表示实型,6表示双精度
      print*,'ndims=',ndims        !变量维数
      print*,'dimids=',dimids      !每一维的ID
      end
!-------------------------------------------------------------------------------------------------
编译运行后会输出以下结果:
 ierr=           0
 vname=HGT       
 xtype=           5
 ndims=           3
 dimids=           3           4           1
其实上面这个程序在运行之前,我们可能并不知道HGT这个变量是几维的,所以不能根据ndims的大小给数定义dimids,我也是先试运行了一下,才定义dimids的维数为3的,但实际上你可以将dimids的维数定义的大一些,例如定义为dimids(10),如果它只是三维的,输出的dimids后面7个数会显示为0.
有了变量的各维的ID之后就可以用nf_inq_dim函数进一步获取其维数的大小了,同时还可以知道其维数的名字。在上面程序的基础上追加以下几行:
      do i=1,ndims
         ierr=nf_inq_dim(ncid,dimids(i),dimname,len)
         print*,'name of',i,'is  ',  dimname,' length=',len
      enddo
别忘了在程序的最前面加上新变量的定义:
      integer i,len
      character*15 dimname
运行这个修改后的程序可以输出如下结果:
 ierr=           0
 vname=HGT       
 xtype=           5
 ndims=           3
 dimids=           3           4           1           0           0
 name of           1 is  west_east       length=         144      
 name of           2 is  south_north     length=          96
 name of           3 is  Time            length=          25
最后三行即为各维数的名字和大小,这样就可以确定HGT这个变量的维数是(144,96,25)了。
 
2. 读取变量的部分内容
如果要读取部分的变量内容,则要用nf_get_vars_real这个函数,注意与读全部变量的函数(nf_get_var_real)几乎相同,就多个s。
下面是只读第一层的HGT,把变量h定义为一个二维数组:
!This program is only for reading the first layer of HGT
      implicit none
      integer, parameter :: LLm=145, MMm=96, Ndas=25
      integer ierr,ncid,varid,len_file,xtype,ndims,dimids(5),natts
      real  h(LLm,MMm)
      character*299 file_wrfout
      integer start(3),counts(3),strid(3)
      data start /1,1,1/                       !   开始位置
      data counts /144,96,1/             !   计数长度
      data strid /1,1,1/                       !   间隔步长
      include 'netcdf.inc'
 
      file_wrfout='/data/zlzang/fly_wrfout/wrfout_d03_2010-06-03_00:00:00'
      len_file=len_trim(file_wrfout)
      ierr=nf_open(trim(file_wrfout),nf_nowrite,ncid)
      ierr=nf_inq_varid (ncid, 'HGT', varid)
      ierr=nf_get_vars_real (ncid,varid,start,counts,strid,h)
      print*,ierr
      print*,h(1,1),h(144,1)
      end
注意在fortran中读netcdf默认的是从(1,1,1)开始的,而在matalab中是从(0,0,0)开始的。如果相读第二层的数据,则只要改一下start为/1,1,2/就可以了。
标签: linux Fortran
扩展知识