Python中将非网格化经纬度CSV转为栅格对象的合适方法

Python中没有R语言terra::rast()那样一行直接转换的极简API，但可以通过rasterio结合pandas、numpy快速实现，以下是两种常见场景的简洁实现方案：

场景1：点数据栅格化（计数/聚合） #

适用于统计每个栅格单元内的点数量，或者对属性值求和、取均值等。

import pandas as pd
import rasterio
from rasterio.transform import from_origin
import numpy as np

# 读取数据
df = pd.read_csv("gbif.csv")

# 定义栅格参数（分辨率可自定义）
res = 0.1  # 栅格单元的度数
x_min, x_max = df.longitude.min(), df.longitude.max()
y_min, y_max = df.latitude.min(), df.latitude.max()

# 计算栅格行列数与转换矩阵
width = int((x_max - x_min) / res)
height = int((y_max - y_min) / res)
transform = from_origin(x_min, y_max, res, res)

# 聚合点数据到栅格（这里以计数为例，可替换为其他聚合逻辑）
raster_data, _, _ = np.histogram2d(
    df.longitude, df.latitude,
    bins=[width, height],
    range=[[x_min, x_max], [y_min, y_max]]
)

# 创建并保存栅格（或返回内存中的rasterio对象）
with rasterio.open(
    "output_count_raster.tif",
    "w",
    driver="GTiff",
    height=height,
    width=width,
    count=1,
    dtype=raster_data.dtype,
    crs="EPSG:4326",
    transform=transform,
) as dst:
    dst.write(raster_data, 1)

场景2：离散点属性插值到栅格 #

如果需要将点的属性值（如密度、观测值等）插值填充到栅格中，可使用scipy的插值工具：

import pandas as pd
import rasterio
from rasterio.transform import from_origin
import numpy as np
from scipy.interpolate import griddata

# 读取数据（假设包含需要插值的属性列，比如occurrenceCount）
df = pd.read_csv("gbif.csv")

# 定义栅格参数
res = 0.1
x_min, x_max = df.longitude.min(), df.longitude.max()
y_min, y_max = df.latitude.min(), df.latitude.max()
width = int((x_max - x_min) / res)
height = int((y_max - y_min) / res)
transform = from_origin(x_min, y_max, res, res)

# 创建栅格网格
x_grid, y_grid = np.meshgrid(
    np.linspace(x_min, x_max, width),
    np.linspace(y_max, y_min, height)  # 匹配rasterio的y轴从上到下的顺序
)

# 插值（可选nearest/linear/cubic方法）
interpolated_data = griddata(
    (df.longitude, df.latitude),
    df["occurrenceCount"],
    (x_grid, y_grid),
    method="nearest"
)

# 保存栅格
with rasterio.open(
    "output_interpolated_raster.tif",
    "w",
    driver="GTiff",
    height=height,
    width=width,
    count=1,
    dtype=interpolated_data.dtype,
    crs="EPSG:4326",
    transform=transform,
) as dst:
    dst.write(interpolated_data, 1)

封装为复用函数（接近R的简洁调用） #

如果需要频繁使用，可以封装成函数，实现类似R的一行式调用体验：

import pandas as pd
import rasterio
from rasterio.transform import from_origin
import numpy as np
from scipy.interpolate import griddata

def df_to_raster(df, x_col="longitude", y_col="latitude", value_col=None, res=0.1, crs="EPSG:4326"):
    # 计算栅格范围与尺寸
    x_min, x_max = df[x_col].min(), df[x_col].max()
    y_min, y_max = df[y_col].min(), df[y_col].max()
    width = int((x_max - x_min) / res)
    height = int((y_max - y_min) / res)
    transform = from_origin(x_min, y_max, res, res)
    
    # 生成栅格数据
    if value_col is None:
        # 默认统计点数量
        data, _, _ = np.histogram2d(
            df[x_col], df[y_col],
            bins=[width, height],
            range=[[x_min, x_max], [y_min, y_max]]
        )
    else:
        # 插值属性值
        x_grid, y_grid = np.meshgrid(
            np.linspace(x_min, x_max, width),
            np.linspace(y_max, y_min, height)
        )
        data = griddata(
            (df[x_col], df[y_col]),
            df[value_col],
            (x_grid, y_grid),
            method="nearest"
        )
    
    # 返回内存中的rasterio Dataset对象
    memfile = rasterio.MemoryFile()
    with memfile.open(
        driver="GTiff", height=height, width=width, count=1,
        dtype=data.dtype, crs=crs, transform=transform
    ) as dataset:
        dataset.write(data, 1)
    return memfile.open()

# 调用示例
df = pd.read_csv("gbif.csv")
raster = df_to_raster(df, value_col="occurrenceCount")

# 查看栅格信息
print(raster.meta)

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