#working with vectors

#characters:

vChr<-c("One","Two","Three","Four")


#numeric vector
vNum<-c(1,2,3,4)


#creating a data frame
df<-data.frame(vChr,vNum)

head(df)
##    vChr vNum
## 1   One    1
## 2   Two    2
## 3 Three    3
## 4  Four    4
#set column names

colnames(df)<-c("Character","Number")

#selecting rows and columns.

#use squared brackets to select x(rows) and y(columns) as in [x,y] 

#select first two rows of df and all columns (a space before or after the comma means select all of the rows or columns respectively)

df2<-df[1:2,]

#select first three rows and first column:

df3.1<-df[1:3,1]



#subsetting numerics using mathmatical operators
vSub<-subset(df,df$Number<2)

vSub
##   Character Number
## 1       One      1
#subset using "less than or equal to"
vSub1<-subset(df,df$Number <=2)

vSub1
##   Character Number
## 1       One      1
## 2       Two      2
#subsetting with characters
vSub2<-subset(df,df$Character=="One")

vSub2
##   Character Number
## 1       One      1
#subsetting with != ("is not")
vSub3<-subset(df,df$Character!="Four")


vSub3
##   Character Number
## 1       One      1
## 2       Two      2
## 3     Three      3
#using boolean operators "&"(AND) and "|"  (OR)

ANDsub<-subset(df,df$Character=="One" & df$Number==1) #AND

ANDsub
##   Character Number
## 1       One      1
ORsub<-subset(df,df$Number==1 | df$Number==2) #OR

ORsub
##   Character Number
## 1       One      1
## 2       Two      2

Saving your workspace to file using the save.image() function - this allows you to load in your results when you start your next R session

save.image(file = "test.RData")

#loading back in
load("test.RData")

Loading most up-to-date version of packages

Sometimes you may want to use the most updated version of an R package for compatibility purposes. By default packages are downloaded from CRAN (the Comprehensive R Archive Network) when using install.packages() in RStudio. Package developers also host development version of their packages elsewhere (usually on GitHub) and these can be installed using the devTools() and/or the remotes() package.

For example, currently it is best to install development (newest) versions of the raster and sp packages (two packages that we use a lot on the course) due to changes around how the proj library (a library for performing conversions between cartographic projections) handles different CRS (coordinate referencing systems). Many packages (including rgdal) rely on the Proj library in their functions so there is currently A LOT of work going into updating these packages to align with these changes. For the raster and sp packages we can install the latest versions (which deal with these recent changes to some degree at least). Before doing so, it is better to remove the related libraries from your R library folder on your machine. To locate this folder you can type .libPaths() into the console.

When you have removed the relevant libraries then, to load development versions of raster and sp, you can then run:

install.packages("devtools")
install.packages("remotes")
devtools::install_version("sp", version = "1.4-1")
remotes::install_github("rspatial/raster")

Then call the libraries as normal.

Basic model formula in R

Most modeling functions (glm, maxent, randomForest, support vector machine) use the same basic R syntax of:

function(responseVariable~ predictor variables,family,data)

Where family is an optional argument that relates to error distributions (e.g. normal, poisson, binomial) and data is the data frame that contains all your variables.

In our SDM work in Week 4 we used “.,” in place of specifying predictor variables but to select and build up a model formula we can also use individual variable names and link them together with the “+” symbol (as we did in Week 6).An example with a glm model would look like:

glmModel<-glm(Pres~temperature+preciptation+elevation, binomial(link='logit'),data = all.cov)

One advantage of building a model this way is that we can incorporate transformations to our data. For example to model the relationship between temperature and species presence as non-linear (curved) in a distribution model we can use the poly() function (for polynomial) as used in Week 6 E.g.

#quadratic non-linear term
glmModel<-glm(Pres~poly(temperature,2)+preciptation+elevation, binomial(link='logit'),data = all.cov)

#cubic non-linear term
glmModel<-glm(Pres~poly(temperature,3)+preciptation+elevation, binomial(link='logit'),data = all.cov)

Converting between vector and raster

Sometimes it is useful to be able to convert between vector and raster formats in our work. The terra() package provides this functionality. A simple example would be to convert species occurrence points to raster. This might be useful if, for example, you wanted to create a raster layer where the value of each cell denotes the number of occurrences (i.e. abundance) per unit area (i.e. the area represented by each cell).

