date()
## [1] "Tue Feb 12 21:57:36 2013"
Use the Columbus, OH crime data from “http://geodacenter.org/downloads/data-files/columbus.zip”.
1.Read the data into R and create a choropleth map of the CRIME variable using the ssplot method.
suppressWarnings(require(maptools, quietly = TRUE))
## Checking rgeos availability: TRUE
tmp = download.file("http://geodacenter.org/downloads/data-files/columbus.zip",
"columbus.zip", mode = "wb")
unzip("columbus.zip")
CO = readShapeSpatial("columbus", proj4string = CRS("+proj=longlat +ellps=clrk66"))
slotNames(CO)
## [1] "data" "polygons" "plotOrder" "bbox" "proj4string"
str(CO, max.level = 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 49 obs. of 20 variables:
## .. ..- attr(*, "data_types")= chr [1:20] "N" "N" "N" "N" ...
## ..@ polygons :List of 49
## ..@ plotOrder : int [1:49] 21 9 5 20 1 40 6 42 2 7 ...
## ..@ bbox : num [1:2, 1:2] 5.87 10.79 11.29 14.74
## .. ..- attr(*, "dimnames")=List of 2
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
head(CO@data)
## AREA PERIMETER COLUMBUS_ COLUMBUS_I POLYID NEIG HOVAL INC CRIME
## 0 0.30944 2.441 2 5 1 5 80.47 19.531 15.73
## 1 0.25933 2.237 3 1 2 1 44.57 21.232 18.80
## 2 0.19247 2.188 4 6 3 6 26.35 15.956 30.63
## 3 0.08384 1.428 5 2 4 2 33.20 4.477 32.39
## 4 0.48889 2.997 6 7 5 7 23.23 11.252 50.73
## 5 0.28308 2.336 7 8 6 8 28.75 16.029 26.07
## OPEN PLUMB DISCBD X Y NSA NSB EW CP THOUS NEIGNO
## 0 2.8507 0.2172 5.03 38.80 44.07 1 1 1 0 1000 1005
## 1 5.2967 0.3206 4.27 35.62 42.38 1 1 0 0 1000 1001
## 2 4.5346 0.3744 3.89 39.82 41.18 1 1 1 0 1000 1006
## 3 0.3944 1.1869 3.70 36.50 40.52 1 1 0 0 1000 1002
## 4 0.4057 0.6246 2.83 40.01 38.00 1 1 1 0 1000 1007
## 5 0.5631 0.2541 3.78 43.75 39.28 1 1 1 0 1000 1008
require(RColorBrewer)
## Loading required package: RColorBrewer
## Warning: package 'RColorBrewer' was built under R version 2.15.2
range(CO$CRIME)
## [1] 0.1783 68.8920
rng = seq(0, 70, 10)
cls = brewer.pal(7, "Reds")
spplot(CO, "CRIME", col.regions = cls, at = rng)
2.Compute Moran's I for the CRIME variable.
suppressWarnings(require(spdep, quietly = TRUE))
## Attaching package: 'boot'
## The following object(s) are masked from 'package:lattice':
##
## melanoma
## deldir 0.0-21
CO.nb = poly2nb(CO)
CO.wts = nb2listw(CO.nb)
m = length(CO$CRIME)
s = Szero(CO.wts)
moran(CO$CRIME, CO.wts, n = m, S0 = s)
## $I
## [1] 0.5002
##
## $K
## [1] 2.226
##
3.Create a Moran's scatter plot using the CRIME variable. Label the axes and include the regression line.
crime = CO$CRIME
Ccrime = lag.listw(CO.wts, crime)
require(ggplot2)
## Loading required package: ggplot2
dat = data.frame(crime = crime, Ccrime = Ccrime)
ggplot(dat, aes(x = crime, y = Ccrime)) + geom_point() + geom_smooth(method = "lm") +
xlab("crime") + ylab("Spatial Lag of crime")
4.Check the value of Moran's I by finding the slope of the regression line.
lm(Ccrime ~ crime, data = dat)
##
## Call:
## lm(formula = Ccrime ~ crime, data = dat)
##
## Coefficients:
## (Intercept) crime
## 17.5 0.5
##