0

I tried to reproject a shapefile based on the projection of a raster file following this example (section "Reproject a Layer").

My code is as follows:

    # open the raster tif and extract crs
    t = "path/to/raster.tif"
    tif = gdal.Open(t)
    crs1 = tif.GetSpatialRef()#.GetProjection()
    tif = None
    # open shapefile
    drv = ogr.GetDriverByName("ESRI Shapefile")
    s = "path/to/shapefile.shp"
    # set name of the new shapefile
    sUTM = "shapefile_UTM.shp"
    LayerName = "shapefile_UTM"
    if os.path.exists(sUTM):
        drv.DeleteDataSource(sUTM)
    shp = drv.Open(s)
    layer = shp.GetLayer()
    crs0 = layer.GetSpatialRef()
    # transformation from crs0 to crs1
    coordTrans = osr.CoordinateTransformation(crs0, crs1)
    # create output shapefile
    outSHP = drv.CreateDataSource(sUTM)
    outLayer = outSHP.CreateLayer(LayerName,
                                  geom_type = ogr.wkbMultiPolygon)
    lyrdef = layer.GetLayerDefn()
    for i in range(0, lyrdef.GetFieldCount()):
        fieldDefn = lyrdef.GetFieldDefn(i)
        outLayer.CreateField(fieldDefn)
    outdef = outLayer.GetLayerDefn()
    feature = layer.GetNextFeature()
    while feature:
        geom = feature.GetGeometryRef()
        geom.Transform(coordTrans)
        outFeature = ogr.Feature(outdef)
        outFeature.SetGeometry(geom)
        for i in range(0, outdef.GetFieldCount()):
            outFeature.SetField(outdef.GetFieldDefn(i).GetNameRef(),
                                feature.GetField(i))
        outLayer.CreateFeature(outFeature)
        outFeature = None
        feature = layer.GetNextFeature()
    layer = None
    shp = None
    outSHP = None

However, when I run this code, it seems to run into some issues, because it ends with Restarting kernel... and all variables are gone.

Moreover, the output shapefile has only .dbf, .shp and .shx files. There is no .prj file. Of course, I could write this file in an additional step, but I think there is some underlying issue somewhere which I would prefer to resolve instead of fixing the outcome.

What is the error in the above shown code?

Edit:

The reference system I want to transform the shapefile to—i.e., what print(crs1) gives—is:

PROJCS["WGS 84 / UTM zone 35S",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",27],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32735"]]
2
  • Please Edit the question to specify the CRS of the raster. It may not be supported.
    – Vince
    May 24, 2021 at 12:22
  • @Vince I want to have the Shapefile in UTM zone 35S. I added the full info that was stored in the crs1 variable to the question. May 24, 2021 at 12:29

1 Answer 1

2

You can greatly simplify the operation using gdal.VectorTranslate

from osgeo import gdal

# Open the raster tif and extract crs
t = "path/to/raster.tif"
tif = gdal.Open(t)
crs1 = tif.GetProjectionRef()
tif = None

# Input shp path
s = "path/to/shapefile.shp"

# Set name of the new shapefile
sUTM = "shapefile_UTM.shp"
if os.path.exists(sUTM):
    drv.DeleteDataSource(sUTM)

# Open shp input, set destination crs output to
# raster crs and write to output shp name
srcDS = gdal.OpenEx(s)
ds = gdal.VectorTranslate(sUTM, srcDS, format='ESRI Shapefile', dstSRS=crs1, reproject=True)

ds = None
srcDS = None

# Write prj
with open(f"{os.path.splitext(sUTM)[0]}.prj", 'w') as f:
    f.write(crs1)

2
  • Thanks a lot, this solved the problem and additionally, made the script shorter and easier to comprehend. I wonder why they have this cumbersome example on the Python GDAL/OGR Cookbook documentation. Perhaps it is outdated. May 28, 2021 at 14:41
  • 1
    Exactly. Although most parts of the website from your sample remain relevant, it has been published on a code repository not updated since 5 years github.com/pcjericks/py-gdalogr-cookbook and for most content since 7/8 years
    – ThomasG77
    May 28, 2021 at 14:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.