Reproyectando y guardando un shapefile tipo polígono con GDAL/OGR

El codigo siguiente corrige los problemas de ejecución que presenta este otro para reproyectar un shapefile con base en el sistema de referencia espacial que es tomado de un ráster.

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

from osgeo import ogr, osr, gdal #tif with projections I want tif = gdal.Open("/home/zeito/pyqgis_data/utah_dem4326.tif") #shapefile with the from projection driver = ogr.GetDriverByName("ESRI Shapefile") dataSource = driver.Open("/home/zeito/pyqgis_data/polygon8.shp", 1) layer = dataSource.GetLayer() #set spatial reference and transformation sourceprj = layer.GetSpatialRef() targetprj = osr.SpatialReference(wkt = tif.GetProjection()) transform = osr.CoordinateTransformation(sourceprj, targetprj) to_fill = ogr.GetDriverByName("Esri Shapefile") ds = to_fill.CreateDataSource("/home/zeito/pyqgis_data/projected.shp") outlayer = ds.CreateLayer('', targetprj, ogr.wkbPolygon) outlayer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger)) #apply transformation i = 0 for feature in layer: transformed = feature.GetGeometryRef() transformed.Transform(transform) geom = ogr.CreateGeometryFromWkb(transformed.ExportToWkb()) defn = outlayer.GetLayerDefn() feat = ogr.Feature(defn) feat.SetField('id', i) feat.SetGeometry(geom) outlayer.CreateFeature(feat) i += 1 feat = None ds = None

Después de la ejecución del código anterior se obtiene el vectorial reproyectado que se exhibe en la imagen siguiente. La tabla atributiva refleja el llenado adecuado del campo id; algo que no haría el código objeto de corrección.