# R scipr to create a risk surface from a Risk Terrain Model calculated using
# the RTMDx software.

# Set this to your current directory:
setwd("/Users/nick/Google Drive/research/CrimeCourse/risk-terrain/results")

# Load the required libraries
library(GISTools)
library(rgdal) # For reading shapefiles
library(OpenStreetMap) # For basemaps
library(classInt) # For creating class levels for maps
library(RColorBrewer) # For creating colour ranges
library(rgeos) # For spatial operations (e.g. intersect)
library(sp) # For creating points in polygons

cell.size <- 250 # Cell size of the output raster (meters)

# Bandwidth used to create the kernel density in the detached houses input
det.bw <- 500

# Distance threshold for the bus stops data
bus.threshold <- 500

# Coefficients of the input layers
det.coef <- 2.1578
bus.coef <- 3.0617

# Other model variables
intercept <- -3.3565

# Read in the input data (bus stops and detached houses)
det.points <- readOGR(dsn = "../rtm_data/", "det_points" ) # Detached houses
bus.points <- readOGR(dsn = "../rtm_data/", "bus_stops" ) # Bus stops


# KDE funtion needs to know number of cells in which to create the raster (e.g. will create n x n grid)
leeds.district <- readOGR(dsn = "../rtm_data/", "leeds_district" ) 
bb <- bbox(leeds.district)
width <- bb[1,2] - bb[1,1]
height <- bb[2,2] - bb[2.1]
length <- max(c(width,height))
num.cells <- ceiling(length / cell.size)

# Create the input layers. The function below takes the input point layers and
# returnd raster and returns the raster risk surface (either as density or proximity types)
calc.risk <- function(x, type, bandwidth) {
  
  if (type == "DENSITY") { # Calculates risk for density inputs. This is 1 for all
    # cells which are more than 2 sd above the mean.
    # (Note, need to convert from SpatialPixelsDataFrame (kde output) to RasterLayer)
    kde.surface = as(kde.points(x, h = bandwidth, n=num.cells, lims=leeds.district), "RasterLayer")
    mean <- cellStats(kde.surface, 'mean')
    sd <- cellStats(kde.surface, 'sd')
    # Each cell should be 1 for areas 2 standard deviations above the mean, 0 otherwise
    y <- calc(kde.surface, function(a) {  if( a > (mean+2*sd) ) 1 else 0 })
    return(y)
  }
  
  else if (type == "PROXIMITY") { # Calculates risk for proximity inputs. This is 1
    # for all cells that are within the threshold distance ('bandwidth')
    # (Note:use kde function as above to make sure that the extents will be the same)
    r = as(kde.points(x, h = bandwidth, n=num.cells, lims=leeds.district), "RasterLayer")
    r[] <- 0
    #r <- raster(res=cell.size,xmn=bb[1,1],ymn=bb[2,1],xmx=bb[1,2],ymx=bb[2,2]) # Empty raster
    r <- rasterize(x, r, field=1) # rasterise the points
    r <-  buffer(r, width=bandwidth) # buffer the points
    # Convert nas to 0
    r[is.na(r)] <- 0
    return(r)
  }
  
}

det.risk <- calc.risk(det.points, type="DENSITY", bandwidth=det.bw)
bus.risk <- calc.risk(bus.points, type="PROXIMITY", bandwidth=bus.threshold)

# Finally combine the risk layers into a Risk Terrain Model
all.risk <- stack(det.risk, bus.risk) # Creates a RasterStack

rtm <- calc(all.risk, fun=function(x) 
  # Here each x is a single raster cell. Each cell is represented as an array with
  # three elements (one for the value of each input layer)
  exp(intercept + det.coef*x[1] + bus.coef*x[2] ) / exp(intercept)
)

# Plot the Risk Terrain Model:

par(mfrow=c(1,1))
plot(rtm, main="Risk Terrain Model")