如何统计与栅格像素相交的矢量要素数量及OSM转权重栅格方案

需求背景 #

本需求旨在寻找更快更准确的方法，将OpenStreetMap小片提取结果栅格化为人口加权栅格。
我想要将小型.pbf文件转换为便于后续空间分析的GeoTiff格式。本次问题的需求仅限定为处理多边形几何，因为我已经找到适用于线要素的高性能解决方案，该方案效果极佳，我甚至考虑将所有多边形转为线要素处理。

数据处理示例 #

以下是我想要转换的一类数据的处理示例：

wget https://download.geofabrik.de/europe/liechtenstein-latest.osm.pbf
osmium tags-filter liechtenstein-latest.osm.pbf landuse=grass -o liechtenstein_grass.pbf
ogr2ogr -t_srs EPSG:3857 liechtenstein_grass.sqlite -dsco SPATIALITE=YES -nln multipolygons -nlt POLYGON -skipfailures liechtenstein_grass.pbf

现有基础与预期目标 #

我找到了一个分区统计脚本可以作为解决该问题的基础，该脚本输入矢量层和栅格层，遍历矢量要素裁剪栅格并执行对应统计。
我不需要常规的分区统计结果，想要统计每个栅格像素相交的矢量要素数量。我已经生成了Int32格式的全球栅格网格，每个像素带有唯一值：

{qgis_process} run native:creategrid -- TYPE=2 EXTENT="-20037760, -8399416, 20037760, 18454624 [EPSG:3857]" HSPACING=1912 VSPACING=1912 HOVERLAY=0 VOVERLAY=0 CRS="EPSG:3857" OUTPUT="grid.gpkg"

sqlite3 land.gpkg
SELECT load_extension("mod_spatialite");
alter table output add column ogcfod int;
update output set ogcfod = fid;

gdal_rasterize -l output -a ogcfod -tap -tr 1912.0 1912.0 -a_nodata 0.0 -ot Int32 -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=2 grid.gpkg grid.tif -te -20037760 -8399416 20037760 18454624

我认为遍历矢量要素的方案性能更高，因为矢量要素数量远少于栅格的8800万+分区。
我想要实现的脚本逻辑为：输入矢量层和栅格层，遍历矢量要素，查询每个要素覆盖的所有像素值，为字典{px_id: qty}中对应像素的计数加1。

问题现象 #

但我调试脚本时始终返回相同几何，反复只输出1个像素ID，代码如下：

import sys
import gdal
import numpy
import ogr
import osr
from rich import inspect, print


def zonal_stats(feat, lyr, raster):
    # Get raster georeference info
    transform = raster.GetGeoTransform()
    xOrigin = transform[0]
    yOrigin = transform[3]
    pixelWidth = transform[1]
    pixelHeight = transform[5]

    # Reproject vector geometry to same projection as raster
    sourceSR = lyr.GetSpatialRef()
    targetSR = osr.SpatialReference()
    targetSR.ImportFromWkt(raster.GetProjectionRef())
    coordTrans = osr.CoordinateTransformation(sourceSR, targetSR)
    feat = lyr.GetNextFeature()
    geom = feat.GetGeometryRef()
    geom.Transform(coordTrans)

    # Get extent of feat
    geom = feat.GetGeometryRef()
    if geom.GetGeometryName() == "MULTIPOLYGON":
        count = 0
        pointsX = []
        pointsY = []
        for polygon in geom:
            geomInner = geom.GetGeometryRef(count)
            ring = geomInner.GetGeometryRef(0)
            numpoints = ring.GetPointCount()
            for p in range(numpoints):
                lon, lat, z = ring.GetPoint(p)
                pointsX.append(lon)
                pointsY.append(lat)
            count += 1
    elif geom.GetGeometryName() == "POLYGON":
        ring = geom.GetGeometryRef(0)
        numpoints = ring.GetPointCount()
        pointsX = []
        pointsY = []
        for p in range(numpoints):
            lon, lat, z = ring.GetPoint(p)
            pointsX.append(lon)
            pointsY.append(lat)

    else:
        sys.exit("ERROR: Geometry needs to be either Polygon or Multipolygon")

    xmin = min(pointsX)
    xmax = max(pointsX)
    ymin = min(pointsY)
    ymax = max(pointsY)

    print(xmin, xmax)
    print(ymin, ymax)

    # Specify offset and rows and columns to read
    xoff = int((xmin - xOrigin) / pixelWidth)
    yoff = int((yOrigin - ymax) / pixelWidth)
    xcount = int((xmax - xmin) / pixelWidth) + 1
    ycount = int((ymax - ymin) / pixelWidth) + 1

    # Create memory target raster
    target_ds = gdal.GetDriverByName("MEM").Create(
        "", xcount, ycount, 1, gdal.GDT_Int32
    )
    target_ds.SetGeoTransform(
        (
            xmin,
            pixelWidth,
            0,
            ymax,
            0,
            pixelHeight,
        )
    )

    # Create for target raster the same projection as for the value raster
    raster_srs = osr.SpatialReference()
    raster_srs.ImportFromWkt(raster.GetProjectionRef())
    target_ds.SetProjection(raster_srs.ExportToWkt())

    # Rasterize zone polygon to raster
    gdal.RasterizeLayer(target_ds, [1], lyr, burn_values=[1])

    # Read raster as arrays
    banddataraster = raster.GetRasterBand(1)
    dataraster = banddataraster.ReadAsArray(xoff, yoff, xcount, ycount)

    bandmask = target_ds.GetRasterBand(1)
    datamask = bandmask.ReadAsArray(0, 0, xcount, ycount)

    print(dataraster)

    # Mask zone of raster
    # zoneraster = numpy.ma.masked_array(dataraster, numpy.logical_not(datamask))

    # print(zoneraster)

    # exit()


def loop_zonal_stats(lyr, raster):
    featList = range(lyr.GetFeatureCount())
    statDict = {}

    for FID in featList:
        feat = lyr.GetFeature(FID)
        meanValue = zonal_stats(feat, lyr, raster)
        statDict[FID] = meanValue
    return statDict


def main(input_zonal_raster, input_value_polygon):
    raster = gdal.Open(input_zonal_raster)
    shp = ogr.Open(input_value_polygon)
    lyr = shp.GetLayer()

    return loop_zonal_stats(lyr, raster)


if __name__ == "__main__":
    if len(sys.argv) != 3:
        print(
            "[ ERROR ] you must supply two arguments: input-zone-raster-name.tif input-value-shapefile-name.shp "
        )
        sys.exit(1)

    main(sys.argv[1], sys.argv[2])

前期调研内容 #

• 统计包含重复项的重叠多边形数量相关方案
• 栅格化叠加计数的Python实现参考
• 若gdal_rasterize支持烧录每个像素相交的多边形总数而非固定值，也可满足我的需求
• OSM汇总热力图相关实现参考
• QGIS统计重叠多边形相关方案
• heatmapkerneldensity效果不佳，可能是我使用方式有误，但输出结果不符合预期。

我测试的热力图生成命令如下：

{qgis_process} run qgis:heatmapkerneldensityestimation -- INPUT="{basen}.json.geojson" RADIUS=2868 RADIUS_FIELD=None PIXEL_SIZE=1912 WEIGHT_FIELD=None KERNEL=4 DECAY=0 OUTPUT_VALUE=0 OUTPUT="{basen}.tif

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