11

I am doing a small-scale workflow in ipython notebook with geopandas and shapely, pulling a bunch of geospatial data, sometimes from shapefiles, sometimes from postgis (where more expensive processing is performed).

Now, I then pull the postgis tables into python using sqlalchemy, transforming the geometry to WKT along the way, getting something like:

sql = """
SELECT ST_AsText(ST_Transform(the_geom,4326)) as newgeom,* 
  FROM public.parcels2010_small limit 5;
 """
parcels = pd.read_sql(sql, engine)
parcels

+----+---------------------------------------------------+---------------------------------------------------+-----------+-------------+
| id | newgeom                                           | the_geom                                          | parcel_id | osm_node_id |
+----+---------------------------------------------------+---------------------------------------------------+-----------+-------------+
| 0  | MULTIPOLYGON(((-122.991093691444 38.4878691106... | 01060000209C0E00000100000001030000000100000097... | 1805792   | 66237       |
+----+---------------------------------------------------+---------------------------------------------------+-----------+-------------+
| 1  | MULTIPOLYGON(((-122.444576448624 37.7346386006... | 01060000209C0E0000010000000103000000010000008A... | 10435     | 123826      |
+----+---------------------------------------------------+---------------------------------------------------+-----------+-------------+
| 2  | MULTIPOLYGON(((-122.796785208193 38.5427593334... | 01060000209C0E0000010000000103000000010000007D... | 1817842   | 313047      |
+----+---------------------------------------------------+---------------------------------------------------+-----------+-------------+
| 3  | MULTIPOLYGON(((-122.695538506163 38.3618570798... | 01060000209C0E0000010000000103000000010000009B... | 1934612   | 63776       |
+----+---------------------------------------------------+---------------------------------------------------+-----------+-------------+
| 4  | MULTIPOLYGON(((-122.223424422869 37.8416019090... | 01060000209C0E00000100000001030000000100000072... | 861785    | 26369       |
+----+---------------------------------------------------+---------------------------------------------------+-----------+-------------+

This looks a lot like the geometry when loaded from a shapefile, but it is not cast as a shapely geometry. I couldn't find the canonical way to do it, either using shapely alone, or perhaps with Descartes.

8

GeoAlchemy understands SQLAlchemy PostGIS and Shapely.

from geoalchemy2.shape import to_shape

for parcel in parcels:
    parcel_shape = to_shape(parcel.the_geom)
12

GeoAlchemy may be the better route in this case, but here's a solution that answers the question more directly and doesn't require any packages beyond what is described in the question:

from shapely import wkb

# ....

sql = """SELECT * FROM public.parcels2010_small LIMIT 5;"""
parcels = pd.read_sql(sql, engine)

for parcel in parcels:
    parcel.the_geom = wkb.loads(parcel.the_geom, hex=True)

if you were to then print the geometry it should look something like this:

print parcels[0].the_geom

<shapely.geometry.multipolygon.MultiPolygon object at ...>

See docs on the shapely.wkb module here.

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.