Rasterizando con el modulo GDAL/OGR de Python

Con el módulo GDAL/OGR es posible acceder a una herramienta para rasterizar capas vectoriales y usarlas como máscaras ráster. El código que se facilita a continuación permite realizar este proceso dándole al rasgo rasterizado el valor constante de uno.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

from osgeo import gdal, ogr vector_layer = "/home/zeito/pyqgis_data/cut_polygon3.shp" raster_layer = "/home/zeito/pyqgis_data/aleatorio.tif" target_layer = "/home/zeito/pyqgis_data/mask.tif" # open the raster layer and get its relevant properties raster_ds = gdal.Open(raster_layer, gdal.GA_ReadOnly) xSize = raster_ds.RasterXSize ySize = raster_ds.RasterYSize geotransform = raster_ds.GetGeoTransform() projection = raster_ds.GetProjection() # create the target layer (1 band) driver = gdal.GetDriverByName('GTiff') target_ds = driver.Create(target_layer, xSize, ySize, 1, gdal.GDT_Byte) target_ds.SetGeoTransform(geotransform) target_ds.SetProjection(projection) source_ds = ogr.Open(vector_layer) source_layer = source_ds.GetLayer() # rasterize the vector layer into the target one ds = gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[1]) target_ds = None

El código anterior fue usado con las capas vectorial y ráster que se presentan a continuación:

y después de su ejecución se obtuvo la imagen rasterizada siguiente:

Sin embargo, cuando se tiene más de un rasgo a ser rasterizado y cada uno de ellos tiene que tener un valor diferencial distinto, es necesario asignarlo a un campo de la tabla atributiva del vectorial. Es lo que realiza el codigo siguiente considerando el campo 'value' como atributo.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

from osgeo import gdal, ogr vector_layer = "/home/zeito/pyqgis_data/cut_polygon4.shp" raster_layer = "/home/zeito/pyqgis_data/aleatorio.tif" target_layer = "/home/zeito/pyqgis_data/mask.tif" # open the raster layer and get its relevant properties raster_ds = gdal.Open(raster_layer, gdal.GA_ReadOnly) xSize = raster_ds.RasterXSize ySize = raster_ds.RasterYSize geotransform = raster_ds.GetGeoTransform() projection = raster_ds.GetProjection() # create the target layer (1 band) driver = gdal.GetDriverByName('GTiff') target_ds = driver.Create(target_layer, xSize, ySize, 1, gdal.GDT_Byte) target_ds.SetGeoTransform(geotransform) target_ds.SetProjection(projection) source_ds = ogr.Open(vector_layer) source_layer = source_ds.GetLayer() # rasterize the vector layer into the target one ds = gdal.RasterizeLayer(target_ds, [1], source_layer, options = ["ATTRIBUTE=value"]) target_ds = None

El código anterior se probó con las layers que se observan a continuación:

y después de su ejecución se obtuvo el ráster de la imagen siguiente. Los valores para cada rasgo son los esperados según el campo value de la tabla atributiva del vectorial.