I have a shapefile I want to reproject based on another projection, and then save the results.
I am trying to do this like so:
from osgeo import ogr, osr, gdal
#tif with projections I want
tif = gdal.Open("C:/path_to_tif/file.tif")
#shapefile with the from projection
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open("C:/path_to_shape/file.shp", 1)
layer = dataSource.GetLayer()
#set spatial reference and transformation
sourceprj = layer.GetSpatialRef()
targetprj = osr.SpatialReference(wkt = tif.GetProjection())
transform = osr.CoordinateTransformation(sourceprj, targetprj)
#apply transformation
for feature in layer:
transformed = feature.GetGeometryRef()
transformed.Transform(transformed)
#save shapefile
to_fill = ogr.GetDriverByName("Esri Shapefile")
ds = to_fill.CreateDataSource("C:/path_to_projected/projected.shp")
outlayer = ds.CreateLayer('', None, ogr.wkbPolygon)
outlayer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
defn = layer.GetLayerDefn()
feat = ogr.Feature(defn)
feat.SetField('id', 123)
geom = ogr.CreateGeometryFromWkb(transformed.ExportToWkb())
feat.SetGeometry(geom)
outlayer.CreateFeature(feat)
but the resulting shapefile is empty and feat.SetGeometry(geom)
returns 0
Answer
Following code fixes your issues (compare with your code) and works adequately. I tried it out with my layers at referred paths (change for your paths).
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
At following image it can be observed reprojected shapefile (EPSG:4326 from EPSG:32612) with its attributes table. It works as expected.
No comments:
Post a Comment