Load Function Libraries

library(sp)
## Warning: package 'sp' was built under R version 3.3.2
library(maptools)
## Warning: package 'maptools' was built under R version 3.3.2
## Checking rgeos availability: TRUE
library(rgeos)
## Warning: package 'rgeos' was built under R version 3.3.2
## rgeos version: 0.3-23, (SVN revision 546)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Linking to sp version: 1.2-4 
##  Polygon checking: TRUE
library(rgdal)
## Warning: package 'rgdal' was built under R version 3.3.2
## rgdal: version: 1.2-7, (SVN revision 660)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-4
library(raster)
library(gstat)
## Warning: package 'gstat' was built under R version 3.3.2
library(ggmap)
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.3.2
## Google Maps API Terms of Service: http://developers.google.com/maps/terms.
## Please cite ggmap if you use it: see citation("ggmap") for details.
library(ggplot2)

Take a look at packages

# vignette(package=c("sp", "maptools", "raster"))
# sp1 <- vignette("intro_sp") # The intro sp vingette
# print(sp1) # The pdf
# edit(sp1) # Just the R code

# Coordinate Reference system
# reference ellipsoid - defines shape of the earth
# datum - reference point(s) on the earth's surface
# projection - how we project a spherical surface of the earth onto 2 dimensions (equirectangular, cylindrical or Mercator's, conical); The meuse coordinates refernce to Netherlands RDH topographical convention
# sp:proj4string() uses rgdal library to specify CRS
# sp:spTransform() will change from one CRS to another

Meuse River Dataset

Investigate the structure of the meuse data set, convert into spatial, and play with aspects of mapping these data.

data(meuse)
class(meuse) # Right now it is a data.frame
## [1] "data.frame"
str(meuse) # $x and $y hold location info
## 'data.frame':    155 obs. of  14 variables:
##  $ x      : num  181072 181025 181165 181298 181307 ...
##  $ y      : num  333611 333558 333537 333484 333330 ...
##  $ cadmium: num  11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
##  $ copper : num  85 81 68 81 48 61 31 29 37 24 ...
##  $ lead   : num  299 277 199 116 117 137 132 150 133 80 ...
##  $ zinc   : num  1022 1141 640 257 269 ...
##  $ elev   : num  7.91 6.98 7.8 7.66 7.48 ...
##  $ dist   : num  0.00136 0.01222 0.10303 0.19009 0.27709 ...
##  $ om     : num  13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
##  $ ffreq  : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ soil   : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
##  $ lime   : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
##  $ landuse: Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
##  $ dist.m : num  50 30 150 270 380 470 240 120 240 420 ...
plot(meuse$x, meuse$y) # Plot data using coordinates, but this is not spatial data yet

coordinates(meuse) <- ~x+y # Coerce meuse into a SpatialPointsDataFrame, "sp" and tell it which columns contain the coordinates.

class(meuse) # Double check this coerced as expected and now we can work with it as spatial data
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
str(meuse)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  155 obs. of  12 variables:
##   .. ..$ cadmium: num [1:155] 11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
##   .. ..$ copper : num [1:155] 85 81 68 81 48 61 31 29 37 24 ...
##   .. ..$ lead   : num [1:155] 299 277 199 116 117 137 132 150 133 80 ...
##   .. ..$ zinc   : num [1:155] 1022 1141 640 257 269 ...
##   .. ..$ elev   : num [1:155] 7.91 6.98 7.8 7.66 7.48 ...
##   .. ..$ dist   : num [1:155] 0.00136 0.01222 0.10303 0.19009 0.27709 ...
##   .. ..$ om     : num [1:155] 13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
##   .. ..$ ffreq  : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..$ soil   : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
##   .. ..$ lime   : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
##   .. ..$ landuse: Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
##   .. ..$ dist.m : num [1:155] 50 30 150 270 380 470 240 120 240 420 ...
##   ..@ coords.nrs : int [1:2] 1 2
##   ..@ coords     : num [1:155, 1:2] 181072 181025 181165 181298 181307 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:155] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "x" "y"
##   ..@ bbox       : num [1:2, 1:2] 178605 329714 181390 333611
##   .. ..- 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 NA
# The @data contains all of the attritutes

# Now, tell R that these are in a geographical projection
proj4string(meuse) <- CRS("+init=epsg:28992")
meuse@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
# SpatialLines object joins up these points in sequence
m.cc <- coordinates(meuse)
m.sl <- SpatialLines(list(Lines(list(Line(m.cc)), "line1")))
plot(m.sl, main="Lines")

