1

I'm generating a shapefile using the GDAL and shapely libraries in Python. There's a problem with the projection, when I import the resulting shapefile in ArcMap, there is not projection associated with the file. How can I correct that, or how can I project the resulting layer ?

Here's my code:

output_shp = "circle.shp"
data_projection = 26919

[...]

srs = osr.SpatialReference()
srs.ImportFromEPSG(data_projection)

# output of srs: <osgeo.osr.SpatialReference; proxy of <Swig Object of type 'OSRSpatialReferenceShadow *' at 0x00000000099FF5A0> >

driver = ogr.GetDriverByName('Esri Shapefile')
ds = driver.CreateDataSource(output_shp)
layer = ds.CreateLayer(layer_name, srs, ogr.wkbPolygon)
layer.CreateField(ogr.FieldDefn('id',ogr.OFTInteger))
defn = layer.GetLayerDefn()  

for x in range(0,length):
    feat = ogr.Feature(defn)
    feat.SetField('id', x+1)

    geome = ogr.CreateGeometryFromWkt(geom_out[x].wkt)
    feat.SetGeometry(geome)   
    layer.CreateFeature(feat)
kaycee
  • 901
  • 1
  • 9
  • 35

1 Answers1

0

Ok. I've found the solution, thanks to MikeT answer in this post.

The problem was that my GDAL_DATA environment variable was not correctly set. Therefore, the SRS I was adding to my layer was empty because it couldn't find it in the directory.

A little checkup to see if the importation fails (thanks to MikeT):

from osgeo import osr
testSR = osr.SpatialReference()
res = testSR.ImportFromEPSG(4326)
if res != 0:
    raise RuntimeError(repr(res) + ': could not import from EPSG')
print testSR.ExportToPrettyWkt()
Community
  • 1
  • 1
kaycee
  • 901
  • 1
  • 9
  • 35