Mapping

require(ggmap)
## Loading required package: ggmap
## Loading required package: ggplot2
require(sp)
## Loading required package: sp
require(maptools)
## Loading required package: 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()
require(rgdal)
## Loading required package: rgdal
## rgdal: version: 1.0-4, (SVN revision 548)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
##  Path to GDAL shared files: /Users/shannoncall/Library/R/3.2/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Users/shannoncall/Library/R/3.2/library/rgdal/proj
##  Linking to sp version: 1.1-1
require(raster)
## Loading required package: raster
citation('ggmap')
## 
## To cite ggmap in publications, please use:
## 
##   D. Kahle and H. Wickham. ggmap: Spatial Visualization with
##   ggplot2. The R Journal, 5(1), 144-161. URL
##   http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf
## 
## A BibTeX entry for LaTeX users is
## 
##   @Article{,
##     author = {David Kahle and Hadley Wickham},
##     title = {ggmap: Spatial Visualization with ggplot2},
##     journal = {The R Journal},
##     year = {2013},
##     volume = {5},
##     number = {1},
##     pages = {144--161},
##     url = {http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf},
##   }
I decided to create maps for my thesis proposal of each site I will be doing a vegetation survey later this summer. First, I created a map of the entire region, West of the Puget Sound on Kitsap Peninsula. Then I created a map of each creek/estuary.Site 1. Beaver Creek; Site 2. Dogfish Creek; and 3. Carpenter Creek; located in Kitsap County, Washington State.

Kitsap County

Kitsap.County<-c(lon=-122.652583,lat=47.749333)
Kitsap.County<-get_map(location=Kitsap.County, source="google",maptype="terrain",crop=FALSE,zoom=10)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=47.749333,-122.652583&zoom=10&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
ggmap(Kitsap.County)

Site 1: Beaver Creek

B.Creek.LatLon <-c(lon=-122.556443,lat=47.570382)
B.Creek<-get_map(location=B.Creek.LatLon, source="google",maptype="terrain", crop=FALSE,zoom=14)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=47.570382,-122.556443&zoom=14&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
ggmap(B.Creek)+
geom_point(aes(x=-122.556443,y=47.570382), data=NULL, alpha=.5,pch=21, col="darkred", bg="orange", size=2)

Site 2:Dogfish Creek

D.Creek.LatLon<-c(lon=-122.652583,lat=47.749333)
D.Creek<-get_map(location=D.Creek.LatLon, source="google",maptype="terrain", crop=FALSE,zoom=14)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=47.749333,-122.652583&zoom=14&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
ggmap(D.Creek)+
geom_point(aes(x=-122.652583,y=47.749333), data=NULL, alpha=.5,pch=21, col="darkred", bg="orange", size=2)

Site 3:Carpenter Creek

C.Creek.LatLon<-c(lon=-122.514859,lat=47.801068)
C.Creek<-get_map(location=C.Creek.LatLon, source="google",maptype="terrain", crop=FALSE,zoom=14)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=47.801068,-122.514859&zoom=14&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
ggmap(C.Creek)+
geom_point(aes(x=-122.514859,y=47.801068), data=NULL, alpha=.5,pch=21, col="darkred", bg="orange", size=2)

Variogram

require(gstat)
## Loading required package: gstat
data(meuse.all)
meuse<-meuse.all

class(meuse)
## [1] "data.frame"
coordinates(meuse)<- ~x+y
proj4string(meuse) <- CRS("+init=epsg:28992")
class(meuse)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
bubble(meuse, xcol="x",ycol="y",zcol="lead", maxsize= 2.5, main = "Lead Concentrations (ppm)")

Lead Variogram, a measure of environmental distance

