寻找 xarray 的质心

Question

我有一个 xarray 表示布尔值（例如，forest/no 森林）地理空间数据，其维度 x、y 表示纬度和经度，我想要质心坐标。

import xarray as xr
import numpy as np
A = xr.DataArray(np.array([[0,0,1,1,0,0],
                           [0,1,1,1,0,0],
                           [0,1,1,1,0,0],
                           [0,1,1,0,0,0],
                           [0,0,0,0,1,0]]),
                 dims=['y','x'],
                 coords={'x': [10,20,30,40,50,60],
                         'y': [50,40,30,20,10]})

我想到了以下解决方案：

centr_x = float(np.sum(A.sum('y')/np.sum(A.sum('y')) * A.x))
centr_y = float(np.sum(A.sum('x')/np.sum(A.sum('x')) * A.y))

只是想知道我是否在 xarray 中遗漏了一个专门做这个的函数？在我看来，这将是一个相当普遍的计算。

感谢您的建议！

Answer 1

xarray 中没有任何内容可以将栅格数据解释为几何或点集合，或计算此类要素的质心。

您的方法很聪明 - 您当然可以使用加权平均来获得简单的质心，例如：

In [8]: A.x.weighted(A).mean()
Out[8]:
<xarray.DataArray 'x' ()>
array(31.81818182)

In [9]: A.y.weighted(A).mean()
Out[9]:
<xarray.DataArray 'y' ()>
array(32.72727273)

您也可以使用 geometries/points 使用 shapely 来做到这一点。为此，您可以将栅格数据转换为一系列点：

In [4]: points = A.where(A).to_series().dropna()

In [5]: points
Out[5]:
y   x
50  30    1.0
    40    1.0
40  20    1.0
    30    1.0
    40    1.0
30  20    1.0
    30    1.0
    40    1.0
20  20    1.0
    30    1.0
10  50    1.0
dtype: float64

In [6]: points = shapely.geometry.MultiPoint(list(zip(
   ...:     points.index.get_level_values('x'),
   ...:     points.index.get_level_values('y'),
   ...: )))

In [7]: points
Out[7]: <shapely.geometry.multipoint.MultiPoint at 0x1050b76a0>

然后你可以使用形状工具计算质心：

In [8]: points.centroid
Out[8]: <shapely.geometry.point.Point at 0x1533fada0>

In [9]: points.centroid.xy
Out[9]: (array('d', [31.818181818181817]), array('d', [32.72727272727273]))

请注意，如果您的数据对应于地理空间坐标，您可能需要在等面积投影内计算质心。类似的 xarray-based 方法需要将 lat/lon 坐标转换为 equal-area 投影：


In [10]: xx, yy = np.meshgrid(A.x, A.y)

In [11]: points = gpd.points_from_xy(
    ...:     xx.ravel(), yy.ravel(), crs='epsg:4326'
    ...: )

In [12]: equal_area_xy = points.to_crs('+proj=cea')

In [13]: A.coords["equal_area_y"] = (
    ...:     ("y", "x"), equal_area_xy.y.reshape(A.shape)
    ...: )

In [14]: A.coords["equal_area_x"] = (
    ...:     ("y", "x"), equal_area_xy.x.reshape(A.shape)
    ...: )

In [15]: A
Out[15]:
<xarray.DataArray (y: 5, x: 6)>
array([[0, 0, 1, 1, 0, 0],
       [0, 1, 1, 1, 0, 0],
       [0, 1, 1, 1, 0, 0],
       [0, 1, 1, 0, 0, 0],
       [0, 0, 0, 0, 1, 0]])
Coordinates:
  * x             (x) int64 10 20 30 40 50 60
  * y             (y) int64 50 40 30 20 10
    equal_area_y  (y, x) float64 4.866e+06 4.866e+06 ... 1.1e+06
    equal_area_x  (y, x) float64 1.113e+06 2.226e+06 ... 6.679e+06

现在可以使用同样的加权平均法计算质心：


In [16]: equal_area_centroid_y = A.equal_area_y.weighted(A).mean()

In [17]: equal_area_centroid_x = A.equal_area_x.weighted(A).mean()

In [18]: centroid = gpd.points_from_xy(
    ...:     [equal_area_centroid_x],
    ...:     [equal_area_centroid_y],
    ...:     crs="+proj=cea",
    ...: ).to_crs("epsg:4326")

In [19]: centroid.x, centroid.y
Out[19]: (array([31.81818182]), array([31.94714021]))

请注意，在我们通过 equal-area 投影 round-tripped 的版本中，纬度明显较低 - 这说明经线之间的距离在较高纬度处较小，因此质心将正确考虑此曲率时会稍微低一些。这当然仍将地球的形状近似为球体 - 更精确的质心将需要更复杂的工作流程。

寻找 xarray 的质心

Finding the centroid of an xarray

python

python-xarray