Monday, 6 March 2017

python - Reprojecting and saving shapefile in gdal


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.


enter image description here


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...