Reprojecting and saving shapefile in gdal

Question

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 1

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.