Document Sample

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

Shared By:

Categories:

Tags:
video sharing, Laboratory Equipment, file transfer, file sharing service, file hosting, free file sharing, file sharing, class definition, free web space, instance variable

Stats:

views: | 19 |

posted: | 6/25/2011 |

language: | English |

pages: | 8 |

OTHER DOCS BY keralaguest

How are you planning on using Docstoc?
BUSINESS
PERSONAL

By registering with docstoc.com you agree to our
privacy policy and
terms of service, and to receive content and offer notifications.

Docstoc is the premier online destination to start and grow small businesses. It hosts the best quality and widest selection of professional documents (over 20 million) and resources including expert videos, articles and productivity tools to make every small business better.

Search or Browse for any specific document or resource you need for your business. Or explore our curated resources for Starting a Business, Growing a Business or for Professional Development.

Feel free to Contact Us with any questions you might have.