前言

這兩天幫一個(gè)朋友處理了些 nc 數(shù)據(jù)，本以為很簡(jiǎn)單的事情，沒(méi)想到里面涉及到了很多的細(xì)節(jié)和坑，無(wú)論是“知難行易”還是“知易行難”都不能充分的說(shuō)明問(wèn)題，還是“知行合一”來(lái)的更靠譜些，既要知道理論又要知道如何實(shí)現(xiàn)，于是經(jīng)過(guò)不太充分的研究后總結(jié)成此文，以記錄如何使用 python 處理 nc 數(shù)據(jù)。

一、nc 數(shù)據(jù)介紹

nc 全稱 netCDF（The Network Common Data Form），可以用來(lái)存儲(chǔ)一系列的數(shù)組，就是這么簡(jiǎn)單（參考https://www.unidata.ucar.edu/software/netcdf/docs/netcdf_introduction.html）。

既然 nc 可以用來(lái)一系列的數(shù)組，所以經(jīng)常被用來(lái)存儲(chǔ)科學(xué)觀測(cè)數(shù)據(jù)，最好還是長(zhǎng)時(shí)間序列的。

試想一下一個(gè)科學(xué)家每隔一分鐘采集一次實(shí)驗(yàn)數(shù)據(jù)并存儲(chǔ)了下來(lái)，如果不用這種格式存儲(chǔ)，時(shí)間長(zhǎng)了可能就需要?jiǎng)?chuàng)建一系列的 csv 或者 txt 等，而采用 nc 一個(gè)文件就可以搞定，是不是很方便。

更方便的是如果這個(gè)科學(xué)實(shí)驗(yàn)與氣象、水文、溫度等地理信息稍微沾點(diǎn)邊的，完全也可以用 nc 進(jìn)行存儲(chǔ)， GeoTiff 頂多能多存幾個(gè)波段（此處波段可以認(rèn)為是氣象、水文等不同信號(hào)），而 nc 可以存儲(chǔ)不同波段的長(zhǎng)時(shí)間觀測(cè)結(jié)果，是不是非常方便。

可以使用 gdal 查看數(shù)據(jù)信息，執(zhí)行：

gdalinfo name.nc

即可得到如下信息：

Driver: netCDF/Network Common Data Format
Files: test.nc
Size is 512, 512
Coordinate System is `'
Subdatasets:
 SUBDATASET_1_NAME=NETCDF:"test.nc":T2
 SUBDATASET_1_DESC=[696x130x120] T2 (32-bit floating-point)
 SUBDATASET_2_NAME=NETCDF:"test.nc":PSFC
 SUBDATASET_2_DESC=[696x130x120] PSFC (32-bit floating-point)
 SUBDATASET_3_NAME=NETCDF:"test.nc":Q2
 SUBDATASET_3_DESC=[696x130x120] Q2 (32-bit floating-point)
 SUBDATASET_4_NAME=NETCDF:"test.nc":U10
 SUBDATASET_4_DESC=[696x130x120] U10 (32-bit floating-point)
 SUBDATASET_5_NAME=NETCDF:"test.nc":V10
 SUBDATASET_5_DESC=[696x130x120] V10 (32-bit floating-point)
 SUBDATASET_6_NAME=NETCDF:"test.nc":RAINC
 SUBDATASET_6_DESC=[696x130x120] RAINC (32-bit floating-point)
 SUBDATASET_7_NAME=NETCDF:"test.nc":SWDOWN
 SUBDATASET_7_DESC=[696x130x120] SWDOWN (32-bit floating-point)
 SUBDATASET_8_NAME=NETCDF:"test.nc":GLW
 SUBDATASET_8_DESC=[696x130x120] GLW (32-bit floating-point)
 SUBDATASET_9_NAME=NETCDF:"test.nc":LAT
 SUBDATASET_9_DESC=[130x120] LAT (32-bit floating-point)
 SUBDATASET_10_NAME=NETCDF:"test.nc":LONG
 SUBDATASET_10_DESC=[130x120] LONG (32-bit floating-point)
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0, 512.0)
Upper Right ( 512.0, 0.0)
Lower Right ( 512.0, 512.0)
Center ( 256.0, 256.0)

每一個(gè) SUBDATASET 表示記錄的是一種格式的數(shù)據(jù)（氣象、水文等等），如果要想查看此 SUBDATASET 的具體信息，可以執(zhí)行：

gdalinfo NETCDF:name.nc:SUBDATASET_NAME

此處的 SUBDATASET_NAME 為上面的 T2、PSFC 等等，可以得到如下信息：

Driver: netCDF/Network Common Data Format
Files: test.nc
Size is 120, 130
Coordinate System is `'
Metadata:
 LAT#description=LATITUDE, SOUTH IS NEGATIVE
 LAT#FieldType=104
 LAT#MemoryOrder=XY
 LAT#stagger=
 LAT#units=degree_north
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0, 130.0)
Upper Right ( 120.0, 0.0)
Lower Right ( 120.0, 130.0)
Center ( 60.0, 65.0)
Band 1 Block=120x1 Type=Float32, ColorInterp=Undefined
 NoData Value=9.96920996838686905e+36
 Unit Type: degree_north
 Metadata:
 description=LATITUDE, SOUTH IS NEGATIVE
 FieldType=104
 MemoryOrder=XY
 NETCDF_VARNAME=LAT
 stagger=
 units=degree_north

