## @knitr setup1 options(width=132) opts_chunk$set(cache=TRUE,cache.path="plumecache/",autodep=TRUE,fig.path="figure/plume-",fig.width=9,fig.height=9) set.seed(310366) library=function(x){suppressPackageStartupMessages(base:::library(x,character.only=TRUE))} ## @knitr setup build_dep() # figure out dependencies automatically library("xtable") ## @knitr packages library("sp") library("raster") library("rgdal") library("rgeos") ## @knitr palette library("RColorBrewer") colours = colorRampPalette(brewer.pal(7,"OrRd")[-1])(50) sp.theme(set = TRUE, regions = list(col = colours)) ## @knitr setdd Datadir = "../Data" Datafile = function(f){file.path(Datadir,f)} ## @knitr plume pfunc=function(d2,ang,a,b,k,phi){ ## plume value for distance-squared d2 a * exp(-(d2*exp(-k*cos(ang-phi))^2/b^2)) } plume <- function(src,dst,a,b,k,phi){ ## plume value for src at dst src = coordinates(src) dst = coordinates(dst) d2 = (dst[,1]-src[,1])^2 + (dst[,2]-src[,2])^2 ang = atan2(dst[,2]-src[,2],dst[,1]-src[,1]) pfunc(d2,ang,a,b,k,phi) } ## @knitr data toOSGB = function(s){spTransform(s,CRS("+init=epsg:27700"))} cumbria = toOSGB(readOGR(Datadir,"cumbria",verbose=FALSE)) nuclear = toOSGB(readOGR(Datadir,"nuclear",verbose=FALSE)) ukgaz = readOGR(Datadir,"ukgaz",verbose=FALSE) plot(cumbria) plot(nuclear,add=TRUE,pch=19,col="red") plot(ukgaz,add=TRUE) title("Base map") ## @knitr computePlume plevel = raster(extent(cumbria),nrows=500,ncols=500) plevel[] = plume(nuclear,SpatialPoints(plevel),12000,400,5,pi/4) projection(plevel)=proj4string(cumbria) plot(plevel,col=colours) plot(cumbria,add=TRUE) plot(nuclear,add=TRUE,pch=19,col="red") plot(ukgaz,add=TRUE) title("Plume Spread") ## @knitr exposure ukgaz$exposure = extract(plevel, ukgaz) over9000 = ukgaz$exposure>9000 dangerTowns = data.frame(ukgaz@data[over9000,c("Name","pop","exposure")]) dangerTowns ## @knitr plumecontour invpfunc <- function(ang,a,b,k,phi,f){ ### inverse plume function - compute distance given f d2 = -log(f/a) * b^2 / exp(-k*cos(ang-phi))^2 sqrt(d2) } contourPlume <- function(src,a,b,k,phi,f){ ### return a SpatialPolygon of a plume function at value f ### by sweeping round 360 degrees and computing the distance ang = seq(0,2*pi,len=360) d = invpfunc(ang,a,b,k,phi,f=f) xy=coordinates(src) polys = SpatialPolygons(list(Polygons(list(Polygon(cbind(xy[,1]+d*sin(ang),xy[,2]+d*cos(ang)))),ID=1))) proj4string(polys) = proj4string(src) polys } ## @knitr loadrivers rivers=toOSGB(readOGR(Datadir,"rivers",verbose=FALSE)) lakes = toOSGB(readOGR(Datadir,"naturalwater",verbose=FALSE)) nineK = contourPlume(nuclear,12000,400,5,pi/4,9000) rivers9k = gIntersection(nineK,rivers,byid=TRUE) plot(rivers9k,col="blue",lwd=3) plot(rivers,col="blue",lty=2,add=TRUE) plot(nineK,add=TRUE) plot(lakes,add=TRUE,col="blue") title("River sections with >9000 exposure") ## @knitr loadwater lakes = toOSGB(readOGR(Datadir,"naturalwater",verbose=FALSE)) plot(plevel,col=colours) plot(cumbria,add=TRUE) plot(lakes,add=TRUE,col="#0000FF") title("Lakes") ## @knitr lakesum lakes$pollution = extract(plevel,lakes,sum,small=TRUE) * prod(res(plevel)) / 1e9 ## @knitr laketable library("lattice") library("RColorBrewer") top10 = order(lakes$pollution,decreasing=TRUE)[1:10] spplot(lakes[top10,],"pollution", panel=function(...){ panel.polygonsplot(...) sp.text(coordinates(lakes[top10,]),lakes$name[top10],adj=0) },main="Lake exposure") lakes@data[top10,] ## @knitr noisyraster library("automap") noisyRaster <- function(plevel,crs,npts=300,idp=4,a=3000){ pts = spsample(SpatialPoints(plevel),npts,type="random") pts = SpatialPointsDataFrame(pts,data.frame(Z=rnorm(npts,0,1))) proj4string(pts)=CRS(proj4string(ukgaz)) fk = idw(Z~1,pts,SpatialPoints(plevel,proj4string=crs),idp=idp) fkp=plevel fkp[]=fk$var1.pred pmax = max(values(plevel)) return(fkp*a*(plevel/pmax)+plevel) } sensorValues = noisyRaster(plevel,CRS(proj4string(ukgaz))) ## @knitr writesamples samples = ukgaz[,"Name"] samples$sample = extract(sensorValues,samples) writeOGR(samples,Datadir,"samples","ESRI Shapefile",overwrite=TRUE) rm(samples) ## @knitr readSamples samples = readOGR(Datadir,"samples",verbose=FALSE) spplot(samples,"sample",main="Samples from towns/villages") ## @knitr idwfit library("automap") rIDW = raster(extent(cumbria),ncol=100,nrow=100,crs=projection(samples)) sampleFitIDW = idw(sample~1,samples,SpatialPoints(rIDW,proj4string=CRS(proj4string(samples)))) rIDW[]=sampleFitIDW$var1.pred plot(rIDW) title("Inverse Distance Weighting") ## @knitr krigefit sampleFit = autoKrige(sample~1,samples) plot(sampleFit) ## @knitr exceedencefunctions pOver = function(thresh, krigeFit){ ### convert kriging output to exceedence ### pnorm computes probabilities from a Normal/Gaussian distribution pValue = pnorm(thresh,krigeFit$krige_output$var1.pred,krigeFit$krige_output$var1.stdev,lower.tail=FALSE) ### create a raster object from the kriging coordinates and the probability pExceed = rasterFromXYZ(cbind(coordinates(krigeFit$krige_output),(pValue>0.05) + (pValue>0.95))) attr(pExceed,"thresh")=thresh return(pExceed) } eColours = c("#4040FF","grey","red") plotExceed = function(pe,overlay,region){ ### plot the exceedence with a couple of additional overlays thresh = attr(pe,"thresh") plot(pe,col=eColours,legend=FALSE) legend("topright",c(paste("Under ",thresh),"Gray Area",paste("Over ",thresh)),fill=eColours,bg="white") plot(overlay,add=TRUE) plot(region,add=TRUE) } ## @knitr exceedencecalc pK = pOver(9000,sampleFit) plotExceed(pK,nineK,cumbria) title("Exceedence Map") ## @knitr extrasamples # grey cells are == 1 greyPts = SpatialPoints(xyFromCell(pK,1:length(pK))[values(pK)==1 & !is.na(values(pK)),]) ch = chull(coordinates(greyPts)) ch = c(ch[length(ch)],ch) nNew = 25 newPts = spsample(Polygon(coordinates(greyPts[ch,])),nNew,"hexagonal") extra = nNew-length(newPts) if(extra > 0){ more = spsample(Polygon(coordinates(greyPts[ch,])),extra,"random") newPts = rbind(newPts,more) } plot(pK,col=eColours,legend=FALSE) points(newPts,pch=19) title("New Sensor Locations") ## @knitr extrasensorvalues newPts2 = SpatialPointsDataFrame(newPts,data.frame(sample=extract(sensorValues,newPts))) proj4string(newPts2)=CRS(proj4string(samples)) writeOGR(newPts2,Datadir,"sensors","ESRI Shapefile",overwrite=TRUE) rm(newPts2) ## @knitr krige2 sensors = readOGR(Datadir,"sensors") newData = rbind(sensors[,"sample"],samples[,"sample"]) sampleFit2 = autoKrige(sample~1,newData) pK2 = pOver(9000,sampleFit2) plotExceed(pK2,nineK,cumbria) title("Exceedence Map") ## @knitr showetable eTable = cbind(table(values(pK)),table(values(pK2))) rownames(eTable)=c("Low","Grey","High") colnames(eTable)=c("Settlements","Plus Sensors") eTable ## @knitr sensorsurface plot(sensorValues)