# Load river data (lines format)
data(meuse.riv)
# Convert to a polygon:
rivm <- list(Polygons(list(Polygon(meuse.riv)), "meuse.riv"))
meuse.riv <- SpatialPolygons(rivm)
class(meuse.riv)
## [1] "SpatialPolygons"
## attr(,"package")
## [1] "sp"
plot(meuse.riv, col="royalblue")

# Give same graphical projection to the river polygon, so they can be mapped together
proj4string(meuse.riv) <- CRS("+init=epsg:28992")

### MAPPING 1: Create maps utilizing the bubble function, indicating which attribute using "zcol=", and add the river polygon created above.
bubble(meuse, zcol="zinc", scales=list(draw=TRUE), col="red", pch=1, maxsize=1.5, sp.layout=list("sp.polygons", meuse.riv, col="royalblue"))

bubble(meuse, zcol="lead", scales=list(draw=TRUE), col="black", pch=1, maxsize=1.5, sp.layout=list("sp.polygons", meuse.riv, col="royalblue"))

# Ensure proper libraries are loaded for this chunck of code
library(rgdal)
library(raster)
library(mapmisc)
## 
## Attaching package: 'mapmisc'
## The following object is masked from 'package:ggmap':
## 
##     geocode
proj4string(meuse) <- CRS("+init=epsg:28992")
meuse.ll <- spTransform(meuse, CRS("+init=epsg:4326"))
data(netherlands)
class(nldElev)
## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"
### MAPPING 2: Create a scale for zinc and elevation levels,w which are then mapped with the map.new() function

# Crop the raster layer nldElev to fit the meuse dataset better
nldElev.m <- crop(nldElev, extend(extent(meuse), 1000))
plot(nldElev.m)
plot(meuse, add=TRUE, pch=20)

# Create scale for zinc levels and for the elevation
zincScale <- colourScale(meuse$zinc, breaks=5, style='equal', opacity=0.8, dec=-1, firstBreak=0)
elevScale <- colourScale(nldElev.m, style='equal', breaks=6, col=terrain.colors, firstBreak=0, dec=-1, opacity=c(0.2, 0.9))

map.new(meuse)
image(nldElev, breaks=elevScale$breaks, col=elevScale$colOpacity, legend=F, add=TRUE)
plot(meuse.riv, col="lightblue", add=TRUE)
plot(meuse, col=zincScale$plot, add=TRUE, pch=16)
legendBreaks("topleft", breaks=zincScale, title="Zinc (ppm)")

### MAPPING 3

# Pull interesting pieces of the raster file nldElev, so we can utilize these in a map
mContour <- rasterToContour(nldElev.m)
mSlope <- terrain(nldElev.m, opt="slope")
mAspect <- terrain(nldElev.m, opt="aspect")
mHill <- hillShade(mSlope, mAspect, 40, 270)

# Attribute
meuse$Elev <- extract(nldElev.m, meuse)
meuse$Slope <- extract(mSlope, meuse)
meuse$Aspect <- extract(mAspect, meuse)
head(meuse@data)
##   cadmium copper lead zinc  elev       dist   om ffreq soil lime landuse
## 1    11.7     85  299 1022 7.909 0.00135803 13.6     1    1    1      Ah
## 2     8.6     81  277 1141 6.983 0.01222430 14.0     1    1    1      Ah
## 3     6.5     68  199  640 7.800 0.10302900 13.0     1    1    1      Ah
## 4     2.6     81  116  257 7.655 0.19009400  8.0     1    2    0      Ga
## 5     2.8     48  117  269 7.480 0.27709000  8.7     1    2    0      Ah
## 6     3.0     61  137  281 7.791 0.36406700  7.8     1    2    0      Ga
##   dist.m     Elev       Slope   Aspect
## 1     50 32.81339 0.015897487 4.239358
## 2     30 32.05648 0.026658793 4.613836
## 3    150 33.18457 0.007027912 1.158795
## 4    270 34.05482 0.023137072 4.608548
## 5    380 31.76153 0.009861425 3.684898
## 6    470 33.02278 0.036534802 4.682915
# Map all of these elements together
plot(mHill, col=grey(0:100/100), legend=F, main="Meuse River",
     axes=F, bty="n", box=F)
