Getting CRS to export shapely polygon to shapefile

Question

I am trying to export a shapely polygon to a shapefile. I have tried this, but the CRS is missing, when I try to import it in to QGIS.

The shapely polygon is from this OSMNX example but edited to work with location.

My (list of two) polygons:

In [68]:

isochrone_polys


Out[68]:
[<shapely.geometry.polygon.Polygon at 0x1d884a52cc0>,
 <shapely.geometry.polygon.Polygon at 0x1d883974f60>]

I tried this using Fiona:

from pyproj import Proj, transform
import fiona
from fiona.crs import from_epsg
Proj(isochrone_polys[0].crs)

But have an attribute error:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-100-f3e93644040a> in <module>
      2 import fiona
      3 from fiona.crs import from_epsg
----> 4 Proj(isochrone_polys[0].crs)

AttributeError: 'Polygon' object has no attribute 'crs'

The edges in and nodes in the street network in the tutorial seem to be in 4326, but when I put this in as the CRS in QGIS, it was incorrect.

edges is a geopandas dataframe, while the isochrone_polys[1] is in shapely.

type(edges)
Out[11]: geopandas.geodataframe.GeoDataFrame

type(isochrone_polys[1])
Out[12]: shapely.geometry.polygon.Polygon

print(edges.crs)
+proj=utm +zone=55 +ellps=WGS84 +datum=WGS84 +units=m +no_defs

print(isochrone_polys[1].crs)
Traceback (most recent call last):

  File "<ipython-input-16-e4c06ab91976>", line 1, in <module>
    print(isochrone_polys[1].crs)

AttributeError: 'Polygon' object has no attribute 'crs'

The polygon itself seems to be in a projected coordinate system (not sure which one):

isochrone_polys[0].bounds
Out[7]: 
(316180.5323058479,
 -4182505.539783961,
 317463.60173865134,
 -4181214.9721103134)

How would I export this properly with the correct CRS?

Answer 1

The link that you tried didn't set CRS so when you import the shapefile into QGIS. It's incorrect.

If you know the CRS is EPSG 4326, then you can set it while creating the shapefile like this:

from osgeo import osr, ogr # import osr, ogr 
sr = osr.SpatialReference()
sr.ImportFromEPSG(4326)
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.CreateDataSource('abc.shp')
lyr = ds.CreateLayer('route', sr, ogr.wkbPolygon)

NEW UPDATE:

from osgeo import ogr, osr
from shapely.geometry import Polygon
import shapely.wkt

# read wkt and use it to create a shapely polygon P
with open('/Users/**/osmnx_test-master/isochrone_polys_0.txt') as f:
    x = f.readline()
P = shapely.wkt.loads(x)

# Now convert it to a shapefile with OGR    
sr = osr.SpatialReference()
# or sr.ImportFromProj4(edges.crs)
sr.ImportFromProj4('+proj=utm +zone=55 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.CreateDataSource('abc.shp')
layer = ds.CreateLayer('route', sr, ogr.wkbPolygon)# Add one attribute
layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
defn = layer.GetLayerDefn()

## If there are multiple geometries, put the "for" loop here

# Create a new feature (attribute and geometry)
feat = ogr.Feature(defn)
feat.SetField('id', 123)

# Make a geometry, from Shapely object
geom = ogr.CreateGeometryFromWkb(P.wkb)
feat.SetGeometry(geom)
layer.CreateFeature(feat)
feat =  None 

# Save and close the data source
ds = None