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

Loading packages

library(spaMM)
## spaMM (Rousset & Ferdy, 2014, version 3.1.2) is loaded.
## Type 'help(spaMM)' for a short introduction,
## 'news(package='spaMM')' for news,
## and 'citation(spaMM)' for proper citation.
library(rgdal)
## Loading required package: sp
## rgdal: version: 1.3-6, (SVN revision 773)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/Francois.rousset/Documents/R/win-library/3.6/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/Francois.rousset/Documents/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.3-1
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 object worldcountries and oceanmask to be loaded are data objects created from public domain shapefiles downloaded from http://www.naturalearth.com on 2015/10/21 and included in spaMM up to version 3.0.0. Higher-resolution files are available on http://www.naturalearthdata.com, from where they can be loaded as shown below. or inlusin into spaMM, worldcountries had to be edited for non-ASCII characters before inclusion in spaMM: worldcountries@data$formal_fr was removed and the Côte d’Ivoire level of some factor variables was renamed.

load(url('http://kimura.univ-montp2.fr/~rousset/spaMM/worldcountries.RData'))
load(url('http://kimura.univ-montp2.fr/~rousset/spaMM/oceanmask.RData'))
  
# All shape files can be found here: http://www.naturalearthdata.com/downloads/
# and unedited versions can be loaded into R by e.g.
# if (requireNamespace("rgdal", quietly = TRUE)) {
#   oceanmask <- readOGR("ne_110m_ocean.shp", layer="ne_110m_ocean")
#   worldcountries <- readOGR("ne_110m_admin_0_countries_lakes.shp",
#                            layer="ne_110m_admin_0_countries_lakes")
# }

The following mask extends beyond the rectangular region shown in later plots. We crop the ocean mask before projecting it, otherwise artefacts will appear.

palette <- spaMM.colors()
# (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