Preparation

In this practical we will work with different types of data: GPS tracking data on individual Oystercatchers, GIS data (bathymetry) of the Balgzand, and sampling data on the benthic food of Oystercatchers. Before we start, first make sure that all packages are installed & loaded.

install.packages(c('sp','gstat','adehabitatHR','raster','maptools','StreamMetabolism','sampling','lattice','latticeExtra'))
library(sp)
library(gstat)
library(adehabitatHR)
library(raster)
library(maptools)
library(StreamMetabolism)
library(sampling)
library(lattice)
library(latticeExtra)

Next, we load the data with which we’ll work in this practical.

setwd('~/Dropbox/')

# download data
download.file("http://adriaandokter.com/SPEE15/Rpractical/Balgzand2009.asc", destfile = "Balgzand2009.asc")
download.file("http://adriaandokter.com/SPEE15/Rpractical/practical_functions.R", destfile = "practical_functions.R")
download.file("http://adriaandokter.com/SPEE15/Rpractical/DataSampling.RData", destfile = "DataSampling.RData")
download.file("http://adriaandokter.com/SPEE15/Rpractical/PatchPolygons.RData", destfile = "PatchPolygons.RData")
download.file("http://adriaandokter.com/SPEE15/Rpractical/GPSData.all.RData", destfile = "GPSData.all.RData")
## load the data
BathymetryMap<-readAsciiGrid('Balgzand2009.asc',colname='bathymetry')
load("DataSampling.RData")
load("PatchPolygons.RData")
load("GPSData.all.RData")
source("practical_functions.R")

Read in the bathymetric map, and plot it

BathymetryMap<-readAsciiGrid('Balgzand2009.asc',colname='bathymetry')
image(BathymetryMap)

The practical is organised in three parts:

Part I: kernel home range estimation

The GPS tracking data on individual Oystercatchers is stored in the GPSData.all dataframe. Have a look at its structure, and plot the GPS points on the base map:

str(GPSData.all)
points(GPSData.all[,c("RDX","RDY")],col='red',pch=".")

First we make a subselection of the data, to exclude the two birds that lived outside the core study area of the western Balgzand. We also convert from the simple data.frame class to a SpatialPointsDataFrame class object, by which we tell R that these are spatial data.

# Select data only in the core area on Balgzand:
GPSData=GPSData.all[GPSData.all$RDX>115000 & GPSData.all$RDX<120000 & GPSData.all$RDY>548000 & GPSData.all$RDY<554000,]
# The above subselection has removed some birds
# the command droplevels removes factor levels (of field ID) that are no longer used
GPSData=droplevels(GPSData)
# convert to SpatialPointsDataFrame object,
# assigning RDX and RDY (Dutch Amersfoort coordinates) as the spatial coordinates
coordinates(GPSData) <- ~ RDX+RDY

First we make a grid with the right resolution and extent, on which we will estimate each animals utilisation distribution. We choose a resolution of 20 meter, the same resolution as the bathymetric basemap.

# Extra info about mkraster (a user function in Practical_Functions.r):
#
# The default resolution (parameter name: 'resolution) is 100 m 
# and the default range (parameter name: 'rng') in mkraster is
# 'full' - this the entire space covering the Balgzand base map.
#
# Other options for the rng-parameter are 
# 'core', 'texel' and 'wadden'
#
# Alternatively you dan also supply a spatial object
# to the function (the parameter for that is 'data'), 
# then the extent of the domain will be the bounding box of 
# the data in the spatial object, e.g.  
#   allgrid = mkraster(data=SpatialObject, resolution=500)

fullgrid = mkraster(resolution=20, rng='core')

Apply the function kernelUD for estimation of the utilisation distribution. We can set various parameters on this function, most importantly the properties of the kernels that we will overlay with the GPS data points:

# manual bandwidth of the kernel
kd.manual = kernelUD(GPSData['ID'],grid=fullgrid,h=100)
# automated estimation of the bandwidth
kd.auto = kernelUD(GPSData['ID'],grid=fullgrid,h="LSCV",hlim = c(0.001, 1.5))
# print the identifiers of the birds in the dataset (as a reminder):
unique(GPSData$ID)
# print the estimated bandwidth parameter h (example for bird 406)
kd.auto$"406"@h$h

