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
python类TermProgress()的实例源码
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