plot(nldElev.m, col=topo.colors(25, alpha=0.35), add=TRUE,
     legend.args=list(text="Elevation (m)"))
plot(mContour, add=T, col="grey10")
plot(meuse.riv, col="lightblue", add=TRUE)
plot(meuse, col=zincScale$plot, add=TRUE, pch=16)
legendBreaks("bottomleft", breaks=zincScale, title="Zinc (ppm)")

# Use a different projection for the meuse data
library(rgdal)
proj4string(meuse) <- CRS("+init=epsg:28992")
meuse.ll <- spTransform(meuse, CRS("+init=epsg:4326"))
meuse.latlong <- spTransform(meuse.ll, CRS("+proj=longlat")) # Transform into lat/long data

map.new(meuse.latlong)

plot(meuse.latlong, axes=T)

subset.lead <- meuse.latlong@data[,3]
plot(subset.lead)

# Use Google maps
m.map <- get_map(location= c(lon=5.742061, lat= 50.971337), source="google", maptype="satellite", crop=FALSE, zoom=13)
## Source : https://maps.googleapis.com/maps/api/staticmap?center=50.971337,5.742061&zoom=13&size=640x640&scale=2&maptype=satellite&language=en-EN
ggmap(m.map)

# I would like to plot onto this map, but I am not sure how to accomplish this with this type of file

# ggmap(m.map) + 
# Create bubbleplot of data
bubble(meuse, xcol="x", ycol="y", zcol="lead", main= "Lead concentrations (ppm)")

Variogram showing lead concentrations (ppm) along the Meuse rive

The variogram shows the distance between locations along the x-axis and the difference of attribute values squared along the y-axis.

# Look at all point values, the "cloud", by indicating cloud=TRUE in the function
leadvarcloud <- variogram(lead~1, meuse, cloud=T)
plot(leadvarcloud, pch=20, cex=0.5, col="black",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", 
     main= "Lead concentrations (ppm)")

## Look at the average model value at fewer distances to get a clearer visual of what's going on
leadvar <- variogram(lead~1, meuse, cloud=FALSE)
plot(leadvar, pch=20, cex=0.5, col="black",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", 
     main= "Lead concentrations (ppm)")

leadvar.log <- variogram(log(lead)~1, meuse)
plot(leadvar.log, pch=20, cex=0.5, col="black",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", 
     main= "Log Transformed Lead concentrations (ppm)") # log transform didn't change the overall distribution

The variogram indicates that lead concentrations are autocorrelated out to a distance of about 1000 meters, because this is when the semivariance begins to level off. Semivariance increases at greater distances, which means lead concentration values are less similar at greater distances. The log transformed variogram plot looks a little cleaner, but shows the same pattern with semivariance increasing until it levels off and actually decreases again by a little less than 1500 meters.

Mystery Data

load("mysteryData.Rdata")

summary(dat)
##        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(dat) <- c("x", "y")
bubble(dat, xcol="x", ycol="y", zcol="z", main= "Mystery Data")

proj4string(dat) <- CRS("+init=epsg:28992")
class(dat)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
str(dat)
## 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
##   .. .. ..$ : NULL
##   .. .. ..$ : 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__
plot(dat, axes=T)

# Plot variogram
zvar.cloud <- variogram(z~1, dat, cloud=T)
plot(zvar.cloud, pch=20, col="royalblue", ylab=expression("Semivariance (" *gamma* ")"), xlab="Distance", main = "Mystery Z Attribute")

zvar <- variogram(z~1, dat, cloud=F)
plot(zvar, pch=20, col="royalblue", ylab=expression("Semivariance (" *gamma* ")"), xlab="Distance", main = "Mystery Z Attribute")

The sill appears to be around 15000 for these data.

Playing around with other mapping ideas

# See for examples
vignette("mapmisc", package="mapmisc")

Practice getting map to utilize

fidalgo.bay <- c(long=-122.579859, lat=48.473768)
fidalgo.bay <- get_map(location= fidalgo.bay, source="google", maptype="terrain", crop=FALSE, zoom=13)
## Source : https://maps.googleapis.com/maps/api/staticmap?center=48.473768,-122.579859&zoom=13&size=640x640&scale=2&maptype=terrain&language=en-EN
ggmap(fidalgo.bay)

# Plot the data locations
map.new(meuse)
plot(nldTiles, add=T)
points(meuse, col="blue", cex=0.5)
scaleBar(proj4string(meuse), pos="topleft", bg="white")