首先聲明的是前一系列并沒(méi)有完結(jié)，只是暫時(shí)新開(kāi)一個(gè)新系列。關(guān)于畫(huà)圖部分目前也沒(méi)有什么太好的想法，所以暫時(shí)擱置，等有想法了會(huì)繼續(xù)寫(xiě)下去的?？紤]到目前大家最感興趣的還是數(shù)據(jù)讀取和處理部分，并且這部分也是python處理數(shù)據(jù)的優(yōu)勢(shì)，因此就專(zhuān)門(mén)開(kāi)一個(gè)系列來(lái)寫(xiě)這部分內(nèi)容。這也是對(duì)“Python氣象數(shù)據(jù)處理與繪圖(1)：數(shù)據(jù)讀取”的深入。按目前的構(gòu)想，這部分應(yīng)該是以幾個(gè)庫(kù)使用為主。首先是數(shù)據(jù)讀取部分，我使用的是Xarray，因此這個(gè)系列也從Xarray的介紹開(kāi)始。主要內(nèi)容還是圍繞Xarray的官網(wǎng)文檔進(jìn)行。
1、Xarray中的數(shù)據(jù)結(jié)構(gòu)是怎樣的？
在Xarray中，數(shù)組是具有結(jié)構(gòu)和標(biāo)簽的，它們分為以下幾種：
1.DataArray:
帶有標(biāo)注或命名維度的多維數(shù)組。DataArray將metadata(例如維名稱(chēng)，坐標(biāo)和屬性)添加到基礎(chǔ)的“未標(biāo)記”的數(shù)據(jù)結(jié)構(gòu)，例如numpy和Dask數(shù)組。
2.Dataset:
具有類(lèi)似字典結(jié)構(gòu)的尺寸對(duì)齊的DataArray對(duì)象的集合。因此，可以在單個(gè)DataArray的維度上執(zhí)行的大多數(shù)操作都可以在數(shù)據(jù)集上執(zhí)行。
3.Variable:
類(lèi)似于NetCDF的變量，由dimensions, data, 和 attributes組成。變量和numpy數(shù)組之間的主要功能區(qū)別在于，對(duì)變量的數(shù)字運(yùn)算可以通過(guò)維名稱(chēng)實(shí)現(xiàn)數(shù)組廣播。
通俗的講Variable< DataArray< Dataset （<指包含于的意思），這樣解釋并不完全準(zhǔn)確，但是對(duì)于初學(xué)者來(lái)說(shuō)這樣理解是沒(méi)有問(wèn)題的。
Xarray中的數(shù)據(jù)結(jié)構(gòu)是怎樣被識(shí)別和標(biāo)記的呢？
前邊提到，可以通過(guò)對(duì)維命名的操作實(shí)現(xiàn)數(shù)據(jù)篩選和處理，實(shí)現(xiàn)數(shù)據(jù)的標(biāo)記和命名是通過(guò)以下幾個(gè)定義實(shí)現(xiàn)的：
1.Dimension:
維，維度的維，比如說(shuō)在一個(gè)二維直角坐標(biāo)系，維就是x和y，在一個(gè)二維圓柱投影地理坐標(biāo)系，維就是lat和lon，我們通常下載的數(shù)據(jù)最多也就是四維，即時(shí)間維，高度維，緯度維和經(jīng)度維。
2.Coordinate：
坐標(biāo)或者說(shuō)刻度。還是舉例子解釋?zhuān)热缯f(shuō)Dimension是緯度維，那么對(duì)應(yīng)的Coordinate就是緯度坐標(biāo)(90°N，89°N，88°N........89°S，90°S)
3.Index:
索引號(hào)，也可以說(shuō)位置標(biāo)號(hào)。a[0]就代表a數(shù)組的第一個(gè)數(shù)，0就是index
這次不舉一個(gè)例子，舉一些例子。
#CN05.1格點(diǎn)資料
f = xr.open_dataset('CN05.1_Tmax_1961_2017_daily_05x05.nc')
print(f)
#<xarray.Dataset>
#Dimensions: (latitude: 82, longitude: 142, time: 20574)
#Coordinates:
# * longitude (longitude) float64 69.75 70.25 70.75 71.25 ... 139.2 139.8 140.2
# * latitude (latitude) float64 14.75 15.25 15.75 16.25 ... 54.25 54.75 55.25
# * time (time) datetime64[ns] 1961-01-01 1961-01-02 ... 2017-04-30
#Data variables:
# tmax (time, latitude, longitude) float32 ...
#Attributes:
# CDI: Climate Data Interface version 1.6.5rc3 (http://code.zmaw.d...
# Conventions: CF-1.4
# history: Thu Aug 23 09:34:52 2018: cdo -r remapcon,grid05x05 daily/0...
# CDO: Climate Data Operators version 1.6.5rc3 (http://code.zmaw.d...
我們可以看到，這個(gè)文件，是一個(gè)dataset，里邊含有變量：Data variables，數(shù)據(jù)集的維度有經(jīng)度緯度和時(shí)間，各自有各自的坐標(biāo)Coordinates，同樣數(shù)據(jù)集還有一些屬性Attributes來(lái)表明數(shù)據(jù)集信息。
我們可以通過(guò):
print(f.variables)
print(f.dims)
print(f.coords)
來(lái)分別查看數(shù)據(jù)集中所含有的變量，維，坐標(biāo)。
再比如NCEP的位勢(shì)高度資料
f = xr.open_dataset('hgt.1948.nc')
print(f)
#<xarray.Dataset>
#Dimensions: (lat: 73, level: 17, lon: 144, time: 366)
#Coordinates:
# * level (level) float32 1000.0 925.0 850.0 700.0 ... 50.0 30.0 20.0 10.0
# * lat (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
# * lon (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
# * time (time) datetime64[ns] 1948-01-01 1948-01-02 ... 1948-12-31
#Data variables:
# hgt (time, level, lat, lon) float32 ...
#Attributes:
# Conventions: COARDS
# title: mean daily NMC Reanalysis (1948)
# description: Data is from NMC initialized reanalysis\n(4x/day). It co...
# platform: Model
# history: created 99/05/11 by Hoop (netCDF2.3)
# References: (http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reana...
# dataset_title: NCEP-NCAR Reanalysis 1
基本同上。
需要說(shuō)明一點(diǎn)的是，ncl中數(shù)據(jù)存在short格式，在讀取時(shí)需要使用short2flt()函數(shù)，但是在python中是不存在short格式的，默認(rèn)均為float，無(wú)需考慮這點(diǎn)。
如何創(chuàng)建一個(gè)DataArray?
有時(shí)候我們可能通過(guò)其他手段讀取了相關(guān)數(shù)據(jù)，但是數(shù)據(jù)是np.array格式的，我們需要將其轉(zhuǎn)換為DataArray，亦或者我們需要輸出一個(gè)NC文件，需要將計(jì)算后的數(shù)組轉(zhuǎn)為DataArray格式，這就用到了創(chuàng)建的方法。
創(chuàng)建一個(gè)DataArray需要什么？
1.Data
廢話，創(chuàng)建數(shù)據(jù)當(dāng)然需要數(shù)據(jù)啦。數(shù)據(jù)可以是numpy ndarray，series，DataFrame等等格式的
2.coords：坐標(biāo)列表或字典
如果是列表，則應(yīng)為元組列表，其中第一個(gè)元素是維名稱(chēng)，第二個(gè)元素是對(duì)應(yīng)的坐標(biāo)array_like對(duì)象。個(gè)人建議用字典結(jié)構(gòu)去定義
3.dims：維名稱(chēng)列表
如果省略，并且coords是元組列表，則維度名稱(chēng)取自coords。
4.attrs：屬性
5.name:變量名
以上除了data以外，都不是必須的。
創(chuàng)建如下：
data = np.array([[1,2,3],[4,5,6]])
level = ['500', '850', '1000']
times = pd.date_range('2000-01-01', periods=2)
foo = xr.DataArray(data, coords=[times, level], dims=['time', 'level'])
print(foo)
#<xarray.DataArray (time: 2, level: 3)>
#array([[1, 2, 3],
# [4, 5, 6]])
#Coordinates:
# * time (time) datetime64[ns] 2000-01-01 2000-01-02
# * level (level) <U4 '500' '850' '1000'
例子基本于官網(wǎng)一致，對(duì)標(biāo)簽略作了修改，level更加便于氣象專(zhuān)業(yè)的理解。
剛剛說(shuō)了，除了data以外，其他都是不必要的
foo = xr.DataArray(data)
print(foo)
#<xarray.DataArray (dim_0: 2, dim_1: 3)>
#array([[1, 2, 3],
# [4, 5, 6]])
#Dimensions without coordinates: dim_0, dim_1
如果我們是從一個(gè)DataFrame數(shù)據(jù)轉(zhuǎn)化為DataArray的話(這種操作通常是為了將Pandas和Xarray聯(lián)合使用)：
df = pd.DataFrame({'x': [0, 1], 'y': [2, 3]}, index=['a', 'b'])
df.index.name = 'abc'
df.columns.name = 'xyz'
print(df)
#xyz x y
#abc
#a 0 2
#b 1 3
print(xr.DataArray(df))
#<xarray.DataArray (abc: 2, xyz: 2)>
#array([[0, 2],
# [1, 3]])
#Coordinates:
# * abc (abc) object 'a' 'b'
# * xyz (xyz) object 'x' 'y'
會(huì)自動(dòng)識(shí)別行列的名稱(chēng)和序號(hào)。
官方文檔還給出了更復(fù)雜的實(shí)例，但是個(gè)人感覺(jué)很少用到，因此不寫(xiě)在這里了，通常用以上的例子可以完成大部分操作了。
在創(chuàng)建了數(shù)據(jù)之后，我們同樣可以使用相關(guān)的操作獲取DataArray的各種信息：
a = foo.values
a = foo.dims
a = foo.coords
a = foo.attrs
如果想對(duì)DataArray的值進(jìn)行修改可以通過(guò)以下兩種方法：
foo.values = foo.values +2
foo = foo +2
目前測(cè)試結(jié)果兩者是等價(jià)的，但是我不知道官方為什么沒(méi)給出第二種方法。
通過(guò)指令
foo.attrs['units'] = 'meters'
賦予屬性信息，比如說(shuō)給一個(gè)單位，備注等等。
通過(guò)指令
foo.name = 'hgt'
賦予名稱(chēng)信息
通過(guò)指令：
foo.rename('temperature')
改名，比如說(shuō)通過(guò)hgt計(jì)算得到了一個(gè)新變量，需要改名，就可以用這個(gè)指令了。
在得到了一個(gè)DataArray之后，用于畫(huà)圖時(shí)，比如說(shuō)我們需要獲取他的經(jīng)度和緯度(在這里，剛剛的例子是時(shí)間和高度)，
那么可以直接通過(guò)
foo.coords['time']
foo['time']
這兩種方式取出坐標(biāo)信息。
要修改或者刪除某坐標(biāo)信息的話，原理和修改數(shù)據(jù)是一樣的：
foo['time'] = pd.date_range('1999-01-02', periods=2)
del foo['time']
如何創(chuàng)建一個(gè)Dataset?
官網(wǎng)給出了一個(gè)以氣候數(shù)據(jù)為例的Dataset結(jié)構(gòu)：
一個(gè)數(shù)據(jù)集，包含了數(shù)據(jù)主體(Temperature, Precipitation)，維度坐標(biāo)(latitude，longitude)
根據(jù)官網(wǎng)的例子，一個(gè)Dataset是這樣創(chuàng)建的，其實(shí)與DataArray類(lèi)似：
temp = 15 + 8 * np.random.randn(2, 2, 3)
precip = 10 * np.random.rand(2, 2, 3)
lon = [[-99.83, -99.32], [-99.79, -99.23]]
lat = [[42.25, 42.21], [42.63, 42.59]]
ds = xr.Dataset({'temperature': (['x', 'y', 'time'], temp),
'precipitation': (['x', 'y', 'time'], precip)},
coords={'lon': (['x', 'y'], lon),
'lat': (['x', 'y'], lat),
'time': pd.date_range('2014-09-06', periods=3),
'reference_time': pd.Timestamp('2014-09-05')})
#<xarray.Dataset>
#Dimensions: (time: 3, x: 2, y: 2)
#Coordinates:
# lon (x, y) float64 -99.83 -99.32 -99.79 -99.23
# lat (x, y) float64 42.25 42.21 42.63 42.59
# * time (time) datetime64[ns] 2014-09-06 2014-09-07 2014-09-08
# reference_time datetime64[ns] 2014-09-05
#Dimensions without coordinates: x, y
#Data variables:
# temperature (x, y, time) float64 15.09 7.656 20.82 ... 2.477 10.53 17.56
# precipitation (x, y, time) float64 3.444 2.694 6.921 ... 7.351 2.099 5.972
實(shí)際上這個(gè)例子與我們通常接觸的不太一樣，因?yàn)榇蟛糠值臄?shù)據(jù)的lat和lon是一個(gè)一維的。
對(duì)Dataset的操作與DataArray基本是一致的，就不重復(fù)了。
