1

I'm trying to learn the use of PySal and Python. As I'm already regular user of PostGIS, I would like to integrate the database as the storage medium for the data to be analysed. I'm using the test data on Saint Louis homicides loaded into PostGIS for testing

                                                              Table "pysal.stl_hom"
   Column   |             Type             |                          Modifiers                          | Storage  | Stats target | Description 
------------+------------------------------+-------------------------------------------------------------+----------+--------------+-------------
 gid        | integer                      | not null default nextval('pysal.stl_hom_gid_seq'::regclass) | plain    |              | 
 poly_id_og | numeric(10,0)                |                                                             | main     |              | 
 name       | character varying(32)        |                                                             | extended |              | 
 state_name | character varying(25)        |                                                             | extended |              | 
 state_fips | character varying(2)         |                                                             | extended |              | 
 cnty_fips  | character varying(3)         |                                                             | extended |              | 
 fips       | character varying(5)         |                                                             | extended |              | 
 fipsno     | double precision             |                                                             | plain    |              | 
 hr7984     | double precision             |                                                             | plain    |              | 
 hr8488     | numeric                      |                                                             | main     |              | 
 hr8893     | double precision             |                                                             | plain    |              | 
 hc7984     | integer                      |                                                             | plain    |              | 
 hc8488     | integer                      |                                                             | plain    |              | 
 hc8893     | integer                      |                                                             | plain    |              | 
 po7984     | integer                      |                                                             | plain    |              | 
 po8488     | integer                      |                                                             | plain    |              | 
 po8893     | integer                      |                                                             | plain    |              | 
 pe77       | double precision             |                                                             | plain    |              | 
 pe82       | double precision             |                                                             | plain    |              | 
 pe87       | double precision             |                                                             | plain    |              | 
 rdac80     | double precision             |                                                             | plain    |              | 
 rdac85     | double precision             |                                                             | plain    |              | 
 rdac90     | double precision             |                                                             | plain    |              | 
 geom       | geometry(MultiPolygon,26913) |                                                             | main     |              | 

I'm still very much a rookie when it comes to Python and I can't seem to find many clues as to how I can use the PySal module with postGIS data. But it seems that shapely can be used to create objects that can be manipulated by PySal.

The method I've tried in the interactive Python interpreter is the following, where I've left out the the connection part and cursor (I can print the rows cursor object):

import shapely.geometry psycopg2 
sql = """select gid,st_asewkb(geom) as geom from pysal.stl_hom;"""
rows = conn.execute(sql)
sh = shapely.geometry.base.geom_from_wkb(row['geom'])

The result is the error below, which I can't figure out how to interpret:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python2.7/dist-packages/shapely/geometry/base.py", line 116, in geom_from_wkb
    return geom_factory(deserialize_wkb(data))
  File "/usr/local/lib/python2.7/dist-packages/shapely/geometry/base.py", line 106, in deserialize_wkb
    geom = lgeos.GEOSGeomFromWKB_buf(c_char_p(data), c_size_t(len(data)))
TypeError: string or integer address expected instead of buffer instance

So my questions are

  1. Can I use Shapely to create an object from PostGIS usable by PySal and how can I write this code
  2. Are there other viable approaches to use PySal with PostGIS without transforming the data into shape files first
  • Your question is on the broad side, there being at least two of them. However, to answer one of them, in the input side, the docs show that you can import data as WKT. So, you could use ST_AsWKT in the above query, to avoid having to read/write from shapes. GeoJSON is another possibility. – John Powell Jul 6 '16 at 18:01
  • I see what you mean, but the docs only show functions that can be used with actual files (and I don't see any geojson examples) – JonasPedersen Jul 6 '16 at 19:10
2

The easiest solution is to get a dictionary from query (psycopg2.extras: Real dictionary cursor or Dictionary-like cursor ) and the GeoJson format (ST_AsGeoJSON)

import psycopg2
conn = psycopg2.connect("dbname='strati 'host='localhost' user='me'")
import psycopg2.extras
dict_cur = conn.cursor(cursor_factory=psycopg2.extras.RealDictCursor)
dict_cur.execute("""SELECT "DIP_DIR","DIP", ST_AsGeoJSON(the_geom) from teststrati;""")
# first element
rec = dict_cur.fetchone()
# now rec is a Python dictionary
print rec 
{'DIP_DIR': 130, 'DIP': 30, 'st_asgeojson': '{"type":"Point","coordinates":[272070.600040999997873,155389.387920000008307]}'}
rec.keys()
['DIP_DIR', 'DIP', 'st_asgeojson']
print rec['DIP_DIR'], rec['DIP']
130 30
# geometry
geom = rec['st_asgeojson']
print geom
{"type":"Point","coordinates":[272070.600040999997873,155389.387920000008307]}

Now convert to a shapely geometry

# here geom is a string and shapely needs a dictionary
import json
from shapely.geometry import shape
shape(json.loads(geom))
<shapely.geometry.point.Point object at 0x114c53550>
print shape(json.loads(geom))
POINT (272070.600041 155389.38792)

You can directly use PySAL: Using PySAL with Shapely for GIS Operations

  • That definitely answered my question. Now I just need to incorporate the other attributes, to be able to calculate stuff like morans I. That is of course a completely separate question. As a little help to others looking at this I transformed the Shapely object into a PySal object using shppy = shape(json.loads(rec[0]['geometry']))and pysalgeom = pysal.cg.asShape(shppy) – JonasPedersen Jul 7 '16 at 12:11
3

The quickest way is to keep the geometry in hex-encoded WKB.

import psycopg2
from shapely import wkb

conn = psycopg2.connect('...')
curs = conn.cursor()

shps = {}  # key: gid, value: Shapely geom
curs.execute('select gid, geom as geom from pysal.stl_hom;')
for gid, geom in curs:
    shps[gid] = wkb.loads(geom, hex=True)

There's no need to use GeoJSON, which is slower and lossy compared to WKB.

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.