自制Raster与Shapefile叠加错位问题修复求助

我在项目中制作巴西南里奥格兰德州的Raster数据，生成过程正常，但与该州Shapefile同轴绘图时，Shapefile显示过小且倒置，无法正确叠加——后续需要基于叠加结果做分区统计，必须解决此问题。

裁剪NetCDF并生成Raster的代码 #

# 定义南里奥格兰德州的经纬度范围
min_lat, max_lat = -37.0, -23.0
min_lon, max_lon = -61.0, -46.0

# 打开数据集并裁剪
ds = xarray.open_dataset("path")
ds_clipped = ds.sel(lat=slice(min_lat, max_lat), lon=slice(min_lon, max_lon))

# 从xarray数据集GDD_arr生成Raster
outpath = "path"
with rasterio.open(
    outpath,
    mode="w",
    driver="GTiff",
    height=GDD_arr.shape[1],
    width=GDD_arr.shape[2],
    count=1,
    dtype='float32',
    crs= CRS.from_epsg(4326),
    transform = from_origin(-60.875, -23.125, 0.25, 0.25),
    nodata = -9999) as new_dataset:
  new_dataset.write(GDD_arr[0], 1)

绘图代码及异常效果 #

# 打开Raster和Shapefile
tiff = rasterio.open("GDD_Count_ssp245_2020_tobacco.tif")
shapefile = gpd.read_file("RGdS_AgData.shp")

# 设置绘图
f, ax = plt.subplots()

# 绘制Raster
rasterplot.show(np.flipud(tiff.read(1)), ax=ax)

# 绘制Shapefile
shapefile.plot(ax=ax, facecolor='w', edgecolor='k')
plt.show()

叠加错位效果 #

我曾尝试将Raster经度范围从0-360转换为-180-180，但问题依旧，Shapefile仍显示过小且倒置。以下是预期效果示意图：

预期效果 #

问题核心在于Raster的地理变换参数错误，以及绘图时破坏了地理坐标对应关系，以下是分步修复方案：

1. 修正Raster生成时的地理变换（Transform） #

手动输入from_origin参数容易出错，建议直接从裁剪后的xarray数据集提取地理信息生成变换矩阵，确保和原始网格对齐：

from rasterio.transform import from_bounds

# 从裁剪后的数据集提取经纬度范围与分辨率
left = ds_clipped.lon.min().values
right = ds_clipped.lon.max().values
bottom = ds_clipped.lat.min().values
top = ds_clipped.lat.max().values
x_res = (right - left) / len(ds_clipped.lon)
y_res = (top - bottom) / len(ds_clipped.lat)

# 若纬度数组是从北到南排列，y分辨率需设为负数（纬度随像素行增加而减小）
if ds_clipped.lat[0] > ds_clipped.lat[-1]:
    y_res = -y_res

# 生成正确的地理变换矩阵
transform = from_bounds(left, bottom, right, top, len(ds_clipped.lon), len(ds_clipped.lat))

# 重新写入Raster
with rasterio.open(
    outpath,
    mode="w",
    driver="GTiff",
    height=ds_clipped.lat.shape[0],
    width=ds_clipped.lon.shape[0],
    count=1,
    dtype='float32',
    crs=CRS.from_epsg(4326),
    transform=transform,
    nodata=-9999
) as new_dataset:
    new_dataset.write(GDD_arr[0], 1)

2. 绘图时对齐地理坐标，取消手动翻转数组 #

之前用np.flipud翻转Raster会破坏坐标对应关系，正确做法是利用rasterio的元数据直接绘图：

tiff = rasterio.open("GDD_Count_ssp245_2020_tobacco.tif")
shapefile = gpd.read_file("RGdS_AgData.shp")

# 确保Shapefile与Raster的CRS一致
if shapefile.crs != tiff.crs:
    shapefile = shapefile.to_crs(tiff.crs)

f, ax = plt.subplots()
# 直接读取Raster并绑定变换矩阵，无需手动翻转
rasterplot.show(tiff.read(1), transform=tiff.transform, ax=ax)
shapefile.plot(ax=ax, facecolor='w', edgecolor='k')
plt.show()

3. 验证Shapefile的CRS与范围 #

检查并统一两者的坐标参考系，同时确认范围匹配：

# 打印CRS信息
print("Shapefile CRS:", shapefile.crs)
print("Raster CRS:", tiff.crs)
# 转换Shapefile到Raster的CRS
shapefile = shapefile.to_crs(CRS.from_epsg(4326))

# 打印范围信息，确认两者覆盖区域匹配
print("Shapefile bounds:", shapefile.total_bounds)
print("Raster bounds:", tiff.bounds)

内容的提问来源于stack exchange，提问作者emmahaggerty