zonalstats.py 文件源码

def dump_raster(self, filename, driver='GTiff', attr=None,
                    pixel_size=1., remove=True):
        """ Output layer to GDAL Rasterfile

        Parameters
        ----------
        filename : string
            path to shape-filename
        driver : string
            GDAL Raster Driver
        attr : string
            attribute to burn into raster
        pixel_size : float
            pixel Size in source units
        remove : bool
            if True removes existing output file

        """
        layer = self.ds.GetLayer()
        layer.ResetReading()

        x_min, x_max, y_min, y_max = layer.GetExtent()

        cols = int((x_max - x_min) / pixel_size)
        rows = int((y_max - y_min) / pixel_size)

        # Todo: at the moment, always writing floats
        ds_out = io.gdal_create_dataset('MEM', '', cols, rows, 1,
                                        gdal_type=gdal.GDT_Float32)

        ds_out.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
        proj = layer.GetSpatialRef()
        if proj is None:
            proj = self._srs
        ds_out.SetProjection(proj.ExportToWkt())

        band = ds_out.GetRasterBand(1)
        band.FlushCache()
        print("Rasterize layers")
        if attr is not None:
            gdal.RasterizeLayer(ds_out, [1], layer, burn_values=[0],
                                options=["ATTRIBUTE={0}".format(attr),
                                         "ALL_TOUCHED=TRUE"],
                                callback=gdal.TermProgress)
        else:
            gdal.RasterizeLayer(ds_out, [1], layer, burn_values=[1],
                                options=["ALL_TOUCHED=TRUE"],
                                callback=gdal.TermProgress)

        io.write_raster_dataset(filename, ds_out, driver, remove=remove)

        del ds_out