Ergebnis der Analyse
Das OpenData-Portal von San Francisco stellt eine detaillierte Kriminalitätsstatistik inklusive Geodaten zur Verfügung. Anhand dieser Daten wird ein Vorgehen vorgestellt mit dem Geodaten anhand der Geokoordinaten gefiltertet werden können. Dieser Filter wird anschließend durch einen Kreis mit dem vorgegebenen Radius auf der Karte verdeutlicht. In diesem Beispiel soll dafür dargestellt werden, wie viele Kriminalitätsfälle bei zwei Hotels in einem bestimmten Umkreis stattgefunden haben.

Das Vorgehen untergliedert sich in folgende Teilschritte:
1. Import der Kriminalitätsstatistik
2. Geocoding mit OpenStreetMap
3. Funktionen für das Geofiltering
4. Grundlage für Erstellen der Kreise auf der Map
5. Zusammenführung und Erstellung der Map

Für dieses Vorgehen werden folgende Packages in R benötigt.

library(RCurl)    #zum Import von Daten aus OpenData SF
library(RJSONIO)  #zum Konvertierung von JSON-Objekte
library(ggplot2)  #zur Datenvisualisierung
library(ggmap)    #zur Integration von Kartenmaterial bei der Datenvisualisierung in 
library(dplyr)    #zur (komfortableren) Verarbeitung von Data-Frames
library(knitr)    #zur Erstellung dieses Blogposts
library(rmarkdown)#zur Erstellung dieses Blogposts

1. Schritt: Download Kriminalitätsstatistik von OpenData SF

Die Daten der Kriminalitätsstatistik für das Jahr 2015 stammen aus Open Data-Portal der Stadt und des County San Franciscos. Über die API des OpenData-Portals können die Daten direkt in R importiert werden.

data.url <- getURL("https://data.sfgov.org/api/views/ritf-b9ki/rows.csv")
crime.data <- read.csv(textConnection(data.url))

Die Kriminalitätsstatistik enthält eine detaillierte Übersicht aller Vorfälle in San Francisco für das Jahr 2015. Insgesamt liegen 154.694 Beobachtungen vor. Zu diesen Beobachtungen enthält das Dataframe insbesondere die in der folgenden Tabelle aufgeführten und später verwendeten Informationen.

Variable Beschreibung
Dates Kategorie des Vorfalls (insg. 38 levels)
Category detaillierte Beschreibung des Vorfalls
DayOfWeek Wochentag des Vorfalls
X Longitude
Y Latitude

Für einen Touristen sind jicht alle Vorfälle gleichermaßen relevant. Mit der Funktion myCrimeDataFilter lässt sich der Datensatz nach der Art der Kriminalität, dem Wochentag (z.B. nur Wochenenden) und dem Monat einschränken.

myCrimeDataFilter <- function(data, cat, day, mon){
    crime.fil <- na.omit(data)
    crime.fil <- filter(crime.fil, Category %in% cat)
    crime.fil <- filter(crime.fil, DayOfWeek %in% day)
    crime.fil <- filter(crime.fil, Month %in% mon)
    return(crime.fil)
}

Der Monat ist noch nicht in dem Datensatz enthalten und muss aus dem Datum extrahiert werden.

crime.data$Month <- as.character(crime.data$Date)
crime.data$Month <- substr(crime.data$Month, 1, 10)
crime.data$Month <- as.POSIXct(crime.data$Month, format="%m/%d/%Y")
crime.data$Month <- months(crime.data$Month)

2. Schritt: Funktion zur Ermittlung der Geodaten für eine vorgegebene Adresse (Geocoding)

Wenn man für eine vorgegebene Adresse (z.B. ein Hotel) die Geokoordinaten ermittelt, wird dieses Vorgehen Geocoding oder Geotagging genannt. In diesem Vorgehen wird hierfür Nominatim verwendet. Nominatim ist ein Werkzeug, um in OpenStreetMap über den Namen und Adresse nach Objekten zu suchen. Der Dienst ist unter nominatim.openstreetmap.org zu finden. Der Zugriff auf Nominatim und damit auf OpenStreetMap erfolgt über die nachfolgend definierte Hilfsfunktion myOSMGeoCode.

