zonalstats.py 文件源码

def _create_dst_features(self, dst, trg, **kwargs):
        """ Create needed OGR.Features in dst OGR.Layer

        Parameters
        ----------
        dst : OGR.Layer
            destination layer
        trg : OGR.Geometry
            target polygon
        """
        # TODO: kwargs necessary?

        # claim and reset source ogr layer
        layer = self.src.ds.GetLayerByName('src')
        layer.ResetReading()

        # if given, we apply a buffer value to the target polygon filter
        trg_index = trg.GetField('index')
        trg = trg.GetGeometryRef()
        trg = trg.Buffer(self._buffer)
        layer.SetSpatialFilter(trg)

        # iterate over layer features
        for ogr_src in layer:
            geom = ogr_src.GetGeometryRef()

            # calculate intersection, if not fully contained
            if not trg.Contains(geom):
                geom = trg.Intersection(geom)

            # checking GeometryCollection, convert to only Polygons,
            #  Multipolygons
            if geom.GetGeometryType() in [7]:
                geocol = georef.ogr_geocol_to_numpy(geom)
                geom = georef.numpy_to_ogr(geocol, 'MultiPolygon')

            # only geometries containing points
            if geom.IsEmpty():
                continue

            if geom.GetGeometryType() in [3, 6, 12]:
                idx = ogr_src.GetField('index')
                georef.ogr_add_geometry(dst, geom, [idx, trg_index])