在 Python 3.6 中从点（gpd.geodataframe 对象）创建栅格

Question

我想使用 geopandas.geodataframe.GeoDataFrame 对象从点文件

创建光栅文件 (.tif)。

我的数据框有两列：[几何] 和 [值]。目标是使用

[Value] 值在 [几何] 点制作 10m 分辨率栅格。

我的数据集是：

geometry |  Value
0  | POINT (520595.000 5720335.000) | 536.678345
1  | POINT (520605.000 5720335.000) | 637.052185
2  | POINT (520615.000 5720335.000) | 1230.553955
3  | POINT (520625.000 5720335.000) | 944.970642
4  | POINT (520635.000 5720335.000) | 1094.613281
5  | POINT (520645.000 5720335.000) | 1123.185181
6  | POINT (520655.000 5720335.000) | 849.37634
7  | POINT (520665.000 5720335.000) | 1333.459839
8  | POINT (520675.000 5720335.000) | 492.866608
9  | POINT (520685.000 5720335.000) | 960.957214
10 | POINT (520695.000 5720335.000) | 539.401978
11 | POINT (520705.000 5720335.000) | 573.015625
12 | POINT (520715.000 5720335.000) | 970.386536
13 | POINT (520725.000 5720335.000) | 390.315094
14 | POINT (520735.000 5720335.000) | 642.036865

我之前尝试过，所以，我知道使用

from geocube.api.core import make_geocube

我可以做到这一点，但由于某些库我有限制，我无法使用

make_geocube

。
有什么想法吗？

Answer 1

分配 x 和 y 列，转换为 xarray，然后使用

rioxarray 导出到 tiff：

# do this before sending to xarray
# to ensure extension is loaded
import rioxarray

# assuming your GeoDataFrame is called `gdf`
gdf["x"] = gdf.x
gdf["y"] = gdf.y

da = (
    gdf.set_index(["y", "x"])
    .Value
    .to_xarray()
)

da.rio.to_raster("myfile.tif")

为了使其发挥作用，这些点必须组成一个完整的规则网格，每个组合的 x 和 y 值都重复。如果这只是转换为 xarray 的任意点的集合，其中 x 和 y 作为垂直维度，将会爆炸你的记忆，结果将几乎完全是 NaN。

Answer 2

如果您不希望单个像素内包含多个点，那么您可以使用

rasterio.features.rasterize

 来实现此目的。但请记住，如果确实有多个像素与单个像素重叠，则 

rasterio.features.rasterize

 只能：

刻录单个硬编码值，
将重叠点值加在一起，或者
只需使用最后看到的点值

在栅格化分类值的情况下，这应该足够了，但对于栅格化连续数据，您可能想要取每个像素的平均值。

为此，您可以执行以下操作：

import numpy as np
from typing import Union


def rasterize_points(
    points: np.ndarray, res: Union[int, float], bbox: tuple
) -> np.ndarray:
    """Rasterize points into a grid with a given resolution.

    Args:
        points (np.ndarray): Points to rasterize, with columns (x, y, value) (for
            geographic coordinates, use (lon, lat, value))
        res (Union[int, float]): Resolution of the grid
        bbox (tuple): Bounding box of the grid

    Returns:
        np.ndarray: Rasterized grid
    """
    width = int((bbox[2] - bbox[0]) / res)
    height = int((bbox[3] - bbox[1]) / res)

    rast = np.zeros((height, width), dtype=np.float32)
    count_array = np.zeros_like(rast)

    for x, y, value in points:
        col = int((x - bbox[0]) / res)
        row = int((bbox[3] - y) / res)
        rast[row, col] += value
        count_array[row, col] += 1

    # Avoid division by zero
    count_array[count_array == 0] = 1

    # Calculate the average
    rast = rast / count_array

    return rast

上面的例子不是很优化，但它完成了工作。请注意，它仅返回一个 2D numpy 数组。您需要使用

rasterio

 或 

rioxarray

 合并相关坐标和其他地理空间信息（例如 CRS 和变换）。