GIS Beyond ArcMap

• One use is to make quick maps
• First step is always to load the needed libraries

library(tidyverse)
library(readxl)
library(rgdal)
library(ggmap)

• Available from NYC.gov; I use 2014 file from NYCLU

nyc_sqf <- read_xlsx("~/Google Drive/Computer Backup/R Workshop/Data/2014_sqf_web.xlsx")

select(nyc_sqf, xcoord, ycoord) %>%
  summary()

     xcoord            ycoord      
 Min.   : 914151   Min.   :121420  
 1st Qu.: 990886   1st Qu.:176014  
 Median :1005210   Median :192836  
 Mean   :1005444   Mean   :198626  
 3rd Qu.:1023191   3rd Qu.:221223  
 Max.   :1066574   Max.   :271882  
 NA's   :1650      NA's   :1650    

• The coordinates in the data are from a UTM projection
• We want latitude/longitude for mapping

lon_lat <- nyc_sqf %>% 
  select(xcoord, ycoord) %>%
  as.matrix() %>%
  project(proj = "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs", inv = TRUE)

nyc_sqf$lon <- lon_lat[,1]
nyc_sqf$lat <- lon_lat[,2]

select(nyc_sqf, lon, lat) %>%
  summary()

      lon              lat       
 Min.   :-74.25   Min.   :40.50  
 1st Qu.:-73.98   1st Qu.:40.65  
 Median :-73.92   Median :40.70  
 Mean   :-73.92   Mean   :40.71  
 3rd Qu.:-73.86   3rd Qu.:40.77  
 Max.   :-73.70   Max.   :40.91  
 NA's   :1650     NA's   :1650   

• ggmap is an extension of ggplot2 library for graphics
• Relies on static web map tiles
• Data still just in a table, not yet spatial

get_map(location = "new york", 
        maptype = "satellite",
        zoom = 11) %>%
  ggmap() +
  coord_cartesian() +
  geom_point(data = nyc_sqf, aes(lon, lat),
             alpha = 0.05, color = "#ce1c1c", na.rm = TRUE)

library(tmap)

data(Europe, land, rivers, metro)

tm_shape(land) + 
    tm_raster("trees", breaks=seq(0, 100, by=20), legend.show = FALSE) +
tm_shape(Europe, is.master = TRUE) +
    tm_borders() +
tm_shape(rivers) +
    tm_lines(lwd="strokelwd", scale=5, legend.lwd.show = FALSE) +
tm_shape(metro) +
    tm_bubbles("pop2010", "red", border.col = "black", border.lwd=1, 
        size.lim = c(0, 11e6), sizes.legend = c(1e6, 2e6, 4e6, 6e6, 10e6), 
        title.size="Metropolitan Population") +
    tm_text("name", size="pop2010", scale=1, root=4, size.lowerbound = .6, 
        bg.color="white", bg.alpha = .75, 
        auto.placement = 1, legend.size.show = FALSE) + 
tm_format_Europe() +
tm_style_natural()

• Can create points file similar to “Display XY Data” in ArcMap

library(sp)

nyc_sqf <- filter(nyc_sqf, is.na(lon) == F)
coordinates(nyc_sqf) <- select(nyc_sqf, lon, lat)

class(nyc_sqf)

[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"

slotNames(nyc_sqf)

[1] "data"        "coords.nrs"  "coords"      "bbox"        "proj4string"

• Let's take a random day and see if there's clustering in stop-and-frisk locations
• spatstat has many tools to analyze point patterns

library(spatstat)
library(maptools)

sample <- nyc_sqf[nyc_sqf$datestop==1092014,] %>% as.ppp()
window <- convexhull(sample) %>% as.owin()

sample[window] %>%
  envelope(fun = Kest) %>%
  plot(main = "Clustering in one day of NYPD stops")

sa <- readOGR("/Users/benjaminbellman/Dropbox (Brown)/Police/Data/Final Files", "SanAntonio_1015")

OGR data source with driver: ESRI Shapefile 
Source: "/Users/benjaminbellman/Dropbox (Brown)/Police/Data/Final Files", layer: "SanAntonio_1015"
with 366 features
It has 81 fields
Integer64 fields read as strings:  REGIONA DIVISIONA STATEA COUNTYA COUSUBA PLACEA TRACTA CONCITA AIANHHA RES_ONLYA TRUSTA AITSCEA ANRCA CBSAA CSAA METDIVA NECTAA CNECTAA NECTADIVA UAA CDCURRA SLDUA SLDLA ZCTA5A SUBMCDA SDELMA SDSECA SDUNIA PUMA5A BTTRA BLKGRPA BTBGA medhhinc totpop white black asian hisp nonwhite Rowid_ ZONE_CODE COUNT 

library(spdep)
library(classInt)

weights <- poly2nb(sa) %>% nb2listw()
locm <- localmoran(sa$gini, weights)
sa$Sgini <- scale(sa$gini)
sa$lag <- lag.listw(weights, sa$Sgini)
sa$pval <- locm[, 5]

sa$quad_sig <- ifelse(sa$Sgini >= 0 & sa$lag >= 0 & sa$pval <= 0.05, 1,
                      ifelse(sa$Sgini <= 0 & sa$lag <= 0 & sa$pval <= 0.05, 2,
                      ifelse(sa$Sgini >= 0 & sa$lag <= 0 & sa$pval <= 0.05, 3,
                      ifelse(sa$Sgini >= 0 & sa$lag <= 0 & sa$pval <= 0.05, 4, 5))))

breaks <- seq(1, 5, 1)
labels <- c("High-High", "Low-Low", "High-Low", "Low-High", "Not Signif.")
np <- findInterval(sa$quad_sig, breaks)
colors <- c("red", "blue", "lightpink", "skyblue2", "white")