Cómo rasterizar rasgos individuales de un shapefile tipo polígono con GDAL/OGR en Python

En este post se incluye un procedimiento para rasterizar rasgos individuales de un shapefile tipo polígono mediante los módulos Python de GDAL/OGR. La rutina introduce una proyección, tanto en las capas vectoriales como en las ráster, e incluye la asignación de los nodata a cero en las capas ráster.

El código completo es el siguiente:

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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98

from osgeo import gdal from osgeo import ogr, osr from osgeo import gdalconst import os raster_file = "/home/zeito/pyqgis_data/mask.tif" shape_file = "/home/zeito/pyqgis_data/test2.shp" # open raster and obtain extent and properties data = gdal.Open(raster_file, gdalconst.GA_ReadOnly) wkt = data.GetProjection() srs = osr.SpatialReference() srs.ImportFromWkt(wkt) spatialRef = osr.SpatialReference() spatialRef.ImportFromWkt(wkt) spatialRef.MorphToESRI() geo_transform = data.GetGeoTransform() x_min = geo_transform[0] y_max = geo_transform[3] x_max = x_min + geo_transform[1] * data.RasterXSize y_min = y_max + geo_transform[5] * data.RasterYSize x_res = data.RasterXSize y_res = data.RasterYSize pixel_width = geo_transform[1] # open shapefile and create layer driver = ogr.GetDriverByName("ESRI Shapefile") dataSource = driver.Open(shape_file, 0) layer = dataSource.GetLayer() nr_features = layer.GetFeatureCount() for feat_id in range(nr_features): feature = layer.GetFeature(feat_id) id_no = feature.GetField("id_no") presence = feature.GetField("presence") outShapefile = os.path.join("/home/zeito/pyqgis_data/", "feature" + str(feat_id)+".shp") outPrj = os.path.join("/home/zeito/pyqgis_data/", "feature" + str(feat_id)+".prj") output = os.path.join("/home/zeito/pyqgis_data/", "rasterized" + str(feat_id)+".tif") outDriver = ogr.GetDriverByName("ESRI Shapefile") # Remove output shapefile if it already exists if os.path.exists(outShapefile): outDriver.DeleteDataSource(outShapefile) # Create the output shapefile outDataSource = outDriver.CreateDataSource(outShapefile) outLayer = outDataSource.CreateLayer("feature" + str(feat_id), geom_type=ogr.wkbPolygon) file = open(outPrj, 'w') file.write(spatialRef.ExportToWkt()) file.close() # Add an ID field idField = ogr.FieldDefn("id_no", ogr.OFTInteger) outLayer.CreateField(idField) # Add a presence field presField = ogr.FieldDefn("presence", ogr.OFTInteger) outLayer.CreateField(presField) # Create the feature and set values featureDefn = outLayer.GetLayerDefn() feat = ogr.Feature(featureDefn) geom = feature.GetGeometryRef() feat.SetGeometry(geom) feat.SetField("id_no", id_no) feat.SetField("presence", presence) outLayer.CreateFeature(feat) feat = None # Save and close DataSource outDataSource = None target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Byte) target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width)) new_dataSource = driver.Open(outShapefile, 0) new_layer = new_dataSource.GetLayer() # setting nodata values of output raster band = target_ds.GetRasterBand(1) band.SetNoDataValue(0) # setting spatial reference of output raster target_ds.SetProjection( srs.ExportToWkt() ) # this works now gdal.RasterizeLayer(target_ds, [1], new_layer, options=["ATTRIBUTE=presence"]) target_ds = None new_dataSource = None dataSource = None

La situación inicial luce como sigue:

Después de ejecutado el script en la Python Console de QGIS 3.4, así se observan los rasgos individualizados al ser cargados en la Map View; señalando que la individualización de rasgos fue ejecutada correctamente.

Los rasgos (features) rasterizados se observan a continuación y son tal como se esperaban.