python中提供了多種方式來處理netcdf文件,這裡主要講一下常用的 netcdf4-python 模塊。
netcdf4-python是 netCDF4 C庫的python模塊。新版本(V4)的 netcdf 中有很多以前版本中沒有的優點,而且新版本是在 HDF5 上建立的。此模塊可以讀寫 netCDF4 及 netCDF3 格式的文件,同時也可創建 HDF5 客戶端只讀的文件。
netCDF4 格式的許多特徵都實現了,比如:多個無限維度,組(groups)及zlib數據壓縮。除此之外,所有新的數據類型(64-bit 和 無符號整型)也已實現。同時也支持複數(struct),變長及枚舉等數據類型,但不支持 opaque 數據類型,同時也不支持複數,變長型和枚舉型數據的組合,比如複數數據中包含枚舉型數據,或是變長類型數據中包含複數數據。
創建,打開和關閉 netCDF 文件
通過調用 Dataset 構造器可以創建 netCDF 文件,同時也可以用來打開已存在的文件。當然也可以執行寫入操作,可以寫入包括維度,組,變量和屬性的數據類型。
netCDF文件包含五種類型的數據,即 NETCDF3_CLASSIC,NETCDF3_64BIT_OFFSET, NETCDF3_64BIT_DATA, NETCDF4_CLASSIC,NETCDF4.
NETCDF3_CLASSIC 是原始的netcdf二進位文件,此格式文件大小不能超過 2 G;V3.6版本庫引入了 NETCDF3_64BIT_OFFSET 格式,從而使其大小可以超過 2 G;NETCDF3_64BIT_DATA 是V4.4.0版本庫引入的,其擴展了NETCDF3_64BIT_OFFSET 格式,從而可以支持無符號64位整型數據及64位維度大小;NETCDF3_64BIT 是 NETCDF3_64BIT_OFFSET 格式的別名;NETCDF4_CLASSIC 使用了V4磁碟格式(HDF5),但是忽略了V3 API 中沒有的特徵。只有當重新連結 netcdf 庫時,才可以通過 netCDF3 客戶端讀取,同時也可以通過HDF5客戶端讀取。
netCDF4 模塊可以讀取和寫入上述格式中的文件。當創建文件時,可以通過 Dataset 構造器的 format 關鍵詞參數指定格式。默認的格式是 NETCDF4。通過查看 data_model 屬性可以確定文件的格式。通過 Dataset 實例的 close 方法可以關閉已經打開的 netcdf 文件。
>>> from netCDF4 import Dataset
>>> rootgrp = Dataset("test.nc", "w", format="NETCDF4")
>>> print rootgrp.data_model
NETCDF4
>>> rootgrp.close()
如果使用 url 代替傳給 Dataset 構造器的文件名,那麼就可以獲取遠程數據集了。前提是支持獲取遠程數據功能。
netcdf 文件中的 Groups
版本4的 netcdf 支持按層級來劃分數據,這類似文件系統中的目錄。Groups 可以包含變量,維度和屬性,同時也可以包含其他 groups。Dataset 回創建一個特殊的組,即 root group,類似unix文件系統中的根目錄。通過 Dataset 或 Group 實例的 createGroup 方法可以創建 Group 實例。createGroup 方法接受一個包含組名的 python 字符串。通過Dataset 實例的 groups 目錄屬性可以獲取包含在 root group 中的新 Group 實例的信息。註:只有 NETCDF4 格式文件支持 Groups,使用其他格式創建Group時會報錯。
>>> rootgrp = Dataset("test.nc", "a")
>>> fcstgrp = rootgrp.createGroup("forecasts")
>>> analgrp = rootgrp.createGroup("analyses")
>>> print rootgrp.groups
OrderedDict([("forecasts",
<netCDF4._netCDF4.Group object at 0x1b4b7b0>),
("analyses",
<netCDF4._netCDF4.Group object at 0x1b4b970>)])
組是可以進行嵌套的,也就是說一個組可以包含其他組,這就像unix文件系統的目錄下可以包含其他目錄一樣。每一個Group都包含一個 groups 屬性字典,其中包含了該組中包含的所有組實例。同樣,每一個組實例都包含了path屬性,指明了group所處的 」路徑「(類似unix文件系統中路徑)。為了簡化創建嵌套組的方式,可以傳遞類似unix文件系統中路徑的方式給 createGroup 方法:
>>> fcstgrp1 = rootgrp.createGroup("/forecasts/model1")
>>> fcstgrp2 = rootgrp.createGroup("/forecasts/model2")
如果路徑中不存在其他中間元素的話,創建組時就像在unix系統中使用 mkdir -p 命令創建路徑一樣。如果你試著創建已存在的組的話,不會導致錯誤,只會返回已存在的組。
下面是一個瀏覽 Dataset 中所有組的例子。函數 walktree 是一個生成器函數,用來遍歷目錄樹。注意所列印出的組信息。
>>> def walktree(top):
>>> values = top.groups.values()
>>> yield values
>>> for value in top.groups.values():
>>> for children in walktree(value):
>>> yield children
>>> print rootgrp
>>> for children in walktree(rootgrp):
>>> for child in children:
>>> print child
<type "netCDF4._netCDF4.Dataset">
root group (NETCDF4 file format):
dimensions:
variables:
groups: forecasts, analyses
netcdf 文件中的維度
netcdf根據維度信息創建所有變量的大小,所以在創建變量之前必須要創建維度信息。對於標量變量來說,不需要維度信息。Dateset 或 Group 實例的 createDimension 方法可以用以創建一個維度,傳遞給此方法的 python 字符串和整數用來表示維度的名稱和大小。如果要創建無限維度(即可以隨時添加數據),可以將大小設置為 None 或 0。
下例中, time 和 level 都是無限維變量。可以包含多個無限維變量是netcdf 的一個新特徵,之前的netcdf僅支持包含一個無限維變量,而且必須要包含在最左邊,即第一個維度。
>>> level = rootgrp.createDimension("level", None)
>>> time = rootgrp.createDimension("time", None)
>>> lat = rootgrp.createDimension("lat", 73)
>>> lon = rootgrp.createDimension("lon", 144)
所有的 Dimension 實例都被存儲在一個 python 字典中。
>>> print rootgrp.dimensions
OrderedDict([("level", <netCDF4._netCDF4.Dimension object at 0x1b48030>),
("time", <netCDF4._netCDF4.Dimension object at 0x1b481c0>),
("lat", <netCDF4._netCDF4.Dimension object at 0x1b480f8>),
("lon", <netCDF4._netCDF4.Dimension object at 0x1b48a08>)])
對 Dimension 實例執行 len 函數可以獲取此維度的大小,Dimension 實例的 isunlimited 方法可以確定維度是不是無限維。
>>> print len(lon)
144
>>> print lon.isunlimited()
False
>>> print time.isunlimited()
True
直接列印 Dimension 對象獲取更多有用信息,比如維度的名稱,大小以及是否是無限維等。
>>> for dimobj in rootgrp.dimensions.values():
>>> print dimobj
<type "netCDF4._netCDF4.Dimension"> (unlimited): name = "level", size = 0
<type "netCDF4._netCDF4.Dimension"> (unlimited): name = "time", size = 0
<type "netCDF4._netCDF4.Dimension">: name = "lat", size = 73
<type "netCDF4._netCDF4.Dimension">: name = "lon", size = 144
<type "netCDF4._netCDF4.Dimension"> (unlimited): name = "time", size = 0
通過 Dataset 或 Group 實例的 netCDF4.Datatset.renameDimension 方法可以改變維度的名稱。
netcdf 中的變量
netcdf 中的變量就像 numpy 模塊中的 python 多維數組。然而,不像 numpy 數組,可以在一個或多個無限維添加netcdf 變量。
使用 Dataset 或 Group 實例的 createVariable 方法可以創建一個 netcdf 變量。createVariable 方法需要至少傳遞兩個參數(變量名,變量數據類型)。通過包含變量名的元組確定變量的維度。如果要創建標量變量,只需要忽略維度關鍵詞即可。
變量的數據類型和 numpy 數據的類型是一致的。你可以指定數據類型為 numpy dtype 對象,或是可以轉換為numpy dtype 對象的任何東西。
有效的數據類型包括:
f4 : 32位浮點數
f8 : 64位浮點數
i/u4 : 32位有/無符號整型
i/u2 : 16位有/無符號整型
i/u8 : 64位有/無符號整型
i/u1 : 8位有/無符號整型
舊的單字符代碼 (f, d, h, s, b, B, c, i, l) 分別表示 (f4, f8, i2, i2, i1, i1, S1, i4, i4),這種設置方式仍然有效。如果文件格式是 NETCDF4, 無符號整型和64位整型可以使用。
維度本身也可以被定義為變量,稱為 坐標變量。
createVariable 方法會返回一個 Variable 類實例,可以通過其方法獲取和設置變量數據及屬性。
>>> times = rootgrp.createVariable("time","f8",("time",))
>>> levels = rootgrp.createVariable("level","i4",("level",))
>>> latitudes = rootgrp.createVariable("lat","f4",("lat",))
>>> longitudes = rootgrp.createVariable("lon","f4",("lon",))
>>>
>>> temp = rootgrp.createVariable("temp","f4",("time","level","lat","lon",))
在交互時可以直接列印變量的所有信息:
>>> print temp
<type "netCDF4._netCDF4.Variable">
float32 temp(time, level, lat, lon)
least_significant_digit: 3
units: K
unlimited dimensions: time, level
current shape = (0, 0, 73, 144)
當然也可以直接在組的層級結構中通過路徑的方式創建變量:
>>> ftemp = rootgrp.createVariable("/forecasts/model1/temp","f4",("time","level","lat","lon",))
如果中間某個組不存在的話,會直接創建,而不會引起錯誤。
你可以通過指定 Group 或 Variable 實例的路徑來查詢 Group 或 Dataset 實例。
>>> print rootgrp["/forecasts/model1"]
<type "netCDF4._netCDF4.Group">
group /forecasts/model1:
dimensions(sizes):
variables(dimensions): float32 temp(time,level,lat,lon)
groups:
>>> print rootgrp["/forecasts/model1/temp"]
<type "netCDF4._netCDF4.Variable">
float32 temp(time, level, lat, lon)
path = /forecasts/model1
unlimited dimensions: time, level
current shape = (0, 0, 73, 144)
filling on, default _FillValue of 9.96920996839e+36 used
Dataset 或 Group 都會存儲在 python 字典中,存儲方式和 維度 相同。
>>> print rootgrp.variables
OrderedDict([("time", <netCDF4.Variable object at 0x1b4ba70>),
("level", <netCDF4.Variable object at 0x1b4bab0>),
("lat", <netCDF4.Variable object at 0x1b4baf0>),
("lon", <netCDF4.Variable object at 0x1b4bb30>),
("temp", <netCDF4.Variable object at 0x1b4bb70>)])
變量名可以通過 Dataset 實例的 renameVariable 方法更改。
netcdf 文件中的屬性
netcdf 文件中包含了兩種類型的屬性:全局屬性和變量屬性。前者提供的是組或整個數據集的信息,後者提供的是組中變量的信息。
指定值給Dataset 或 Group 實例變量可以設置 全局屬性,而位 Variable 實例變量賦值可以設置變量屬性。屬性可以是字符串,數字或序列。
>>> import time
>>> rootgrp.description = "bogus example script"
>>> rootgrp.history = "Created " + time.ctime(time.time())
>>> rootgrp.source = "netCDF4 python module tutorial"
>>> latitudes.units = "degrees north"
>>> longitudes.units = "degrees east"
>>> levels.units = "hPa"
>>> temp.units = "K"
>>> times.units = "hours since 0001-01-01 00:00:00.0"
>>> times.calendar = "gregorian"
Dateset,Group,Variable的 ncattrs 方法可以獲取所有 netcdf 屬性的名稱。使用python 內置的 dir 函數可以返回一些列私有方法和屬性(用戶不能或不應該更改)。
>>> for name in rootgrp.ncattrs():
>>> print "Global attr", name, "=", getattr(rootgrp,name)
Global attr description = bogus example script
Global attr history = Created Mon Nov 7 10.30:56 2005
Global attr source = netCDF4 python module tutorial
Dataset,Group,Variable 實例的 __dict__ 屬性將所有的 netcdf 屬性名/值對存儲在python字典中。
>>> print rootgrp.__dict__
OrderedDict([(u"description", u"bogus example script"),
(u"history", u"Created Thu Mar 3 19:30:33 2011"),
(u"source", u"netCDF4 python module tutorial")])
使用 del 命令可以刪除 Dataset,Group,Variable 的屬性。
寫或讀取netcdf變量數據
現在創建了Variable 實例,那麼如何寫入數據呢?你可以將其視為一個數組,然後傳遞數據給一個切片即可。
>>> import numpy
>>> lats = numpy.arange(-90,91,2.5)
>>> lons = numpy.arange(-180,180,2.5)
>>> latitudes[:] = lats
>>> longitudes[:] = lons
與 numpy 數組對象不同的是,如果你在當前定義的索引範圍外賦值,那麼無限維變量將會沿著這些維度自動擴展。
>>>
>>> nlats = len(rootgrp.dimensions["lat"])
>>> nlons = len(rootgrp.dimensions["lon"])
>>> print "temp shape before adding data = ",temp.shape
temp shape before adding data = (0, 0, 73, 144)
>>>
>>> from numpy.random import uniform
>>> temp[0:5,0:10,:,:] = uniform(size=(5,10,nlats,nlons))
>>> print "temp shape after adding data = ",temp.shape
temp shape after adding data = (6, 10, 73, 144)
>>>
>>>
>>> print "levels shape after adding pressure data = ",levels.shape
levels shape after adding pressure data = (10,)
注意:當沿著 temp 變量的 level 維度增加數據時,level 變量的大小將會增加。
>>>
>>> levels[:] = [1000.,850.,700.,500.,300.,250.,200.,150.,100.,50.]
然而,numpy 和 netcdf 變量在切片規則上稍有不同。通常是使用 start:stop:step 的形式進行切片。對 netcdf 變量而言,布爾數組和整型序列索引的行為與 numpy 數組是不同的。這些索引在每一個維度是單獨作用的(類似 fortran 中的向量下標法)。這意味著:
>>> temp[0, 0, [0,1,2,3], [0,1,2,3]]
切片 netcdf 變量時返回形狀為 (4, 4) 的數組, 但是對 numpy 數組而言,將返回形狀為 (4, )的數組。類似的,用[0, array([True, False, True]), array([False, True, True, True]), :]檢索一個形為 (2, 3, 4, 5)的netcdf數組,將返回一個 (2, 3, 5) 數組,但對於 numpy 來說,這將引起錯誤,因為這相當於 [0, [0,1], [1,2,3], :]。當使用整數序列切片時,索引不需要排序而且其可能包含重複值。對於 numpy 的 'fancy indexing' 可能引起一些迷惑。通過使對維度數組執行邏輯操作來創建切片,可以提取多維 netcdf 變量數據。
比如:
>>> tempdat = temp[::2, [1,3,6], lats>0, lons>0]
將提取時間維的 0, 2, 4層,壓力層的 1, 3, 6層,北半球和東半球交叉部分。返回數組的大小為 (3, 3, 36, 71)。
>>> print "shape of fancy temp slice = ",tempdat.shape
shape of fancy temp slice = (3, 3, 36, 71)
注意:針對標量變量,為了提取標量變量的數據,可以使用 np.asarray(v) 或 v[...],結果是一個 numpy 標量數組。
處理時間坐標
大部分元數據標準(比如CF)指出:時間的測量應該是使用固定的日曆並且相對於一個固定的日期來測量,其單位應該類似於 YY:MM:DD hh-mm-ss。
如果沒有專門的轉換工具的話,這種單位通常比較難處理。此模塊提供了 num2date,date2num函數來處理。
>>>
>>> from datetime import datetime, timedelta
>>> from netCDF4 import num2date, date2num
>>> dates = [datetime(2001,3,1)+n*timedelta(hours=12) for n in range(temp.shape[0])]
>>> times[:] = date2num(dates,units=times.units,calendar=times.calendar)
>>> print "time values (in units %s): " % times.units+"\n",times[:]
time values (in units hours since January 1, 0001):
[ 17533056. 17533068. 17533080. 17533092. 17533104.]
>>> dates = num2date(times[:],units=times.units,calendar=times.calendar)
>>> print "dates corresponding to time values:\n",dates
dates corresponding to time values:
[2001-03-01 00:00:00 2001-03-01 12:00:00 2001-03-02 00:00:00
2001-03-02 12:00:00 2001-03-03 00:00:00]
num2date 轉換指定了 units 和 calendar 的數值為 datetime 對象,而date2num執行相反的操作。當前定義的calendars均是基於 CF元數據轉換標準。date2index函數返回和一系列 datetime 實例對應的netcdf時間變量的索引。
從多個netcdf數據集中獲取數據
如果你想從多個文件中獲取一個變量的數據,可以使用 MFDataset 類進行數據獲取。相比使用單個文件名創建一個 Dataset 實例,MFDataset 實例可以通過一系列文件名或含有通配符的字符串從多個文件中獲取數據。當然前提是這些文件必須具有相同的無限維數據,然後才能執行切片操作。
目前只有 NETCDF3_64BIT_OFFSET, NETCDF3_64BIT_DATA, NETCDF3_CLASSIC orNETCDF4_CLASSIC 格式文件支持上述操作,NETCDF4格式暫不支持多文件獲取。
>>> from netCDF4 import MFDataset
>>> f = MFDataset("mftest*nc")
>>> print f.variables["x"][:]
注意:只能用來讀取數據,不能用於寫入數據。
有效壓縮 netcdf 變量
存儲在 netcdf4 對象中的數據可以執行壓縮和解壓縮操作。通過createVariable 方法的 zlib,complevel,shuffle關鍵詞參數可以確定執行壓縮操作的參數。
要執行壓縮操作,zlib = True;complevel 表示壓縮速度和效率(1表示速度最快,但壓縮比最小,9表示速度最慢,但壓縮比最高);設置 shuffle = False,將關閉 HDF5 的 shuffle 過濾,在重新排序字節進行壓縮之前會對數據塊進行逐行掃描。shuffle 過濾器能顯著改善壓縮比,默認打開此選項。
設置 createVariable 方法的 fletcher32 關鍵詞參數為 True(默認為 False),將利用 Fletcher32 校驗算法進行錯誤檢查,而關鍵詞參數 chunksizes 和 endian 可以用來設置存儲在 HDF5 文件中的 HDF 塊(chunking)參數和二進位字節序(endian-ness)。這些關鍵詞參數僅和 NETCDF4 和 NETCDF4_CLASSIC 文件相關,其他格式文件(NETCDF3_CLASSIC, NETCDF3_64BIT_OFFSET, NETCDF3_64BIT_DATA)會自動忽略這些關鍵詞參數。
如果你的數據有固定精度的話(比如,溫度數據,通常其精度為 0.1 度),可以通過 creatVariable 函數的 least_significant_digit 關鍵詞參數進行數據量化,從而有效改善zlib壓縮。
比如:如果數據的精度為 0.1,設置 least_significant_digit = 1將會使數據被量化為 numpy.around(scale*data)/scale, scale = 2**bits,此例中 bits = 4,這會導致壓縮有損而不是無損,當然這是為了磁碟空間而必須做出的犧牲。
此例中,可以分別執行以下語句看一下文件的大小變化:
>>> temp = rootgrp.createVariable("temp","f4",("time","level","lat","lon",))
>>> temp = dataset.createVariable("temp","f4",("time","level","lat","lon",),zlib=True)
>>> temp = dataset.createVariable("temp","f4",("time","level","lat","lon",),zlib=True,least_significant_digit=3)
如果想要了解更多關於之前提到的三種數據類型可查閱官方文檔。