如何基于同心椭圆数据生成带填充的二维热力图并导出TIFF？

核心思路 #

直接用绘图库绘制椭圆再转TIFF易丢失地理精度，建议采用栅格化赋值思路：先创建对应地理范围的空栅格，再为每个椭圆区域内的像素赋予对应的arf_w权重，最后导出带地理参考的标准TIFF。

具体实现步骤 #

1. 明确地理参数 #

提前确认椭圆所在区域的地理坐标范围（如左上角/右下角投影坐标）、像素分辨率，这是GIS TIFF识别地理信息的核心。

2. 栅格化椭圆区域 #

用numpy创建空栅格，结合matplotlib.patches.Ellipse判断像素是否在椭圆内，批量赋值权重：

import numpy as np
from matplotlib.patches import Ellipse

# 假设已有的椭圆参数列表：(中心x, 中心y, 长半轴, 短半轴, 旋转角, arf_w)
ellipses = [(x0, y0, a, b, angle, w) for ...]

# 定义栅格参数（替换为你的实际地理范围）
x_min, x_max = 100.0, 101.0
y_min, y_max = 30.0, 31.0
pixel_size = 0.001  # 每个像素代表的地理单位
width = int((x_max - x_min) / pixel_size)
height = int((y_max - y_min) / pixel_size)

# 初始化空栅格
grid = np.zeros((height, width), dtype=np.float32)

# 遍历椭圆填充栅格
for x0, y0, a, b, angle, w in ellipses:
    ellipse = Ellipse((x0, y0), 2*a, 2*b, angle=angle)
    # 遍历每个像素中心坐标
    for i in range(height):
        y = y_max - i * pixel_size - pixel_size/2  # TIFF行从上到下对应y递减
        for j in range(width):
            x = x_min + j * pixel_size + pixel_size/2
            if ellipse.contains_point((x, y)):
                # 重叠区域取权重最大值（可按需改为累加）
                if w > grid[i, j]:
                    grid[i, j] = w

3. 导出GIS标准TIFF #

用gdal库添加地理参考并导出，确保GIS软件可识别：

from osgeo import gdal, osr

# 创建TIFF驱动
driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create('ellipse_heatmap.tif', width, height, 1, gdal.GDT_Float32)

# 设置地理变换参数（左上角x, 像素宽度, 旋转, 左上角y, 旋转, 像素高度）
geotransform = (x_min, pixel_size, 0, y_max, 0, -pixel_size)
dataset.SetGeoTransform(geotransform)

# 设置投影（示例为WGS84，替换为你的实际投影）
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
dataset.SetProjection(srs.ExportToWkt())

# 写入数据并设置无数据值
band = dataset.GetRasterBand(1)
band.WriteArray(grid)
band.SetNoDataValue(0)

# 释放资源
dataset = None

性能优化方案 #

若栅格尺寸过大，遍历像素效率低，可改用shapely+rasterio批量栅格化：

from shapely.geometry import Point
from shapely.affinity import rotate, scale
from rasterio.features import rasterize
from rasterio.transform import from_origin

# 将椭圆转为shapely多边形
ellipse_polygons = []
values = []
for x0, y0, a, b, angle, w in ellipses:
    circle = Point(x0, y0).buffer(1)
    ellipse = scale(circle, xfact=a, yfact=b)
    ellipse = rotate(ellipse, angle, origin=(x0, y0))
    ellipse_polygons.append(ellipse)
    values.append(w)

# 创建rasterio变换对象
transform = from_origin(x_min, y_max, pixel_size, pixel_size)

# 批量栅格化
grid = rasterize(
    [(poly, val) for poly, val in zip(ellipse_polygons, values)],
    out_shape=(height, width),
    transform=transform,
    fill=0,
    dtype=np.float32
)

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