1

I use the script in this link to reproject a shapefile and save it. The target projection is derived from a georeferenced raster. The problem is while the target coordinate system is UTM 33N (epsg: 32633) the output shapefile has some other coordinate system (I expected to get Central Europe and Africa and I get Australia and New Zealand). I did notice that the output shapefile didn't have a .prj file, just the .shp, .shx and .dbf files. What could be the problem? Here's the script I use:

def reproject_shapefile(coast, epsg, prj):
import os

"""epsg: target srs
prj: target projection"""
driver = ogr.GetDriverByName('ESRI Shapefile')

# get the input layer
inDataSet = driver.Open(coast)
inLayer = inDataSet.GetLayer()

InSpatialReference = inLayer.GetSpatialRef()
epsg_src = InSpatialReference.GetAttrValue('AUTHORITY', 1)
if not epsg == epsg_src:
    # inSpatialRef = osr.SpatialReference()

    # output SpatialReference
    outSpatialRef = osr.SpatialReference(wkt=prj)

    # create the CoordinateTransformation
    coordTrans = osr.CoordinateTransformation(InSpatialReference, outSpatialRef)

    # create the output layer
    outputShapefile = coast.replace('.shp', '_%s.shp' % epsg)
    if os.path.exists(outputShapefile):
        driver.DeleteDataSource(outputShapefile)
    outDataSet = driver.CreateDataSource(outputShapefile)
    outLayer = outDataSet.CreateLayer('coast_%s' % epsg, geom_type=ogr.wkbMultiPolygon)

    # add fields
    inLayerDefn = inLayer.GetLayerDefn()
    for i in range(0, inLayerDefn.GetFieldCount()):
        fieldDefn = inLayerDefn.GetFieldDefn(i)
        outLayer.CreateField(fieldDefn)

    # get the output layer's feature definition
    outLayerDefn = outLayer.GetLayerDefn()

    # loop through the input features
    inFeature = inLayer.GetNextFeature()
    while inFeature:
        # get the input geometry
        geom = inFeature.GetGeometryRef()
        # reproject the geometry
        geom.Transform(coordTrans)
        # create a new feature
        outFeature = ogr.Feature(outLayerDefn)
        # set the geometry and attribute
        outFeature.SetGeometry(geom)
        for i in range(0, outLayerDefn.GetFieldCount()):
            outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
        # add the feature to the shapefile
        outLayer.CreateFeature(outFeature)
        # dereference the features and get the next input feature
        outFeature = None
        inFeature = inLayer.GetNextFeature()

    # Save and close the shapefiles
    inDataSet = None
    outDataSet = None

    coast = outputShapefile

return coast
  • please add in your actual code so we can see what you are doing. – Ian Turton Jun 9 '18 at 15:57
  • I added the code for the script I use – Nikos Jun 10 '18 at 17:21
  • @Luke still no change. – Nikos Jun 10 '18 at 22:19
  • 1
    If your shapefile has no .prj file, you have to specify the input CRS explicitly. Otherwise you get no reprojection. – AndreJ Jun 11 '18 at 6:01
  • The input shapefile does have a .prj file. But even though when I specify the CRS with the *.ImportFromEPSG() function, nothing happens. – Nikos Jun 11 '18 at 6:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.