R-tutorial-examples.R

setwd("~/Desktop/R-tutorial/")

####
# Mapping example
####

# load packages
require(rgdal)
require(raster)
require(leaflet)
require(htmlwidgets)
require(spatstat)

# read in a shapefile
theevents = readOGR(".","events-burglary-residential-2014")

plot(theevents)

# make a raster

boundary = readOGR(".", "city_limits")
boundary.bbox = bbox(boundary)
boundary.width = boundary.bbox[1,2]-boundary.bbox[1,1]
boundary.height = boundary.bbox[2,2]-boundary.bbox[2,1]

# set raster width
kCellSize = 250

raster.cols = ceiling(boundary.width / kCellSize)
raster.width = raster.cols * kCellSize
raster.rows = ceiling(boundary.height / kCellSize)
raster.height = raster.rows * kCellSize
raster.extent = extent(c(xmin=boundary.bbox[1,1],xmax=boundary.bbox[1,1]+raster.width,
                         ymin=boundary.bbox[2,1],ymax=boundary.bbox[2,1]+raster.height))
# make raster
raster.overall = raster(extent(raster.extent), 
                        nrow=raster.rows, ncol=raster.cols,
                        crs=CRS(proj4string(boundary)))

# rasterize the points
theevents.justthelocations = as(theevents, "SpatialPoints")
points.raster = rasterize(theevents.justthelocations, raster.overall, fun="count")

plot(points.raster)

# density

points.ppp = ppp(coordinates(theevents.justthelocations)[,1], coordinates(theevents.justthelocations)[,2],
                 window=owin(xrange=bbox(raster.overall)[1,], yrange=bbox(raster.overall)[2,]))

points.raster.density = raster(density(points.ppp, dimyx=c(nrow(raster.overall), ncol(raster.overall)), sigma=1000))
projection(points.raster.density) = projection(raster.overall)

plot(points.raster.density)

# make a web map of it

# save a Leaflet map of sig forecasts
raster.palette = colorNumeric(c("transparent",  "yellow", "orange", "red"), values(points.raster.density), na.color = "transparent")

raster.leaflet = leaflet() %>% addProviderTiles("Stamen.Toner") %>%
  addRasterImage(points.raster.density, colors = raster.palette, opacity = 0.75) %>%
  addLegend(position = "bottomright", pal = raster.palette, values = values(points.raster.density), title = "Density of Res Burgs")

# saving a widget as a self-contained HTML file requires the pandoc executable to be installed.
saveWidget(raster.leaflet, "densitymap.html", selfcontained=TRUE)



#####
# Model example
#####

# load packages
require(rpart)

warrants = read.csv("Warrants.csv")

str(warrants)


model = rpart(Resistance ~ ., data=warrants)

plot(model)
text(model, use.n = TRUE)