PbVariCloud<- variogram(lead~1, meuse, cloud = TRUE)
plot(PbVariCloud, pch=21, col="darkred", 
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", main="Lead (Pb) Concentrations (ppm)")

So, as distance between points increases, semivariance increases, meaning they are less similar. The x-axis represents the distance between points, and they y-axis represents the difference of Pb squared values. Each point represents two locations. I don’t find this particularly helpful.
PbVari <- variogram(lead~1, meuse, cloud = TRUE, cutoff=15000)
plot(PbVari,pch=20, cex=1.5,col="deeppink3",
     ylab=expression("Semivariance ("*gamma*")"),xlab="Distance (m)", main = "Lead (Pb) Concentrations (ppm)")

plot(variogram(log(lead)~1, meuse, cutoff=1000)) #Looks like a log transform didn't change the overall distribution between. Data are normally distributed.

The longer the cutoff, the less data show up because it is looking at distances further away from the site. At some point, there are physical constraints that don’t allow the transfer of metals, like 4,000 meters away from the deposition site. The particulates drop out of the water long before they reach any distance like 4,000 meters. Cutoff: This is a spatial separation distance. The length of the diagonal between two points is divided by three. Width: this fuction manipulates the intervals of data pair groups.

Point Pattern - I can’t get the code to knit with this in.. I’ve done it so many times…

Mystery Data

setwd("~/Desktop/ESCI597A/Week05/Part02")
load("mysteryData.Rdata")

mystery<-dat
summary(mystery)
##        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.5846
coordinates(mystery)<- ~x+y

bubble(mystery, xcol="x",ycol="y",zcol="z",maxsize=2.5,main="Mystery Data")

Looks like there is a large concentration toward the center, with lower values near the outside.
proj4string(mystery) <- CRS("+init=epsg:28992")
class(mystery)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
str(mystery)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  100 obs. of  1 variable:
##   .. ..$ z: num [1:100] 292 486 78 168 73 ...
##   ..@ coords.nrs : int [1:2] 1 2
##   ..@ coords     : num [1:100, 1:2] 10.12 -3.7 -3.7 5.76 17.4 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:100] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "x" "y"
##   ..@ bbox       : num [1:2, 1:2] -16.1 -17.6 18.9 12.9
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "x" "y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+init=epsg:28992 +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bess"| __truncated__
mystery@proj4string
## CRS arguments:
##  +init=epsg:28992 +proj=sterea +lat_0=52.15616055555555
## +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000
## +ellps=bessel
## +towgs84=565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725
## +units=m +no_defs
plot(mystery,axes=TRUE)

zVarCloud <-variogram(z~1, mystery, cloud=TRUE)
plot(zVarCloud, pch=21, col="darkred", ylab=expression("Semivariance ("*gamma*")"), xlab="Distance (m)",main="Mystery Z")

zVar <- variogram(z~1, mystery, cloud=FALSE)
plot(zVar, pch=20, cex=1.5,col="darkred", ylab=expression(("Semivariance ("*gamma*")"), xlab="Distance (m)",main="Mystery Z"))

It seems these data are autocorrelated throughout because of the increasing semivariance. However, it seems the range of this variogram is between ~2,500 and ~125,000. The sill (the variance) is probably around just less than 150,000.

Part 02

require(ggplot2)
require(spatstat) 
## Loading required package: 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, shift

Tupelo

setwd("/Users/shannoncall/Desktop/ESCI597A/Week05/Part02")

tupelo <- read.table("tupelo.txt",header=T)
str(tupelo)
## 'data.frame':    751 obs. of  4 variables:
##  $ plot: int  1 1 1 1 1 1 1 1 1 1 ...
##  $ sex : Factor w/ 3 levels "f","j","m": 1 1 1 1 1 1 1 1 1 1 ...
##  $ x   : num  12.76 16.27 3.82 5.74 4.84 ...
##  $ y   : num  11.42 6.09 19.44 2.05 10.53 ...
ggplot(tupelo) +geom_point(aes(x,y,colour=sex)) + facet_grid(.~plot)

coordinates(tupelo) <- c("x","y")
tupeloPPP <- as(tupelo,"ppp")
plot2<- subset(tupeloPPP, plot==2, select=sex, drop=TRUE)
summary(plot2)
## Marked planar point pattern:  257 points
## Average intensity 0.09453631 points per square unit
## 
## Coordinates are given to 4 decimal places
## 
## Multitype:
##   frequency proportion  intensity
## f       104  0.4046693 0.03825594
## j        32  0.1245136 0.01177106
## m       121  0.4708171 0.04450931
## 
## Window: rectangle = [-0.8779, 51.2377] x [0.2889, 52.4524] units
## Window area = 2718.53 square units
plot(plot2)

Male Tupelo

plot2Males <- subset(plot2, marks=="m", drop=TRUE)
summary(plot2Males)
## Marked planar point pattern:  121 points
## Average intensity 0.04450931 points per square unit
## 
## Coordinates are given to 4 decimal places
## 
## Multitype:
##   frequency proportion  intensity
## m       121          1 0.04450931
## 
## Window: rectangle = [-0.8779, 51.2377] x [0.2889, 52.4524] units
## Window area = 2718.53 square units
plot(plot2Males)

The distribution of the males seems fairly random. Let’s see if they are.
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.
plot(plot2Males.k)

Distribution distance is completely CSR. The black line is completely contained within the grey interval.

Female Tupelo

plot2Females <- subset(plot2, marks=="f", drop=TRUE)
summary(plot2Females)
## Marked planar point pattern:  104 points
## Average intensity 0.03825594 points per square unit
## 
## Coordinates are given to 4 decimal places
## 
## Multitype:
##   frequency proportion  intensity
## f       104          1 0.03825594
## 
## Window: rectangle = [-0.8779, 51.2377] x [0.2889, 52.4524] units
## Window area = 2718.53 square units
plot(plot2Females)

Looks pretty random, let’s see what these individuals’ distribution actually is…
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.
plot(plot2Females.k)

Distribution distance is completely CSR (random), nearly clumped between distance 4 and 8.

Juveniles

plot2juvies <- subset(plot2, marks=="j", drop=TRUE)
summary(plot2juvies)
## Marked planar point pattern:  32 points
## Average intensity 0.01177106 points per square unit
## 
## Coordinates are given to 4 decimal places
## 
## Multitype:
##   frequency proportion  intensity
## j        32          1 0.01177106
## 
## Window: rectangle = [-0.8779, 51.2377] x [0.2889, 52.4524] units
## Window area = 2718.53 square units
plot(plot2juvies)

These individuals look clumped. Makes sense since seed dispersal is often a result of a particular biological pattern, like scat. Let’s see if they really are clumped.
plot2juvies.k<- envelope(plot2juvies, 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.
plot(plot2juvies.k)

Juvenile’s are mostly clumped after distance ~3 and before distance ~11.5. We might be able to explain this from a biological factor.

Finnish Hyytiala Data

library(spatstat)
head(hyytiala)
## Marked planar point pattern: 6 points
## Multitype, with levels = aspen, birch, pine, rowan 
## window: rectangle = [0, 20] x [0, 20] metres
data(hyytiala)
str(hyytiala)
## List of 6
##  $ window    :List of 4
##   ..$ type  : chr "rectangle"
##   ..$ xrange: num [1:2] 0 20
##   ..$ yrange: num [1:2] 0 20
##   ..$ units :List of 3
##   .. ..$ singular  : chr "metre"
##   .. ..$ plural    : chr "metres"
##   .. ..$ multiplier: num 1
##   .. ..- attr(*, "class")= chr "units"
##   ..- attr(*, "class")= chr "owin"
##  $ n         : int 168
##  $ x         : num [1:168] 1.912 1.146 0.476 3.539 2.726 ...
##  $ y         : num [1:168] 19.5 18.3 16.8 19.1 17.8 ...
##  $ markformat: chr "vector"
##  $ marks     : Factor w/ 4 levels "aspen","birch",..: 3 3 3 3 3 3 3 3 3 3 ...
##  - attr(*, "class")= chr "ppp"
plot(hyytiala)

plot(split(hyytiala))

plot(density(hyytiala, 10))

I’m going to plot each species separately and look at their distribution

aspen<- subset(hyytiala, marks=="aspen",drop=TRUE)
summary(aspen)
## Marked planar point pattern:  1 points
## Average intensity 0.0025 points per square metre
## 
## Coordinates are given to 5 decimal places
## 
## Multitype:
##       frequency proportion intensity
## aspen         1          1    0.0025
## 
## Window: rectangle = [0, 20] x [0, 20] metres
## Window area = 400 square metres
## Unit of length: 1 metre
plot(aspen) #1 aspen

#Not looking at distribution- there is only one individual


birch<- subset(hyytiala,marks=="birch",drop=TRUE)
summary(birch)
## Marked planar point pattern:  17 points
## Average intensity 0.0425 points per square metre
## 
## Coordinates are given to 7 decimal places
## 
## Multitype:
##       frequency proportion intensity
## birch        17          1    0.0425
## 
## Window: rectangle = [0, 20] x [0, 20] metres
## Window area = 400 square metres
## Unit of length: 1 metre
plot(birch) #17 birch trees

plot(density.ppp(birch))

plot(envelope(birch, Kest, nsim = 1000)) #these look somewhat uniform in distribution, with one clumped area
## Generating 1000 simulations of CSR  ...
## 1, 2, 3, ......10.........20.........30.........40.........50
## .........60.........70.........80.........90.........100
## .........110.........120.........130.........140.........150
## .........160.........170.........180.........190.........200
## .........210.........220.........230.........240.........250
## .........260.........270.........280.........290.........300
## .........310.........320.........330.........340.........350
## .........360.........370.........380.........390.........400
## .........410.........420.........430.........440.........450
## .........460.........470.........480.........490.........500
## .........510.........520.........530.........540.........550
## .........560.........570.........580.........590.........600
## .........610.........620.........630.........640.........650
## .........660.........670.........680.........690.........700
## .........710.........720.........730.........740.........750
## .........760.........770.........780.........790.........800
## .........810.........820.........830.........840.........850
## .........860.........870.........880.........890.........900
## .........910.........920.........930.........940.........950
## .........960.........970.........980.........990......... 1000.
## 
## Done.

pine<- subset(hyytiala, marks=="pine",drop=TRUE)
summary(pine)
## Marked planar point pattern:  128 points
## Average intensity 0.32 points per square metre
## 
## Coordinates are given to 7 decimal places
## 
## Multitype:
##      frequency proportion intensity
## pine       128          1      0.32
## 
## Window: rectangle = [0, 20] x [0, 20] metres
## Window area = 400 square metres
## Unit of length: 1 metre
plot(pine) # 128 pine trees

plot(density.ppp(pine))

plot(envelope(pine, 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.

rowan<- subset(hyytiala, marks=="rowan",drop=TRUE)
summary(rowan)
## Marked planar point pattern:  22 points
## Average intensity 0.055 points per square metre
## 
## Coordinates are given to 7 decimal places
## 
## Multitype:
##       frequency proportion intensity
## rowan        22          1     0.055
## 
## Window: rectangle = [0, 20] x [0, 20] metres
## Window area = 400 square metres
## Unit of length: 1 metre
plot(rowan) #22 rowan trees; average intensity is 0.05 trees per square meter

plot(density.ppp(rowan))

plot(envelope(rowan, 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.

plot(Kest(hyytiala))

#Kpois = theoretical Poisson K(r)
#Kbord = border-corrected estimate of K(r)
#Ktrans = translation-corrected estimate of K(r)
#Kiso = isotropic correction
plot(envelope(hyytiala,Kest, nsim=200)) 
## Generating 200 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24
## .26.28.30.32.34.36.38.40.42.44.46.48
## .50.52.54.56.58.60.62.64.66.68.70.72
## .74.76.78.80.82.84.86.88.90.92.94.96
## .98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144
## .146.148.150.152.154.156.158.160.162.164.166.168
## .170.172.174.176.178.180.182.184.186.188.190.192
## .194.196.198. 200.
## 
## Done.

Mixed forest in Hyytiala, Finland. Tree distribution seems to be mostly randomly distributed, with a small amount of clumping at distanced greater than 4 meters.