zonalstats.py 文件源码

def _create_dst_datasource(self, **kwargs):
        """ Create destination target gdal.Dataset

        Creates one layer for each target polygon, consisting of
        the needed source data attributed with index and weights fields

        Returns
        -------
        ds_mem : gdal.Dataset object
        """

        # TODO: kwargs necessary?

        # create intermediate mem dataset
        ds_mem = io.gdal_create_dataset('Memory', 'dst',
                                        gdal_type=gdal.OF_VECTOR)

        # get src geometry layer
        src_lyr = self.src.ds.GetLayerByName('src')
        src_lyr.ResetReading()
        src_lyr.SetSpatialFilter(None)
        geom_type = src_lyr.GetGeomType()

        # create temp Buffer layer (time consuming)
        ds_tmp = io.gdal_create_dataset('Memory', 'tmp',
                                        gdal_type=gdal.OF_VECTOR)
        georef.ogr_copy_layer(self.trg.ds, 0, ds_tmp)
        tmp_trg_lyr = ds_tmp.GetLayer()

        for i in range(tmp_trg_lyr.GetFeatureCount()):
            feat = tmp_trg_lyr.GetFeature(i)
            feat.SetGeometryDirectly(feat.GetGeometryRef().
                                     Buffer(self._buffer))
            tmp_trg_lyr.SetFeature(feat)

        # get target layer, iterate over polygons and calculate intersections
        tmp_trg_lyr.ResetReading()

        self.tmp_lyr = georef.ogr_create_layer(ds_mem, 'dst', srs=self._srs,
                                               geom_type=geom_type)

        print("Calculate Intersection source/target-layers")
        try:
            tmp_trg_lyr.Intersection(src_lyr, self.tmp_lyr,
                                     options=['SKIP_FAILURES=YES',
                                              'INPUT_PREFIX=trg_',
                                              'METHOD_PREFIX=src_',
                                              'PROMOTE_TO_MULTI=YES',
                                              'PRETEST_CONTAINMENT=YES'],
                                     callback=gdal.TermProgress)
        except RuntimeError:
            # Catch RuntimeError that was reported on gdal 1.11.1
            # on Windows systems
            tmp_trg_lyr.Intersection(src_lyr, self.tmp_lyr,
                                     options=['SKIP_FAILURES=YES',
                                              'INPUT_PREFIX=trg_',
                                              'METHOD_PREFIX=src_',
                                              'PROMOTE_TO_MULTI=YES',
                                              'PRETEST_CONTAINMENT=YES'])

        return ds_mem