myOSMGeoCode <- function(street, zip, city, country){
  
  url <-paste0(
    "http://nominatim.openstreetmap.org/search.php?q=", 
    street, ",+", 
    zip, "+",
    city, ",+", 
    country, 
    "&limit=1&format=json")
 
  osm.json <- fromJSON(url,simplify=FALSE)
  if(length(osm.json) > 0){
    r1 <- osm.json[[1]]$lon
    r2 <- osm.json[[1]]$lat
    r3 <- osm.json[[1]]$display_name
    r4 <- "OK"
    geocode <- c(r1,r2,r3,r4)
  } else{
    url <-paste0("http://nominatim.openstreetmap.org/search.php?q=", zip, "+", city, ",+", country, "&limit=1&format=json")
    osm.json <- fromJSON(url,simplify=FALSE)
    if(length(osm.json) > 0){
      r1 <- osm.json[[1]]$lon
      r2 <- osm.json[[1]]$lat
      r3 <- osm.json[[1]]$display_name
      r4 <- "ZIP"
      geocode <- c(r1,r2,r3,r4)
    } else{
      geocode <- c(NA,NA,NA,"NOK")
    }
  }
  return(geocode)
}

3. Schritt: Funktionen zum Filtern der Daten auf Basis von Geodaten (Geofiltering)

Für einen Touristen sind nicht die Kriminalitätsfälle in dem gesamten Stadtgebiet relevant, sondern in einem bestimmten Stadtteil. In diesem Beispiel steht der Tourist vor der Wahl, welches Hotel er auswählen soll. Ein Kriterium für ihn ist dabei, ob in der näheren Umgebung des Hotel viele Straftaten passieren. Die Kriminalitätsstatistik muss somit in Abhängigkeit der Geodaten gefiltert werden.

Als Basis für die Filterung anhand von Geodaten muss die Entfernung zwischen zwei Standorten ermittelt werden können. Die Orthodrome ist die kürzeste Verbindung zweier Punkte auf einer Kugeloberfläche. Unterstellt man, dass die Erde eine idealen Kugel mit einem Radius von 6378.388 km ist, lässt sich hiermit mit der Formel für die Orthodrome die Entfernung zwischen zwei Adressen ermitteln (“Luftline”). Dieses Vorgehen ist zum Beispiel bei http://www.kompf.de/trekka/distance.php umgesetzt und dient als Grundlage für die nachfolgende Funktion myEarthDistance

myEarthDistance <- function (lon1, lat1, lon2, lat2){
  
  rad <- pi/180
  b1 <- lat1 * rad
  l1 <- lon1 * rad
  b2 <- lat2 * rad
  l2 <- lon2 * rad

  dist = 6378.388 * acos(sin(b1) * sin(b2) + cos(b1) * cos(b2) * cos(l2 - l1))
  return(dist)
}

Die Funktion myGeoGilter verwendet die oben definierte Funktion myEarthDistance. Zu jedem Punkt der Kriminalitätsstatistik wird der Abstand zu einer vorgegebener Geokoordinate (z.B. dem Standort des Hotels) ermittelt. Anschließend wird das Dataframe auf diejenigen Datensätze eingeschränkt, deren Abstand innerhalb eines vorgegebenen Radius liegen.

myGeoFilter <- function(data, lon, lat, radius){
    
  data$Distance <- 0
  data$WithinRadius <- FALSE
    
  for(i in 1:dim(data)[1]){
      data[i,"Distance"] <- myEarthDistance(lon1=lon, lat1=lat, lon2=data[i,"X"], lat2=data[i,"Y"])
  }
    
  data$WithinRadius <- ifelse(data$Distance <= radius, TRUE, FALSE)
  data <- filter(data, WithinRadius == TRUE)
  return(data)
}

4. Schritt: Funktionen zum Plotten der Daten

Der eingesetzte Geofilter soll durch einen Kreis mit dem vorgegebenen Radius auf der Map verdeutlicht werden. Für den Kreis auf der Map ist ein Dataframe zu erzeugen, dass die entsprechenden Datenpunkte in dem vorgegebenen Abstand von der Geokoordinate enthält. Es ist somit eine Umkehrung der Funktion myEarthDistance.

createCircleData <- function(lon, lat, radius){
    rad <- pi/180
    b1.rad <- lat * rad
    l1.rad <- lon * rad

    kw.rad <- seq(1:360)*rad
    b2.rad <- asin( sin(b1.rad) * cos(radius/6378.388) + cos(b1.rad) * sin(radius/6378.388) * cos(kw.rad))
  
    l2.rad <- l1.rad + atan2( sin(kw.rad) * sin(radius/6378.388) * cos(b1.rad) ,  cos(radius/6378.388) - sin(b1.rad)*sin(b2.rad))
  
    b2 <- b2.rad*(1/rad)
    l2 <- l2.rad*(1/rad)
    
    return(data.frame(Latitude = b2,
                      Longitude = l2,
                      ID = rep(as.character(lon),360)))
}

Das nachfolgend definierte Theme für ggplot2 ist die Grundlage für die Karte und enthält keinerlei Achsenbeschriftungen mehr.

