xarray - MissingSpatialDimensionError（将坐标指定为维度）

Question

A 有一个我用 xarray 打开的 netCDF 文件。我想使用 shapefile 剪辑 xarray 数据集；但是，我不知道如何正确设置我的空间维度。

我有以下数据集：

print(ds.keys())

Dimensions:                          (sample: 86401, ddm: 4, delay: 17,
                                      doppler: 11)
Coordinates:
  * sample                           (sample) int32 0 1 2 ... 86398 86399 86400
  * ddm                              (ddm) int8 0 1 2 3
    ddm_timestamp_utc                (sample) datetime64[ns] ...
    sp_lat                           (sample, ddm) float32 ...
    sp_lon                           (sample, ddm) float32 ...
Dimensions without coordinates: delay, doppler
Data variables: (12/126)
    spacecraft_id                    int16 ...
    spacecraft_num                   int8 ...
    ddm_source                       int8 ...
    ddm_time_type_selector           int8 ...
    delay_resolution                 float32 ...

来自： print(ds.dims)

Frozen({'sample': 86401, 'ddm': 4, 'delay': 17, 'doppler': 11})

我尝试扩展维度以包括 sp_lat 和 sp_lon 以及

ds.expand_dims(['x', 'y']
ds.rename_vars({'sp_lon': 'x', 'sp_lat': 'x'})

我也试过了

ds.rename({'sp_lon': 'x', 'sp_lat': 'y'})
ds.rio.set_spatial_dims('x', 'y', inplace=True)

我什至尝试过多重索引。我怎样才能使用我的坐标作为我的空间维度，以便我可以用 ds.rio.clip(...)

我在 python 3.9.

中使用 xarray 和 rioxarray

Answer 1

rioxarray 希望您的数据位于规则的网格上。您的数据似乎是观测数据，在那个时间步为位置（可能是给定的航天器？）给出了纬度和经度值。您将无法扩展数组的维度以包含 (lat, lon)。

相反，您可以直接使用数组 sp_lat、sp_lon 构建一个数组，指示给定的 (sample, ddm) 点是否包含在 shapefile 中。

国家分配示例

例如，如果您有国家/地区的 shapefile：

In [7]: countries = gpd.read_file('https://naturalearth.s3.amazonaws.com/110m_cultural/ne_110m_admin_0_countries.zip')

和以下数据集：

In [13]: sample = np.arange(100)
    ...: ddm = np.arange(4)
    ...: timestep = pd.date_range('2020-01-01', periods=len(sample), freq='H')
    ...: sp_lat = np.random.random(size=(len(sample), len(ddm))) * 180 - 90
    ...: sp_lon = np.random.random(size=(len(sample), len(ddm))) * 360 - 180
    ...:
    ...: ds = xr.Dataset(
    ...:     {},
    ...:     coords={
    ...:         'sample': sample,
    ...:         'ddm': ddm,
    ...:         'ddm_timestamp_utc': (('sample', ), timestep),
    ...:         'sp_lat': (('sample', 'ddm'), sp_lat),
    ...:         'sp_lon': (('sample', 'ddm'), sp_lon),
    ...:     },
    ...: )

In [14]: ds
Out[14]:
<xarray.Dataset>
Dimensions:            (sample: 100, ddm: 4)
Coordinates:
  * sample             (sample) int64 0 1 2 3 4 5 6 7 ... 93 94 95 96 97 98 99
  * ddm                (ddm) int64 0 1 2 3
    ddm_timestamp_utc  (sample) datetime64[ns] 2020-01-01 ... 2020-01-05T03:0...
    sp_lat             (sample, ddm) float64 28.11 -88.63 15.52 ... 70.92 -51.87
    sp_lon             (sample, ddm) float64 -46.9 132.9 ... -70.12 -161.3
Data variables:
    *empty*

您可以将 x、y 观测值转换为 geopandas GeoDataFrame：

In [15]: x_flat = ds.sp_lon.values.ravel()
    ...: y_flat = ds.sp_lat.values.ravel()

In [19]: xy_point_array = gpd.GeoDataFrame(
    ...:     geometry=gpd.points_from_xy(x_flat, y_flat, crs='epsg:4326')
    ...: )

然后，使用 sjoin 将 shapefile 中的所有国家分配给一个点（对于不属于国家的点返回 NaN）：

In [20]: countries_by_point = xy_point_array.sjoin(countries, how='left')

然后可以将结果重新整形为点数组的维度并返回到 xarray:

In [24]: ds.coords['country'] = (
    ...:     ('sample', 'ddm'),
    ...:     countries_by_point.ADM0_A3.values.reshape(sp_lat.shape),
    ...: )

In [25]: ds
Out[25]:
<xarray.Dataset>
Dimensions:            (sample: 100, ddm: 4)
Coordinates:
  * sample             (sample) int64 0 1 2 3 4 5 6 7 ... 93 94 95 96 97 98 99
  * ddm                (ddm) int64 0 1 2 3
    ddm_timestamp_utc  (sample) datetime64[ns] 2020-01-01 ... 2020-01-05T03:0...
    sp_lat             (sample, ddm) float64 28.11 -88.63 15.52 ... 70.92 -51.87
    sp_lon             (sample, ddm) float64 -46.9 132.9 ... -70.12 -161.3
    country            (sample, ddm) object nan 'ATA' nan nan ... 'ATA' nan nan
Data variables:
    *empty*

xarray - MissingSpatialDimensionError（将坐标指定为维度）

xarray - MissingSpatialDimensionError (assign coordinates as dimensions)

python

netcdf

python-xarray

netcdf4

国家分配示例