Am I Near Public Art? How about Crime?

When you are new to a city, it is hard to walk around and carry a map of all the sites. Wouldn’t it be nice if your phone could alert you when you were near something of interest? In this post, I will show you how to create an alert system.

The Application

The application will build off of our earlier post using Albuquerque Crime Data. In that application, the user clicked the map and was presented with all crimes in the last 180 days within a half mile of the click.

This application will use the same logic but replace the click with the users location.

Building the Application

We will start with a Leaflet.js Mobile map. Place the following HTML in the head of your document.

<meta name=”viewport” content=”width=device-width, initial-scale=1.0, maximum-scale=1.0, user-scalable=no” />

<style>
body {
padding: 0;
margin: 0;
}
html, body, #map {
height: 100%;
}
</style>

This will make the map look nice on mobile – and desktop as well. Then create the map.

var map = L.map(‘map’);

L.tileLayer(‘http://{s}.tile.osm.org/{z}/{x}/{y}.png’).addTo(map);

map.locate({setView: true, maxZoom: 16,enableHighAccuracy:true,watch:true});

With mobile maps, we do not pass the setView to the map, we add that parameter to the locate() function. Now we need two functions: one for found and one for not found.

map.on(‘locationfound’, gotcha);

function gotcha(e){ALL OUR CODE GOES HERE}

map.on(‘locationerror’, onLocationError);

function onLocationError(e) {
alert(e.message);}

When our location is found, we want to setup the buffer and prep for calling the Albuquerque Data for Public Art or Incidents.

a=L.marker(e.latlng);
b=a.toGeoJSON();
L.circle(e.latlng, 20).addTo(map);
buffered = turf.buffer(b,.05,”miles”);
var result = turf.featurecollection(buffered.features.concat(b));
g='{“rings”:’+JSON.stringify(buffered.features[0].geometry.coordinates)+’}’;

We create a marker, convert it to GeoJSON to use with Turf.js, then buffer and get the rings to use in the query.

Now we can make our standard AJAX call to an ESRI REST Endpoint. I have commented out the URL for Public Art and the marker. You can uncomment them and comment out the incident url and marker to switch the application.

var params=”f=json&outSR=4326&outFields=*&geometryType=esriGeometryPolygon&spatialRel=esriSpatialRelIntersects&inSR=4326&geometry=”+g;

//var url = “http://coagisweb.cabq.gov/arcgis/rest/services/public/PublicArt/MapServer/0/query&#8221;;
var url = “http://coagisweb.cabq.gov/arcgis/rest/services/public/APD_Incidents/MapServer/0/query&#8221;;
http=new XMLHttpRequest();
http.open(“POST”, url, true);
http.setRequestHeader(“Content-type”, “application/x-www-form-urlencoded”);
http.onreadystatechange = function() {//Call a function when the state changes.
if(http.readyState == 4 && http.status == 200) {
var result= JSON.parse(http.responseText);
…….CODE GOES HERE……..

}}
http.send(params);

Now we get back all the incidents, or art, within 300 feet (.05 miles) of our location. We need to take the results and add them to the map and notify the user they are near something of interest.

for(x=0;x<Object.keys(result.features).length;x++){
xy=result.features[x].geometry;
pattern.push(100);
pattern.push(100);
L.marker([result.features[x].geometry.y,result.features[x].geometry.x]).bindPopup(“<h3>”+result.features[x].attributes.CVINC_TYPE+”</h3>”).addTo(map);}
//L.marker([result.features[x].geometry.y,result.features[x].geometry.x]).bindPopup(“<h3>”+result.features[x].attributes.TITLE+”</h3>”).addTo(map);}

navigator.vibrate(pattern);
pattern=[];

As we move, the service is queried again, constantly watching us and letting us know when we are near something of interest. The vibration pattern is 100ms for each feature we are near and pause for 100ms between each feature. If we are near a lot of features -say crime incidents – our phone will vibrate many times. Possibly signaling a dangerous part of town. If we are using public art, we know we are near a particularly artistic part of the city and maybe we should explore more.

To use the same map for multiple features, you could switch up the vibration patterns, say short fast vibrations for art and longer for crime. Maybe vibrate once for crime but the duration is a function of the number of incidents – 10 incidents results in a 10 second vibration.

Just an idea on how to use open data to alert new comers to the features of the city that they may not be aware of.

Advertisements

Crime Density in Albuquerque

I posted earlier on Air Pollution Density in Albuquerque. I used pollution sources because the crime data was missing x,y values. In this post, I will add a line to the code to allow us to use Crime Data for the last 180 days.

The Application

The application maps the density of crime over the last 180 days. The result is in the image below.

To change the data from pollution to crime, just switch out the url using the link to the Crime Data.

url=”http://coagisweb.cabq.gov/arcgis/rest/services/public/APD_Incidents/MapServer/0/query&#8221;

Because the data is missing (x,y) values, our application will crash. We can add an if statement to check if the coordinates are null and if so, don’t map it.

if(result.features[x].geometry===undefined){}
else{var t = L.marker([result.features[x].geometry.y,result.features[x].geometry.x]);//.addTo(map);
p.push(t.toGeoJSON());}

Lastly, you will need to change the color scale to account for the change in values for each polygon.

function getColor(d) {
return d > 140 ? ‘#800026’ :
d > 120 ? ‘#BD0026’ :
d > 80 ? ‘#E31A1C’ :
d > 60 ? ‘#FC4E2A’ :
d > 40 ? ‘#FD8D3C’ :
d > 20 ? ‘#FEB24C’ :
d > 5 ? ‘#FED976’ :
‘#FFFFFF’;
}

Data does not come perfectly clean. It is our job to try to deal with problems in the data in our applications – or at least make it fail gracefully.

Try changing the grid size or the color values:
JS Bin

 

Improving Albuquerque Crime Incidents Query

In my last post, I showed how to query the Albuquerque Crime Incident. This post will improve on it.

In my original post, you had to do the following:

1. Catch a point on click
2. Draw a buffer
3. Get the envelope of the buffer
4. Query

There was nothing wrong with this, but it does bring in incidents beyond a half mile – the envelope is a square representation of the buffer. But we used the envelope parameter in the REST endpoint so we had no choice.

Well, let’s use the polygon parameter and only grab points in our buffer!

In the original code, after the line:

buffered = turf.buffer(b,.5,”miles”);

add the line below.

g='{“rings”:’+JSON.stringify(buffered.features[0].geometry.coordinates)+’}’;

You can remove the envelope line from the original code.

Now, just change your parameters to reflect the geometryType=Polygon and the new geometry=g.

var params=”f=json&outSR=4326&outFields=*&geometryType=esriGeometryPolygon&spatialRel=esriSpatialRelIntersects&inSR=4326&geometry=”+g;

You will now have incidents only in your buffer.

Query Albuquerque Open Data with Turf.js

I have posted on the Albuquerque Crime Incident data showing how to use Python and for converting the Date field. I have also demostrated how to use Turf.js for analysis. In this post, I will show you how to combine the two to limit your query of the Albuquerque Crime Incident Data.

The Idea

The goal of this post is to allow a user to click on a map and get back all the crime data within a half mile of the click. The end result will look like the image below.

 

Create the Geometry

To make this work, we need to catch where the user clicks on the map and convert it to GeoJSON so we can use it in our Turf.js operations without having to convert it to a Turf.point.

map.on(“click”,function(e){

a=L.marker(e.latlng);
b=a.toGeoJSON();

});

Now that we have a point, we need to create the buffer and then grab the envelope.

var buffered = turf.buffer(b,.5,”miles”);
var result = turf.featurecollection(buffered.features.concat(b));
enveloped = turf.envelope(result);

Query the Service

With the envelope, we can query the REST endpoint and pass the lower left and upper right coordinates to get the Crime Incidents in the bounding box. I have done several posts on making an AJAX request to ArcServer and will not go in to detail here. Just note that we will have two parameters that I do not use often.

geometryType=esriGeometryEnvelope

spatialRel=esriSpatialRelIntersects

And a geometry parameter that will take the corner points of the envelope. The code below shows how to make the query.

var params=”f=json&outSR=4326&outFields=*&geometryType=esriGeometryEnvelope&spatialRel=esriSpatialRelIntersects&inSR=4326&geometry=”+enveloped.geometry.coordinates[0][0][0]+”,”+enveloped.geometry.coordinates[0][0][1]+”,”+enveloped.geometry.coordinates[0][2][0]+”,”+enveloped.geometry.coordinates[0][2][1];
var url = “http://coagisweb.cabq.gov/arcgis/rest/services/public/APD_Incidents/MapServer/0/query&#8221;;
console.log(params);
http=new XMLHttpRequest();
http.open(“POST”, url, true);
http.setRequestHeader(“Content-type”, “application/x-www-form-urlencoded”);
http.onreadystatechange = function() {//Call a function when the state changes.
if(http.readyState == 4 && http.status == 200) {
console.log(http.responseText);
var result= JSON.parse(http.responseText);
for(x=0;x<Object.keys(result.features).length;x++){
xy=result.features[x].geometry;
console.log(result.features[x].geometry);
adate=new Date(result.features[x].attributes.date);
L.marker([result.features[x].geometry.y,result.features[x].geometry.x]).bindPopup(“<h3>”+result.features[x].attributes.CVINC_TYPE+”</h3><p>”+adate.toString()).addTo(map);}
}}

http.send(params);

I added a popup showing the crime type and then formatted the epoch data to a human readable string. Now if you click on the map, you will get back a series of incidents. To improve on this, you could pass a specific date range but the data itself is for the last 180 days (rolling, meaning tomorrow one day will be dropped off).

Public Safety GIS 2: Python, MongoDB, C# and ArcGIS Desktop Explorer

I wrote a plugin for ArcGIS Desktop explorer that called a python script and passed the map center as an argument. The python script queried a MongoDB and wrote out a csv of closest points. The button then loads the csv as points. I had to do this because I didn’t know how to do it in C#. Well, I figured it out.

Now C# does all the work querying the database and drawing the points. I still have a python script that reads data from a csv export from sharepoint, geocodes it and creates the mongoDB.

Here is my code:

PYTHON TO CREATE THE DB:

from pymongo import Connection,GEO2D
import urllib, urllib2, simplejson
import csv

data = csv.reader(open(‘C:\Documents and Settings…\myfile.csv’))
fields=data.next()
db=Connection().rtcc
db.cnau.create_index([(“loc”,GEO2D)])

for row in data:
try:
address=row[0]
combinedaddress=address+”,albuquerque, NM”

param = {‘q’: combinedaddress,’output’: ‘json’, ‘oe’: ‘utf8’}
url = ‘http://maps.google.com/maps/geo?&#8217; + urllib.urlencode(param)
rawreply = urllib2.urlopen(url).read()
reply = simplejson.loads(rawreply)
print reply[‘Placemark’][0][‘Point’][‘coordinates’][0]
long=reply[‘Placemark’][0][‘Point’][‘coordinates’][0]
lat=reply[‘Placemark’][0][‘Point’][‘coordinates’][1]
z={“loc”:[long,lat],”Activity”:row[1]}

db.cnau.insert(z)
except KeyError:
pass

 

C# TO PLUGIN ARCGIS DESKTOP EXLPORER:

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using ESRI.ArcGISExplorer;
using ESRI.ArcGISExplorer.Application;
using ESRI.ArcGISExplorer.Mapping;
using ESRI.ArcGISExplorer.Geometry;
using ESRI.ArcGISExplorer.Data;
using ESRI.ArcGISExplorer.Threading;

using MongoDB;
using MongoDB.Bson;
using MongoDB.Driver.Linq;
using MongoDB.Driver;
using MongoDB.Driver.Builders;

namespace MongoC
{
public class Button : ESRI.ArcGISExplorer.Application.Button
{

public override void OnClick()
{

Folder folderMapItem = new Folder();
MapDisplay mdisp = ESRI.ArcGISExplorer.Application.Application.ActiveMapDisplay;

//use this to change the projections 35,-106
FeatureLayer proj = FeatureLayer.OpenShapefile(“C:\\Documents and Settings…\\forProjection.shp”);
// Path for Windows 7: FeatureLayer proj = FeatureLayer.OpenShapefile(“C:\\Users\\user\\…\\forProjection.shp”);
mdisp.CoordinateSystem2D = proj.CoordinateSystem;
string center = mdisp.Center.ToString();
folderMapItem.Name = “Close To Me”;
string[] c = center.Split(null);
String lon = c[3];
String lat = c[6];

MongoClient client = new MongoClient(); // connect to localhost. Deploy at: “mongodb://ip.of.host:27017”
MongoServer server = client.GetServer();
MongoDatabase db = server.GetDatabase(“rtcc”);

var collection = db.GetCollection<BsonDocument>(“cnau”);

var query = Query.Near(“loc”, double.Parse(c[3]), double.Parse(c[6]));
//var cursor = collection.Find(query); Put this in foreach so dont need it anymore.
foreach (BsonDocument item in collection.Find(query).SetLimit(5))
{
BsonElement loc = item.GetElement(“loc”);
string g = loc.Value.ToString();
string x = g.Trim(new Char[] { ‘[‘, ‘]’ });
String[] a = x.Split(‘,’);
// BsonElement name = item.GetElement(“name”); Work for all indexes in MongoDB
// Pass the other data to the notes popups below

ESRI.ArcGISExplorer.Geometry.Point mypoint = new ESRI.ArcGISExplorer.Geometry.Point(double.Parse(a[0]), double.Parse(a[1]));

// Create a Note using the Point geometry.
Note mypointNote = new Note(“MyNote”, mypoint, Symbol.Marker.Pushpin.Orange);

// Populate the Note popup with information about the Point geometry using the ToString method.
StringBuilder info = new StringBuilder(@”<p><a href=””http://google.com””>Google</a></p>&#8221;);
info.AppendLine(@”<p>” + mypoint.ToString() + @”</p>”);

//MapDisplay disp = ESRI.ArcGISExplorer.Application.Application.ActiveMapDisplay;
mypointNote.Popup.Content = info.ToString();

// Add the Note to the current Map.
//mdisp.Map.ChildItems.Add(mypointNote);
folderMapItem.ChildItems.Add(mypointNote);
}

mdisp.Map.ChildItems.Add(folderMapItem);
}
}
}

Python and SQLite3

I need to generate a report every month. This report needs to have 30 names on it and I need to make sure than none of the people have appeared in previous months. I used Python and SQLite to store my items and to verify them. Here is what I did:

I wrote a script to empty all the data from an Excel sheet in to a SQLite Database Table. Then I wrote my program. The if __name__ == “__main__”: takes my command line arguments (load or scan and the data file.xls) and calls the appropriate function. If there is an error, it prints out the correct usage.

if __name__ == “__main__”:

if (len(sys.argv) < 3):

print “USAGE: warrant.py load/scan excelfile.xls”

elif sys.argv[1]==”load”:

load(sys.argv[2])

elif sys.argv[1]==”scan”:

scan(sys.argv[2])

else: print “Error…..Unknown Argument. \n\nUSAGE: warrant.py load/scan excel.xls”

The function scan() runs through an excel sheet and compares the names to the database. Currently, if the same last and first appear it prints out the info for the person. I will fix this by using ss# as another filter.

The load() function takes an excel file and uploads them to the database. No error checking. It will crash if the name already exists, but the list should have been cleaned with scan() first. I need to fix this anyway.

I hate error checking. My programs are usually for me so I don’t care if they fail because I can read the error and track it down. But now that I need anyone to be able to use it,  I cant have this – when the excel file doesn’t exist:

Error when Python can’t find the Excel file.

Instead, I need to catch the errors and fail gracefully – like when the commands entered are wrong:

Catching an error in usage – gracefully.

Fixed the error if the file is not found by using from os.path import exists then updating the if name=main to include:

elif sys.argv[1]==”load”:
           if exists(sys.argv[2]):
                  load(sys.argv[2])
           else:
                  print “file not found”

So here is my code thus far:

import sqlite3 as lite
import sys
from os.path import exists
from xlrd import open_workbook

def load(data):
    wb=open_workbook(data)
    sheet=wb.sheet_by_index(0)
    con = lite.connect(‘warrant.db’)
    with con:
        cur = con.cursor()
        for rownum in range(sheet.nrows):
            fname=sheet.cell(rownum,0).value
            lname=sheet.cell(rownum,1).value
            dob=sheet.cell(rownum,2).value
            ss=sheet.cell(rownum,3).value
            cur.execute(“INSERT INTO previous VALUES(fname,lname,dob,ss)”)
            print “Added: “+fname+”, “+lname+”, “+dob+”, “+ss

def scan(data):
    wb=open_workbook(data)
    sheet=wb.sheet_by_index(0)
    con = lite.connect(‘warrant.db’)
    with con:
        cur = con.cursor()
        for rownum in range(sheet.nrows):
            lastname=sheet.cell(rownum,1).value
            firstname=sheet.cell(rownum,0).value
            query=”SELECT * FROM previous WHERE lname='”+lastname+”‘ AND fname='”+firstname+”‘”
            cur.execute(query)    
            while True:
                row=cur.fetchone()
                if row==None:
                    break
                print row[0], row[1], row[2], row[3]

if __name__ == “__main__”:
    if (len(sys.argv) < 3):
        print “\n*************************\n\n”+”USAGE: warrant.py load/scan excelfile.xls\n\n”+”*************************”
    elif sys.argv[1]==”load”:
        if exists(sys.argv[2]):        
            load(sys.argv[2])
        else:
            print “\n***********************\nFile: “+ sys.argv[2] +” doesn’t exist.\n***********************\n”
    elif sys.argv[1]==”scan”:
        if exists(sys.argv[2]):
            scan(sys.argv[2])
        else:
            print “\n***********************\nFile: “+ sys.argv[2] +” doesn’t exist.\n***********************\n”
    else:
        print “\n*************************\n\n”+”Error…..Unknown Argument. \n\nUSAGE: warrant.py load/scan excel.xls\n\n”+”*************************”

 

Social Network in Python Using NetworkX

My First Network Graph using NetworkX and Gephi.

I recently stumbled upon the PyData NYC 2012 Channel on Vimeo. This in turn led me to a talk by Gilad Lotan on NetworkX and Gephi. His talk is awesome, please check it out. I have played with NetworkX a little, mostly for a demonstration on Python and mapping I gave, so it was just an overview and I didn’t get to spend much time on it. After watching Gilad’s talk, I thought I would revisit NetworkX. I want to show a simple example of how to create a graph using NetworkX. I will save it as a PNG and a .GEXF to be edited in Gephi.

I hard coded a series of Nodes and Edges in my Python code. The nodes are the Names of people in the circles on the above image. The edges are the lines between them. I then used matplotlib to show the graph and to save it out as a .gexf. Here is the code:

import networkx as nx
import matplotlib.pyplot as plt

g=nx.Graph()

g.add_node(“Paul”)
g.add_node(“Joe”)
g.add_node(“Bill”)
g.add_node(“John”)
g.add_node(“Sarah”)
g.add_node(“Mary”)
g.add_node(“Karen”)
g.add_node(“Betsy”)
g.add_edge(“Betsy”,”Karen”)
g.add_edge(“Betsy”,”Paul”)
g.add_edge(“Betsy”,”Sarah”)
g.add_edge(“Paul”,”Sarah”)
g.add_edge(“Karen”,”Sarah”)
g.add_edge(“Sarah”,”Mary”)
g.add_edge(“Paul”,”Mary”)
g.add_edge(“Bill”,”Mary”)
g.add_edge(“John”,”Mary”)

nx.draw(g)
plt.savefig(“Something.png”)
nx.write_gexf(g, “Something.gexf”)

Import the libraries, create a graph, add nodes, add edges, then save. Simple. Now you can pull it in to Gephi and make it pretty- since this is the output I got from Matplot:

Output from Matplotlib – plt.show()

I won’t go in to detail on how to use Gephi – because honesty, I just used it for the first time – watch Gilad’s talk and learn from him – he knows what he’s doing!