This file is part of the additional documentation for the spaMM package.

Loading packages

library(spaMM)
## spaMM (version 1.11.1) is loaded.
## Type 'help(spaMM)' for a short introduction,
## and news(package='spaMM') for news.
library(rgdal)
## Loading required package: sp
## rgdal: version: 1.1-7, (SVN revision 612)
##  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/francois/Documents/R/win-library/3.4/rgdal/gdal
##  GDAL does not use iconv for recoding strings.
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: C:/Users/francois/Documents/R/win-library/3.4/rgdal/proj
##  Linking to sp version: 1.2-2
library(rasterVis)
## Loading required package: raster
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer

Toy example of fitting a spatial model

This uses the ‘Bonne’ map projection for computing distances for the fun of showing how various projections can be used (but we would instead advise using the control.dist argument to select a geodesic distance):

data(blackcap)
projection <- "+proj=bonne +lat_1=30n +lon_0=10e +units=km"
bonne <- as.data.frame(project(as.matrix(blackcap[,c("longitude","latitude")]), projection))
colnames(bonne) <- c("bonnex","bonney")
bc <- cbind(blackcap,bonne)
bfit <- corrHLfit(migStatus ~ 1 + Matern(1|bonnex+bonney), data=bc, HLmethod="ML")

Predictions from the fitted model

xrange <- range(bc$bonnex)
yrange <- range(bc$bonney)
gridSteps <- 101
border <- 20
xGrid <- seq(xrange[1]-border, xrange[2]+border, length.out = gridSteps)
yGrid <- seq(yrange[1]-border, yrange[2]+border, length.out = gridSteps)
newdata <- expand.grid(bonnex=xGrid, bonney=yGrid)
gridpred <- predict(bfit, newdata = newdata)

Raster for predictions

CreateRaster <- function(
    long,
    lat,
    values,
    proj="+proj=longlat +datum=WGS84",
    save.spatial.files=FALSE,
    filename="data_raster",
    overwrite.spatial.files=TRUE
  ) {
    # Args:
    #   long: a vector of the longitudes of the raster cells
    #   lat: a vector of the latitudes of the raster cells
    #   values: a vector of the values of the raster cells
    #   proj: the projection system for the raster
    #   save.spatial.files: logical indicating if 
    #          an hard copy of the raster should be saved (as ascii)
    #   filename: name of the file for the hard copy
    #   overwrite.spatial.files: logical indicating if 
    #          an existing hard copy should be overwritten or not
    #
    # Returns:
    #   The raster.
    #
    data <- data.frame(long=long, lat=lat, values=values)  # a dataframe 
    #                 with longitudes, lattitudes, and values is being created
    coordinates(data) <- ~long+lat  # coordinates are being set for the raster
    proj4string(data) <- CRS(proj)  # projection is being set for the raster
    gridded(data) <- TRUE  # a gridded structure is being set for the raster
    data.raster <- raster(data)  # the raster is being created
    if(save.spatial.files) writeRaster(
      data.raster,
      filename=paste(filename, ".asc", sep=""),
      overwrite=overwrite.spatial.files
    )  # if save=TRUE the raster is exported as an ascii file
    return(data.raster) # the raster is being returned
  }

raster.predict <- CreateRaster(
  long=newdata$bonnex,
  lat=newdata$bonney,
  values=gridpred,
  proj=projection)

Defining the ocean mask consistent with the map projection

The following mask extends beyond the rectangular region shown in later plots.

palette <- spaMM.colors()
  
data(worldcountries)
data(oceanmask)
# All shape files can be found here: http://www.naturalearthdata.com/downloads/
# and loaded into R by e.g.
# oceanmask <- readOGR("ne_110m_ocean.shp", layer="ne_110m_ocean")

## crop the ocean mask before projecting it, otherwise artefacts will appear 
# (1) project world
worldcountries <- spTransform(worldcountries, projection)
# (2) define extent of region to be cropped in projected world
extent.crop <- extent(projectExtent(raster.predict, oceanmask@proj4string))+20
# (3) crop oceanmask according to this extent
oceanmask <- crop(oceanmask, extent.crop)
## Loading required namespace: rgeos
oceanmask <- spTransform(oceanmask, projection)
plot(oceanmask)

Level plot from raster

p <- levelplot( ## filled contour with legend bar
  raster.predict, ## levelplot from rasterVis handles any raster 
  maxpixels=4e6,
  margin=FALSE,
  cuts=length(palette)-1,
  col.regions=palette
)

country.layer <- layer(
  sp.polygons(worldcountries, fill=fill),
  data=list(sp.polygons=sp.polygons, worldcountries=worldcountries, 
            fill="transparent") 
)

p+country.layer

Adding points, countries, and ocean mask

CreateSpatialPoints <-
  function(
    long,
    lat,
    values=-9999,
    proj="+proj=longlat +datum=WGS84",
    save.spatial.files=FALSE,
    filename="my_points",
    overwrite.spatial.files=TRUE
  ) {
    data.sp <- data.frame(long=long, lat=lat, values=values)
    coordinates(data.sp) <- ~long+lat
    proj4string(data.sp) <- CRS(proj)
    if(save.spatial.files)   writeOGR(
      data.sp,
      dsn="./",
      layer=filename,
      morphToESRI=TRUE,
      driver="ESRI Shapefile",
      overwrite_layer=overwrite.spatial.files
    )
    return(data.sp)
  }

my.points <- CreateSpatialPoints( ## userdefined
  long=bc$bonnex,
  lat=bc$bonney,
  proj=projection
)

points.layer <- layer(
  sp.points(my.points, col=col, cex=cex, pch=pch, lwd=lwd),
  data=list(sp.points=sp.points, my.points=my.points, 
            col="white", cex=2, pch=15, lwd=2)
)

ocean.layer <- layer(
  sp.polygons(oceanmask, fill=fill, col=col),
  data=list(sp.polygons=sp.polygons, oceanmask=oceanmask, 
            fill="black", col="transparent")
)
p+ocean.layer+country.layer+points.layer