此處只有一個(gè) Band ，每一個(gè) Band 記錄了一個(gè)時(shí)間點(diǎn)（或者其他區(qū)分形式）的一條記錄，這個(gè)記錄是一個(gè)數(shù)組。

所以看到這里，各位應(yīng)該已經(jīng)明白了，可以直接使用 GDAL 處理 nc 數(shù)據(jù)，比如直接使用 gdalwarp 將某個(gè) SUBDATASET 轉(zhuǎn)成 GeoTiff 等等，此處暫且不表，各位只需要查閱一下 gdalwarp 手冊(cè)即可知道如何處理。

明白了以上信息基本也就清楚了如何處理此數(shù)據(jù)。

二、數(shù)據(jù)處理

python 是運(yùn)用非常廣泛，自然其下各種類庫(kù)非常豐富，專業(yè)一點(diǎn)的說(shuō)法就叫生態(tài)豐富。

2.1 netCDF4

此框架可以直接將 nc 讀取成數(shù)組（詳細(xì)信息參考https://github.com/Unidata/netcdf4-python （本地下載））。讀取方式如下：

dataset = netCDF4.Dataset('name.nc') # open the dataset

這樣即可讀出整個(gè) nc 中的數(shù)據(jù)信息，如果需要獲取某個(gè) SUBDATASET 只需要使用 dataset[SUBDATASET_NAME] 即可，返回的是一個(gè)三維數(shù)組，表示不同時(shí)間段（或其他區(qū)分方式下）的數(shù)據(jù)信息。

我們可以對(duì)此數(shù)組做各種操作，如求平均值、方差等等，又讓我想起了大學(xué)里的那一堆枯燥但又讓人很有興趣的實(shí)驗(yàn)課程。當(dāng)然，此處如果使用 numpy 框架進(jìn)行處理，會(huì)起到事半功倍的效果，如求長(zhǎng)時(shí)間序列下的平均值：

np_arr = np.asarray(dataset[SUBDATASET_NAME])
average_arr = np.average(np_arr, axis=0)

到這里跟地信有關(guān)的同志都會(huì)看出一個(gè)問(wèn)題，此框架只能對(duì)數(shù)據(jù)進(jìn)行處理，而不能進(jìn)行與位置有關(guān)的操作，這就導(dǎo)致數(shù)據(jù)無(wú)法變成直白的地圖可視化效果。其實(shí)任何數(shù)據(jù)都是相通的，我們可以采用此種方式處理完后轉(zhuǎn)為 GeoTiff 等，當(dāng)然我們也可以直接采用 GeoTiff 的處理流程來(lái)進(jìn)行處理。