Now that we have estimated the utilisation distribution, we can calculate home ranges (at any percentage level you want) by intersecting it:

# calculate the kernel home range (in ha) for each individual, at the 50% and 80% percentage levels
kd.manual.area <- kernel.area(kd.manual, percent=c(50,80))
# calculate the contours of the homeranges for plotting
khr80=getverticeshr(kd.manual, percent=80)

And we may plot the home range polygons on the base map, together with the points on which it was estimated:

# plot the basemap
image(BathymetryMap,xlim=c(115000,120000),ylim=c(547000,555000))
# plot the homerange (example again for bird 406)
plot(khr80["406",],add=T)
# plot the underlying data points
points(GPSData[which(GPSData$ID=="406"),]@coords,col='red',pch='.')
  1. Estimate 80% kernel home ranges with a bandwidth parameter of 200 meter and 20 meter. What are the differences, and which bandwidth do you think is more appropriate given the data?

  2. Explore how home range sizes varies at different percentage levels. Is there a cut-off point when the size suddenly becomes much larger? If so, why would this happen?

The Uva-Bits GPS tags that were used also record whether the birds are moving or not, using an accelerometer that measures body acceleration. We can thus remotely detect whether a bird is moving because it is searching for food / handling prey, or that it is standing still, e.g. during high-tide roosting. In the dataset the activity is stored in the ACCFractionActive field of the GPS data frame.

  1. Inspect the birds activity values and define a threshold to classify birds as active or inactive, and explore the point patterns on the basemap. Are the locations where birds are active or inactive obviously different?

  2. Split the utilisation distribution into two components, by estimating separately a distribution for the times when the birds were active, and the times when birds were standing still. Estimate the kernel bandwidth for each distribution by least-squares cross-validation. Are estimated bandwidths are very different in the two cases? Explain whether you can expect such a difference in estimated bandwith, given the differences in the Oystercatchers movement pattern between the two cases.

  3. Is the kernel home range of birds during low tide (when birds are foraging) different in size than during high tide (when birds cannot forage)? Statistically test this difference.

Oystercatchers can be highly territorial, and will defend a preferred foraging area against competitors. These are often local breeding birds. Wintering birds and juveniles arrive from inland breeding areas in the Netherlands, Germany and Scandivinavia often do not have established territories, and may be more mobile.

  1. Based on your estimation of the home ranges sizes, which birds would you classify as territorial, and which birds not? Is a sample size of 8 birds sufficient to study individual differences in movement patterns?

Besides home range size, other measures of similarity between home ranges can be defined. A common measure to quantify home range overlap is Bhattacharyya’s affinity. This statistic is defined as follows: \[BA=\int_{-\infty}^\infty \int_{-\infty}^\infty \sqrt{UD_1(x,y)} \sqrt{UD_2(x,y)} dx dy\]

  1. Given that utilisation distributions are normalised probability density functions, and the above definition, which range of values can BA assume? Which value indicates complete overlap?

In R you may calculate the overlap between all pairs of estimated utilisation distributions as follows:

ko=kerneloverlaphr(kd.manual,method="BA")
  1. Based on the overlap between all pairs, could you group individuals in different categories with a different use of the Balgzand?

Because intertidal birds like Oystercatchers can only forage during low tide, their daily foraging routines are much more in sync with the alternation of the tides, than with the alternation between day and night. Nonetheless the foraging conditions change dramatically at night (e.g. in terms of vision, or predators), which may have consequences for the space use of the birds.

Sunrise and sunset may be calculated at Balgzand with the following function, which is already loaded in the header:

# plot the basemap
SunRiseSet = function(date,lat=52.927336,lon=4.849632){
  sunrise.set(lat,lon,strftime(date,format="%Y/%m/%d"))
}
  1. (Extra) Explore the differences in the point patterns of the individuals across day and night

  2. (Extra) Are the home range sizes different between day and night.

Part II: the food landscape of the Oystercatchers

To understand the movement patterns in the intertidal, knowledge on the spatial distribution of food is critically important. The Balgzand was therefore visited in the autumn of 2012 to sample the benthic bivalves that Oystercatchers use: Cockles (Cerastoderma edule) and - as we found out - the spat (new recruits) of the razor shell clam Ensis directus. Out of time constraints not the full Balgzand could be surveyed, and sampling was confined to several patches that overlayed with the kernel home ranges of the Oystercatchers. Patches were chosen to extend beyond the homeranges, to also include sufficient area in the sampling program that the birds were not using.

