0

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?

There is some code here, the python file might be better than the Jupyter notebook.

  • More information added. – william3031 Oct 17 at 3:53
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

  • 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 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 Wang Oct 16 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 at 23:13
  • Can you add your code? – Zac Wang Oct 17 at 0:18
  • 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 Wang Oct 17 at 1:52

Your Answer

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

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