geoprofiling_w_map.R

geoprofiling_w_map.R

# Geographic profiling
# Methodology developed by Kim Rossmo
# R program by Jeremy Coleman

library(maps)
library(grDevices)

### Define the geoprofiling function

geo.prob <- function(xi, yi, xn, yn, B, phi, f, g){
  
  if(g<1 || f<1){stop("f and g must be greater than 1")}
  
  if(phi>=1 || phi<=0){stop("phi must be between 0 and 1")}
  
  n <- length(xn)
  
  prob.1 <- (phi/((abs(xi - xn) + abs(yi - yn))^f))
  
  prob.2 <- (((1 - phi)*(B^(g-f)))/((2*B - abs(xi - xn) - (yi - yn))^g))
  
  same <- which(is.infinite(prob.1)); prob.1[same] <- 0
  same <- which(is.infinite(prob.2)); prob.2[same] <- 0
  
  prob <- (1/n)*sum(prob.1 + prob.2)
  
  return(prob)
  
}

### Take a map and build an empty probability grid

create.grid <- function(mymap, by=0.001, precision=3){
  
  max.x <- round(mymap$range[2], precision); min.x <- round(mymap$range[1], precision)
  max.y <- round(mymap$range[4], precision); min.y <- round(mymap$range[3], precision)
  
  x.range <- seq(min.x, max.x, by=by)
  
  y.range <- seq(min.y, max.y, by=by)
  
  grid <- matrix(nrow=length(y.range), ncol=length(x.range))
  
  row.names(grid) <- y.range
  
  colnames(grid) <- x.range
  
  return(grid)
  
}

### Create color grid using filled in probability grid

create.colorgrid <- function(grid, probs = seq(0,1,length=100), col=NULL){
  
  # Define color ramp, with transparency for mapping
  if(is.null(col)){
    colr <- colorRamp(c("white", "darkgreen", "yellow", "red"))
  }else{colr <- colorRamp(col)}
  
  colr2 <- colr(probs)
  
  colr3 <- rgb(colr2[,1], colr2[,2], colr2[,3], maxColorValue=255)
  
  quantile.grid <- as.vector(grid)
  colorgrid <- grid
  
  quants <- quantile(quantile.grid, probs=probs, names=F)
  
  for(i in length(quants):1){
    colorgrid[grid<=quants[i]] <- colr3[i]
  }
  
  return(colorgrid)
}

##### Test the program with a few random crimes

burke <- map('county', 'north carolina,burke', plot=F)

grid <- create.grid(burke)

B <- 30/(70)
phi <- 0.5
f <- 1
g <- 1

set.seed(10)
xn <- sample(1:ncol(grid), 8); xn <- as.numeric(colnames(grid)[xn])
yn <- sample(1:nrow(grid), 8); yn <- as.numeric(row.names(grid)[yn])

lats <- as.numeric(row.names(grid))
longs <- as.numeric(colnames(grid))

for(j in 1:ncol(grid)){
  for(i in 1:nrow(grid)){
    
    grid[i,j] <- geo.prob(longs[j], lats[i], xn, yn, B, phi, f, g)
    
  }
}

colorgrid <- create.colorgrid(grid)

### Map county with probability grid overlaid

map(burke, mar=c(0,0,0,0))
for(j in 1:ncol(grid)){
  for(i in 1:nrow(grid)){
    xleft <- longs[j] - 0.0005
    xright <- longs[j] + 0.0005
    ybottom <- lats[i] - 0.0005
    ytop <- lats[i] + 0.0005
    rect(xleft, ybottom, xright, ytop, col=colorgrid[i,j], border=NA)
  }
}

points(x=xn, y=yn, pch=19, cex=2)
map('county', 'north carolina',  add=T, lwd=3)

### Try again, using some specific crime locations

xn <- c(-81.81418, -81.74384, -81.75856, -81.79946, -81.91070, -81.91233, -81.71603)
yn <- c(35.78089, 35.79287, 35.73427, 35.72495, 35.76091, 35.79154, 35.71562)

grid <- create.grid(burke)

B <- (10/70)
phi <- 0.5
f <- 2
g <- 1

lats <- as.numeric(row.names(grid))
longs <- as.numeric(colnames(grid))

for(j in 1:ncol(grid)){
  for(i in 1:nrow(grid)){
    
    grid[i,j] <- geo.prob(longs[j], lats[i], xn, yn, B, phi, f, g)
    
  }
}

colorgrid <- create.colorgrid(grid)

### Map county with probability grid overlaid

map(burke, mar=c(0,0,0,0))
for(j in 1:ncol(grid)){
  for(i in 1:nrow(grid)){
    xleft <- longs[j] - 0.0005
    xright <- longs[j] + 0.0005
    ybottom <- lats[i] - 0.0005
    ytop <- lats[i] + 0.0005
    rect(xleft, ybottom, xright, ytop, col=colorgrid[i,j], border=NA)
  }
}

points(x=xn, y=yn, pch=19, cex=2)
map('county', 'north carolina',  add=T, lwd=3)