The outline of the sampled patches can be drawn on the basemap using the PatchPolygons object, that was loaded with the other data:

## food landscape has been sampled to include as much of the Oystercatcher home ranges as possible:
# plot the basemap
image(BathymetryMap,xlim=c(115000,120000),ylim=c(547000,555000))
# plot the outline of the sampled patches (or 'beds', in dutch 'bank')
plot(PatchPolygons,add=T)
# overlay the GPS points:
points(GPSData@coords,col='red',pch=".")

We will use regression kriging to make an interpolated map of cockle density (CockleNDensity) and cockle size (CockleLengthMean).

Step 1: fit a variogram model (example for the cockle density data):

## convert DataSampling, the data frame containing the information collected in the 
## sampling program on Balgzand to a SpatialPointsDataFrame
coordinates(DataSampling) <- ~ RDX+RDY
## initialize a gstat object
g <- gstat(id="CockleNDensity", formula=CockleNDensity ~ 1, data=DataSampling)
## define a sensible maximum cutoff distance up to which you want to plot/fit the variogram
## (side question: why would you want to limit this range in our arrangement with multiple patches?)
mycutoff=500
## calculate the variogram
v<-variogram(g, cutoff=mycutoff, width=25)
## plot the variogram
plot(v)
## lets fit a spherical variogram model to the variogram (nugget is a dummy value here):
v.fit <- fit.variogram(v, model=vgm(model='Sph',range=mycutoff,nugget=500))
## update the gstat object, by adding the fit
g <- gstat(g, id="CockleNDensity", model=v.fit )

Step 2: inspect the fit, and adjust the variogram model if necessary (back to Step 1)

## plot the variogram with fit
## calculate the variogram asymptote (question!)
v.asymptote=var(DataSampling$CockleNDensity)
# plot the fit
p1=plot(v, model=v.fit, as.table=TRUE,ylim=c(0,1.1*v.asymptote),xlab="Distance [m]",ylab="Semivariance Cockle density")    
# plot the asymptote
p2=xyplot(rep(v.asymptote,2)~c(0,1000),type="l",col="red")
# both together
p1+p2

Step 3: Use regression kriging to interpolate on a grid, under assumption of the estimated variogram model of step 2:

## create a grid onto which we will interpolate:
grid = mkraster(data=DataSampling,resolution=25)
## interpolate using the estimated variogram model (may take a minute or so):
g.predict <- predict(g, newdata=grid)

Step 4: Plot the results:

## some settings for the color scale:
trellis.par.set(sp.theme(set = FALSE,frame.plot=F,regions = list(col = rev(bpy.colors(n=102,cutoff.tails=0)))))
## let's plot the results for bed/patch/bank 8 as an example:
bank=8
# extract the bounding box of bed 8:
bbox=DataSampling[which(DataSampling$Bank==bank),]@bbox
# also plot the location of the sampling stations
samplingstations <- list("sp.points", DataSampling, pch = 4, col = "black", cex=0.5)
# plot it
spplot(g.predict, zcol="CockleNDensity.pred", at = seq(-10,250,10), sp.layout=list(samplingstations), contour=TRUE, labels=FALSE, pretty=TRUE, col='brown', main='CockleNDensity - Kriging prediction',xlim=bbox["RDX",],ylim=bbox["RDY",])
  1. In this study we had a maximum of weeks to sample the food landscape, because over longer time scales the food landscape starts to change as a result of predation / depletion. If it takes approximately two minutes to collect a sample (once a boat has brought you to the right mudflat), how many samples do you think you can (roughly) collect per person in two weeks? How many people do you estimate were in the team that collected the sampling data in this study?

  2. Most of the sampling was performed on a 50 meter grid. Given the shape of the variogram for cockle density, was this a good choice? Would you have preferred sampling on a grid with larger spacing (such that a larger area of the Balgzand and the homerange of the Oystercatchers could have been covered?)

  3. Explain why the sample variance approaches the asymptote of the variogram

