以下文章来源于深雨露 ,作者嵊雨饧
封图:Photo by Eiliv-Sonas Aceron on Unsplash
“东风夜放花千树,更吹落,星如雨。 ——辛弃疾(宋)《青玉案·元夕》 ”
在长期的生产实践过程中,人类掌握了四季的规律,懂得何时、何地进行田间管理。近代以来,特别是近百年以来,科学技术的革命快速影响人类理解天气、气候的深度和广度。各类现代化观测站点积累的地球科学数据为理解天气、气候提供了极好的理论依据。
python
语言作为一种高级语言提供了一个与这类地球科学数据提供了一个良好的交互环境基础,而由python语言编写的xarray包[1]则为该类数据的处理提供了良好的平台。
多维数组(Multi-dimensional, N-dimensional, ND Arrays,Tensors)在计算科学、物理学、天文学、地球科学、生物信息学、工程学、金融等领域应用颇为广泛。
如果你之前学习过python语言,那你一定不会对NumPy包[2]陌生。那有个问题便可以提出,为什么不直接使用numpy数组读取地球科学数据,而要使用xarray提供的读取方式?
下面举一个例子
对于地球上的各个地区的温度(Temperature)和湿度(Precipitation)而言。如果你生活在北半球,就平均而言,南方地区常比北方地区更暖、更湿,所以说不同地方的温度、湿度不是一样的;同时地球上的位置通常是用经纬度来表示的,比如北京中心的经纬度为116°20′E, 39°56′N;因而为了描述温度、湿度变量在平面(二维空间)不同位置的变化,我们引入了两个维度(Dimension)进行描述,当然也就可以命名维度名称分别为x
和y
。
那如何将现实生活中的数据存储在计算机中。计算机比较愚蠢,只能类似矩阵一样储存信息。无论是一维数组、二维数组、三维数组下标只能从0开始。
“C语言、python 数组下标从0开始; FORTRAN、MATLAB、R、Julia 数组下标从1开始。 ”
这就造成了一个麻烦的事情。对于实际生活中的经纬度,纬度一般是从90°S至90°N(-90°N~90°N),经度一般是从0°~360°。这两个数组并不吻合。怎么办?应用数学中映射
的思想,将Python中的数组和现实生活中的坐标联系起来。
比如将实际位置(0°,-90°N)即(0°,90°S)映射为Python中的数组(0,0)。其中实际位置的位置(latitude, longitude),我们称为坐标(Coordinate);Python的数组的位置(x, y),我们称为维度(Dimension)。同时定义Python数组的x方向的1个步长等于实际位置坐标90°,y方向的1个步长等于实际位置坐标45°,那么可知(1, 1)<->(90°E, -45°N)
平面映射
还有一个问题?不同时间的温度、降水量也是不同的。一般而言,夏季的降水量、温度都是大于冬季的。那么我们就必须在温度、湿度变量上再引入一个维度——时间t
进行描述。这时候你就可以知道任意时间、任意地点的温度、湿度大小了。与位置的描述有同样的问题,维度t
在Python数组中是从0开始的,不是一个现实生活中的时间。基于上面同样的思想,我们可以定义t=0
时,映射的时间坐标time = 2021-01-01
。假设Python数组1个步长对应时间坐标的一天,那么当时间维度t=1
时,时间坐标为time = 2021-01-02
。
时间映射
所以说维度是针对计算机内数值储存位置,而坐标表示的是实际时空中的位置信息。这两者通过一定的映射关系得以联系到一起。
Xarray包提供了两种数据储存结构:DataArray类
和Dataset类
。 DataArray将维度名称,坐标和属性添加到多维数组,而Dataset则是多个DataArray数组的集合。
“什么叫做类,什么叫做实例? 苹果是水果的组成部分,那么苹果称为水果的类。类表示了一类事物的集合。而针对苹果中的一个苹果,我们称为这是一个苹果(类)中的一个实例。实例是针对一个类中具体的事物。水果摊上的A苹果,B苹果分别都是苹果的实例,它们都属于苹果这一类。 ”
这两个类的创建一般是通过读取数据得到的。为了理解上述概念,首先来看一下如何写代码来创建。
写代码如同造房子一样,python安装后只是完成了地基。为了能让造的房子实现特定的功能,比如住宅、写字楼、商业用房,那么相应的图纸也就不同。这里我们是用来处理地理数据的,那么我们就需要引入处理地理数据的专用包xarray以及配套需要的numpy和pandas包。
import numpy as np
import xarray as xr
import pandas as pd
# 给一个随机数种子,使得每次运行得到的随机数是相同的
rng = np.random.default_rng(seed=0)
xr.DataArray是一个容器,能够装一个变量的信息,比如可以把温度在时空的变化放在里面。
假如温度数据在时空上全为1,那么有
da = xr.DataArray(
# 温度数据
np.ones((3, 4, 2)),
# 维度名称
dims=("x", "y", "t"),
# 数组名称
name="Temperature Data"
)
可以注意到DataArray括号内的三行内容:第一行,利用Numpy包的np.ones
函数创建一个3x4x2 numpy类的全1数组;第二行,我们把三个维度分别命名为x, y, z;第三行,取了一个名字Temperature Data。
假如只有一个维度的情况下可以在写维度名称的时候偷一下懒,不用写括号了
# 一个维度的情形
xr.DataArray([1, 1],
dims="x"
)
运行结果示例
维度名称(dims)以及数组的名称(name)只要满足python set类
的要求即可,但为了查看数据的方面,最好还是命名为字符串的形式(就如同上面的例子写的那样)
只有维度也不行啊,如何将数据与实际时空关联在一起呢?那就要加上坐标信息了
da = xr.DataArray(
# 温度数据
np.ones((3, 4, 2)),
# 维度名称
dims=("x", "y", "t"),
# 数组名称
name="Temperature Data",
# 坐标数据
coords={"longtitude":("x",[0,90,180]),"latitude":("y",[-90,-45,45,90]),"time":("t",['2021-01-01','2021-01-02'])}
)
对坐标数据而言,是用coords={}
大括号框定的区域。"longtitude":("x",[0,90,180])
中的"longtitude"
是坐标名称,"x"
是坐标名称对应的维度,[0,90,180]
是维度(x方向有三个数据,那么对应就有三个坐标值)映射对应坐标数据。
当然有时还需要对数据加上一些说明,比如说明这个数据的来源、获得方法等。这时候便可以加上属性(attributes)
da = xr.DataArray(
# 温度数据
np.ones((3, 4, 2)),
# 维度名称
dims=("x", "y", "t"),
# 数组名称
name="Temperature Data",
# 坐标数据
coords={"longtitude":("x",[0,90,180]),
"latitude":("y",[-90,-45,45,90]),
"time":("t",['2021-01-01','2021-01-02'])
},
# 属性
attrs={"Autor": "Shen yulu"}
)
注意在加上属性的时候,坐标属性项的右大括号要加上英文逗号,不然会报错。
现在我们有了DataArray
类的一个实例da
,如何对其进行查看呢?
xarray对于数据的显示有两种显示形式:
html形式
(仅在Jupyter 笔记本中可用,Jupyter 笔记本中默认以html形式
显示);text形式
。这两种显示形式的选择需指定display_style
选项。首先来看一下text形式
的表示形式:
with xr.set_options(display_style="text"):
display(da)
运行结果
上述输出包含了
DataArray
数组名称(Name):'a'
。如果没有名称,则不会显示。(x:3,y:4,z:2)
:第一个维度叫x
且大小为3
,第二个维度叫y
且大小为4,第三个维度叫z
且大小为2。*
。html形式
与text形式
类似
with xr.set_options(display_style="html"):
display(da)
运行结果
数据列表项(如“Coordinates, Attributes”)左侧的三角符号将详细信息折叠为单行。多个盘状垛堞的图标可以查看对应变量的部分数据。如果坐标名称与维度名称重名,则用粗体
标记维度名称,而非text形式
的*
。默认情况下,若在笔记本中直接查看某个xarray对象,直接写对象名称即可。这时候展示的是html形式
da
运行结果
如果我只想要da
中其中的数据怎么办呢?
da
数据(data)da.data
运行结果
da
维度名称(dimensions)da.dims
提取维度名称
da
坐标信息(coordinates)da.coords
提取坐标信息
da
属性(attributes)da.attrs
提取属性
有时候坐标参数(如纬度)需要附加一些信息,例如可以附加纬度坐标变化的步长(step)为多少,这时候就需要对坐标参数添加属性了
da = xr.DataArray(
# 温度数据
np.ones((3, 4, 2)),
# 维度名称
dims=("x", "y", "t"),
# 数组名称
name="Temperature Data",
# 坐标数据
coords={"longtitude":("x",[0,90,180]),
"latitude":("y",[-90,-45,45,90],{"step": "5 degree"}),
"time":("t",['2021-01-01','2021-01-02'])
},
# 属性
attrs={"Autor": "Shen yulu"}
)
da
为坐标参数添加属性
注意到coords中的"latitude":("y",[-90,-45,45,90],{"step": "5 degree"})
,添加坐标参数属性的方法即在坐标值的后面添加上一个大括号括好的数据。冒号的左侧为属性名称,右侧为对应的值。对于多个属性的添加,需用逗号进行间隔,如{"step": "5 degree","first value":1}
da = xr.DataArray(
# 温度数据
np.ones((3, 4, 2)),
# 维度名称
dims=("x", "y", "t"),
# 数组名称
name="Temperature Data",
# 坐标数据
coords={"longtitude":("x",[0,90,180]),
"latitude":("y",[-90,-45,45,90],{"step": "5 degree","first value":1}),
"time":("t",['2021-01-01','2021-01-02'])
},
# 属性
attrs={"Autor": "Shen yulu"}
)
da
为坐标参数添加多个属性
现已通过随机函数常见了一个名为height
伪数据
height
为基础数据创建一个DataArray
对象。下面提供了部分代码。“提示:均分函数np.linspace(起点, 终点, 分隔数) ”
height = rng.random((180, 360)) * 400
# 在这里写你的代码
Dataset
的对象可以将多个变量放在一起。可以定义每一个都有相应不同维度。
Dataset
由下列三个部分组成
data_vars
:类似于python字典从名称至值的映射关系。对于每一个变量都必须要提供维度名称 和DataArray
对象或元组语法。coords
: 与 DataArray
类似attrs
: 与 DataArray
类似作为例子,下面我们来创建一个有两个变量的Dataset
对象 :
ds = xr.Dataset(
data_vars={
"a": (("x", "y"), np.ones((3, 4))),
"b": ("t", np.full(8, 3), {"b atrri": "b value"}),
},
coords={
"longtitude":("x",[-1, 0, 1]),
"latitude":("y",[-1, 0, 1, 2]),
"time":("t",["2021-01-01", "2021-01-02", "2021-01-03", "2021-01-04", "2021-01-05", "2021-01-06", "2021-01-07", "2021-01-08"])
},
attrs={"Dataset attr": "Num1"}
)
ds
运行结果示例
data_vars
的大括号中包含了a
和b
两个变量,变量名后以小括号扩住这个变量的信息。小括号的信息包含下列信息x
和y
。np.ones
函数构建了一个3×4的全1矩阵,利用np.full
函数构建项数为8全为3的一维数组。:
+属性值字符串。不同属性之间用逗号间隔,类似于{"atrri1": "First", "atrri2": "Second"}
data_vars
末尾的右大括号添加逗号。coords
的大括号包含了具体的坐标参数。与 DataArray
类似。要确保维度名称代表的数据的数目和坐标参数的数目相一致。attrs
的大括号包含了属性参数。与 DataArray
类似。类似于DataArray
,Dataset
也有text形式
输出和html形式
输出。
文本输出情况
with xr.set_options(display_style="text"):
display(ds)
文本输出
html输出情况
with xr.set_options(display_style="html"):
display(ds)
html输出
如果要创建一些在同一个坐标(Coordinates)上但有不同值的变量,我们不能采用简化的语法。相反,我们需要使用到DataArray
对象
x_a = np.arange(1, 4)
x_b = np.arange(-1, 3)
a = xr.DataArray(np.linspace(0, 1, 3), dims="x", coords={"x": x_a})
b = xr.DataArray(np.zeros(4), dims="x", coords={"x": x_b})
xr.Dataset(data_vars={"a": a, "b": b})
运行结果示例
x_a
利用np.arange
函数创建了一个线性序列数组,以1开始,4结束(最终得到的数组不包含本身),步长为1(np.arange
函数未分配参数默认步长为1)。
“x_a = array([1, 2, 3]) ”
x_b
创建了一个线性序列数组,以-1开始,3结束,步长为1。
“x_b = array([-1, 0, 1, 2]) ”
变量a和变量b的维度都在x
上,但对这个维度我们针对变量a和b分别建立了两个坐标
和
,那么
与
的和
与
的对应关系可为下表所示
(坐标) | 1 | 2 | 3 | ||
---|---|---|---|---|---|
(变量1) | / | / | 0 | 0.5 | 1 |
(坐标) | -1 | 0 | 1 | 2 | |
(变量2) | 0 | 0 | 0 | 0 | / |
(坐标) 123
(变量1)//00.51
(坐标)-1012
(变量2)0000/
在这种情况下,综合了两个变量的坐标,并且将缺失的数据补足了nan
值(已经转换为float
类型)。例如b
没有一个对应x==3
的值,所以就用nan
代替了。
DataArray
创建一个具有height
和 gravity_anomaly
两个变量且具有x
和 y
两个维度的Dataset
height = rng.random((360, 180)) * 400
gravity_anomaly = rng.random((360, 180)) * 400 - 200
height_da = xr.DataArray(height,
dims=("x","y"))
gravity_anomaly_da = xr.DataArray(gravity_anomaly,
dims=("x","y"))
# 在这里写你的代码
latitude
和 longitude
两个坐标longitude
: 从 -180 至 180,步长为1(采用np.linspace
函数)latitude
: 从 -90 至 90 ,步长为1xr.Dataset(
data_vars={'height':height_da,
'gravity_anomaly':gravity_anomaly_da},
coords={ # 在这里写你的代码
}
)
latitude
: "type": "geodetic"longitude
: "prime_meridian": "greenwich"height
: "ellipsoid": "wgs84"gravity_anomaly
: "ellipsoid": "grs80"height_da = xr.DataArray(height,
dims=("x","y"),
# 在这里写你的代码)
gravity_anomaly_da = xr.DataArray(gravity_anomaly,
dims=("x","y"),
# 在这里写你的代码)
xr.Dataset(
data_vars={'height':height_da,
'gravity_anomaly':gravity_anomaly_da},
coords={'longtitude':('x',np.linspace(-180,180,360), # 在这里写你的代码),
'latitude':('y',np.linspace(-90,90,180),# 在这里写你的代码)
}
)
height = rng.random((180, 360)) * 400
xr.DataArray(
height,
dims={"x", "y"},
name="My random array",
coords={"longtitude":("x", np.linspace(-180,180,360)),
"latitude":("y", np.linspace(-90,90,180))
}
)
练习1
# 第一问
height = rng.random((360, 180)) * 400
gravity_anomaly = rng.random((360, 180)) * 400 - 200
height_da = xr.DataArray(height,
dims=("x","y"))
gravity_anomaly_da = xr.DataArray(gravity_anomaly,
dims=("x","y"))
xr.Dataset(data_vars={'height':height_da,
'gravity_anomaly':gravity_anomaly_da}
)
练习2第一问
# 第二问
xr.Dataset(
data_vars={'height':height_da,
'gravity_anomaly':gravity_anomaly_da},
coords={'longtitude':('x',np.linspace(-180,180,360)),
'latitude':('y',np.linspace(-90,90,180))
}
)
练习2第二问
# 第三问
height_da = xr.DataArray(height,
dims=("x","y"),
attrs={"ellipsoid": "wgs84"})
gravity_anomaly_da = xr.DataArray(gravity_anomaly,
dims=("x","y"),
attrs={"ellipsoid": "grs80"})
xr.Dataset(
data_vars={'height':height_da,
'gravity_anomaly':gravity_anomaly_da},
coords={'longtitude':('x',np.linspace(-180,180,360),{"prime_meridian": "greenwich"}),
'latitude':('y',np.linspace(-90,90,180),{"type": "geodetic"})
}
)
练习2第三问
本文部分内容来自xarray-contrib/xarray-tutorial[3]
[1]
xarray包: http://xarray.pydata.org/
[2]
NumPy包: https://numpy.org/
[3]
xarray-contrib/xarray-tutorial: https://github.com/xarray-contrib/xarray-tutorial