Toni — Apr 1, 2013, 2:00 PM
library(spatstat)
Loading required package: mgcv
This is mgcv 1.7-22. For overview type 'help("mgcv-package")'.
Loading required package: deldir
deldir 0.0-21
spatstat 1.31-0 Type 'help(spatstat)' for an overview of spatstat
'latest.news()' for news on latest version 'licence.polygons()' for
licence information on polygon calculations
library(sp)
library(rgdal)
rgdal: version: 0.8-4, (SVN revision 431) Geospatial Data Abstraction
Library extensions to R successfully loaded Loaded GDAL runtime: GDAL
1.9.2, released 2012/10/08 Path to GDAL shared files:
Q:/RStu97.64b/R/R-2.15.2/library/rgdal/gdal GDAL does not use iconv for
recoding strings. Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009,
[PJ_VERSION: 470] Path to PROJ.4 shared files:
Q:/RStu97.64b/R/R-2.15.2/library/rgdal/proj
library(maptools)
Loading required package: foreign
Loading required package: grid
Loading required package: lattice
Checking rgeos availability: TRUE
library(ggmap)
Loading required package: ggplot2
id <- read.csv(file.choose())
Error: more columns than column names
summary(id)
Error: error in evaluating the argument 'object' in selecting a method for
function 'summary': Error: object 'id' not found
id$Type <- sub(id$Type, pattern = "CRIMINAL EVENT", replacement = "Criminal Event")
Error: object 'id' not found
id$Type <- sub(id$Type, pattern = "criminal event", replacement = "Criminal Event")
Error: object 'id' not found
id$Type <- sub(id$Type, pattern = "EXPLOSIVE HAZARD", replacement = "Explosive Hazard")
Error: object 'id' not found
#In order to be sure the old categories are removed from the list of possible event types we have to run:
id$Type <- factor(id$Type)
Error: object 'id' not found
levels(id$Type)
Error: object 'id' not found
#The spatstat library is for the analysis and manipulation of spatial points. In a spatstat point file the attributes of each location are stored in âmarksâ. Each point on the map can have zero, one, or many associated marks. A file with too many marks is cumbersome, select a subset of attributes to serve as potential marks:
id.marks.simple <- id[, 4]
Error: object 'id' not found
id.marks.full <- id[, c(4, 7, 9:18)]
Error: object 'id' not found
#If you looked carefully at summary(id) you noticed that latitude and longitude are attributes in the table, we'll use these to map out the events:
id.spp <- SpatialPointsDataFrame(coords = id[, 19:20], proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"),
data = id.marks.full)
Error: object 'id' not found
###EXPLORATORY MAPS#### To make maps we'll use a very new package called ggmaps that allows us to send a query to online mapping services like Google maps. It returns a georeferenced image of the specified type and scale. We can then plot on top of that image. The code below uses a package called ggplot2 which we'll look at more carefully later in the semester. The first step is to get the base map.
map <- get_map(location = "Iraq", maptype = "hybrid", zoom = 6) #get map
ggmap(map, extent = "panel")
#overlay wikileaks war data:
ggmap(map, extent = "panel") + geom_point(data = id, aes(x = id$Longitude, y = id$Latitude),
alpha = 0.1, color = "red") + ggtitle("WikiLeaks Iraq Data \n Incidents involving Death(s)")
Error: object 'id' not found
#Alter the map the so that the size of the dots are linked to the total number of deaths in each events (note the use of size inside of aes()):
ggmap(map, extent = "panel") + geom_point(data = id, aes(x = id$Longitude, y = id$Latitude,
size = Total.deaths), alpha = 0.3, color = "red") + ggtitle("WikiLeaks Iraq Data \n Incidents involving Death(s)")
Error: object 'id' not found
#Map a subset of events, for example, those involving detentions:
#this plots 1 layer which is the Google map and then overlays a set of points, the data for the points is the data.frame id. We then map specific columns in id to visual variables, such as size and position using aes(). I'm a big fan of ggplot2's functionality for visualization (even though it kinda sucks for mapping).
ggmap(map, extent = "panel") + geom_point(data = id[id$Enemy.detained > 0, ],
aes(x = Longitude, y = Latitude, size = Enemy.detained), alpha = 0.5, color = "green") +
ggtitle("WikiLeaks Iraq Data \n Incidents involving Dententions(s)")
Error: object 'id' not found
######EXPLORATORY SPATIAL ANALYSIS#############
#The logic of spatstat is a little odd from a geographic perspective, but once you understand how it works it is a really useful and robust library.
#Many of the analytic tools in spatstat are based on distance. Working with decimal degrees is a pain, they're difficult to interpret as distance units because 1 degree is a different distance on different parts of the globe. Re-project the data to a coordinate system that uses meters instead of degrees of latitude/longitude. Knowing how many deaths occurred per decimal degree isn't that useful.
## Convert to a different coordinate system. So that coordinates are in
## meters not degrees
id.spp.rp <- spTransform(id.spp, CRS("+init=epsg:3839"))
Error: error in evaluating the argument 'x' in selecting a method for
function 'spTransform': Error: object 'id.spp' not found
id.spp.rp$Easting <- coordinates(id.spp.rp)[, 1]
Error: error in evaluating the argument 'obj' in selecting a method for
function 'coordinates': Error: object 'id.spp.rp' not found
id.spp.rp$Northing <- coordinates(id.spp.rp)[, 2]
Error: error in evaluating the argument 'obj' in selecting a method for
function 'coordinates': Error: object 'id.spp.rp' not found
## Remove all events outside of Baghdad and convert back to a table
id.df.rp.bag <- id.spp.rp@data[id.spp.rp@data$Region == "MND-BAGHDAD", ]
Error: object 'id.spp.rp' not found
#Notice my use of @data I have pulled out a subset of the map and saved it as a non-spatial data.frame. The spatial information is preserved in the columns Easting and Northing. We can still make a âmapâ of the points:
plot(x = id.df.rp.bag$Easting, y = id.df.rp.bag$Northing)
Error: error in evaluating the argument 'x' in selecting a method for
function 'plot': Error: object 'id.df.rp.bag' not found
#remove outliers
id.df.rp.bag <- id.df.rp.bag[id.df.rp.bag$Northing < 4960000 & id.df.rp.bag$Easting >9960000, ]
Error: object 'id.df.rp.bag' not found
plot(x = id.df.rp.bag$Easting, y = id.df.rp.bag$Northing)
Error: error in evaluating the argument 'x' in selecting a method for
function 'plot': Error: object 'id.df.rp.bag' not found
#zoom into the region with dense events:
id.df.rp.bag <- id.df.rp.bag[(id.df.rp.bag$Northing < 4950000 & id.df.rp.bag$Northing > 4940000) & (id.df.rp.bag$Easting > 1e+07 & id.df.rp.bag$Easting < 10010000), ]
Error: object 'id.df.rp.bag' not found
#spatstat objects:
#ppp - set of points located in space. Points can have attributes called "marks" or the ppp object can be "unmarked" in which case the object stores only locations. Marks are used to differentiate types of points or to store a variable that describes some continuous variable for each point.
#owin - the WINDOW which defines the spatial extent of the ppp object. all ppp objects must have an owin
#im - a pixel image like a raster layer in GIS often derived from the ppp object.
#step 1 for ppp spatial analysis is define the window.
#Define the window using convexhull.xy() which traces the outermost points of a scatter plot. Give the function the x and y coordinates of the points, and it returns an owin object â which is necessary for making a point pattern objectâ¦
id.chull <- convexhull.xy(x = id.df.rp.bag$Easting, y = id.df.rp.bag$Northing)
Error: object 'id.df.rp.bag' not found
#draw window manually - must draw COUNTER CLOCKWISE
plot(x = id.df.rp.bag$Easting, y = id.df.rp.bag$Northing)
Error: error in evaluating the argument 'x' in selecting a method for
function 'plot': Error: object 'id.df.rp.bag' not found
bdry <- locator()
Error: plot.new has not been called yet
id.win <- owin(poly = bdry)
Error: object 'bdry' not found
#create window
id.ppp <- ppp(x = id.df.rp.bag$Easting, y = id.df.rp.bag$Northing, window = id.chull,
marks = id.df.rp.bag[, c(1, 3:12)])
Error: object 'id.chull' not found
unitname(id.ppp) <- c("meter", "meters") ##Assign names to units
Error: object 'id.ppp' not found
# To change the window of an existing ppp object
id.ppp <- id.ppp[id.win]
Error: object 'id.ppp' not found
#visualize:
p1 <- plot.ppp(id.ppp)
Error: object 'id.ppp' not found
legend(max(coords(id.ppp))[1] + 1000, mean(coords(id.ppp))[2], pch = p1, legend = names(p1),
cex = 0.5)
Error: object 'p1' not found
#BY DEFAULT, the first "mark" is used to label the points.
# This mark was the âTypeâ column in the original data, it describes the type of incident. The data are so dense that the resulting graph is impossible to read. There are alternative visualizations:
## Plot the density of points using a gaussian (normal) kernel with a
## 1000m standard deviation
plot(density(id.ppp, sigma = 1000))
Error: error in evaluating the argument 'x' in selecting a method for
function 'plot': Error in density(id.ppp, sigma = 1000) : object 'id.ppp'
not found
#This density plot places a kernel over the center of each pixel and counts up the number of points per unit area. The points further away from the center of each pixel count less toward the density of the pixel than points that are close to the center. This is the standard kernel density estimator that is used in most GIS programs.
#The density command creates a pixel image or im object. These objects allow us to do cool stuff when combined with points, as we will see laterâ¦
#We can make a contour map of the density:
contour(density(id.ppp, 1000), axes = F)
Error: object 'id.ppp' not found
#3D
layout(matrix(c(1, 1, 1, 2, 3, 4), 2, 3, byrow = TRUE))
persp(density(id.ppp, 1000))
Error: object 'id.ppp' not found
persp(density(id.ppp, 1000), theta = 150, phi = 15) #rotate by adjusting theta
Error: object 'id.ppp' not found
persp(density(id.ppp, 1000), theta = 120, phi = 15) #rotate by adjusting theta
Error: object 'id.ppp' not found
persp(density(id.ppp, 1000), theta = 100, phi = 15) #rotate by adjusting theta
Error: object 'id.ppp' not found
dev.off() #clearlayout
null device
1
#visualize each type of event
plot(split(id.ppp))
Error: error in evaluating the argument 'x' in selecting a method for
function 'plot': Error in split(id.ppp) : error in evaluating the argument
'x' in selecting a method for function 'split': Error: object 'id.ppp' not
found
plot(density(split(id.ppp)))
Error: error in evaluating the argument 'x' in selecting a method for
function 'plot': Error in split(id.ppp) : error in evaluating the argument
'x' in selecting a method for function 'split': Error: object 'id.ppp' not
found Calls: density -> split
contour(density(split(id.ppp)), axes = F)
Error: error in evaluating the argument 'x' in selecting a method for
function 'split': Error: object 'id.ppp' not found
#map marks other than "type"
plot(id.ppp, which.marks = "Total.deaths", maxsize = 3000, main = "Total Deaths")
Error: error in evaluating the argument 'x' in selecting a method for
function 'plot': Error: object 'id.ppp' not found
#############EXPLORATORY ANALYSIS OF POINTS###############
#SPATSTAT can use pixel images like density maps in the analysis of points.
#here's a weighted density map that reflects the total number of deaths at each location in the window.
#start with ppp object.
id.ppp.simp <- unmark(id.ppp) #remove all data from the points
Error: object 'id.ppp' not found
deaths <- marks(id.ppp)["Total.deaths"] #create a list of # of deaths
Error: object 'id.ppp' not found
id.ppp.simp.death <- id.ppp.simp %mark% deaths #make a point pattern (ppp) with a single variable
Error: object 'deaths' not found
#The density function has a lot of options, for example we can weight each point by the total deaths using the âweightâ option, we can specify the type of kernel (using kernel), and the size of the pixels (using eps):
id.death.im <- density(id.ppp.simp.death, sigma = 1000, weights = marks(id.ppp.simp.death),
kernel = "rectangular", eps = 100) #eps sets pixel size in map units
Error: object 'id.ppp.simp.death' not found
plot(id.death.im, main = "Weighted Density Map of \nWikiLeaks Deaths in Baghdad")
Error: error in evaluating the argument 'x' in selecting a method for
function 'plot': Error: object 'id.death.im' not found
# We can change the colors by making a custom color ramp
grRd <- colorRampPalette(c("grey60", "red"))
co <- colourmap(grRd(5), range = c(3.37e-06, 0.00036)) #range from summary(id.death.im)
plot(id.death.im, main = "Weighted Density Map of \nWikiLeaks Deaths in Baghdad",
col = co)
Error: error in evaluating the argument 'x' in selecting a method for
function 'plot': Error: object 'id.death.im' not found
#map algebra is possible. also, export of images to a table, such taht each pixel becomes a cell in the table.
#The Death Density image can be used to divide the study region into areas of high, medium, and low density.
# the function tess is used to create a special type of spatstat object called a tessellation
########tesselations can be used to group points.
death.map.breaks <- quantile(id.death.im, probs = (0:3)/3) #define density cut points
Error: object 'id.death.im' not found
death.map.cuts <- cut(id.death.im, breaks = death.map.breaks, labels = c("Low DD",
"Medium DD", "High DD")) #cut the image
Error: object 'id.death.im' not found
death.map.quartile <- tess(image = death.map.cuts)
Error: object 'death.map.cuts' not found
plot(death.map.quartile)
Error: error in evaluating the argument 'x' in selecting a method for
function 'plot': Error: object 'death.map.quartile' not found
#We can also use the tessellation to count the number of events in each category on the map:
qCount.im <- quadratcount(X = id.ppp, tess = death.map.quartile)
Error: object 'id.ppp' not found
plot(qCount.im)
Error: error in evaluating the argument 'x' in selecting a method for
function 'plot': Error: object 'qCount.im' not found
#We could also divide the window up into a n by n grid and count the number of events in each grid cell:
qCount <- quadratcount(X = id.ppp, nx = 4, ny = 4) #returns the count of points in each quadrant
Error: object 'id.ppp' not found
plot(qCount)
Error: error in evaluating the argument 'x' in selecting a method for
function 'plot': Error: object 'qCount' not found
#We can also use the density map to look at how the intensity of events (number of points per unit area) relates to the density of death image:
plot(rhohat(unmark(id.ppp), id.death.im))
Error: error in evaluating the argument 'x' in selecting a method for
function 'plot': Error in unmark(id.ppp) : object 'id.ppp' not found
Calls: rhohat -> unmark
#The above plot basically shows the density of points as a function of the values on the weighted death density map.
#The yellow HighDD regions on the map contain not only more points, but the points are higher density than parts of the death density of image that have lower numbers of death.
#By "points" here, I think Seth means "EVENTS" again.
id.ppp.types <- split(id.ppp, which.marks = "Type")
Error: error in evaluating the argument 'x' in selecting a method for
function 'split': Error: object 'id.ppp' not found
id.ppp.types <- unmark(id.ppp.types)
Error: object 'id.ppp.types' not found
#the object id.ppp.types now holds a bunch of unmarked point patterns, each type of events is in its own ppp object.
id.ppp.types$"Criminal Event"
Error: object 'id.ppp.types' not found
# filter the significant clusters (high density areas) from the background rate (noise). We'll do this using an approach developed by Byers and Raferty (1998) called Nearest Neighbor Cleaning. We'll discuss the theory in class, the basic idea is to estimate the parameters of two point generating processes (one for noise and one for signal). Once the parameters of each process are estimated points are assigned to one of the two processes.
nnc <-nnclean(X = id.ppp.types$"Criminal Event", k=20)
Error: object 'id.ppp.types' not found
plot(split(nnc))
Error: error in evaluating the argument 'x' in selecting a method for
function 'plot': Error in split(nnc) : error in evaluating the argument
'x' in selecting a method for function 'split': Error: object 'nnc' not
found
#The map above shows intense clusters of criminal activity against the background ânoise.â In the next exercise we'll do hypothesis tests to understand the processes that generated the observed point patterns.