In the example below we first aggregate our template raster (the raster used to set the dimensions and resolution of the output raster) to an appropriate size (we’ll use 1000m^2 as an example) and then rasterize our points according to this template. We set the “field” argument in the rasterize() function to 1 so that each converted point is given a value of 1 and fun=“sum” so that the values of the resulting raster cells represent the sum (count) of points inside the area of each cell.

library(terra)
## terra 1.7.3
#read in the data
sciurus<- read.csv("Sciurus.csv")

#subset the data to only include points with complete coordinates.
sciurus<-subset(sciurus,!is.na(sciurus$Latitude))

#remove all points with uncertainty > 100m
sciurus<-subset(sciurus, sciurus$Coordinate.uncertainty_m<100)

#note that you can also use braces to do thsi:

sciurus<-sciurus[sciurus$Coordinate.uncertainty_m<100,]

#make vector layer

#create crs object

sciurus.latlong<-data.frame(x=sciurus$Longitude,y=sciurus$Latitude)



#Use coordinates object to create our spatial points object

sciurus.sp<-vect(sciurus.latlong,geom=c("x","y"))

#check that the points now have our desired crs. 
crs(sciurus.sp)<-"epsg:4326"


#Crop and aggregate the raster layer

studyExtent<-c(-4.2,-2.7,56.5,57.5) #Set the extent to something workable

#now crop points to studyExtent

C<-crop(sciurus.sp,studyExtent)


LCM<-rast("LCMUK.tif") #we'll use this raster as the template in the rasterize() function. 

LCM<-LCM$LCMUK_1

#set CRS to that of LCM


sciurusFin<-terra::project(C,crs(LCM))

#get coordinates
sciurusCoords<-crds(sciurusFin)

#set cropping extent
x.min <- min(sciurusCoords[,1]) - 5000
x.max <- max(sciurusCoords[,1]) + 5000
y.min <- min(sciurusCoords[,2]) - 5000
y.max <- max(sciurusCoords[,2]) + 5000

extent.new <- ext(x.min, x.max, y.min, y.max)


LCM <- crop(LCM, extent.new)

plot(LCM$LCMUK_1)

LCM<-terra::aggregate(LCM,fact=40,fun="modal")#aggregate to 1km^2 cells using band 1. Note fun=modal for a categorical raster. 

sciurusR<-rasterize(sciurusFin,LCM,field=1,fun="sum",background=0) #create raster using LCM as a template and where the value of each cell is the count of points in the cell area.

One reason why converting to raster can be useful when working in R is that we can create distance rasters from the results, useful for a range of ecological investigations. The distance() function in the raster package can do this but I suggest you use the gridDistance() function (from the same package) as it works much faster!

sciurusR_Pres<-sciurusR #create copy of raster
  
sciurusR_Pres[sciurusR_Pres>0]=1 # truncate raster to presence/absence (1-0)


sciurusDist<-gridDistance(sciurusR_Pres,target=1)
## Warning: [gridDistance] 'gridDistance' was renamed to 'gridDist'.
## 'gridDistance' will be removed in a future version
plot(sciurusDist)

plot(sciurusFin,add=TRUE)

Working with the focal() function to create buffers

The focal function can be used to create rasters that describe a “neighbourhood” around each cell of an input raster. The focalMat() function is used to set the template raster (to determine the resolution, extent etc.), the radius of the neighbourhood and the type of buffer (“circle” for a standard buffer, “Gauss” for a Gaussian filter). The default for the focalMat function is that all values within the focal neighbourhood sum to one, so for a binary raster (like our woodland layer in week 2) if all cells in the neighbourhood around the focal cell contain woodland then the focal cell receives a value of 1 (i.e. 100% cover). If we want to return the sum of the raster values (rather than a proportion between 0 and 1) we can truncate the focalMatrix (i.e. the neighbourhood) so that each cell in the neighbourhood gets a values for 1 if it is in the neighbourhood, and 0 otherwise. This can be useful if working with rasters that represent counts (for example, of species occurrence). The code below uses the ifelse() function to set all values >0 to one, else it returns 0. The resulting raster values created by the focal function then represents a count (of red squirrels in this case) for the neighbourhood around each cell.