The bivalve food stocks in the intertidal go through a very interesting reproductive cycle. Only once every few years there is major recruitment (spat fall) of one (or several) bivalve species. In those years many new Cockle beds establish, while in other years spat fall is almost absent. This highly pulsed recruitment - which is not properly understood yet - leads to strong spatial structures in the age of the bivalves (Cockle spat is smaller than 1 cm in length, 2y Cockles typically around 2 cm, and >2y Cockles 3cm or more).

  1. Also make a kriging interpolation for the size of the cockles, i.e. CockleLengthMean. Compare the variogram for Cockle density and Cockle size. What are the differences, and can you hint at a cause for these differences, given the recruitment cycle of the Cockle?

  2. (extra) Because predation on Cockles is stronger on patches that are long exposed by the tides, you expect a strong trend with altitude in the abundance of Cockles. How are the cockles distributed in altitude? Explore whether the technique of universal kriging provides benefits, as compared to ordinary kriging.

Part III: Predicting Oystercatcher distributions using optimality theory

A common functional form for a Type II functional response is the so-called Holling disk equation (it is called the disc equation because of foraging experiments performed by Holling, where he had a blindfolded person (i.e. his secretary) forage on hypothetical resources, consisting of round paper discs). Almost all functional responses for (shore)birds are described in the literature by different parametrisations of this equation.

\[f(n)=\frac{a n}{1+a h n}\] with n the bird density, a the so-called attack rate (an is called the encounter rate) , and h the handling time of the prey. Handling time is defined as the time spent opening and consuming a prey item after encountering it.

# example to plot an analytical function:
f <- function(x,a,h) (a*x)/(1+a*h*x)
curve(f(x,1,1),from=0,to=10)
curve(f(x,1,2),from=0,to=10,add=T)
  1. Plot the equation and verify it has the shape of a Type II functional response. Why does the functional response reach a plateau? What variable determines the level of the plateau?

  2. Holling’s disk equation assumes that attack rates and handling times are independent of prey density. Are these reasonable assumptions, under which circumstances?

Although Oystercatchers are specialized in eating large bivalve prey like Cockles, which they open with their powerful bill. Based on both field observations and experimental studies, the following expression for the capture rate (i.e. the rate at which Cockle preys are swallowed) of Cockles has been determined

# capture rate (in cockles captured per second)
fCockle <- function(density,length) 0.000860*density/(1 + 0.0001897*density*length^1.792)
  1. plot the capture rate that Oystercatchers can realize for several realistic size and density ranges

A parameter highly relevant to a foraging Oystercatcher is the intake rate it can realize on different prey items. Intake rate is often considered the central currency in optimal foraging theory. The energetic intake for the Oystercatcher depends not only on the capture rate, but also on the energetic contents of the prey. The energetic content is usually expressed in AFDM (ash-free dry mass, i.e. the dry weight of the flesh, corrected for components that cannot be burnt, like remains of sand). The following equation converts the length of a Cockle to AFDM:

# the body mass index of Cockles in autumn, as measured on Balgzand
BMICOCKLE=9
# conversion of Cockle length [mm] to AFDM (ash-free dry mass) [mg] given a BMI value
Length2AFDM <- function(length) BMICOCKLE*(length^3)/(10^3) 

Capture rate and flesh content can be combined into an intake rate in the absence of conspecific competition. This intake rate we call the interference-free intake rate (IFIR):

IFIR <- function(density,length) fCockle(density,length)*Length2AFDM(length)
  1. plot the intake rate in AFDM/s that Oystercatchers can realize for several realistic size and density ranges. Is the longer handling time for large Cockles compensated by a higher flesh content, or not?

We are now ready to convert the landscape prey densities, as reconstructed in part II, into a landscape of intake rate, as perceived by Oystercatchers. Let us focus - for the sake of simplicity - on one of the patches in the study area, patch number 4.