myMapTheme <- theme(legend.position="bottom",
                    panel.grid.major=element_blank(),
                    panel.grid.minor=element_blank(),
                    axis.title=element_blank(),
                    axis.ticks=element_blank(),
                    axis.text=element_blank(),
                    panel.background = element_rect(fill = "white"))

5. Schritt: Zusammenführung und Erstellung der Map

Im letzten Schritt werden die obigen Schritte zusammengeführt. Es sollen die beiden Hotels Best Western Tuscan Inn und W San Francisco miteinader verglichen werden. Die Adressen der beiden Hotels sind in einer csv-Datei gespeichert.

hotel.data <- read.csv("Data/HotelDataSF.csv", sep=";", stringsAsFactors=FALSE)
Name Street Disctric City Zip Country
W San Francisco 181 3rd Street South of Market (SOMA) San Francisco CA 94124 USA
Best Western Tuscan Inn 425 North Point Street Russian Hill San Francisco CA 94133 USA

Zu beiden Hotels müssen nun die Geokoordinaten mit OpenStreetMap und der erstellten Funktion myOSMGeoCode ermittelt werden.

for(i in 1:dim(hotel.data)[1]){
    
  street <- hotel.data[i, "Street"]
  zip <- hotel.data[i, "ZIP"]
  city <- hotel.data[i, "City"]
  country <- hotel.data[i, "Country"]
    
  geo.data <- myOSMGeoCode(street, zip, city, country)
    
  hotel.data[i,"Lon"] <- as.numeric(geo.data[1])
  hotel.data[i, "Lat"] <- as.numeric(geo.data[2])
  hotel.data[i, "Status"] <- geo.data[4]
}

In diesem Beispiel sollen lediglich Raubüberfalle und Diebstähle betrachtet werden. Die anstehenden Monate Mai bis August sind von besonderem Interesse und entsprechend ist die Kriminalitätsstatistik mit der Funktion myCrimeDataFilter gemäß diesen Kriterien einzuschränken.

crime.category <- c("STOLEN PROPERTY", "ROBBERY")
crime.weekdays <- c("Monday", "Tuesday", "Wedensday", "Thursday", "Friday","Saturday","Sunday")
crime.months <- c("Mai", "Juni", "Juli", "August")

crime.filtered <- myCrimeDataFilter(data=crime.data, cat = crime.category, day = crime.weekdays, mon = crime.months)

Im zweiten Schritt wird die Funktion myGeoFilter angewendet. Es sind nur noch diejenigen Kriminalitätsfälle von Interesse, die innerhalb eines Radius von 1km um die beiden Hotels beobachtet wurden.

crime.filtered.1 <- myGeoFilter(data=crime.filtered, lon = hotel.data[1,"Lon"], lat = hotel.data[1,"Lat"], radius = 1)
crime.filtered.2 <- myGeoFilter(data=crime.filtered, lon = hotel.data[2,"Lon"], lat = hotel.data[2,"Lat"], radius = 1)
crime.filtered <- rbind(crime.filtered.1, crime.filtered.2)

Für die Visualisierung auf der Karte sind nun noch die Daten für die Kreise auf der Karte zu ermitteln. Dies erfolgt mit der Funktion createCircleData.

  circle.data1 <- createCircleData(lon = hotel.data[1,"Lon"], lat = hotel.data[1,"Lat"], radius = 1)
  circle.data2 <- createCircleData(lon = hotel.data[2,"Lon"], lat = hotel.data[2,"Lat"], radius = 1) 
  circle.data <- rbind(circle.data1, circle.data2)

Als Hintergrund der Karte soll der Stadtplan von San Francisco dienen. Dieser wird von OpenStreetMap mit der Funktion get_map aus dem Package ggmap importiert.

map.sf <- get_map(location = c(-122.44, 37.77, -122.38, 37.82), source = "osm", color = "bw")

Zum Abschluss werden die gesamten Daten zu einer gesammelten Map zusammengeführt und das Ergebnis gespeichert.

g <- ggmap(map.sf)
g <- g + geom_point(data=crime.filtered, aes(x=X, y=Y, color=Category), alpha=1, size=1.5)
g <- g + geom_point(data=hotel.data, aes(x=Lon, y=Lat), color="black", alpha=1, size=1.5)
g <- g + geom_point(data=circle.data, aes(x=Longitude, y=Latitude), color="black", alpha=0.5, size=1)
g <- g + geom_text(data=hotel.data, aes(x=Lon, y=Lat, label=Name), hjust=0, vjust=0, fontface="bold", color="black", size=3)
g <- g + myMapTheme
g

ggsave("SF_Crime_Map.jpg")