fwRed<-focalMat(sciurusR,2000,type= 'circle') #set the neighbourhood
head(fwRed) #inspect
##            [,1]       [,2]       [,3]       [,4]       [,5]
## [1,] 0.00000000 0.00000000 0.07692308 0.00000000 0.00000000
## [2,] 0.00000000 0.07692308 0.07692308 0.07692308 0.00000000
## [3,] 0.07692308 0.07692308 0.07692308 0.07692308 0.07692308
## [4,] 0.00000000 0.07692308 0.07692308 0.07692308 0.00000000
## [5,] 0.00000000 0.00000000 0.07692308 0.00000000 0.00000000
#re-scale weight matrix to 0,1
fwSciurusBuffer <- ifelse(fwRed > 0, 1, 0)
head(fwSciurusBuffer)#inspect binary neighbourhood
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0    1    0    0
## [2,]    0    1    1    1    0
## [3,]    1    1    1    1    1
## [4,]    0    1    1    1    0
## [5,]    0    0    1    0    0
focalRed<-focal(sciurusR,fwSciurusBuffer,sum)#create focal raster

plot(focalRed)#plot map

Converting from raster to vector is as easy as using the rasterToPoints() function. We set “spatial =TRUE” so that a SpatialPointsDataFrame is returned. Other functions Note also that to avoid returning a point for every raster cell, we can specify a function that will return only cells with positive values (i.e. containing at least one occurrence point):

sciurusNA<-sciurusR #make copy

sciurusNA[sciurusNA<1]=NA #remove NAs so only cells with values are returned

redPoints<-as.points(sciurusNA, values=FALSE,na.rm=TRUE, na.all=FALSE)

plot(LCM)

plot(redPoints,add=TRUE)

Extending options in the adehabitat kernel density estimation function.

We can extend some of the functionality of the kernelUD approach to density estimation used in Week 3 in order to produce raster objects of the resulting density estimate. For example we can create use an existing layer to set the resolution of an output density raster. First, create a “SpatialPixelsDataFrame” (similar to a raster object but compatible with the adehabitat package), from a raster of the extent and resolution we want to work with, then use this as a template for the kernelUD function before converting the output back to raster (for use with other functions).

library(adehabitatHR)
## Loading required package: sp
## Loading required package: deldir
## deldir 1.0-6      Nickname: "Mendacious Cosmonaut"
## 
##      The syntax of deldir() has had an important change. 
##      The arguments have been re-ordered (the first three 
##      are now "x, y, z") and some arguments have been 
##      eliminated.  The handling of the z ("tags") 
##      argument has been improved.
##  
##      The "dummy points" facility has been removed. 
##      This facility was a historical artefact, was really 
##      of no use to anyone, and had hung around much too 
##      long.  Since there are no longer any "dummy points", 
##      the structure of the value returned by deldir() has 
##      changed slightly.  The arguments of plot.deldir() 
##      have been adjusted accordingly; e.g. the character 
##      string "wpoints" ("which points") has been 
##      replaced by the logical scalar "showpoints". 
##      The user should consult the help files.
## Loading required package: ade4
## Loading required package: adehabitatMA
## Registered S3 methods overwritten by 'adehabitatMA':
##   method                       from
##   print.SpatialPixelsDataFrame sp  
##   print.SpatialPixels          sp
## 
## Attaching package: 'adehabitatMA'
## The following object is masked from 'package:terra':
## 
##     buffer
## Loading required package: adehabitatLT
## Loading required package: CircStats
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:terra':
## 
##     area
## Loading required package: boot
library(raster)
## 
## Attaching package: 'raster'
## The following object is masked from 'package:MASS':
## 
##     select
sciurusR<-raster(sciurusR)#adehabitat works with raster library

pixelSciurus<-as(sciurusR,"SpatialPixelsDataFrame") # create SpatialPixelsDataFrame from raster

sciurusFinSP<-as(sciurusFin,"Spatial")

redDen<-kernelUD(sciurusFinSP, h = "href", kern ="bivnorm",grid=pixelSciurus) # carry out the kernel density function setting "grid = pixelSciurus" to set the resolution to that of the original sciurusR object.

rasDen<-raster(redDen) # convert back to raster