# get the bounding box of patch 4 from the Polygon object (consists of two polygons, 04 and 04B)
bbox=PatchPolygons[c("04","04B")]@bbox
# make a subset of the predicted cockle density map, and store it in patch4.pred object
tmp=g.predict["CockleNDensity.pred"]
patch4.pred=tmp[tmp$x>bbox["x","min"] & tmp$x<bbox["x","max"] & tmp$y>bbox["y","min"] & tmp$y<bbox["y","max"],]
# remove negative Cockle density, which are non-sensicle, by setting them to zero.
# (side question: how is it possible that kriging produces negative values?)
patch4.pred@data$CockleNDensity.pred=pmax(patch4.pred@data$CockleNDensity.pred,0)
# specify the range of Cockle density values to display
breaks=seq(-10,200,10)
# plot the image
image(patch4.pred,xlim=bbox["x",],ylim=bbox["y",],col=terrain.colors(length(breaks)-1),breaks=breaks)
# overlay the GPS data points
points(GPSData@coords,col='red',pch=".")

Since the functional response depends on cockle size, this information needs to be added:

# For now, let's assume all cockles have a length of 30 mm
# Use your results from Part II to fill in the real values
patch4.pred@data$CockleLengthMean.pred=30
  1. Use your results from part II to calculate the landscape of interference-free intake rates
#  the area of a single grid cell
PatchSize=patch4.pred@grid@cellsize[1]*patch4.pred@grid@cellsize[2]

We can combine all the results to calculate the ideal-free distribution. The intake rate of Oystercatchers does not only depend on the food landscape, but also on the number of competitors present. This effect of competition (also called interference) has been shown to follow an exponential relationship, as shown by both experimental and field studies: \[I(n)=\exp(-k n)\] with n the bird density and k the exponential interference constant. The intake rate including interference effects (IR) then equals: \[IR(n,s)=IFIR(n,s)\times I(n)\] with s as before is the size of the prey.

  1. When there is no competition or costs of interference (i.e. k approaches zero), how do you expect that identical, intake-rate maximizing birds will distribute?

  2. Assume we have a system with 2 patches, both of size 1 (we do not worry about dimensions for now). Further assume the exponential interference constant k=.01, and the interference free intake rate in the patches are 3 and 2. If 100 birds are released in the system, how will they distribute? Note: since all birds are identical, they will distribute such that they all have the same, maximized intake rate. Hint: express the bird density in terms of the interference-free intake rate. To find the maximized intake rate, invoke the condition that the bird densities multiplied the patch areas should add up to the total number of birds present (100).

Let’s now calculate the IFD numerically. We use the function IFD (this function is coded in practical_functions.R), which calculates the ideal free distribution for a set of patches, given their interference-free intake rate, area, total number of birds, and interference constant k. The bird densities that IFD returns are in [birds/m^2], patch sizes should be in [m^2] and the interference constant k in [m^2/bird]

CockleLenghts=patch4.pred@data$CockleLengthMean.pred
CockleDensities=patch4.pred@data$CockleNDensity.pred
CockleIFIR = IFIR(CockleDensities,CockleLenghts)
Areas=rep(PatchSize,length(CockleDensities))
# calculate the ideal-free distribution, assuming 400 birds on the patch, and an
# exponential interference constant of k=100 m^2/bird
IFD.pred=IFD(CockleIFIR,Areas,400,100)
# add the IFD prediction to the patch4.pred object:
patch4.pred@data$IFD.pred=IFD.pred$animaldensity
# plot the IFD
breaks=seq(0,max(IFD.pred$animaldensity),max(IFD.pred$animaldensity)/100)
image(patch4.pred["IFD.pred"],xlim=bbox["x",],ylim=bbox["y",],col=terrain.colors(length(breaks)-1),breaks=breaks)
# overlay the GPS points
points(GPSData@coords,col='red',pch=".")

Supporting field observations in patch 4 have indicated that during low tide around 400 Oystercatchers are foraging in this patch. Also, experimental studies in aviaries indicate that the interference constant as a result of direct competition equals k=5.

  1. Predict the distribution of the 400 Oystercatchers in patch 4 for these parameters. What are the most striking differences when compared with the GPS data?

  2. Adjust the interference constant such that it best matches the GPS data (rough estimate required only). If your estimate and the model were to be true, at which distance would Oystercatchers start to experience density dependent competition effects? Is it reasonable to assume this density dependence is the result of direct interactions, like kleptoparasitism?

  3. The assumption of the ideal-free distribution is that birds are both ‘ideal’ and ‘free’. Are these good assumptions for an Oystercatcher? How do you expect the Oystercatcher distribution will change if birds do not behave ideal? And if they do not behave free?

Answers

Answers to the exercises are now available.