Geoprocessing with Geodatabases

Daniel Fewell

In [48]:
import psycopg2
import requests
from shapely.geometry import Point,Polygon,MultiPolygon, mapping
import datetime

from shapely.wkb import loads
from shapely.wkt import dumps, loads
import json
import folium

from shapely.geometry import mapping

import ipywidgets as widgets

Importing libraries.

In [49]:
connection = psycopg2.connect(database="myspatialdb",user="useroh",password="geog482")
cursor = connection.cursor()

Establishing connection to DB and creating cursor to read data in DB.

In [50]:
center = [35.106196,-106.629515]
zoom = 10
p = Point([-106.578677,35.062485])

Setting up map options and creating point for buffer and finding incidents.

Finding the Number of Incidents in the Past 150 Days within a 500 Meter Buffer from the Old Town (-106.6717236, 35.0959965)

In [51]:
map_q1 = folium.Map(location=center, zoom_start=zoom)

Creating map object.

In [52]:
cursor.execute("SELECT ST_AsGeoJSON(ST_Buffer(ST_GeomFromText('{}')::geography,500));".format(p.wkt))
buff=cursor.fetchall()
buffer=json.loads(buff[0][0])
folium.GeoJson(data=buffer).add_to(map_q1)
Out[52]:
<folium.features.GeoJson at 0x7f43f46c3fd0>

In this cell I am creating a 500m buffer around the Point object, saving the buffer geometry to a variable, and adding it to the map.

In [53]:
cursor.execute("SELECT ST_AsText(ST_Buffer(ST_GeomFromText('{}')::geography,500));".format(p.wkt))
bufferwkt=cursor.fetchall()
b=loads(bufferwkt[0][0])

In this cell I am creating a 500m buffer around the Point object and saving the buffer as a variable to extract incidents in the next step.

In [54]:
cursor.execute("SELECT ST_AsGeoJSON(incidents.geom) FROM incidents where ST_Intersects(ST_GeomFromText('{}'), incidents.geom) and date >= NOW() - interval '150 day';".format(b.wkt))
crime=cursor.fetchall()
for x in crime:
    layer=json.loads(x[0])
    folium.GeoJson(data=layer).add_to(map_q1)

map_q1
Out[54]:
Make this Notebook Trusted to load map: File -> Trust Notebook

In this cell I am using the buffer to select all incidents within its area for the last 150 days. After retreiving the data from the DB, I add all of these points onto the map.

Finding the Area with the Highest Number of Incidents in the Past 150 Days Using the Areacommand Widget

In [55]:
cursor.execute("SELECT name FROM areacommand;")
arean=cursor.fetchall()
arean
Out[55]:
[('FOOTHILLS',),
 ('NORTHEAST',),
 ('NORTHWEST',),
 ('SOUTHEAST',),
 ('SOUTHWEST',),
 ('VALLEY',)]

In this cell I select and view all of the different area names.

In [56]:
map_q2 = folium.Map(location=center, zoom_start=zoom)

Creating map object.

In [57]:
@widgets.interact(x="None")
def areacommand(x):
    
    if x:
        cursor.execute("SELECT ST_AsGeoJSON(i.geom) FROM incidents i JOIN areacommand acmd ON ST_Intersects(acmd.geom, i.geom) WHERE acmd.name like'{}' and date >= NOW() - interval '150day';".format(x))
        c=cursor.fetchall()
    
        for x in c:
            layer=json.loads(x[0])
            folium.GeoJson(data=layer).add_to(map_q2)
        return c
    
    else:
        pass

Here is the interactive areacommand widget. It selects all data within 150 days and then adds them to the map based on the name of the area typed into the interactive box.

In [66]:
map_q2
Out[66]:
Make this Notebook Trusted to load map: File -> Trust Notebook

the plotted areacommand data

In the cells below, I selected all incidents by area and got the length of the returned object.

In [59]:
# FOOTHILLS
cursor.execute("SELECT ST_AsGeoJSON(i.geom) FROM incidents i JOIN areacommand acmd ON ST_Intersects(acmd.geom, i.geom) WHERE acmd.name like'{}' and date >= NOW() - interval '150day';".format('FOOTHILLS'))
len(cursor.fetchall())
Out[59]:
881
In [60]:
# NORTHEAST
cursor.execute("SELECT ST_AsGeoJSON(i.geom) FROM incidents i JOIN areacommand acmd ON ST_Intersects(acmd.geom, i.geom) WHERE acmd.name like'{}' and date >= NOW() - interval '150day';".format('NORTHEAST'))
len(cursor.fetchall())
Out[60]:
1199
In [61]:
# NORTHWEST
cursor.execute("SELECT ST_AsGeoJSON(i.geom) FROM incidents i JOIN areacommand acmd ON ST_Intersects(acmd.geom, i.geom) WHERE acmd.name like'{}' and date >= NOW() - interval '150day';".format('NORTHWEST'))
len(cursor.fetchall())
Out[61]:
721
In [62]:
# SOUTHEAST
cursor.execute("SELECT ST_AsGeoJSON(i.geom) FROM incidents i JOIN areacommand acmd ON ST_Intersects(acmd.geom, i.geom) WHERE acmd.name like'{}' and date >= NOW() - interval '150day';".format('SOUTHEAST'))
len(cursor.fetchall())
Out[62]:
1806
In [63]:
# SOUTHWEST
cursor.execute("SELECT ST_AsGeoJSON(i.geom) FROM incidents i JOIN areacommand acmd ON ST_Intersects(acmd.geom, i.geom) WHERE acmd.name like'{}' and date >= NOW() - interval '150day';".format('SOUTHWEST'))
len(cursor.fetchall())
Out[63]:
745
In [64]:
# VALLEY
cursor.execute("SELECT ST_AsGeoJSON(i.geom) FROM incidents i JOIN areacommand acmd ON ST_Intersects(acmd.geom, i.geom) WHERE acmd.name like'{}' and date >= NOW() - interval '150day';".format('VALLEY'))
len(cursor.fetchall())
Out[64]:
1234

The area with the highest number of incidents is 'SOUTHEAST' at 1806 incidents.