res(rasDen) # check resolution
## [1] 1000 1000
plot(rasDen) # plot

Density in Spatstat

Another option for creating density surfaces is, when working with the Spatstat package, to use the density.ppp() function:

library(spatstat)
## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## spatstat.geom 3.0-6
## 
## Attaching package: 'spatstat.geom'
## The following objects are masked from 'package:raster':
## 
##     area, rotate, shift
## The following object is masked from 'package:MASS':
## 
##     area
## The following object is masked from 'package:ade4':
## 
##     disc
## The following objects are masked from 'package:terra':
## 
##     area, delaunay, rescale, rotate, shift, where.max, where.min
## Loading required package: spatstat.random
## spatstat.random 3.1-3
## Loading required package: spatstat.explore
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:raster':
## 
##     getData
## spatstat.explore 3.0-6
## 
## Attaching package: 'spatstat.explore'
## The following object is masked from 'package:boot':
## 
##     envelope
## Loading required package: spatstat.model
## Loading required package: rpart
## spatstat.model 3.2-1
## Loading required package: spatstat.linnet
## spatstat.linnet 3.0-6
## 
## spatstat 3.0-3 
## For an introduction to spatstat, type 'beginner'
library(maptools)
## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired during 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.
sciurusWin<-as.owin(as.im.RasterLayer(sciurusR))

pppSciurus<-ppp(sciurusCoords[,1],sciurusCoords[,2], window = sciurusWin)
## Warning: data contain duplicated points
plot(pppSciurus)

redDensity<-density.ppp(pppSciurus)

plot(redDensity)

Plotting the results of iGraph objects onto a map

Usually it is helpful to plot the results of patch-based graphs onto another map for context. This can be done using Add=TRUE as for other spatial object you have worked with.

Note that, in the latest version of the iGraph package, coordinates, if provided, are automatically scaled to an arbitrary range of -1 to 1. I am not sure why this change has been implemented (perhaps it is more intuitive for some kinds of graphs such as those related to social networks for which geographical location is typically not needed). So, to get our graphs onto a map that you can produce with correct coordinates, we need to add an extra command at AFTER “add=TRUE” and that is “rescale=FALSE” (see below).

First I’ll plot the national parks and conifer patches for context:

library(gdistance)
## Loading required package: igraph
## 
## Attaching package: 'igraph'
## The following objects are masked from 'package:spatstat.geom':
## 
##     diameter, edges, is.connected, vertices
## The following object is masked from 'package:raster':
## 
##     union
## The following objects are masked from 'package:terra':
## 
##     blocks, union
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
## Loading required package: Matrix
## 
## Attaching package: 'gdistance'
## The following object is masked from 'package:igraph':
## 
##     normalize
## The following object is masked from 'package:terra':
## 
##     costDistance
library(igraph)

#read in land cover data
clc<-rast("CLC_UK.tif")



library(terra)

#read in park boundaries
Cairn<-vect("SG_CairngormsNationalPark_2010.shp")

Lomond<-vect("SG_LochLomondTrossachsNationalPark_2002.shp")

#read in confier patches
Conifer<-vect("Conifer_Patch.shp")

# combine park spatvectors
Parks<-union(Cairn,Lomond)

###############################################################################################
#############################Converting between SP and Spatvector Objects#######################
################################################################################################


#Converting between libraries for creating spatialobjects can be important as some analytical libraries only work with one of them (e.g. the gdistance library for creating least cost paths currently only works with RasterLayer objects created using the raster package)

#if you need to convert between SpatVector and SpatialPolygonsDataFrame you can do so with the as() function:
CairnSP<-as(Cairn,"Spatial")

LomondSP<-as(Lomond,"Spatial")

##To convert the other way:

Cairn<-vect(CairnSP)

#With raster data it's even easier:

#from SpatRaster to Raster:
clcRaster<-raster(clc)

#from Raster to SpatRaster:
clc<-rast(clcRaster)

# Next crop the CLC dataset to some suitable extent (you choose, below is just an example)
plot(ScotCLC)

plot(Cairn,add=TRUE,border="Red") # plot the park boundary and give it a red border


plot(Lomond,add=TRUE,border="Blue") # plot park boundary and give it a blue border

plot(Conifer,add=TRUE)

