Docstoc

lab3

Document Sample
lab3 Powered By Docstoc
					1


                                         Lab 3
                                     (1 March-2011)
                                        Maps in R
                  Spatial Statistics & Knowledge Discovery (DT876)


In this lab we will cover:
    1) Loading packages
    2) Displaying a maps
    3) Calculating Moran’s I
    4) Producing a scatter plot
    5) Moran’s I using points (ape package)
    6) Using fields package

Please refer to handout from the book Spatial Data Analysis: by LLoyd

   1) Loading packages
Check loaded package:
(.packages())
For example if the maptools library is not installed then install it with the command.
install.packages("maptools", dependencies=TRUE)

In general, make sure your machine is connected to the Internet
    Start R
    Type the following at the R prompt
    install.packages('requiredPackage', dependencies=TRUE))
    You should be prompted with a list of sites that have a copy of the package.




    For example the above screen shot shows a request to install fields package.
    Select Ireland and the package should be downloaded and installed automatically.
    You should wait until the R prompt (>) returns.

    2.1 Displaying a maps
    Display a map of Ireland. Use help to undersrand the data and function e.g.
    The eire data set is in the maptools package.
    library(spdep)
2


    help(eire)
    The following function is in the base package
    help(findInterval).

    We now display a map of Ireland. We load in a shape file which should be stored in
    your local library.
    library(spdep)
    eireMap <- readShapePoly("C:\\Documents and Settings\\pbrowne.CS\\My
    Documents\\R\\win-library\\2.12\\spdep\\etc\\shapes\\eire.shp"[1],ID="names",
    proj4string=CRS("+proj=utm +zone=30 +units=km"))
    names(eireMap)
    factor(eireMap $pale)
    eireMap $names

    summary(eireMap$A)
    res <- eireMap$A
    brks <- round(fivenum(eireMap$A), digits=2)
    cols <- rev(heat.colors(4))
    plot(eireMap, col=cols[findInterval(res, brks, all.inside=TRUE)])
    title(main="Regression residuals")
    legend(x=c(-20, 70), y=c(6120, 6050), legend=leglabs(brks), fill=cols, bty="n")
    text(coordinates(eireMap), labels=as.character(eireMap$names), cex=0.4)


2.2 Displaying more maps
Copy all files are stored in offaly     student distrib, MY-R-Dir.
to your own machine:
C:\MY-R-Dir\
Start R and load in the county file
library(spdep)
county <- readShapePoly("C:\\MY-R-Dir\\county_REGION.shp")
plot(county)
text(coordinates(county), labels=as.character(county$NAME), cex=0.2)
Note that readShapePoly is from maptools package.


3) Calculating Moran’s I
The following calculates the Moran’s I statistic for two girds.
Please see section 4.8 of handout.
Note that readShapePoly is from maptools package.

setwd("C:\\My-R-Dir")
library(spdep)
getinfo.shape("lab3a.shp")
getinfo.shape("lab3b.shp")
lab3a <- readShapePoly("C:\\MY-R-Dir\\lab3a.shp")
lab3b <- readShapePoly("C:\\MY-R-Dir\\lab3b.shp")

#Examine the components of lab3a
names(lab3a)
is(lab3a$rainfall)
lab3a$rainfall
coordinates(lab3a)
getClass(lab3a)
dim(lab3a)
slotNames(lab3a)
lab3a@data
lab3a@polygons
3

#Examine the statistical properties of labd3a
summary(lab3a$rainfall)
summary(lab3b$rainfall)
sd(lab3a$rainfall)
sd(lab3b$rainfall)
stem(lab3a$rainfall)

# View both maps
plot(lab3a)
text(coordinates(lab3a), labels=as.character(lab3a$rainfall), cex=1.0)
title("Lab 3a")
plot(lab3b)
text(coordinates(lab3b), labels=as.character(lab3b$rainfall), cex=1.0)
title("Lab 3b")

Both maps have the same labelling




# First use rooks neighbours
lab3a.nb <- poly2nb(lab3a,queen = FALSE)
lab3b.nb <- poly2nb(lab3b,queen = FALSE)
# Look at the connectivity
plot(lab3a.nb, coordinates(lab3a), add=TRUE, col="red")
# Note the results in both cases
# and check against values in section 4.8 of handout
moran.test(lab3a$rainfall, nb2listw(lab3a.nb))
moran.test(lab3b$rainfall, nb2listw(lab3b.nb))

# Now use queens neighbours
lab3a.nb <- poly2nb(lab3a,queen = TRUE)
lab3b.nb <- poly2nb(lab3b,queen = TRUE)
# Look at the connectivity
plot(lab3a.nb, coordinates(lab3a), add=TRUE, col="red")
# Note the results in both cases
# and check against values in section 4.8 of handout
moran.test(lab3a$rainfall, nb2listw(lab3a.nb))
moran.test(lab3b$rainfall, nb2listw(lab3b.nb))
4


