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

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

或者：寻找一种更快、更准确的方法将小型 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

所以我在想，如果我们仍然可以对矢量特征进行迭代（因为这些特征要少得多，而且栅格网格中有 8800 多万个区域），它的性能可能会更高。

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

然而，当试图让这个脚本工作时，它一直给我相同的几何形状......它只一遍又一遍地向我显示像素 ID 中的 1 个

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,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,ymax,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],burn_values=[1])

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

    bandmask = target_ds.GetRasterBand(1)
    datamask = bandmask.ReadAsArray(0,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,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

解决方法

这似乎有效。输出是 CSV 列 =["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,ymax,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],burn_values=[1])

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

    bandmask = target_ds.GetRasterBand(1)
    datamask = bandmask.ReadAsArray(0,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,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,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,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,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,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,raster)
            return px_ids

    return Parallel(n_jobs=1)(
        delayed(processFID)(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,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