如何使用Python将地球球形图像展平为适配Leaflet的平面投影图像

完全可以通过Python实现和ArcGIS效果一致的投影变换，输出适配Leaflet的无偏移叠加图，具体实现方案如下：

依赖安装 #

你需要安装地理空间投影处理和图像处理的相关库：
pip install rasterio pyproj numpy pillow

核心实现逻辑 #

你之前用ArcGIS做的“展平”操作本质是坐标系投影转换：把原始NASA图像的球坐标系（一般是EPSG:4326 WGS84地理坐标系，或者对应卫星的自定义投影坐标系）转换为Leaflet默认使用的EPSG:3857 Web墨卡托投影坐标系，转换过程会自动计算目标图像的尺寸、重采样像素值，不是简单的resize。

具体实现步骤 #

1. 确认原始图像的空间元数据
   你需要先获取两个核心参数，这两个参数在你用ArcGIS打开HDF5文件时会自动读取，代码中可以直接从原始HDF5文件提取，也可以根据NASA产品的公开参数手动填写：

• 原始图像的坐标系（CRS），比如大部分NASA公开地球影像默认是EPSG:4326
• 原始图像对应的地理边界：最小经度、最大经度、最小纬度、最大纬度，比如你合并后的4800x2400图像如果是覆盖全球的话就是[-180, 180, -90, 90]

2. 坐标转换与图像重采样
   示例代码如下，你可以根据自己的参数修改对应字段：

import numpy as np
from PIL import Image
import rasterio
from rasterio.transform import from_bounds
from rasterio.warp import calculate_default_transform, reproject, Resampling
from pyproj import CRS

# -------------------------- 替换为你自己的参数 --------------------------
# 原始图像参数
raw_img = Image.open("你的4800x2400热力图路径.png")
raw_arr = np.array(raw_img).transpose(2, 0, 1)  # 把(H,W,C)转为(C,H,W)适配rasterio格式
raw_crs = CRS.from_epsg(4326)  # 原始图像的坐标系，这里以WGS84为例
raw_bounds = [-180, 180, -90, 90]  # 原始图像的经纬度范围[min_lon, max_lon, min_lat, max_lat]
# 目标投影参数（Leaflet默认用EPSG:3857）
dst_crs = CRS.from_epsg(3857)
resampling_method = Resampling.bilinear  # 热力图用双线性重采样即可，和ArcGIS效果一致
# -----------------------------------------------------------------------

# 计算原始图像的变换矩阵
raw_height, raw_width = raw_arr.shape[1], raw_arr.shape[2]
raw_transform = from_bounds(*raw_bounds, width=raw_width, height=raw_height)

# 计算投影后的目标图像参数：变换矩阵、宽、高
dst_transform, dst_width, dst_height = calculate_default_transform(
    raw_crs, dst_crs, raw_width, raw_height, *raw_bounds
)

# 创建目标数组
dst_arr = np.zeros((raw_arr.shape[0], dst_height, dst_width), dtype=raw_arr.dtype)

# 执行投影重采样
reproject(
    source=raw_arr,
    destination=dst_arr,
    src_transform=raw_transform,
    src_crs=raw_crs,
    dst_transform=dst_transform,
    dst_crs=dst_crs,
    resampling=resampling_method
)

# 转为PIL Image并保存
dst_img = Image.fromarray(dst_arr.transpose(1,2,0))
dst_img.save("适配Leaflet的投影后图像.png")

# 输出投影后的地理边界，Leaflet叠加时直接用这个做bounds即可
dst_bounds = rasterio.transform.array_bounds(dst_height, dst_width, dst_transform)
print(f"Leaflet叠加用边界：[ [{dst_bounds[1]}, {dst_bounds[0]}], [{dst_bounds[3]}, {dst_bounds[2]}] ]")

注意事项 #

• 如果你不确定原始HDF5对应的坐标系和边界，可以用gdalinfo命令直接读取HDF5文件的元数据，和ArcGIS读取到的参数完全一致。
• 重采样方法可以根据你的图像类型调整：分类数据用最近邻，连续的热力图、影像用双线性/三次卷积即可。

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