Python气象数据处理进阶之Xarray(1)：Xarray的
首先声明的是前一系列并没有完结，只是暂时新开一个新系列。关于画图部分目前也没有什么太好的想法，所以暂时搁置，等有想法了会继续写下去的。考虑到目前大家最感兴趣的还是数据读取和处理部分，并且这部分也是python处理数据的优势，因此就专门开一个系列来写这部分内容。这也是对“Python气象数据处理与绘图(1)：数据读取”的深入。按目前的构想，这部分应该是以几个库使用为主。首先是数据读取部分，我使用的是Xarray，因此这个系列也从Xarray的介绍开始。主要内容还是围绕Xarray的官网文档进行。
1、Xarray中的数据结构是怎样的？
在Xarray中，数组是具有结构和标签的，它们分为以下几种：
1.DataArray:
带有标注或命名维度的多维数组。DataArray将metadata(例如维名称，坐标和属性)添加到基础的“未标记”的数据结构，例如numpy和Dask数组。
2.Dataset:
具有类似字典结构的尺寸对齐的DataArray对象的集合。因此，可以在单个DataArray的维度上执行的大多数操作都可以在数据集上执行。
3.Variable:
类似于NetCDF的变量，由dimensions, data, 和 attributes组成。变量和numpy数组之间的主要功能区别在于，对变量的数字运算可以通过维名称实现数组广播。
通俗的讲Variable< DataArray< Dataset （<指包含于的意思），这样解释并不完全准确，但是对于初学者来说这样理解是没有问题的。
Xarray中的数据结构是怎样被识别和标记的呢？
前边提到，可以通过对维命名的操作实现数据筛选和处理，实现数据的标记和命名是通过以下几个定义实现的：
1.Dimension:
维，维度的维，比如说在一个二维直角坐标系，维就是x和y，在一个二维圆柱投影地理坐标系，维就是lat和lon，我们通常下载的数据最多也就是四维，即时间维，高度维，纬度维和经度维。
2.Coordinate：
坐标或者说刻度。还是举例子解释，比如说Dimension是纬度维，那么对应的Coordinate就是纬度坐标(90°N，89°N，88°N........89°S，90°S)
3.Index:
索引号，也可以说位置标号。a[0]就代表a数组的第一个数，0就是index
这次不举一个例子，举一些例子。
#CN05.1格点资料
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...
我们可以看到，这个文件，是一个dataset，里边含有变量：Data variables，数据集的维度有经度纬度和时间，各自有各自的坐标Coordinates，同样数据集还有一些属性Attributes来表明数据集信息。
我们可以通过:
print(f.variables)
print(f.dims)
print(f.coords)
来分别查看数据集中所含有的变量，维，坐标。
再比如NCEP的位势高度资料
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
基本同上。
需要说明一点的是，ncl中数据存在short格式，在读取时需要使用short2flt()函数，但是在python中是不存在short格式的，默认均为float，无需考虑这点。
如何创建一个DataArray?
有时候我们可能通过其他手段读取了相关数据，但是数据是np.array格式的，我们需要将其转换为DataArray，亦或者我们需要输出一个NC文件，需要将计算后的数组转为DataArray格式，这就用到了创建的方法。
创建一个DataArray需要什么？
1.Data
废话，创建数据当然需要数据啦。数据可以是numpy ndarray，series，DataFrame等等格式的
2.coords：坐标列表或字典
如果是列表，则应为元组列表，其中第一个元素是维名称，第二个元素是对应的坐标array_like对象。个人建议用字典结构去定义
3.dims：维名称列表
如果省略，并且coords是元组列表，则维度名称取自coords。
4.attrs：属性
5.name:变量名
以上除了data以外，都不是必须的。
创建如下：
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'
例子基本于官网一致，对标签略作了修改，level更加便于气象专业的理解。
刚刚说了，除了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
如果我们是从一个DataFrame数据转化为DataArray的话(这种操作通常是为了将Pandas和Xarray联合使用)：
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'
会自动识别行列的名称和序号。
官方文档还给出了更复杂的实例，但是个人感觉很少用到，因此不写在这里了，通常用以上的例子可以完成大部分操作了。
在创建了数据之后，我们同样可以使用相关的操作获取DataArray的各种信息：
a = foo.values
a = foo.dims
a = foo.coords
a = foo.attrs
如果想对DataArray的值进行修改可以通过以下两种方法：
foo.values = foo.values +2
foo = foo +2
目前测试结果两者是等价的，但是我不知道官方为什么没给出第二种方法。
通过指令
foo.attrs['units'] = 'meters'
赋予属性信息，比如说给一个单位，备注等等。
通过指令
foo.name = 'hgt'
赋予名称信息
通过指令：
foo.rename('temperature')
改名，比如说通过hgt计算得到了一个新变量，需要改名，就可以用这个指令了。
在得到了一个DataArray之后，用于画图时，比如说我们需要获取他的经度和纬度(在这里，刚刚的例子是时间和高度)，
那么可以直接通过
foo.coords['time']
foo['time']
这两种方式取出坐标信息。
要修改或者删除某坐标信息的话，原理和修改数据是一样的：
foo['time'] = pd.date_range('1999-01-02', periods=2)
del foo['time']
如何创建一个Dataset?
官网给出了一个以气候数据为例的Dataset结构：
Dataset结构
一个数据集，包含了数据主体(Temperature, Precipitation)，维度坐标(latitude，longitude)
根据官网的例子，一个Dataset是这样创建的，其实与DataArray类似：
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
实际上这个例子与我们通常接触的不太一样，因为大部分的数据的lat和lon是一个一维的。
对Dataset的操作与DataArray基本是一致的，就不重复了。
