Python在Google Maps/Mapbox底图叠加GeoTIFF栅格图失败求助

问题背景 #

现有EPSG:4326投影的rainfall_clipped.tif栅格文件，尝试用Contextily包在Mapbox底图上叠加显示。已尝试将栅格转换为EPSG:3857投影，但仍无效。Contextily能识别投影并生成对应范围的底图，栅格数值读取正常（颜色条显示正确），但栅格图层未出现在地图上，怀疑是图层顺序或其他问题。

现有代码 #

data = rasterio.open("rainfall_clipped.tif")

# Read the bounds
left, bottom, right, top = data.bounds

# Create a figure and axes
fig, ax = plt.subplots(figsize=(10, 10))

# Add the raster data to the plot using imshow
im = ax.imshow(data.read(1), extent=[left, right, bottom, top], cmap="Blues")

# Add a colorbar
fig.colorbar(im, ax=ax)

# Add a basemap using contextily
ctx.add_basemap(ax, crs=data.crs, source=ctx.providers.MapBox(accessToken="my_key", id="mapbox/satellite-v9"))

# Show the plot
plt.show()

效果说明 #

• 当前输出：仅显示Mapbox底图和颜色条，栅格图层未显示
• 预期效果：栅格图层叠加在底图上方，清晰呈现降雨分布

解决建议 #

• 调整图层堆叠顺序：Contextily的add_basemap默认会将底图画在最上层，覆盖之前的栅格。通过zorder参数明确层级，确保栅格在底图之上：

  # 给栅格设置更高的zorder值
  im = ax.imshow(data.read(1), extent=[left, right, bottom, top], cmap="Blues", zorder=1)
  # 底图设置更低的zorder值
  ctx.add_basemap(ax, crs=data.crs, source=ctx.providers.MapBox(accessToken="my_key", id="mapbox/satellite-v9"), zorder=0)
  

• 明确栅格显示范围与透明度：若栅格数值范围和显示设置不匹配，可能导致图层不可见。可手动指定显示范围并调整透明度：

  raster_data = data.read(1)
  im = ax.imshow(
      raster_data, 
      extent=[left, right, bottom, top], 
      cmap="Blues", 
      zorder=1,
      vmin=raster_data.min(), 
      vmax=raster_data.max(),
      alpha=0.7  # 调整透明度，避免完全遮挡底图
  )
  

• 强制统一投影：即使文档说明支持双投影，手动将栅格转换为EPSG:3857后再绘制，避免潜在的投影不匹配问题：

  from rasterio.warp import reproject, Resampling, calculate_default_transform
  
  dst_crs = 'EPSG:3857'
  # 计算转换参数
  transform, width, height = calculate_default_transform(
      data.crs, dst_crs, data.width, data.height, *data.bounds)
  kwargs = data.meta.copy()
  kwargs.update({
      'crs': dst_crs,
      'transform': transform,
      'width': width,
      'height': height
  })
  
  # 保存转换后的栅格
  with rasterio.open('rainfall_3857.tif', 'w', **kwargs) as dst:
      for i in range(1, data.count + 1):
          reproject(
              source=rasterio.band(data, i),
              destination=rasterio.band(dst, i),
              src_transform=data.transform,
              src_crs=data.crs,
              dst_transform=transform,
              dst_crs=dst_crs,
              resampling=Resampling.nearest)
  
  # 使用转换后的栅格重新绘制
  data_3857 = rasterio.open('rainfall_3857.tif')
  left, bottom, right, top = data_3857.bounds
  im = ax.imshow(data_3857.read(1), extent=[left, right, bottom, top], cmap="Blues", zorder=1)
  ctx.add_basemap(ax, crs=dst_crs, source=ctx.providers.MapBox(accessToken="my_key", id="mapbox/satellite-v9"), zorder=0)
  

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