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)

        feat_cnt = layer.GetFeatureCount()

        if feat_cnt:
            [georef.ogr_add_geometry(dst, ogr_src.GetGeometryRef(),
                                     [ogr_src.GetField('index'), trg_index])
             for ogr_src in layer]
        else:
            layer.SetSpatialFilter(None)
            src_pts = np.array([ogr_src.GetGeometryRef().GetPoint_2D(0)
                                for ogr_src in layer])
            centroid = georef.get_centroid(trg)
            tree = cKDTree(src_pts)
            distnext, ixnext = tree.query([centroid[0], centroid[1]], k=1)
            feat = layer.GetFeature(ixnext)
            georef.ogr_add_geometry(dst, feat.GetGeometryRef(),
                                    [feat.GetField('index'), trg_index])