Now add the graph (I have not shown the code for making the i.graph objects - that’s your job!).

Note that setting vertex.label.colour and edge.color to “transparent” effectively removes these elements which gives me a less cluttered output. Here the size of the nodes (i.e. patches) is is set to a fixed value for simplicity, in your actual work you might use the respective degree centrality in the network to set the size (e.g. how many links each node shares with other nodes in the network of patches).

plot(ScotCLC)

plot(Cairn,add=TRUE)


plot(Lomond,add=TRUE)

plot(graph.Amean, layout=coords,
     vertex.label.dist=0.5,vertex.label.color="transparent",vertex.size=50^3,edge.color="transparent", add=TRUE,rescale=FALSE)

Saving results to table

Often we may want to tabulate and save our data as a .csv file or similar. We can save results by creating a data frame from them and then writing to file using write.csv() (see the Week 1 practical for the snippet on how to save vector and raster objects to shapefiles and geotiffs). Also note that the values below for betweenness may differ from your own (I just arbitrarily chose a dispersal distance for this example).

Amean.between <- betweenness(graph.Amean)

betweenData<-data.frame(site=1:nrow(coords),coords=coords,betweeness=Amean.between)

head(betweenData)
##   site coords.x coords.y  betweeness
## 1    1 291174.8 729602.5    6.469444
## 2    2 292037.0 729524.3    6.469444
## 3    3 313418.0 729747.1  303.220262
## 4    4 302258.6 730956.7 3331.614536
## 5    5 322447.8 731188.0 5768.702975
## 6    6 312677.1 731391.8 1474.829375
write.csv(betweenData,"betweenessData.csv")

Creating a map legend in R

legendclc<-read.csv("C:/SpatialEcology/clc_legend.csv")#read in spreadsheet with legend key



legendclc<-legendclc[,c(1,4)] #select the relevent columns (ID and class names) 

head(legendclc)
##   CLC_CODE                         LABEL3
## 1      111        Continuous urban fabric
## 2      112     Discontinuous urban fabric
## 3      121 Industrial or commercial units
## 4      122     Road and rail network land
## 5      123                     Port areas
## 6      124                       Airports
ScotCLC<-as.factor(ScotCLC) #convert to factor to access "level" (classes)
landCoverclc<-levels(ScotCLC)[[1]]#create object of landcover classes

legendScot <- legendclc[legendclc$CLC_CODE %in% landCoverclc$ID,]#use %in% command to match classes in the raster to those in the legend


par(mar=c(5, 5, 5, 5)) #set margins (you need to experiment a bit here)

spNodes<-vect(coords)# create vector object from the patch coordinates

plot(ScotCLC,legend = FALSE,col = rainbow(33)) #plot the raster using native colour palette "rainbow" for 33 classes (number of) landcovers in the raster

###Using RcolourBrewer Package to set the colour palette#########.

library(RColorBrewer) #load library

##Use colourRampPalette to create custom palette
coul <- colorRampPalette(brewer.pal(8, "Accent"))(33) #select 8 colours from the "Accent" palette and create 33 combinations

plot(ScotCLC,legend = FALSE,col = coul)#plot again with the new palette




plot(spNodes,pch=21,cex=Amean.between^0.1,col="black", bg="white", lwd=1,add=TRUE)#add the patch nodes, set pch to 21 for a filled point with outline then set fill to white and the outline to black and lwd (line width) = 1



####Add a legend using the legend function. The first argument sets the coordinate for the position (in the map units), the second argument specifies the words that go in the legend (here I use the paste function to simply copy from the legend object made earlier). Cex sets the font size, fill sets the colour of the icons and bty="n" removes the box around the legend.
legend(342000,790000,paste(legendScot[,2]), fill = coul,
       cex = 0.55,bty="n")


#Finally add an entry to symbolize the betweenness layer.
legend(x=183000,y=740000,legend=c("Patch betweenness"),pch=21,col="black",bg="white",cex=0.9,bty="n")

If you want to get more into colour palettes in R, there is a nice cheat sheet with some useful links here: https://www.nceas.ucsb.edu/sites/default/files/2020-04/colorPaletteCheatsheet.pdf

For a guide to the various point character and other graphical parameters see here: https://www.statmethods.net/advgraphs/parameters.html