在y维上切片xarray DataArray时出现意外行为

Question

我刚接触xarray，想弄明白sel()是怎么工作的。

我有一个数据数组tif_xr，它是通过rioxarray打开一个有10个波段的GeoTIFF文件得到的。

tif_xr = rioxarray.open_rasterio('path/to/my/file.tif', chunks={'x':100, 'y':100})

我可以在x维度上切片，这样得到的结果是我预期的。

tif_xr.sel(x=slice(500505,500520))

结果是：

<xarray.DataArray (band: 10, y: 10900, x: 2)> Size: 436kB
dask.array<getitem, shape=(10, 10900, 2), dtype=uint16, chunksize=(10, 100, 2), chunktype=numpy.ndarray>
Coordinates:
  * band         (band) int64 80B 1 2 3 4 5 6 7 8 9 10
  * x            (x) float64 16B 5.005e+05 5.005e+05
  * y            (y) float64 87kB 4.5e+06 4.5e+06 ... 4.391e+06 4.391e+06
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:  Area
    _FillValue:     65535
    scale_factor:   1.0
    add_offset:     0.0

但是，当我用同样的切片逻辑去处理y维度时，却没有得到任何值。

tif_xr.sel(y=slice(4444550,4444560))

结果是：

<xarray.DataArray (band: 10, y: 0, x: 10924)> Size: 0B
dask.array<getitem, shape=(10, 0, 10924), dtype=uint16, chunksize=(10, 0, 100), chunktype=numpy.ndarray>
Coordinates:
  * band         (band) int64 80B 1 2 3 4 5 6 7 8 9 10
  * x            (x) float64 87kB 5.005e+05 5.005e+05 ... 6.097e+05 6.097e+05
  * y            (y) float64 0B 
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:  Area
    _FillValue:     65535
    scale_factor:   1.0
    add_offset:     0.0

这很奇怪，因为当我用sel()并且用等于号而不是切片时（tif_xr.sel(y=4444555)），似乎是可以的：

<xarray.DataArray (band: 10, x: 10924)> Size: 218kB
dask.array<getitem, shape=(10, 10924), dtype=uint16, chunksize=(10, 100), chunktype=numpy.ndarray>
Coordinates:
  * band         (band) int64 80B 1 2 3 4 5 6 7 8 9 10
  * x            (x) float64 87kB 5.005e+05 5.005e+05 ... 6.097e+05 6.097e+05
    y            float64 8B 4.445e+06
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:  Area
    _FillValue:     65535
    scale_factor:   1.0
    add_offset:     0.0

我尝试了不同的切片范围，但都没有成功。

我的GeoTiff文件是用投影坐标系（crs）表示的，所以我在这个例子中使用的x和y值实际上是以米为单位的坐标。

这可能是什么情况呢？我是不是漏掉了什么？

谢谢。

Answer 1

仔细看看你数据集中的y坐标：

其实这些值是逐渐减小的。

所以，当你切片的时候，得先给出上限再给下限：

tif_xr.sel(y=slice(4444560, 4444550))

这样应该能得到你想要的结果。

另外，你也可以指定一个负的步长：

tif_xr.sel(y=slice(4444550, 4444560, -1))

这样也应该可以。

我建议你在加载数据集后，直接把y轴反转，这样就能完全避免这个问题：

# Invert the y-axis to have increasing coordinates
tif_xr = tif_xr.isel(y=slice(None, None, -1))

这样，你所有的 sel() 命令就应该能正常工作了。