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?

Zac · Accepted Answer · 2019-10-17T11:33:59.840

2

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

edited Oct 17 '19 at 11:33

answered Oct 16 '19 at 11:36

Zac

517
3
10

Thanks for the response. I couldn't get it to work though. Is this standalone? I could only get it to work by replacing route with isochrone_polys[0] but I imagine that isn't correct and produces an empty shapefile anyway. Adding it to the original code (replacing the other bits) also didn't seem to work well. – william3031 Oct 16 '19 at 21:24
1

The route in lyr = ds.CreateLayer('route', sr, ogr.wkbPolygon) is just a new layer named in Shapefile, you can name it whatever you want. And this method works with the GDAL/OGR, which you don't need to use Fiona. The basic idea is to create a new shapefile with the CRS you want and add shapely geometry to each feature in the shapefile by geom = ogr.CreateGeometryFromWkb(isochrone_polys[0].wkb) following the link you provided. – Zac Oct 16 '19 at 22:14
I have tried that and still can't get it to work. It exports the abc.shp that is empty when I open it in QGIS. – william3031 Oct 16 '19 at 23:13
Can you add your code? – Zac Oct 17 '19 at 00:18
In here: https://github.com/william3031/osmnx_test Thanks. – william3031 Oct 17 '19 at 00:47
I have added a python script which has pretty much the same information as the notebook. Thanks for looking at this. – william3031 Oct 17 '19 at 01:30
1

I checked your code. First, you should delete the parts of test1 and test2, just use test3. Because I am not familiar with networkx, therefore I created a polygon in shapely to simulate the isochrone_polys[0]. And it works, check my update in my answer where I added the code. – Zac Oct 17 '19 at 01:52
Thanks for checking. Yes, it does work with Polygon([(0, 0), (1, 100), (100, 0)]) as the coordinates fit in to a WGS84 structure but it doesn't with isochrone_polys[0]. The coordinates of isochrone_polys[0] seem to be projected, but not sure which projection. However it gets picked up by default in QGIS as WGS84, then is in the wrong place. I have updated the question with this information. – william3031 Oct 17 '19 at 03:03
Thanks for the update. I haven't used the data from the bbox section in this. Everything is run from the place which is just a centroid of the suburb. Whlie the edges has a CRS, the polygon doesn't. AttributeError: 'Polygon' object has no attribute 'crs' – william3031 Oct 17 '19 at 03:45
Shapely geometry doesn't include CRS. +proj=utm +zone=55 +ellps=WGS84 +datum=WGS84 +units=m +no_defs is a proj4 format CRS, so when you create shapefile, change 4326 to 32755, and check the result. – Zac Oct 17 '19 at 04:04
I exported it with the crs of 32755. When I imported it in to QGIS, it didn't pick it up and defaulted in to WGS84. I manually changed it in to 32755, but that isn't the right one either. – william3031 Oct 17 '19 at 04:20
Could you export isochrone_polys[0] as a wkt format to a txt file, so I may try to import it to a shapefile and check whether it can work or not? isochrone_polys[0].wkt – Zac Oct 17 '19 at 05:24
I have added those to the github repo. – william3031 Oct 17 '19 at 06:17
@william3031, check the new update code and screenshot. Hope it works this time. – Zac Oct 17 '19 at 10:36
Yes, that works. Thanks heaps for the work you've put in to answering the question. – william3031 Oct 17 '19 at 10:52

Getting CRS to export shapely polygon to shapefile

1 Answers1