天天看點

cdf檔案 mysql,在linux下用fortran讀取netcdf檔案(以WRF模式輸出的資料為例)

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/就可以了。