4) Moran’s I scatter plot
Assuming section 3 has been completed we can now construct a Morans scatter plot.
The spatial lag of the variable on the vertical axis and the original variable on the
horizontal axis.
The spatial lag refers to the values of a location's neighbours, see Lecture 3.
moran.plot(lab3a$rainfall, listw=nb2listw(lab3a.nb,style="C"))
moran.plot(lab3b$rainfall, listw=nb2listw(lab3b.nb,style="C"))




Notice that there is practically no autocorrelation present in map lab3a.

The points with greater influence can be labelled and plotted as follows:
moran.plot(lab3a$rainfall, labels=lab3a$name,
                    listw=nb2listw(lab3a.nb,style="C"))
moran.plot(lab3b$rainfall, labels=lab3b$name,
                   listw=nb2listw(lab3b.nb,style="C"))
Run these commands and study the results.


3) Calculating Local Moran’s I
Assuming section 3 has been completed we can now calculate the local Morans for each
the value at each location.

resI <- localmoran(lab3a$rainfall, nb2listw(lab3a.nb))

The components of resI are
# Ii local moran statistic for each location i
# E.Ii expectation of local moran statistic
# Var.Ii variance of local moran statistic
# Z.Ii standard deviate of local moran statistic
# Pr() p-value of local moran statistic

Lets examine Z.Ii which is a “standard score”.
It is the difference between an actual value and the mean divided by the standard
deviation.
5




where
  x is the Moran I for each location
  μ is the mean of all Moran I’s;
  σ is the standard deviation of all Moran I’s.

As an example we calculate the Z value of the first point

(0.87886598+0.125)/sqrt(0.3788857)


These values were drawn from the first row of resI
The result is 1.630879 which agrees with Z.Ii below.

> View(resI)




4) Generating random trees (ape package)
library(ape)
layout(matrix(1:9, 3, 3))
# Nine random trees:
for (i in 1:9) plot(rtree(20))
# Nine random cladograms:
for (i in 1:9) plot(rtree(20, FALSE), type = "c")
#generate 4 random trees of bird orders:
data(bird.orders)
layout(matrix(1:4, 2, 2))
for (i in 1:4)
  plot(rcoal(23, tip.label = bird.orders$tip.label), no.margin = TRUE)
layout(matrix(1))



5) Moran’s I using points (ape package)

library(ape)
ozone.dists <- as.matrix(dist(cbind(ozone$Lon, ozone$Lat)))
ozone.dists.inv <- 1/ozone.dists
diag(ozone.dists.inv) <- 0
ozone.dists.inv[1:5, 1:5]
6

Moran.I(ozone$Av8top, ozone.dists.inv)

Another example using trees
tr <- rtree(30)
x <- rnorm(30)
# weights w[i,j] = 1/d[i,j]:
w <- 1/cophenetic(tr)
# set the diagonal w[i,i] = 0 (instead of Inf...):
diag(w) <- 0
Moran.I(x, w)
Moran.I(x, w, alt = "l")
Moran.I(x, w, alt = "g")
Moran.I(x, w, scaled = TRUE) # usualy the same


6. Running field example
At the R command prompt type:
library(fields)
help(fields)
This will open your defualt browser and display the help page for the fields package.
We will use the ozone data set, which is supplied with the field package, there is no need
to load the data. For further details type:
 help(ozone2)
When a plot is created it can be viewed along side the R main screen by selecting Tile
Verically as below.




library(fields)
# some air quality data, daily surface ozone measurements for the Midwest USA
data(ozone2)

# The response is 8-hour average (surface) ozone ( from 9AM-4PM) measured in parts per
# billion (PPB) for 153 sites in the midwestern US over the period June 3,1987 through
# August 31, 1987, 89 days. This season of high ozone corresponds with a large modeling
7

#   experiment using the EPA Regional Oxidant Model.
#   Ground-level ozone, though less concentrated than athmospheric ozone,
#   is more of a problem because of its health effects.
#   Typical peak concentrations in American cities is 0.15 to 0.51 ppm


x<-ozone2$lon.lat
y<- ozone2$y[16,] # June 18, 1987
# pixel plot of spatial data
quilt.plot( x,y)
# add map
US( add=TRUE) # add US map
# fits a thin plate smoothing spline surface (tps) to ozone measurements.
fit<- Tps(x,y)
summary(fit) #diagnostic summary of the fit
set.panel(2,2)
plot(fit) # four diagnostic plots of fit and residuals.
set.panel()
surface(fit) # contour/image plot of the fitted surface
US( add=TRUE, col="magenta", lwd=2) # US map overlaid
title("Daily max 8 hour ozone in PPB, June 18th, 1987")




install.packages('GeoXp',dependencies=TRUE)
library(GeoXp)
moranplotmap(baltimore.spdf, "PRICE", W.listw ,flower=TRUE,
locmoran=TRUE,criteria=(baltimore.spdf$AC==1), identify=TRUE)
8

				
DOCUMENT INFO