使用 python 處理 nc 資料

標籤：floating   res   資料   nod   update   lower   可視化   就是   總結   

前言

這兩天幫一個朋友處理了些 nc 資料，本以為很簡單的事情，沒想到裡面涉及到了很多的細節和坑，無論是“知難行易”還是“知易行難”都不能充分的說明問題，還是“知行合一”來的更靠譜些，既要知道理論又要知道如何?，於是經過不太充分的研究後總結成此文，以記錄如何使用 python 處理 nc 資料。

一、nc 資料介紹

nc 全稱 netCDF（The Network Common Data Form），可以用來儲存一系列的數組，就是這麼簡單（參考https://www.unidata.ucar.edu/software/netcdf/docs/netcdf_introduction.html）。

既然 nc 可以用來一系列的數組，所以經常被用來儲存科學觀測資料，最好還是長時間序列的。

試想一下一個科學家每隔一分鐘採集一次實驗資料並儲存了下來，如果不用這種格式儲存，時間長了可能就需要建立一系列的 csv 或者 txt 等，而採用 nc 一個檔案就可以搞定，是不是很方便。

更方便的是如果這個科學實驗與氣象、水文、溫度等地理資訊稍微沾點邊的，完全也可以用 nc 進行儲存， GeoTiff 頂多能多存幾個波段（此處波段可以認為是氣象、水文等不同訊號），而 nc 可以儲存不同波段的長時間觀測結果，是不是非常方便。

可以使用 gdal 查看資料資訊，執行：

gdalinfo name.nc

即可得到如下資訊：

Driver: netCDF/Network Common Data FormatFiles: test.ncSize is 512, 512Coordinate 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)

每一個 SUBDATASET 表示記錄的是一種格式的資料（氣象、水文等等），如果要想查看此 SUBDATASET 的具體資訊，可以執行：

gdalinfo NETCDF:name.nc:SUBDATASET_NAME

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

Driver: netCDF/Network Common Data FormatFiles: test.ncSize is 120, 130Coordinate System is `'Metadata:  LAT#description=LATITUDE, SOUTH IS NEGATIVE  LAT#FieldType=104  LAT#MemoryOrder=XY  LAT#stagger=  LAT#units=degree_northCorner 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

此處只有一個 Band ，每一個 Band 記錄了一個時間點（或者其他區分形式）的一條記錄，這個記錄是一個數組。

所以看到這裡，各位應該已經明白了，可以直接使用 GDAL 處理 nc 資料，比如直接使用 gdalwarp 將某個 SUBDATASET 轉成 GeoTiff 等等，此處暫且不表，各位只需要查閱一下 gdalwarp 手冊即可知道如何處理。

明白了以上資訊基本也就清楚了如何處理此資料。

二、資料處理

python 是運用非常廣泛，自然其下各種類庫非常豐富，專業一點的說法就叫生態豐富。

2.1 netCDF4

此架構可以直接將 nc 讀取成數組（詳細資料參考https://github.com/Unidata/netcdf4-python）。讀取方式如下：

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

這樣即可讀出整個 nc 中的資料資訊，如果需要擷取某個 SUBDATASET 只需要使用 dataset[SUBDATASET_NAME] 即可，返回的是一個三維數組，表示不同時間段（或其他區分方式下）的資料資訊。

我們可以對此數組做各種操作，如求平均值、方差等等，又讓我想起了大學裡的那一堆枯燥但又讓人很有興趣的實驗課程。當然，此處如果使用 numpy 架構進行處理，會起到事半功倍的效果，如求長時間序列下的平均值：

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

到這裡跟地信有關的同志都會看出一個問題，此架構只能對資料進行處理，而不能進行與位置有關的操作，這就導致資料無法變成直白的地圖視覺效果。其實任何資料都是相通的，我們可以採用此種方式處理完後轉為 GeoTiff 等，當然我們也可以直接採用 GeoTiff 的處理流程來進行處理。

2.2 rasterio

rasterio 是 Mapbox 開源的空間資料處理架構，功能非常強大，此處不細說，只表如何處理我們的 nc 資料。

當然第一種方式就是使用 netCDF4 處理完之後，使用此架構寫入 GeoTiff，但是這樣不太優雅，而且使用了兩個架構，明顯過於麻煩，我們直接使用此架構從讀資料開始處理。

此處讀的時候就有技巧了，要像採用 gdalinfo 讀取 SUBDATASET 一樣來直接讀取此 SUBDATASET 資料，如下：

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 函數傳入 NETCDF:name.nc:SUBDATASET_NAME，採用 src.read(range(1, dim + 1)) 可以直接讀出此範圍內所有 Band （時間點）的資訊，範圍可以自己設定，注意從 0 開始，當然也可以僅讀取某個 Band 的資訊。

src.meta 記錄了此 SUBDATASET 的中繼資料資訊，與 gdalinfo 看到的基本相同。

這樣我們就可以繼續將此資料使用 numpy 等架構進行處理，處理完之後更重要的是要寫入 GeoTiff 中（直白的說就是增加空間資訊）。

也很簡單，如下即可：

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

newfile 為儲存路徑，res_arr 為計算結果數組，注意尺寸不要發生變化（width*height），out_meta 為目標檔案的中繼資料描述資訊，可以直接將上面 src.meta 進行簡單處理即可。

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 表示目標資料空間投影資訊，transform 表示目標檔案 空間範圍資訊，可以通過經緯度資訊和映像尺寸等計算得到。

dst.write_band 將資料寫入對應波段，當然此處也可以寫入多個波段，根據計算結果而定，同樣從 1 開始。

三、總結

本文簡單介紹了 nc 資料的特點及如何使用 python 處理 nc 資料。每個目標都有多條路可以達到，重要的是找到那條自己喜歡的和適合自己的路，然而話又說回來，即使走的不是想要的那條路，不是一樣可以達到目標嘛！所以關鍵是要找到自己的目標。

使用 python 處理 nc 資料