Mastering Geospatial Analysis with Python

Incidents in a polygon

Our crime database has a polygon area—area commands and beats—as well as incident points. To build a crime dashboard, we want to be able to map incidents within a specific area command or beat. We can do that by using JOIN and ST_Intersects. The following code shows you how:

cursor.execute("SELECT ST_AsGeoJSON(i.geom) FROM incidents i JOIN areacommand acmd ON ST_Intersects(acmd.geom, i.geom) WHERE acmd.name like'FOOTHILLS' and date >= NOW() - interval '10 day';")crime=cursor.fetchall()for x in crime:    layer=json.loads(x[0])    layergeojson=GeoJSON(data=layer)    map.add_layer(layergeojson)

The previous code selects the geometry from incidents as GeoJSON (ST_AsGeoJSON(i.geom) from incidents), where the incident ST_Intersects ...