【xarray库(一) 】创建xarray对象

自学气象人

发布于 2022-11-02 10:14:05

5.1K0

发布于 2022-11-02 10:14:05

文章被收录于专栏：自学气象人

以下文章来源于深雨露，作者嵊雨饧

封图：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类

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。
部分数组数据
坐标列表，每一行的项都是坐标的一项。每一项包含坐标名称，括号代表与之关联的维度（可有多个），坐标值数据类型，部分坐标值。另外，如果坐标名称与维度名称重名，则将标有*。
按字母顺序排序的维度名称列表（不含坐标名称）。
无序属性（attributes）的列表。

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

为坐标参数添加多个属性

练习1

现已通过随机函数常见了一个名为height伪数据

不考虑创建的其他要素，试着以height为基础数据创建一个DataArray对象。下面提供了部分代码。
添加维度名字“x”和“y”。
添加数组名称“My random array”
分别对“x”和“y”维度添加坐标“longtitude”（经度）和“latitude”（纬度）。；经度从-180至180，步长为1；纬度从-90至90，步长为1。

“提示：均分函数np.linspace(起点, 终点, 分隔数) ”

height = rng.random((180, 360)) * 400
# 在这里写你的代码

xr.Dataset类

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分别建立了两个坐标

x_a

和

x_b

，那么

与

x_a

的和

与

x_b

的对应关系可为下表所示

（坐标）			1	2	3
（变量1）	/	/	0	0.5	1
（坐标）	-1	0	1	2
（变量2）	0	0	0	0	/

x_a

（坐标） 123

（变量1）//00.51

x_b

（坐标）-1012

（变量2）0000/

在这种情况下，综合了两个变量的坐标，并且将缺失的数据补足了nan值（已经转换为float类型）。例如b没有一个对应x==3的值，所以就用nan代替了。

练习2

利用两个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 ，步长为1

xr.Dataset(
    data_vars={'height':height_da,
                'gravity_anomaly':gravity_anomaly_da},
    coords={  # 在这里写你的代码
            }
)

添加属性到坐标（coordinates）和变量（variables）：

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),# 在这里写你的代码)
            }
)

参考答案

练习1参考答案

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

练习2参考答案

# 第一问
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

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2021-11-22，如有侵权请联系 cloudcommunity@tencent.com 删除

html

编程算法

python

本文分享自自学气象人微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

html

编程算法

python

登录后参与评论

0 条评论

热度