2.2 rasterio

rasterio 是 Mapbox 開(kāi)源的空間數(shù)據(jù)處理框架，功能非常強(qiáng)大，此處不細(xì)說(shuō)，只表如何處理我們的 nc 數(shù)據(jù)。

當(dāng)然第一種方式就是使用 netCDF4 處理完之后，使用此框架寫(xiě)入 GeoTiff，但是這樣不太優(yōu)雅，而且使用了兩個(gè)框架，明顯過(guò)于麻煩，我們直接使用此框架從讀數(shù)據(jù)開(kāi)始處理。

此處讀的時(shí)候就有技巧了，要像采用 gdalinfo 讀取 SUBDATASET 一樣來(lái)直接讀取此 SUBDATASET 數(shù)據(jù)，如下：

with rio.open('NETCDF:name.nc:SUBDATASET_NAME') as src:
 print(src.meta)
 dim = int(src.meta['count'])
 src.read(range(1, dim + 1))

即給 open 函數(shù)傳入 NETCDF:name.nc:SUBDATASET_NAME，采用 src.read(range(1, dim + 1)) 可以直接讀出此范圍內(nèi)所有 Band （時(shí)間點(diǎn)）的信息，范圍可以自己設(shè)定，注意從 0 開(kāi)始，當(dāng)然也可以僅讀取某個(gè) Band 的信息。

src.meta 記錄了此 SUBDATASET 的元數(shù)據(jù)信息，與 gdalinfo 看到的基本相同。

這樣我們就可以繼續(xù)將此數(shù)據(jù)使用 numpy 等框架進(jìn)行處理，處理完之后更重要的是要寫(xiě)入 GeoTiff 中（直白的說(shuō)就是添加空間信息）。

也很簡(jiǎn)單，如下即可：

with rio.open(newfile, 'w', **out_meta) as dst:
 dst.write_band(1, res_arr)

newfile 為存儲(chǔ)路徑，res_arr 為計(jì)算結(jié)果數(shù)組，注意尺寸不要發(fā)生變化（width*height），out_meta 為目標(biāo)文件的元數(shù)據(jù)描述信息，可以直接將上面 src.meta 進(jìn)行簡(jiǎn)單處理即可。

out_meta = 
 meta.update({"driver": "GTiff",
   "dtype": "float32",
   'count': 1,
   'crs': 'Proj4: +proj=longlat +datum=WGS84 +no_defs',
   'transform': rasterio.transform.from_bounds(west, south, east, north, width, height)
  })

crs 表示目標(biāo)數(shù)據(jù)空間投影信息，transform 表示目標(biāo)文件 空間范圍信息，可以通過(guò)經(jīng)緯度信息和圖像尺寸等計(jì)算得到。

dst.write_band 將數(shù)據(jù)寫(xiě)入對(duì)應(yīng)波段，當(dāng)然此處也可以寫(xiě)入多個(gè)波段，根據(jù)計(jì)算結(jié)果而定，同樣從 1 開(kāi)始。

三、總結(jié)

本文簡(jiǎn)單介紹了 nc 數(shù)據(jù)的特點(diǎn)及如何使用 python 處理 nc 數(shù)據(jù)。每個(gè)目標(biāo)都有多條路可以達(dá)到，重要的是找到那條自己喜歡的和適合自己的路，然而話又說(shuō)回來(lái)，即使走的不是想要的那條路，不是一樣可以達(dá)到目標(biāo)嘛！所以關(guān)鍵是要找到自己的目標(biāo)。

好了，以上就是這篇文章的全部?jī)?nèi)容了，希望本文的內(nèi)容對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值，如果有疑問(wèn)大家可以留言交流，謝謝大家對(duì)腳本之家的支持。

最新評(píng)論