你遇到的错误:
```
rasterio.errors.WindowError: Bounds and transform are inconsistent
```
这是一个 **常见的 `rasterio` 几何操作错误**,通常出现在使用 `geometry_mask`、重投影或裁剪栅格时。它表示:**你提供的地理范围(bounds)与仿射变换矩阵(transform)之间存在不一致**,即坐标系统对不上。
---
## ✅ 错误原因详解
在你的代码中,很可能是这一段出错:
```python
from rasterio.features import geometry_mask
affine_transform = rasterio.Affine.translation(
TARGET_LONS[0, 0] - target_resolution / 2,
TARGET_LATS[0, 0] - target_resolution / 2
) * rasterio.Affine.scale(target_resolution, target_resolution)
mask = geometry_mask(
geometries=shapes,
out_shape=spei_interp.shape,
transform=affine_transform,
invert=True
)
```
### ❌ 问题出在哪?
- `TARGET_LONS[0, 0]` 是经纬度坐标(如 `104.8`)
- 但 `rasterio.Affine` 中的 `translation(xoff, yoff)` 要求的是左上角像元中心的地理坐标
- 如果仿射变换构造不当(比如偏移量没算准),就会导致:
- 输出 shape 和 transform 不匹配 bounds
- 几何体落在 transform 定义的网格之外
- 抛出 `WindowError: Bounds and transform are inconsistent`
---
## ✅ 正确解决方案
我们应使用更可靠的方式构造 `transform`,推荐通过 `rasterio.transform.from_origin` 或直接利用 `rioxarray` 提供的 `.rio.transform()` 方法。
---
### ✅ 推荐修复方式(修改原代码中的 `geometry_mask` 构造部分)
替换原来的 `affine_transform` 构造方式为以下**科学且健壮的方法**:
```python
import rasterio
from rasterio.transform import from_origin
from rasterio.features import geometry_mask
import numpy as np
# --- 正确构造 transform ---
# 左上角坐标(注意:是左上角外边缘,不是第一个像素中心)
left = extent[0] # W
top = extent[3] # N
res = target_resolution
# 使用 rasterio 的 from_origin 创建正确仿射变换
# from_origin(left, top, x_size, y_size)
transform = from_origin(west=left, north=top, xsize=res, ysize=res)
# --- 验证 transform 是否与目标网格一致 ---
height, width = spei_interp.shape # 插值后的高分辨率网格
# 可选:打印验证
print(f"Expected shape: ({height}, {width})")
print(f"Transform:\n{transform}")
# --- 创建掩膜 ---
mask = geometry_mask(
geometries=shapes,
out_shape=(height, width),
transform=transform,
invert=True, # 内部为 True
all_touched=False # 只有完全包含的像素才设为 True
)
```
---
### 🔍 关键说明
| 概念 | 说明 |
|------|------|
| `from_origin(west, north, xsize, ysize)` | 自动构建仿射变换,确保 bounds 和 transform 一致 |
| `west=extent[0], north=extent[3]` | 网格左上角地理坐标 |
| `xsize=ysize=target_resolution` | 像元大小(0.01°) |
| `out_shape=(height, width)` | 必须与你要生成的数组形状完全一致 |
> ⚠️ 特别注意:`geometry_mask` 中的 `out_shape` 必须是 `(rows, cols)` 即 `(lat_len, lon_len)`,顺序不能反!
---
## ✅ 更进一步:使用 `rioxarray` 替代手动构造(推荐!)
既然你已经在用 `rioxarray` 处理 DEM,不如也用它来处理 mask,避免手动管理 transform。
### ✅ 改进方案:将插值后的 SPEI 转为 `xarray.DataArray` 并绑定 CRS 和 transform
```python
import xarray as xr
# 将插值后的 SPEI 包装成带坐标的 DataArray
spei_da = xr.DataArray(
data=spei_interp,
coords={'y': target_lats[::-1], 'x': target_lons}, # 注意:target_lats 应该是从北到南
dims=['y', 'x']
).rio.write_crs("EPSG:4326")
# 直接使用 rioxarray 的 clip 方法(自动处理 transform)
spei_masked_rio = spei_da.rio.clip(gdf.geometry, gdf.crs, drop=False, invert=False)
# 获取掩膜数组
spei_masked = spei_masked_rio.values
```
👉 这种方法完全绕过 `rasterio.features.geometry_mask` 和手动 transform 构造,**从根本上避免了 `WindowError`**
---
## ✅ 最终建议流程(推荐采用)
```python
# 1. 构造高分辨率网格
target_lons = np.arange(extent[0], extent[1], target_resolution)
target_lats = np.arange(extent[3], extent[2], -target_resolution) # 从北往南
# 2. 插值得到 spei_interp (shape: len(lats), len(lons))
# 3. 转为 xarray 并写入 CRS
spei_da = xr.DataArray(
spei_interp,
coords={'y': target_lats, 'x': target_lons},
dims=['y', 'x']
).rio.write_crs("EPSG:4326")
# 4. 使用 rioxarray 自动裁剪
spei_clipped = spei_da.rio.clip(gdf.geometry, gdf.crs, drop=False)
# 5. 后续绘图直接使用 spei_clipped.values 和 .coords
```
✅ 优点:
- 不用手动构造 `transform`
- 自动处理边界、投影、分辨率
- 兼容性好,适合科研自动化
---
## ✅ 总结:如何解决 `WindowError: Bounds and transform are inconsistent`
| 方法 | 是否推荐 | 说明 |
|------|----------|------|
| ✅ 使用 `rasterio.transform.from_origin` | ✔️ 推荐 | 手动构造中最安全的方式 |
| ✅ 使用 `rioxarray.rio.clip()` | ✅✅ 强烈推荐 | 完全避开底层细节,适合科研项目 |
| ❌ 手动拼 `Affine.translation + scale` | ❌ 不推荐 | 容易因浮点误差导致不一致 |
---