Week5_spat
Sara Cendejas-Zarelli
April 28, 2016
Part 1–Make an interesting map
library(sp)
library(maptools)## Checking rgeos availability: FALSE
## Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
## which has a restricted licence. It is disabled by default;
## to enable gpclib, type gpclibPermit()library(rgdal)## rgdal: version: 1.1-8, (SVN revision 616)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
## Path to GDAL shared files: C:/Users/Sara/Documents/R/win-library/3.3/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
## Path to PROJ.4 shared files: C:/Users/Sara/Documents/R/win-library/3.3/rgdal/proj
## Linking to sp version: 1.2-3library(raster)
library(gstat)
library(ggmap)## Loading required package: ggplot2Messing about with the tree spatial data…
mwaTrees <-read.csv("mwa_trees_xy.csv")
head(mwaTrees)## MSMNT_ID LONG LAT
## 1 MWA014 -114.3149 38.91267
## 2 MWA015 -114.3148 38.91268
## 3 MWA016 -114.3149 38.91257
## 4 MWA017 -114.3147 38.91266
## 5 MWA019 -114.3147 38.91253
## 6 MWA020 -114.3143 38.91242coordinates(mwaTrees)<- c("LONG", "LAT")proj4string(mwaTrees) <- CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0')
mwaTrees@proj4string## CRS arguments:
## +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0plot(mwaTrees,axes=TRUE,pch=21,col="green")
Make a pretty map (Mt. Washington)…
mwaTrees <-read.csv("mwa_trees_xy.csv")mylocation<-c(lon=-114.300,lat=38.915)
myMap<-get_map(location=mylocation, source="stamen",maptype="terrain",crop=TRUE,zoom=14)## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=38.915,-114.3&zoom=14&size=640x640&scale=2&maptype=terrain&sensor=false## Map from URL : http://tile.stamen.com/terrain/14/2988/6265.jpg## Map from URL : http://tile.stamen.com/terrain/14/2989/6265.jpg## Map from URL : http://tile.stamen.com/terrain/14/2990/6265.jpg## Map from URL : http://tile.stamen.com/terrain/14/2991/6265.jpg## Map from URL : http://tile.stamen.com/terrain/14/2988/6266.jpg## Map from URL : http://tile.stamen.com/terrain/14/2989/6266.jpg## Map from URL : http://tile.stamen.com/terrain/14/2990/6266.jpg## Map from URL : http://tile.stamen.com/terrain/14/2991/6266.jpg## Map from URL : http://tile.stamen.com/terrain/14/2988/6267.jpg## Map from URL : http://tile.stamen.com/terrain/14/2989/6267.jpg## Map from URL : http://tile.stamen.com/terrain/14/2990/6267.jpg## Map from URL : http://tile.stamen.com/terrain/14/2991/6267.jpgggmap(myMap)+
geom_point(aes(x=LONG,y=LAT),
data=mwaTrees,alpha = 0.5, color="darkred",size = 2)
Part 2–Interpret a variogram
data(meuse)
head(meuse)## x y cadmium copper lead zinc elev dist om ffreq soil
## 1 181072 333611 11.7 85 299 1022 7.909 0.00135803 13.6 1 1
## 2 181025 333558 8.6 81 277 1141 6.983 0.01222430 14.0 1 1
## 3 181165 333537 6.5 68 199 640 7.800 0.10302900 13.0 1 1
## 4 181298 333484 2.6 81 116 257 7.655 0.19009400 8.0 1 2
## 5 181307 333330 2.8 48 117 269 7.480 0.27709000 8.7 1 2
## 6 181390 333260 3.0 61 137 281 7.791 0.36406700 7.8 1 2
## lime landuse dist.m
## 1 1 Ah 50
## 2 1 Ah 30
## 3 1 Ah 150
## 4 0 Ga 270
## 5 0 Ah 380
## 6 0 Ga 470Turning the meuse data into spatial data…
coordinates(meuse) <- ~x+y
proj4string(meuse) <- CRS("+init=epsg:28992")bubble(meuse, xcol="x",ycol="y",zcol="lead", maxsize = 2.5, main = "Lead concentrations (ppm)")
The Bubble plot shows that high concentrations appear primarily along the river channel.
Here are some Variograms looking at Log transformed Lead concentrations…
With Cloud on
leadVarCloud <- variogram(log(lead)~1, meuse, cloud = TRUE)
plot(leadVarCloud,pch=20,cex=.9,col="gray36",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)",main = "Log Lead (ppm),default cutoff and width")
It looks like the data points are a little skewed towards the right, somewhere above the 1,100 meter point. This would mean that at a distance of 1,100 meters lead samples reach a sill, where lead samples become independent of each other. We can turn cloud off to get a better idea of where this cutoff point is…
With Cloud off
leadVar <- variogram(log(lead)~1, meuse, cloud = FALSE)
plot(leadVar,pch=20,cex=.9,col="gray36",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)",main = "Log Lead (ppm), default cutoff and width")
Yes, you can see that the apex of the curve is around 1,100 meters.
With Cloud off, I’ll mess around with the cutoffs and widths
leadVar<- variogram(log(lead)~1, meuse, cutoff="500",width="0.5", cloud = FALSE)
plot(leadVar,pch=20,cex=.9,col="gray36",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)", main = "Log Lead (ppm), cutoff=500, width=0.5")
leadVar<- variogram(log(lead)~1, meuse, cutoff="500",width="1.5", cloud = FALSE)
plot(leadVar,pch=20,cex=.9,col="gray36",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)", main = "Log Lead (ppm), cutoff=500, width=1.5")
Changing the cutoff distance sets the distance(m) on the x-axis. Changing the width adjusts the spread of the data points. A wider width means a greater spread. Let’s compare cutoffs at 1,000 meters and two different widths…
leadVar<- variogram(log(lead)~1, meuse, cutoff="1000",width="1.5", cloud = FALSE)
plot(leadVar,pch=20,cex=.9,col="gray36",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)", main = "Log Lead (ppm), cutoff=1000, width=1.5")
leadVar<- variogram(log(lead)~1, meuse, cutoff="1000",width="0.5", cloud = FALSE)
plot(leadVar,pch=20,cex=.9,col="gray36",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)", main = "Log Lead (ppm), cutoff=1000, width=0.5")
Yup, looks like that’s the case. Knowing that, I think you would be most likely to adjust those numbers if there was something inherent in your study system that you knew would affect your distance.
Mystery Data–Mine will be the number of Caddisfly larvae emerging from the Sol Duc River.
load("D:/timeseries/mysteryData.Rdata")Caddis<-dat
summary(Caddis)## x y z
## Min. :-16.063 Min. :-17.587 Min. : 0.0002
## 1st Qu.: -8.061 1st Qu.:-10.744 1st Qu.: 30.0933
## Median : 1.032 Median : -3.900 Median : 77.4185
## Mean : 1.155 Mean : -3.795 Mean :116.1836
## 3rd Qu.: 10.852 3rd Qu.: 2.321 3rd Qu.:184.0281
## Max. : 18.854 Max. : 12.897 Max. :485.5846Let’s turn Caddisflies into a spatial data set shall we?
coordinates(Caddis)<- ~x+y
class(Caddis)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"bubble(Caddis, xcol="x",ycol="y",zcol="z", maxsize = 2.5, main = "Caddisfly abundance")
Interesting plot of Caddisfly numbers and their spatial relationships. Looks like Caddisfly abundance could be autocorrelated in a cluster-type fashion, with points near the center being very similar. Let’s see what the variogram shows…
Variogram
cadVarCloud <- variogram(z~1, Caddis, cloud = TRUE)
plot(cadVarCloud,pch=20,cex=1.5,col="black",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)",main = "Caddisfly larvae abundance")
cadVar <- variogram(z~1, Caddis, cloud = FALSE)
plot(cadVar,pch=20,cex=1.5,col="black",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)", main = "Caddisfly larvae abundance")
I think I need to transform my data, It’s looking like a pretty linear pattern with no obvious sill for the semivariance.
par(mfrow=c(1,3))
hist(dat$z)
hist(abs(dat$z))
hist(sqrt(dat$z))
A square-root transform does a better job of normalizing the data.
cadVar <- variogram(sqrt(z)~1, Caddis, cloud = FALSE)
plot(cadVar,pch=20,cex=1.5,col="black",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)", main = "Sqrt of Caddisfly larvae abundance")
The transform reveals that there is a sill around 11 meters. This suggests that Caddisfly larvae are autocorrelated out to a distance of ~11 meters.
Part 2 Point Patterns
library(spatstat)## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:raster':
##
## getData##
## spatstat 1.45-0 (nickname: 'One Lakh')
## For an introduction to spatstat, type 'beginner'##
## Attaching package: 'spatstat'## The following object is masked from 'package:gstat':
##
## idw## The following objects are masked from 'package:raster':
##
## area, rotate, shiftlibrary(nlme)
library(ggplot2)I’ll try for a CSR pattern…
foo <- clickppp(n = 30) plot(density(foo)) points(foo,pch=20) env <- envelope(foo, Kest, nsim = 100) plot(env)
For the most part the points follow a random pattern here, except for a small portion between ~0.1 and ~0.25, where they follow a uniform pattern.
I’ll try for a clustered pattern next…
foo <- clickppp(n = 30,type=“1”) plot(density(foo)) points(foo,pch=20) env <- envelope(foo, Kest, nsim = 100) plot(c)
Now, a uniform pattern…
uni <- clickppp(n = 30) plot(density(uni)) points(uni,pch=20) env.uni <- envelope(uni, Kest, nsim = 100) plot(env.uni)
It was surprisingly difficult to get a uniform pattern. I kept getting uniform out to a distance of 0.15 (or 1/7th of the total area)
Lastly, clumped and uniform…
uni <- clickppp(n = 30) plot(density(uni)) points(uni,pch=20) env.uni <- envelope(uni, Kest, nsim = 100) plot(env.uni)
Ripley’s K shows that the pattern is clumped to a distance of ~0.17 then it heads over to the uniform pattern at 0.25.
Tupelo Trees
tupelo <- read.csv("tupelo.csv",header=T)
summary(tupelo)## plot sex x y
## Min. :1.000 f:357 Min. :-0.88 Min. : 0.29
## 1st Qu.:1.000 j: 41 1st Qu.:12.35 1st Qu.:12.91
## Median :2.000 m:353 Median :23.83 Median :26.69
## Mean :2.059 Mean :24.65 Mean :26.26
## 3rd Qu.:3.000 3rd Qu.:37.69 3rd Qu.:38.38
## Max. :3.000 Max. :51.24 Max. :52.45class(tupelo)## [1] "data.frame"ggplot(tupelo) + geom_point(aes(x,y,colour=sex)) + facet_grid(.~plot)
coordinates(tupelo)<- c("x","y")
tupeloPPP <- as(tupelo,"ppp")
class(tupelo)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"Plotting Male Tupelo tree patterns:
plot1 <- subset(tupeloPPP, plot == 1, select=sex, drop=TRUE)
summary(plot1)## Marked planar point pattern: 225 points
## Average intensity 0.08276382 points per square unit
##
## *Pattern contains duplicated points*
##
## Coordinates are given to 2 decimal places
## i.e. rounded to the nearest multiple of 0.01 units
##
## Multitype:
## frequency proportion intensity
## f 123 0.54666670 0.045244220
## j 5 0.02222222 0.001839196
## m 97 0.43111110 0.035680400
##
## Window: rectangle = [-0.88, 51.24] x [0.29, 52.45] units
## Window area = 2718.58 square unitsplot1Males <- subset(plot1, marks=="m", drop=TRUE)
plot1Males.k <- envelope(plot1Males, Kest, nsim = 100)## Generating 100 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
## 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
## 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
## 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
## 91, 92, 93, 94, 95, 96, 97, 98, 99, 100.
##
## Done.plot2 <- subset(tupeloPPP, plot == 2, select=sex, drop=TRUE)
summary(plot1)## Marked planar point pattern: 225 points
## Average intensity 0.08276382 points per square unit
##
## *Pattern contains duplicated points*
##
## Coordinates are given to 2 decimal places
## i.e. rounded to the nearest multiple of 0.01 units
##
## Multitype:
## frequency proportion intensity
## f 123 0.54666670 0.045244220
## j 5 0.02222222 0.001839196
## m 97 0.43111110 0.035680400
##
## Window: rectangle = [-0.88, 51.24] x [0.29, 52.45] units
## Window area = 2718.58 square unitsplot2Males <- subset(plot2, marks=="m", drop=TRUE)
plot2Males.k <- envelope(plot2Males, Kest, nsim = 100)## Generating 100 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
## 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
## 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
## 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
## 91, 92, 93, 94, 95, 96, 97, 98, 99, 100.
##
## Done.plot3 <- subset(tupeloPPP, plot == 3, select=sex, drop=TRUE)
summary(plot1)## Marked planar point pattern: 225 points
## Average intensity 0.08276382 points per square unit
##
## *Pattern contains duplicated points*
##
## Coordinates are given to 2 decimal places
## i.e. rounded to the nearest multiple of 0.01 units
##
## Multitype:
## frequency proportion intensity
## f 123 0.54666670 0.045244220
## j 5 0.02222222 0.001839196
## m 97 0.43111110 0.035680400
##
## Window: rectangle = [-0.88, 51.24] x [0.29, 52.45] units
## Window area = 2718.58 square unitsplot3Males <- subset(plot3, marks=="m", drop=TRUE)
plot3Males.k <- envelope(plot3Males, Kest, nsim = 100)## Generating 100 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
## 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
## 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
## 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
## 91, 92, 93, 94, 95, 96, 97, 98, 99, 100.
##
## Done.par(mfrow=c(3,3))
plot(plot1Males)
plot(plot2Males)
plot(plot3Males)
plot(plot1Males.k)
plot(plot2Males.k)
plot(plot3Males.k)
Summarized findings: Male Tupelo trees in plot 1 are in a clustered pattern, whereas the males in plots 2 and 3 follow are random (CSR) pattern.
Let’s look at patterns for the female Tupelo trees
plot1females <- subset(plot1, marks=="f", drop=TRUE)
plot1females.k <- envelope(plot1females, Kest, nsim = 100)## Generating 100 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
## 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
## 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
## 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
## 91, 92, 93, 94, 95, 96, 97, 98, 99, 100.
##
## Done.plot2females <- subset(plot2, marks=="f", drop=TRUE)
plot2females.k <- envelope(plot2females, Kest, nsim = 100)## Generating 100 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
## 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
## 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
## 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
## 91, 92, 93, 94, 95, 96, 97, 98, 99, 100.
##
## Done.plot3females <- subset(plot3, marks=="f", drop=TRUE)
plot3females.k <- envelope(plot3females, Kest, nsim = 100)## Generating 100 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
## 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
## 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
## 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
## 91, 92, 93, 94, 95, 96, 97, 98, 99, 100.
##
## Done.par(mfrow=c(3,3))
plot(plot1females)
plot(plot2females)
plot(plot3females)
plot(plot1females.k)
plot(plot2females.k)
plot(plot3females.k)
Summarized findings: It looks like female Tupelo trees in all 3 plots follow a CSR distribution pattern.
Last, let’s compare the juvenile plots…
plot1juveniles <- subset(plot1, marks=="j", drop=TRUE)
plot1juveniles.k <- envelope(plot1juveniles, Kest, nsim = 100)## Generating 100 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
## 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
## 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
## 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
## 91, 92, 93, 94, 95, 96, 97, 98, 99, 100.
##
## Done.plot2juveniles <- subset(plot2, marks=="j", drop=TRUE)
plot2juveniles.k <- envelope(plot2juveniles, Kest, nsim = 100)## Generating 100 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
## 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
## 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
## 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
## 91, 92, 93, 94, 95, 96, 97, 98, 99, 100.
##
## Done.plot3juveniles <- subset(plot3, marks=="j", drop=TRUE)
plot3juveniles.k <- envelope(plot3juveniles, Kest, nsim = 100)## Generating 100 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
## 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
## 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
## 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
## 91, 92, 93, 94, 95, 96, 97, 98, 99, 100.
##
## Done.par(mfrow=c(3,3))
plot(plot1juveniles)
plot(plot2juveniles)
plot(plot3juveniles)
plot(plot1juveniles.k)
plot(plot2juveniles.k)
plot(plot3juveniles.k)
Summarized findings: Due to the small amount of data for juveniles in plots 1 and 3, the K funtions are unreliable (no confidence intervals exist). Juvenile Tupelo trees in plot 2, however, show clustering up until ~10.5 of the total area of the plot.
The overall pattern of distribution was random. This pattern is largely supported in the literature–as Bawa and Opler (1977) indicate–a random distribution of male and female trees is best for the pollen dispersal.