Light-Level geolocator analysis using the GeoLight package
Michael T. Hallworth
Smithsonian Conservation Biology Institute
Migratory Bird Center
Introduction to Archival light-level geolocators
The following document outlines the steps for analyzing data from archival light-level geolocators (hereafter geolocators). Geolocators have been used to track individuals since the early 1990s but were restricted to large organisms because of their large size. Recently, with the miniturization of geolocators, researchers are now able to deploy geolocators on smaller and smaller species. Geolocators are devices that record ambient light levels every 2, 5, or 10 min depending on the model. Geolocators are attached to individuals which then migrate with the device while it records ambient light-levels throughout the year. Once recovered, the data are downloaded and analyzed to determine the time of sunrise and sunset. From the sunrise/sunset data the time of relative noon and midnight are determine. Geographical cooridnates are then derived from the relative ‘noon’ and ‘midnight’ times to give an approximate location of where the individual was throughout the year.
This tutorial uses geolocator data from a male Ovenbird (Seiurus aurocapilla) breeding at Hubbard Brook Experimental Forest, NH and a Wood Thrush (Hylocichla mustelina) breeding in Indiana and is part of an ongoing study modeling regional source-sink dynamics of a migratory songbird. Click here for more information regarding the Wood Thrush project. The geolocators used for the two projects were purchased from British Antarctic Survey (BAS) and Lotek (LightBug). These models operate in the same manner but format the data differently. In order to use GeoLight to analyze the data, the data needed to be formatted correctly so GeoLight can read the data.
Analysis of geolocator data using the GeoLight package in R
The following R packages are needed to conduct the tutorial
library(GeoLight)
library(raster)
library(ks)The data format of geolocators purchased from different vendors are slightly different and need to be converted into a file format that is recognized by GeoLight
OvenRaw<-read.csv("data/2391_68559_000.lig",sep=",",skip=1,header=FALSE,col.names=c("Valid","Date","Julian","Light"))
WothRaw<-read.table("data/0309_RAW.txt",sep="\t",header=TRUE)Ovenbird Raw data from British Antarctic Survey (.lig file)
## Valid Date Julian Light
## 1 ok 10/06/11 22:37:30 40704.94 64
## 2 ok 10/06/11 22:39:30 40704.94 64
## 3 ok 10/06/11 22:41:30 40704.95 64
## 4 ok 10/06/11 22:43:30 40704.95 64
## 5 ok 10/06/11 22:45:30 40704.95 64
## 6 ok 10/06/11 22:47:30 40704.95 64Wood Thrush Raw data from Lotek LightBug (.txt file)
## time date light
## 1 00:07 27/07/11 236
## 2 00:14 27/07/11 236
## 3 00:21 27/07/11 236
## 4 00:28 27/07/11 236
## 5 00:35 27/07/11 236
## 6 00:43 27/07/11 234The following functions were created to convert the different raw data formats into the format needed to process light data in GeoLight. The functions need to be added to your R console to use the function. The source function is used to read in a function stored as an R file in your working directory.
source("source/sourceRead_lig.R") # function written by Simon Wotherspoon accessed from GitHUB
source("source/sourceRead_Lightbug.R") # read_lig function modified by M.T.Hallworth to read in LightBug dataImport .lig file using read.lig fucntion to convert dates that GeoLight recognizes
GL_68559<-read.lig("data/2391_68559_000.lig")
head(GL_68559)## Valid Date Julian Light
## 1 ok 2011-06-10 22:37:30 40704.94 64
## 2 ok 2011-06-10 22:39:30 40704.94 64
## 3 ok 2011-06-10 22:41:30 40704.95 64
## 4 ok 2011-06-10 22:43:30 40704.95 64
## 5 ok 2011-06-10 22:45:30 40704.95 64
## 6 ok 2011-06-10 22:47:30 40704.95 64Now that the data are formatted - you can use GeoLight to determine transitions (sunrise/sunset)
In this example - a threshold of 5 was used - a larger or smaller value can be used but it will decrease/increase the number of transitions that need to be scored.
LightThreshold - determines light levels over 5 as “sun has risen/set” and asks you to accept/reject them.
Note - determining the transitions in this file took approx. 45mins - 1hr
The following code produces an interactive plot which asks the user to either accept or reject each light transition that passes over the threshold specified in the code. This document does not support interactive plots but the plot you will see should look similar to the one below. The blue line identifies the threshold level set in the code.
GL_68559_transitions<-twilightCalc(datetime = GL_68559[,2],
light= GL_68559[,4],
LightThreshold=1, # Here is where you set the threshold level
ask=FALSE) # Here you can go through every twilight
Exercise: Go through and assign a few twilights.
* How did you decide to either accept or reject twilight events?
* How do you think accepting / rejecting twilight events would influence the end result? * Change the threshold to a higher value - how does that change the number of shading events?
Once you have gone through the process of accepting or rejecting the transition events the data will look like this. tFirst and tSecond correspond to the time of the transitions and type illustrates whether the location will be dervied from relative ‘noon’ or ‘midnight’ locations.
head(GL_68559_transitions)## tFirst tSecond type
## 1 2011-06-11 00:19:30 2011-06-11 08:51:30 2
## 2 2011-06-11 08:51:30 2011-06-11 23:59:30 1
## 3 2011-06-11 23:59:30 2011-06-12 08:53:30 2
## 4 2011-06-12 08:53:30 2011-06-13 00:27:30 1
## 5 2011-06-13 00:27:30 2011-06-13 08:45:43 2
## 6 2011-06-13 08:45:43 2011-06-14 00:13:18 1Sun-elevation angle
The next step is to calculate the sun-elevation angle of a known capture location. The sun-elevation angle is the angle of the sun with respect to the horizon at the time the geolocator light data passed the threshold set by the user. Thus, the sun-elevation angle is unique to the threshold used in the analysis.
Here I chose the dates between deployment of the geolocator and July 31 to ensure that only transitions when the bird was at the capture location were used to calculate the sun-elevation angle.
The coordinates also need to be entered - (X,Y) in that order
This Ovenbird was captured at Hubbard Brook Experimental Forest, NH (-71.45,43.945)
SunElev<-getElevation(tFirst= GL_68559_transitions[1:102,1],
tSecond= GL_68559_transitions[1:102,2],
type=GL_68559_transitions[1:102,3],
known.coord=c(-71.45,43.945),
plot=TRUE)
SunElev## [1] -2.55Location estimates assuming no change in sun elevation angle throughout the year
GL_68559Locations<-coord(tFirst= GL_68559_transitions[,1],
tSecond= GL_68559_transitions[,2],
type=GL_68559_transitions[,3],
degElevation=SunElev)## Note: Out of 768 twilight pairs, the calculation of 52 latitudes failed (6 %)head(GL_68559Locations)## [,1] [,2]
## [1,] -68.86315154 41.33162060
## [2,] -66.33819736 38.33935842
## [3,] -66.56264864 38.00401399
## [4,] -70.03643597 41.99894846
## [5,] -69.01068030 43.02228910
## [6,] -67.23486139 41.10153675Plot the location data
Note - We accepted all twilights with the original .lig file. That file includes all transtions the geolocator recorded. Therefore, it may include transitions before it was attached to the bird, or after it was taken off.
plot(GL_68559Locations,
pch="*",
col="red",
xlab="Longitude",
ylab="Latitude")
plot(world,add=TRUE)
Exercise:
How would you subset the data to not include the extreme outliers? i.e. points in the old world.
How could you check to see if those points were either at the beginning of the file, the end or both?
Instead of removing the points - could you set the plot region to focus in on where the bird was?
hint take a look at what the function drawExtent() does in the raster package.
Create Kernel Density Estimates (KDE) around the stationary periods
The following packages were used to create KDE of stationary periods
library(ks)
library(raster)
library(RColorBrewer)Breeding locations were determined using location data described earlier in determining the sun-elevation angle.The non-breeding period determined as 1 November - 3 March (the start of spring Equinox period). See Hallworth et al. 2015 for details.
breeding_68559<-data.frame(GL_68559Locations[1:102,1],GL_68559Locations[1:102,2])
NB_68559<-data.frame(GL_68559Locations[288:534,1],GL_68559Locations[288:534,2])Determine bandwidth for the Kernel density estimate - the bandwidth parameter sets the ‘smoothness’ of the KDE. The bandwidth was estimated using least-square cross validation.
Bwidth<-Hlscv(breeding_68559)
NBwidth<-Hlscv(NB_68559)The following script creates the KDE and converts the KDE to a raster.
Breeding_KDE<-raster(kde(x=breeding_68559,H=Bwidth))
NonBreeding_KDE<-raster(kde(x=NB_68559,H=NBwidth))Plot the results - not elegant but you can dress it up from here anyway you want
Exercise:
Plot the mean and median location for the non-breeding period (Nov 1 - March 31).
The locations are in the object GL_68559Locations which has two columns without headers.
The first column (GL_68559Locations[,1]) is the longitude column.
The corresponding rows in the data for Nov 1 and March 31 are rows 287 & 590.
- Do your conclusions about where the bird wintered change based on how you plot the data (mean, median, kernel density estimate)?
- How might your conclusions change if the bird wintered in South America and not in the Caribbean. Would the way you present the data change your conclusions about where the bird wintered?
Using different sun-elevation angles for different periods of the annual cycle
Notice in the above figure how the non-breeding KDE is exclusively over the Caribbean and does not fall over land. The sun-elevation angle (described above) can make a big difference in the latitude of the locations. The sun-elevation angle can be influenced by a variatey of factors such as habitat type, topography, weather, and bird behavior. Thus, using multiple sun-elevation angles for different portions of the year may be justified. The following code demonstrates how the sun-elevation angle for different portions of the year can be determined. Few studies have in-habitat calibrations for both breeding and non-breeding periods of the year (but see Hallworth et al. 2015, Stanley et al. 2014, McKinnon et al. 2013). Thus the sun-elevation angle for different portions of the year need to be estimated. In GeoLight there is a function called HillEkstromCalib which determines the sun-elevation angle based on the transition events in the data.
In the following example the different periods of the year are determined manually based on the natural history of Ovenbirds.
Create empty column to later fill with numeric values for period -
GL_68559_transitions$period<-rep(NA,dim(GL_68559_transitions)[1])
head(GL_68559_transitions)## tFirst tSecond type period
## 1 2011-06-11 00:19:30 2011-06-11 08:51:30 2 NA
## 2 2011-06-11 08:51:30 2011-06-11 23:59:30 1 NA
## 3 2011-06-11 23:59:30 2011-06-12 08:53:30 2 NA
## 4 2011-06-12 08:53:30 2011-06-13 00:27:30 1 NA
## 5 2011-06-13 00:27:30 2011-06-13 08:45:43 2 NA
## 6 2011-06-13 08:45:43 2011-06-14 00:13:18 1 NADefine Breeding based on natural history of species.
Ideally this would be done based on the biology (breeding season) and so on.
Here we will define the breeding period as from the time of capture to July 31.
We will define the non-breeding period as the dates between 1 November - 31 March.
see Hallworth et al. 2015 for details
Assign a stationary period to the transition data
0= Migratory Period, 1=Breeding, 2=Non-breeding
names(GL_68559_transitions)## [1] "tFirst" "tSecond" "type" "period"GL_68559_transitions[c(1:102,673:768),4]<-1
GL_68559_transitions[103:287,4]<-0
GL_68559_transitions[288:534,4]<-2
GL_68559_transitions[535:672,4]<-0Generate sun-elevation angles for different periods of the year using Ekstrom-Hill calibration. The Ekstrom-Hill calibrartion generates a sun-elevation angle for all non-zero periods.
HillEkstromCalib(tFirst= GL_68559_transitions[,1],
tSecond= GL_68559_transitions[,2],
type=GL_68559_transitions[,3],
site=GL_68559_transitions[,4],
start.angle=c(0,-3), #angle to start process c(breeding,nonbreeding)
plot=TRUE)
## [1] NA -2.1Define a column to store the sun elevation angles that were just generated.
GL_68559_transitions$sun<-rep(NA,dim(GL_68559_transitions)[1])Fill the new column using the sun-elevation angle for the known location (breeding = -2.55) and Ekstrom-Hill calibration for unknown locations (non-breeding = -2.1)
names(GL_68559_transitions)## [1] "tFirst" "tSecond" "type" "period" "sun"GL_68559_transitions[c(1:102,673:768),5]<-(-2.55)
GL_68559_transitions[103:287,5]<-(0) # Migration values are zero in this example
GL_68559_transitions[288:534,5]<-(-2.1)
GL_68559_transitions[535:672,5]<-(0) # Migration values are zero in this exampleNote - One could treat migration points differently - for example, the midpoint between breeding and non-breeding sun elevations could be used as a transition period (this would need to be justified in your manuscript) or you could use either the breeding or non-breeding sun angle but also would need to be justified.
Generate new location data based on newly derived sun-elevation angles
GL_68559LocationsEHC<-coord(tFirst= GL_68559_transitions[,1],
tSecond= GL_68559_transitions[,2],
type=GL_68559_transitions[,3],
degElevation=GL_68559_transitions$sun,
sites=GL_68559_transitions[,4])## Note: Out of 768 twilight pairs, the calculation of 0 latitudes failed (0 %)Plot the new location data
Left = Single sun-elevation angle (-2.55)
Middle = Non-breeding sun-elevation angle determined using Esktrom-Hill calibration (-2.1)
Right = Non-breeding sun-elevation angle determined using in-habitat calibration (-3.4) (not shown) Hallworth et al. 2015
Excercise:
How could using a single sun-elevation angle throughout the year influence your interpretation of the data?
Change the sun-elevation angles to see how the locations change.
* what happens to the non-breeding locations? * what happens to the breeding locations? * What factors may influence the sun-elevation angle?
* Does it make sense to use a single or multiple sun-elevation angles throughout the year?
* How could you design a project (locations to deploy geolocators) to get a better idea of how sun-elevation angles change throughout the year?
Advanced Exercise:
- Here we ran through the entire analysis for the Ovenbird but we have a Wood Thrush example as well.
- Run through the entire analysis to determine where the Wood Thrush spent the non-breeding season.
