Estimating Exposure Example
r-INLA SPDE
date()## [1] "Mon Oct 31 13:02:16 2016"Options
library(knitr)
opts_knit$set(verbose = FALSE)Random Points
Generate fake observations and convert to SpatialPointsDataFrame. A value of “1” (one) is recorded for these fake “presence” locations. Plot the result.
library(sp)
Points.df = as.data.frame(
cbind(X = runif(50),
Y = runif(50)))
Points.df$OBS = 1 #Observation locations given a count of "1" (one).
Points.df$Exposure = 0 #Exposure for Observation locations is set to "0" (zero).
LL = Points.df[, c("X","Y")]
Observations = SpatialPointsDataFrame(LL, Points.df)
plot(Observations, pch=19, col="red")
Fake Study Domain
A bounded area to serve as a fake study domain. We’re not using projected data here, so we get a warning making us aware that the data isn’t projected.
library(rgeos)## rgeos version: 0.3-19, (SVN revision 524)
## GEOS runtime version: 3.5.0-CAPI-1.9.0 r4084
## Linking to sp version: 1.2-3
## Polygon checking: TRUEDomain = gConvexHull(Observations)
plot(Domain, col = "grey", main = "Fake Study Domain")
plot(Observations, pch=19, col="red", add=TRUE)
Construct a mesh
The mesh is constructed with an exaggerated “offset” buffer for visualization purposes. Note that all points are captured within the mesh; but, not all points fall on a mesh vertex.
library(INLA)## Loading required package: Matrix## Loading required package: splines## This is INLA 0.0-1463562937, dated 2016-05-18 (11:01:03+0200).
## See www.r-inla.org/contact-us for how to get help.Mesh0 = inla.mesh.2d(loc = LL, #Adjust mesh vertex to occurence loactions
cutoff = 0.1, #Smoothing edge
max.edge = 0.8, #Eliminate small triangles
offset = 0.50) #Extend mesh beyond points - exaggerated
plot(Mesh0, sub = "Fake Observation Points")
plot(Observations, pch=19, col="red", add=TRUE)
Mesh Verticies
Mesh vertices are selected as integration points. Here we collect the coordinates of these vertices and convert them to a SpatialPointsDataFrame.
Mesh0$n #Number of vertices## [1] 52Integration = as.data.frame(cbind(Mesh0$loc[,1],
Mesh0$loc[,2]))
names(Integration) = c("X", "Y")
LL = Integration[, c("X","Y")]
Integration = SpatialPointsDataFrame(LL, Integration)
Integration$OBS = 0 #Integration locations given a count of "1" (one).
plot(Mesh0, sub = "Integration Points (blue) - at verticies")
plot(Integration, pch=19, col="blue", add=TRUE)
Estimate Exposure
To estimate exposure at the Integration locations we don’t use r-INLA’s “mesh creation” functions; instead, we tessellate using a different function. Here’s the function:
###Fucntion for Voronoi Tessilation
voronoipolygons = function(layer) {
require(deldir)
crds = layer@coords
z = deldir(crds[,1], crds[,2])
w = tile.list(z)
polys = vector(mode='list', length=length(w))
require(sp)
for (i in seq(along=polys)) {
pcrds = cbind(w[[i]]$x, w[[i]]$y)
pcrds = rbind(pcrds, pcrds[1,])
polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i))
}
SP = SpatialPolygons(polys)
voronoi = SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1],
y=crds[,2], row.names=sapply(slot(SP, 'polygons'),
function(x) slot(x, 'ID'))))
}Tessilation around Integration Points
Each bounded polygon encloses all “theoretical points” on the plane that share the integration point (at center) as their closest point.
Int.poly = voronoipolygons(Integration)## Loading required package: deldir## deldir 0.1-12##
## PLEASE NOTE: The components "delsgs" and "summary" of the
## object returned by deldir() are now DATA FRAMES rather than
## matrices (as they were prior to release 0.0-18).
## See help("deldir").
##
## PLEASE NOTE: The process that deldir() uses for determining
## duplicated points has changed from that used in version
## 0.0-9 of this package (and previously). See help("deldir").plot(Int.poly, main = "Tessilation AROUND Integration Points")
plot(Integration, pch=19, col="blue", add=TRUE)
Crop Tessilated Area
Because we specified an “offset” when constructing the mesh, some of the Integration Points fall outside of our study domain. That is, the tessellated area covers an extent greater than our study domain.
Here, we crop the tessilation to coincide with our study domain and then calculate the area for each resulting polygon.
Note in the plot that several Integration Points are outside of the Domain, these are later set to “0” (zero).
library(raster)
New.poly = crop(Int.poly, Domain)
New.poly$Area = gArea(New.poly, byid = TRUE)
plot(Integration, col = "blue", pch=19, main = "Tessilation cropped to fit Domain")
plot(New.poly, add=TRUE)
Compare Intergration Tessilation to Mesh
The “Duel Mesh.” Overlaying the mesh for comparison.
plot(Mesh0, main = " ")
plot(New.poly, add=TRUE, border = "red")
plot(Integration, pch=19, col="blue", add=TRUE)
Extract Area to Integration Points
Exposure is set as the Area for Integration Points within the Domain and “0” (zero) to those outside the Domain.
Area = extract(New.poly, Integration)[,"Area"]
Integration$Exposure = ifelse(is.na(Area), 0, Area)Combine Integration Points with Observation Points
library(maptools)## Checking rgeos availability: TRUEMyPoints = spRbind(Observations, Integration)SpatialPointsDataFrame with Exposure
Combined Observations (“presence,” OBS == 1, Exposure == 0) are shown as red circles and Integration Points (OBS == 0) are shown as blue squares.
For Integration Points inside of the domain Exposure == Area. If the Integration Point is outside the domain (e.g. it’s in the mesh “offset” buffer) Exposure == 0.
rbPal = c('blue','red')
Col = rbPal[as.numeric(cut(MyPoints$OBS, breaks = 2))]
plot(Mesh0, main = " ")
plot(MyPoints, col = Col, add = TRUE, pch =MyPoints$OBS)
MyPoints@data## X Y OBS Exposure
## 1 0.62992989 6.128393e-01 1 0.000000000
## 2 0.36069708 4.672791e-02 1 0.000000000
## 3 0.12872211 5.863965e-01 1 0.000000000
## 4 0.37872857 5.965120e-01 1 0.000000000
## 5 0.65691308 7.712343e-01 1 0.000000000
## 6 0.44368413 9.524700e-01 1 0.000000000
## 7 0.79031549 6.265696e-01 1 0.000000000
## 8 0.67902729 4.780218e-01 1 0.000000000
## 9 0.51018136 2.815514e-01 1 0.000000000
## 10 0.63291023 7.249613e-01 1 0.000000000
## 11 0.04810094 6.623291e-01 1 0.000000000
## 12 0.62920417 9.374384e-01 1 0.000000000
## 13 0.29438176 9.113665e-01 1 0.000000000
## 14 0.05085453 4.235256e-01 1 0.000000000
## 15 0.75326092 9.065964e-01 1 0.000000000
## 16 0.02684669 1.943160e-01 1 0.000000000
## 17 0.69795217 7.060831e-01 1 0.000000000
## 18 0.09225859 9.109489e-02 1 0.000000000
## 19 0.76103936 3.036834e-01 1 0.000000000
## 20 0.02588893 7.903673e-02 1 0.000000000
## 21 0.29241514 5.294816e-01 1 0.000000000
## 22 0.76926444 3.924260e-01 1 0.000000000
## 23 0.86998219 9.549788e-01 1 0.000000000
## 24 0.78364701 5.877028e-01 1 0.000000000
## 25 0.65753734 4.818738e-01 1 0.000000000
## 26 0.33251400 8.177225e-01 1 0.000000000
## 27 0.09859819 4.070855e-02 1 0.000000000
## 28 0.87155409 2.396707e-01 1 0.000000000
## 29 0.70993704 5.915011e-01 1 0.000000000
## 30 0.21091496 3.423274e-01 1 0.000000000
## 31 0.31473859 7.487643e-01 1 0.000000000
## 32 0.25064200 5.120360e-01 1 0.000000000
## 33 0.89781492 9.770123e-02 1 0.000000000
## 34 0.97934395 5.281196e-01 1 0.000000000
## 35 0.89233196 3.317584e-01 1 0.000000000
## 36 0.76988273 3.881219e-01 1 0.000000000
## 37 0.65231647 1.359174e-01 1 0.000000000
## 38 0.94976662 7.588337e-01 1 0.000000000
## 39 0.46188762 8.293674e-01 1 0.000000000
## 40 0.94520252 4.889183e-01 1 0.000000000
## 41 0.18090591 4.516507e-01 1 0.000000000
## 42 0.99708695 7.686296e-01 1 0.000000000
## 43 0.91459499 2.202009e-01 1 0.000000000
## 44 0.27844558 7.381295e-01 1 0.000000000
## 45 0.18927426 8.961655e-02 1 0.000000000
## 46 0.33046872 1.706446e-02 1 0.000000000
## 47 0.53818140 7.839864e-01 1 0.000000000
## 48 0.62493991 2.511208e-01 1 0.000000000
## 49 0.14961174 7.678555e-01 1 0.000000000
## 50 0.41505438 2.537247e-01 1 0.000000000
## 51 -0.07048121 -4.532721e-01 0 0.000000000
## 52 1.05394837 -4.532721e-01 0 0.000000000
## 53 1.47934395 -2.787652e-02 0 0.000000000
## 54 1.47934395 1.052724e+00 0 0.000000000
## 55 1.07708897 1.454979e+00 0 0.000000000
## 56 0.12962827 1.454979e+00 0 0.000000000
## 57 -0.47315331 8.521972e-01 0 0.000000000
## 58 -0.47315331 -5.059999e-02 0 0.000000000
## 59 0.62992989 6.128393e-01 0 0.033287621
## 60 0.36069708 4.672791e-02 0 0.042547351
## 61 0.12872211 5.863965e-01 0 0.023670899
## 62 0.37872857 5.965120e-01 0 0.035699319
## 63 0.65691308 7.712343e-01 0 0.031279308
## 64 0.44368413 9.524700e-01 0 0.009067665
## 65 0.79031549 6.265696e-01 0 0.036765608
## 66 0.67902729 4.780218e-01 0 0.038054567
## 67 0.51018136 2.815514e-01 0 0.040986388
## 68 0.04810094 6.623291e-01 0 0.000000000
## 69 0.62920417 9.374384e-01 0 0.014760337
## 70 0.29438176 9.113665e-01 0 0.000000000
## 71 0.05085453 4.235256e-01 0 0.015775514
## 72 0.75326092 9.065964e-01 0 0.021162206
## 73 0.02684669 1.943160e-01 0 0.000000000
## 74 0.09225859 9.109489e-02 0 0.027154441
## 75 0.76103936 3.036834e-01 0 0.031601858
## 76 0.29241514 5.294816e-01 0 0.023420498
## 77 0.86998219 9.549788e-01 0 0.000000000
## 78 0.33251400 8.177225e-01 0 0.024992797
## 79 0.87155409 2.396707e-01 0 0.025358897
## 80 0.21091496 3.423274e-01 0 0.045397959
## 81 0.89781492 9.770123e-02 0 0.000000000
## 82 0.97934395 5.281196e-01 0 0.000000000
## 83 0.65231647 1.359174e-01 0 0.033573644
## 84 0.94976662 7.588337e-01 0 0.027996407
## 85 0.46188762 8.293674e-01 0 0.029330791
## 86 0.18090591 4.516507e-01 0 0.017427900
## 87 0.62493991 2.511208e-01 0 0.021361459
## 88 0.14961174 7.678555e-01 0 0.000000000
## 89 0.49173358 -4.532721e-01 0 0.000000000
## 90 1.47934395 5.124237e-01 0 0.000000000
## 91 1.47934395 2.422736e-01 0 0.000000000
## 92 0.60335862 1.454979e+00 0 0.000000000
## 93 -0.27181726 -2.519360e-01 0 0.000000000
## 94 -0.17176252 1.153588e+00 0 0.000000000
## 95 -0.16558405 3.475108e-05 0 0.000000000
## 96 1.17552819 2.224848e-01 0 0.000000000
## 97 0.29309566 1.207749e+00 0 0.000000000
## 98 -0.08899920 8.957993e-01 0 0.000000000
## 99 -0.47315331 4.007986e-01 0 0.000000000
## 100 0.73773052 1.109048e+00 0 0.000000000
## 101 0.43492967 4.350562e-01 0 0.038050529
## 102 1.18057170 9.670157e-01 0 0.000000000