您好,登錄后才能下訂單哦!
這篇文章給大家分享的是有關(guān)利用python怎么處理nc數(shù)據(jù)的內(nèi)容。小編覺(jué)得挺實(shí)用的,因此分享給大家做個(gè)參考,一起跟隨小編過(guò)來(lái)看看吧。
如何使用 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 處理完之后,使用此框架寫入 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)行處理,處理完之后更重要的是要寫入 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ù)寫入對(duì)應(yīng)波段,當(dāng)然此處也可以寫入多個(gè)波段,根據(jù)計(jì)算結(jié)果而定,同樣從 1 開(kāi)始。
感謝各位的閱讀!關(guān)于“利用python怎么處理nc數(shù)據(jù)”這篇文章就分享到這里了,希望以上內(nèi)容可以對(duì)大家有一定的幫助,讓大家可以學(xué)到更多知識(shí),如果覺(jué)得文章不錯(cuò),可以把它分享出去讓更多的人看到吧!
免責(zé)聲明:本站發(fā)布的內(nèi)容(圖片、視頻和文字)以原創(chuàng)、轉(zhuǎn)載和分享為主,文章觀點(diǎn)不代表本網(wǎng)站立場(chǎng),如果涉及侵權(quán)請(qǐng)聯(lián)系站長(zhǎng)郵箱:is@yisu.com進(jìn)行舉報(bào),并提供相關(guān)證據(jù),一經(jīng)查實(shí),將立刻刪除涉嫌侵權(quán)內(nèi)容。