使用 python 處理 nc 數據

前言

這兩天幫一個朋友處理了些 nc 數據，本覺得很簡單的事情，沒想到裏面涉及到了不少的細節和坑，不管是「知難行易」仍是「知易行難」都不能充分的說明問題，仍是「知行合一」來的更靠譜些，既要知道理論又要知道如何實現，因而通過不太充分的研究後總結成此文，以記錄如何使用 python 處理 nc 數據。

1、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 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)

每個 SUBDATASET 表示記錄的是一種格式的數據（氣象、水文等等），若是要想查看此 SUBDATASET 的具體信息，能夠執行：

gdalinfo NETCDF:name.nc:SUBDATASET_NAME

此處的 SUBDATASET_NAME 爲上面的 T二、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

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

因此看到這裏，各位應該已經明白了，能夠直接使用 GDAL 處理 nc 數據，好比直接使用 gdalwarp 將某個 SUBDATASET 轉成 GeoTiff 等等，此處暫且不表，各位只須要查閱一下 gdalwarp 手冊便可知道如何處理。

明白了以上信息基本也就清楚瞭如何處理此數據。

2、數據處理

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 開始。

3、總結

本文簡單介紹了 nc 數據的特色及如何使用 python 處理 nc 數據。每一個目標都有多條路能夠達到，重要的是找到那條本身喜歡的和適合本身的路，然而話又說回來，即便走的不是想要的那條路，不是同樣能夠達到目標嘛！因此關鍵是要找到本身的目標。