如何计算与栅格单元相交的矢量几何图形的数量？

Question

或者：寻找一种更快、更准确的方法将小型 OpenStreetMap 提取物栅格化为人口加权栅格。

我想将一个小的 .pbf 文件转换为 GeoTiff，这样可以更轻松地进行进一步的空间分析。出于这个问题的目的，我将限制处理多边形几何的要求，因为我已经找到了a solution that works very well for lines。它工作得很好，我正在考虑将我所有的多边形转换成线。

举一个我想转换的数据类型的例子：

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

我在这里找到了一个区域统计脚本，我们或许可以构建它来解决这个问题：http://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#calculate-zonal-statistics

上面的脚本需要一个矢量图层和一个栅格层，通过裁剪栅格并对其进行一些统计来迭代矢量特征。

我想计算与每个光栅像素相交的矢量特征的数量，而不是普通的区域统计。我有一个全局光栅网格 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

所以我在想，如果我们仍然可以迭代矢量特征（因为矢量特征的数量要少得多，而且栅格网格中有 88m+ 区域），它的性能可能会好得多。

我们想要一个脚本脚本，它采用矢量图层和光栅层，迭代矢量特征，查找特征覆盖的所有像素的值，然后将一个添加到字典中：{px_id: qty}

但是，当尝试使此脚本正常工作时，它一直给我相同的几何图形...它只会一遍又一遍地向我显示 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])

先前的研究：

• https://gis.stackexchange.com/questions/177738/count-overlapping-polygons-including-duplicates
• https://stackoverflow.com/a/47443399/697964
• 如果 gdal_rasterize 可以计算所有与每个像素相交的多边形（而不是固定值），这可能会满足我的需求。
• https://github.com/rory/osm-summary-heatmap/blob/main/Makefile
• https://old.reddit.com/r/gis/comments/4n2q5v/count_overlapping_polygons_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

Answer 1

这似乎有效。输出为 CSV columns=["px​​_id", "qty"]

"""
python rasterlayerzonalvectorcounts.py grid.tif liechtenstein_grass.sqlite

MIT License
Based on https://github.com/pcjericks/py-gdalogr-cookbook/blob/master/raster_layers.rst#calculate-zonal-statistics
"""

import sys
import osr
import os
import ogr
import numpy
import gdal
import pandas
from joblib import Parallel, delayed
from collections import Counter


def zonal_stats(FID, lyr, raster):
    feat = lyr.GetFeature(FID)

    # 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)
    geom = feat.GetGeometryRef()
    geom.Transform(coordTrans)

    # Get extent of feat
    geom = feat.GetGeometryRef()
    xmin, xmax, ymin, ymax = geom.GetEnvelope()

    # 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)

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

    return numpy.array(zoneraster).tolist()


def loop_zonal_stats(input_value_polygon, input_zonal_raster):
    shp = ogr.Open(input_value_polygon)
    lyr = shp.GetLayer()
    print("Processing", lyr.GetFeatureCount(), "features")
    featList = range(lyr.GetFeatureCount())

    def processFID(input_value_polygon, input_zonal_raster, FID):
        shp = ogr.Open(input_value_polygon)
        raster = gdal.Open(input_zonal_raster)
        lyr = shp.GetLayer()
        if FID:
            px_ids = zonal_stats(FID, lyr, raster)
            # print(px_ids)
            px_ids = [item for sublist in px_ids for item in sublist]
            return px_ids

    return Parallel(n_jobs=8)(
        delayed(processFID)(input_value_polygon, input_zonal_raster, FID)
        for FID in featList
    )


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)

    input_zonal_raster = sys.argv[1]
    input_value_polygon = sys.argv[2]

    counts = list(
        filter(None, loop_zonal_stats(input_value_polygon, input_zonal_raster))
    )
    counts = Counter([item for sublist in counts for item in sublist])

    pandas.DataFrame.from_dict(data=counts, orient="index").to_csv(
        os.path.splitext(input_value_polygon)[0] + ".csv", header=False
    )

这将创建一个与源区域 GeoTiff 具有相同网格系统的输出 GeoTiff

我想知道使用What is the purpose of meshgrid in Python / NumPy?是否可以加快速度

"""
python rasterlayerzonalvectorcounts.py grid.tif liechtenstein_grass.sqlite

MIT License
Based on https://github.com/pcjericks/py-gdalogr-cookbook/blob/master/raster_layers.rst#calculate-zonal-statistics
"""

import sys
import osr
import os
import ogr
import numpy
import gdal
import pandas
from joblib import Parallel, delayed
from collections import Counter
from rich import print, inspect


def zonal_stats(FID, lyr, raster):
    feat = lyr.GetFeature(FID)

    # Get raster georeference info
    transform = raster.GetGeoTransform()
    xOrigin = transform[0]
    yOrigin = transform[3]
    pixelWidth = transform[1]
    pixelHeight = transform[5]

    # Get extent of feat
    geom = feat.GetGeometryRef()
    xmin, xmax, ymin, ymax = geom.GetEnvelope()

    # 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

    feat_arr = []
    # if xcount != 1 or ycount != 1:
    # print(xoff, yoff, xcount, ycount)
    for x in range(xcount):
        for y in range(ycount):
            # print(xoff + x, yoff + y)
            feat_arr.append((xoff + x, yoff + y))

    return feat_arr


def loop_zonal_stats(input_value_polygon, input_zonal_raster):
    shp = ogr.Open(input_value_polygon)
    lyr = shp.GetLayer()
    print("Processing", lyr.GetFeatureCount(), "features")
    featList = range(lyr.GetFeatureCount())

    def processFID(input_value_polygon, input_zonal_raster, FID):
        shp = ogr.Open(input_value_polygon)
        raster = gdal.Open(input_zonal_raster)
        lyr = shp.GetLayer()
        if FID:
            px_ids = zonal_stats(FID, lyr, raster)
            return px_ids

    return Parallel(n_jobs=1)(
        delayed(processFID)(input_value_polygon, input_zonal_raster, FID)
        for FID in featList
    )


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)

    input_zonal_raster = sys.argv[1]
    input_value_polygon = sys.argv[2]

    counts = list(
        filter(None, loop_zonal_stats(input_value_polygon, input_zonal_raster))
    )
    counts = Counter([item for sublist in counts for item in sublist])

    raster_srs = osr.SpatialReference()
    raster_srs.ImportFromWkt(gdal.Open(input_zonal_raster).GetProjectionRef())

    raster_arr = numpy.empty((14045, 20960), numpy.int32)
    for px in counts.items():
        # print(px)
        raster_arr[px[0][1]][px[0][0]] = px[1]

    target_ds = gdal.GetDriverByName("GTiff").Create(
        os.path.splitext(input_value_polygon)[0] + ".tif",
        20960,
        14045,
        1,
        gdal.GDT_Int32,
        options=["COMPRESS=LZW"],
    )
    target_ds.SetGeoTransform(gdal.Open(input_zonal_raster).GetGeoTransform())
    target_ds.SetProjection(raster_srs.ExportToWkt())
    target_ds.GetRasterBand(1).WriteArray(raster_arr)
    target_ds.GetRasterBand(1).SetNoDataValue(0)
    target_ds.GetRasterBand(1).GetStatistics(0, 1)
    target_ds.FlushCache()
    target_ds = None