A hidden Markov model for reconstructing animal paths from solar geolocation loggers using templates for light intensity

The ability to track animal movements across long distances has revolutionized our
understanding of animal ecology and has been helpful to conservation 1], 2]. Until recently, our ability to record this information was limited to larger animals
that could carry satellite transmitters. However, recent technological advances have
developed miniaturized devices that extend our ability to track much smaller animals,
especially migratory songbirds 3]. Solar geolocation data loggers (or geolocators), are simple animal tracking devices
that record ambient light levels for the purpose of estimating the latitude and longitude
of an animal that wears the device. One of advantage of geolocators is that they typically
can record data for a year or longer, i.e. cover an annual migration cycle. Despite an ongoing miniaturization of GPS and other
satellite-linked tracking devices, geolocators remain useful because their low mass
(currently ca. 0.35 g) broadens the range of species that may be tagged 3]. Furthermore, current GPS tracking devices for small birds are limited to a relatively
small (e.g. 8–10) number of locations that can be recorded, thus it is difficult to use them
to generate information on departure and arrival dates. In addition, the simplicity
of geolocators design makes them inexpensive, which opens the prospect of affordable
population-wide studies of migration and migratory connectivity.

While the theory of estimating latitude and longitude from the elevation of the sun
is not new, the process of reconstructing animal movements by using light intensity
levels recorded with geolocators presents several challenges. From the perspective
of the hidden Markov modelling framework 4], sketched in Fig. 1, positioning by light level requires at least two parts: (1) an observational model
and (2) a process model. In the case of animal positioning we refer to these parts
as physical and movement models, respectively. The most common approach to solar geolocation
5] uses a simplified physical model. This model, which is referred to as the ‘threshold
method’, requires the definition of each twilight event in a dataset as the time point
corresponding to the moment when solar irradiance reaches some arbitrary, but constant,
threshold level (6], 7], see first column in Table 1). Latitude is then estimated by the duration of time between consequent pairs of
twilights and the longitude by the time of solar noon or midnight. This threshold
approach is still widely used, but it is plagued by many well-known problems such
as biased estimates 8], unrealistic assumptions of constant shading, and a null assumption of no movement
9]. Aside from its general simplicity and accessibility, the advantage of the threshold
method is that it needs only one data point per twilight period, which makes it well
suited to data storage by tiny and simple tags that log a very narrow band of light
intensities and have very limited data storage capacity.

Fig. 1. Sketch of solar geolocation principles and our method of analysis. A solar geolocator
(shown in black on bird’s back) records light levels and times (example raw data in lower panel). When the animal moves, its position is unobserved (hidden), but it can be estimated by the pattern of light changes measured at twilight. We
combine a physical (observation) model about how light levels change with position and time with some basic knowledge of
the patterns of movement between twilights (movement model) along with all previous and subsequent positions in a hidden Markov model framework. Then, using the particle filter, we arrive at the most likely position and movement for each twilight

Table 1. Review of the differences of the existing methods together with rationale for the
current contribution

Early attempts to improve upon the results of threshold analyses involved the application
of a specially developed state-space models on locations estimated with the threshold
method 10]. These models served as post hoc smoothers and generally improved the estimates, but they did not erase biases because
of the original biases in the threshold-based location estimates 11]. In response to the problems of the threshold method stemming from atmospheric properties
implicit in the observational model, Ekstrom (7], 12]; column 2 in Table 1) developed the template fit method. Based on the physics of the atmosphere he derived
the relationship between solar angle (angle above the horizon) and near-surface light
intensity. Ekstrom 12] showed that the template fit model was robust to the effects of shading (see equations
4–7 in the current contribution). Despite the great potential of the template fit
method, it was not used in the next generation of methods nor, to our knowledge, in
any published analyses of empirical data. The next generation of methods employed
specially developed state space models that not only had a better observational model,
but also incorporated assumption of animal movement 13], 14]. These models provide more accurate and precise position estimates and a measure
of their uncertainty (see columns 3 and 4 in Table 1).

All the early approaches to solar geolocation were developed for the tracking of marine
animals. This specialty called for relatively simple movement models that assumed
somewhat constant movement and also incorporated geographic masking under the assumption
that marine animals could only occur in the ocean 14] and are not able to move over the land masses. For many terrestrial organisms (e.g., migratory birds) more complex models are needed that can account for prolonged sedentary
behaviour interspersed with-long distance movements. Here we describe an approach
that attempts to fill the needs associated primarily with tracking small birds by
implementing the following features into a state space model: (1) accommodating the
narrow band light level data that are recorded by many tags used on birds; (2) allowing
for systematic changes in habitats throughout the annual cycle to make the observational
model maximally independent from shading regimes; (3) defining behavioural state parameters
that describe a bird as either sedentary or migrating; (4) modelling total distance
travelled between twilights as opposed to modelling average velocity; (5) masking
that allows movements (migratory state) over unsuitable habitats but not settling
(sedentary state).

The new analytic framework (implemented in FLightR R package 15]) provides the flexibility to incorporate all of these characteristics. We demonstrate
the approach by analysing both simulated and real tracks of small migratory birds.

Data and models

The example data and a sketch of our approach of inferring tag location are shown
at Fig. 1.

Data

Although the method we outline here is applicable to light level data collected by
a wide variety of geolocator tags, we focus our example on data collected by the widely
used Mk geolocators developed by British Antarctic Survey (BAS, Cambridge, United
Kingdom). Although, BAS has produced several models, they generally save a maximum
light level over time intervals ranging from 2 to 10 min (Fig. 1). These tags are optimized for recording low light levels within a narrow band of
light intensities around twilight, which results in 10 points collected during each
twilight. Thus, a typical dataset will be dominated by minimal values (=0) in response
to darkness and maximal values (=64) for measurements in daylight.

To demonstrate our method, we used both real data from BAS tags and simulations of
BAS tag data. The real data consist of two tracks. The first is a from a tree swallow
Tachycineta bicolor tagged at the breeding grounds in southern Canada and flying to wintering grounds
in Cuba by the North American East Coast. The second track is from a golden-crowned
sparrow Zonotrichia atricapilla tagged on its wintering grounds in California and flying to breed to Alaska via the
West Coast. The first track features a 2-min and the second employs a ten-minute light
logging interval. The simulated data are representative of stationary (not moving)
tags, and error associated with weather-related and behavioural shading was derived
from real data collected by a tag attached to the above-mentioned tree swallow while
it was stationary at a known location.

All existing analytical approaches, including the one presented here, focus primarily
on light-level transitions that occur around twilight, so these twilight periods have
to be extracted from the data in an automated 13] or semi-automated 5], 16] way. Light-level data collected with geolocators on birds are noisy because birds
frequently change their light environment by moving in and out of nest boxes, natural
cavities, or dense vegetation, so a completely automated method of detection of twilight
events produces many false positives. For this reason we consider a user-controlled,
semi-automated process for identifying twilight events in a dataset to be most effective.
We used an online interface called TAGS (Totally Awesome Geolocator Service; http://tags.animalmigration.org) to visually inspect data and generate a vector of twilight event periods, both morning
and evening, for each day of the year.

Hidden Markov model

Hidden Markov models are currently the most widely used framework for estimating animal
positions and behaviour through time 17], 18]. For solar geolocation the application of hidden Markov models is intuitive as animal
movements result in unobserved positions that must be estimated probabilistically
(Fig. 1). We developed a model with two parts: the hidden process model of animal movement
and the observational model of light measurements from the tag. For easy reference
we will refer to these models taken together as the FLightR model.

Physical observation model

The observational model matches recorded light levels to theoretical expectations
at different locations on Earth. During each twilight period i the tag records several measurements j. Estimation of the theoretical expectation consists of three steps. The first step
transforms the coordinates of potential twilight location ?_k
at time ? _ij
to angle of sun relative to horizon (solar angle, ? _ij
) with standard astronomical equations 19].

(1)

This step is shared by previous analytical approaches – the Trackit 13] and TripEstimation models 14], 16] and the current contribution.

At the next step, the expected light measurements (ELM) are calculated from solar angle. This step varies significantly across analytical
approaches, all of which try to account for natural variation in the relationship
between solar angle and measured light levels. Sources of this variation are scattering
(which includes cloud cover, shading of the tag by landforms and plant leaves; see
9], 20]) and natural variation in refraction 7].

Nielsen and Siebert 13] assumed the function f relating solar angle to light measurements was unknown and estimated it using cubic
splines with autocorrelation of spline parameters between consecutive twilights

(2)

The potential problem of this approach is that errors affecting light intensity are
multiplicative, not additive 12]. This means that if light intensity errors are significantly larger than zero the
result will be biased. Another drawback is that many light level measurements in each
twilight period are required to fit the splines well, and for the widely used BAS
geolocators, transitions from the maximum to minimum light-level values usually occur
in less than five time steps.

Sumner et al.14] accounted for the multiplicative error by taking the natural log of ELM and fitting
the following equation:

(3)

The function f?’ is nonlinear and is estimated for each tag with penalized splines on a calibration
dataset. The variable K _i
is attenuation (or cloudiness) and is also estimated. This approach is less dependent
on the shading bias but still requires many light level measurements during each twilight
event.

Ekstrom 12] derived the following deterministic relationship (i.e. ‘template’) between light at surface (solar irradiance at Earth surface, SL) and solar angle ?.

(4)

where erf is the error function and erfc is the complementary error function. This equation can be rewritten as follows:

(5)

Note that log (Cloudiness) plays the same role as the parameter K in equation 3.

In the FLightR observational model, we define

Thus

(6)

The function f?’’ is known, in contrast to the unknown function f?’ in equation 3.

Equations 4–6 deal with the relationship between the solar angle and light reaching
surface of Earth but not to the light measured by a tag. To do that, we need to account
for the properties of the tag. Assuming linearity (on a log scale) of the tag measurements
we define

(7)

Where tag intercept I and tag slope Z are tag specific and do not depend on position or time.

The complete equations for the ELM are then:

(8)

Here, the attenuation and noise terms as K’ _j
and ?’
_ij
apply directly to log (ELM _ij
) are not the same as the attenuation and noise terms in eqn. 6. The observed light measurement (OLM) depends on the attenuation K’ _i
and slope Z that both vary from day to day. The dependency of K’ _i
on the animal’s behaviour and habitat led us to focus only on Z for estimating the likelihood that the observed light measurements were produced
at location ?_k
. We assumed the observations come from a model with tag slope Z _i
that is a random variable independent of the animals’ behaviour. This allows us to
estimate the distribution of Z _i
from the calibration dataset. In all the analysed calibration datasets Z _i
had a lognormal distribution:

(9)

Parameters Z _calib
and are assumed to be constant and are estimated from the calibration at the known true sites, where model 8 was fitted.

The complete equations for the observed light measurement are then

(10)

In fact, for given location ?_k
and twilight i, eqn.10 is standard linear model of the form

(11)

with known b _ij
, unknown a _i
, unknown Z _i
, and random error c _j
assumed to be normally distributed with a mean of zero. Thus, for any location ?_k
we can apply a standard least squares procedure to estimate the mean ? _i
and standard deviation of the unknown Z _i
and then use the following probability density as a surrogate for the hypothesis that
the data were observed at location ?_k
:

(12)

Thus we can integrate product of density from eqn. 9 and eqn. 12 over all Z and estimate the required probability density and the likelihood of data at ?_k
.

Parameters Z _calib
and ? _calib
must be estimated for every tag from twilights recorded at a known position. This
means that the tag has to be calibrated in a known position for at least 5 days before
the animal leaves the area or after it arrives back. The calibration does not need
to be a clear-sky calibration and preferably should be done on the bird. The casing
around the light sensor may discolour as the tag ages, such that the calibration generated
when the tag is deployed will not match the calibration when the tag is recovered.
If known locations are available for the beginning and end of the tag use, it is possible
to account for ageing by assuming a linear change in the calibrated slope throughout
the deployment period.

Movement model

Previous movement models for solar geolocation were developed for marine animals that
are assumed to travel at a relatively constant velocity. For birds, we know that there
are long periods (up to 6 months) when birds are primarily sedentary (moving less
than several km a day), punctuated by short periods (1–2 weeks) when birds migrate
rapidly between breeding and non-breeding locations, but during these periods may
spend several days at a single site to refuel before migrating again 21]. For these highly variable movement patterns we needed a more flexible movement model.
We used a simplified “Double” model 22], 23] with just two behavioural states: “Sedentary” and “Migrating”. In the context of
these two states, we defined the position at the twilight ? _i+1
as:

(13)

(14)

Where an offset d _i
follows an uncorrelated random walk with a distribution reflecting bird behaviour
and described by a set of parameters:

(15)

At each moment i, d _i
is distributed as a mixture of a zero increment (no movement) with probability p _i
, and non-zero increment in direction ? _i
and step length S _i
. We assigned distributions ? _i
and S _i
as follows:

(16)

(17)

For all runs in our models, we assumed that the truncation points were a?=?45 and b?=?1000 km. That is, when birds initiated a movement, they could not fly less than
45 km or more than 1500 km in a single between-twilight interval. Parameters ? and
k in the von Mises distribution reflect direction of migration and its concentration,
and the parameters of the truncated normal distribution for distance shape the distribution
of inter-twilight flight distances. R packages ‘circular’ 24] and ‘CircStats’ 25] were used to obtain random draws from von Mises distributions given the parameters
and ‘truncnorm’ 26] for the truncated normal distribution. All movement was modelled as being uncorrelated
in time, assuming that distance flown, direction and behavioural switches during bird
migration may be highly independent one day to the next.

Finally, to allow for migration across water bodies, we introduced a spatial behavioural
mask that prevents birds from entering a sedentary state in locations corresponding
with large water bodies. Hence, birds may fly over water, but cannot switch to sedentary
mode in this habitat.

Bayesian analysis

Many techniques may be used to estimate a hidden Markov model. Our uncorrelated random
walk model has five unobserved variables at each of many time steps. Bayesian methods
are particularly useful for computing posterior distribution over these variables.
The Kalman filter 27] or unscented Kalman filter 28] are other approaches that can be applied to the problem of inference in a long state
space model 13], but these do not apply to highly nonlinear posterior distribution patterns caused
by non-Gaussian light error distribution with spatially explicit masks. To accommodate
these aspects of geolocator data analysis we decided to (1) discretize space and (2)
use the particle filter 29] for approximating the posterior distribution. Spatial discretization with creation
of a regular grid accelerates computational workflow by minimizing amount of possible
between node transitions and also puts no constrains on the spatial error distribution
18], 30], 31]. We used the function regularCoordinates() from the geosphere R package 32] to discretize a state-space into a regular grid with a distance of 50 km between
nodes. Choice of the distance of 50 km between grid nodes was arbitrary, and followed
the idea that it should be small enough to have posterior probabilities distributed
at several grid cells at any time, but at the same time not too small as it will make
grid larger and estimation slower. All the movements were estimated between these
nodes and the observational model was estimated at these nodes. The approach allows
us to “obtain a numerical non-parametric representation of the probability distribution
of the animal’s position. This probability distribution illustrates the uncertainty
of the estimated movement with a high degree of detail. Finally, we can draw inference
about parameters in a likelihood framework, compute the most probable track of the
animal, and sample random tracks that the animal may have traveled” 31].

Following Patterson et al.17], we used a particle filter (PF) or ‘sequential importance resampling’. Our implementation
of the PF was based on the algorithm from Doucet et al.33] as expanded by Andersen et al.34]. The logic of the PF is as follows. At the initialisation phase, our PF creates a
sample of 10
⁶
particles (points with positions) at the actual release point. Then, for each particle,
it generates a new position at i?+?1 from the process model. All new particle locations are then resampled proportionally
to their weight (product of the previous likelihoods and the current likelihood, estimated
by the observational model and a priori constraints as explained below). The resampled set of particles proceeds to the next
step and the process repeats. In other words, PF creates 10
⁶
possible paths that develop from the release point according to the rules for the
movement model. At every twilight all 10
⁶
particles are compared to data passed through the observation model. After each check
all unlikely particles are replaced by likely ones with the probability of their relative
likelihood. Once all iterations are completed, the histories of the remaining particles
render a distribution that approximates the spatial probability distribution for the
bird at each twilight.

Because of the degeneracy problem common for long PF runs (see e.g.35]), we used block sampling 36] in our PF. assuming that at n steps before current step i the particles have reached their global optimum. We selected n?=?90 and then estimated the current weight of the particle as the mean of the logarithms
of particle weights at steps [i-n]:n.

To include information on where the bird was recaptured when the geolocator was removed,
we used the approach described by Andersen et al.34]. Here the last n states of particle histories were resampled with weights proportional to their probability
density function from a normal distribution with a mean of the recapture point coordinates
and standard deviation (SD) equal to some measure of precision (because we had 50 km
between grid points, we used a SD of 25 km).

A priori constraints in the model can be defined according to any biological assumptions.
These constraints are used at the resampling stage as weights, and they are not optimised.
Here, we present the model with two constraints: a general spatial mask and a behavioural
spatial mask. The general spatial mask works through exclusion of selected grid nodes
and is binary. Excluding land, for example, may be useful for modelling fish, marine
mammal, or pelagic bird movements. The behavioural mask is described above and lowers
probability for animals to adopt a stationary state in particular nodes. Value of
zero for this mask will completely prevent sedentary mode at the node.

The PF offers a relatively fast optimization technique, but it is unstable to outliers
37]. Outliers in twilight data can happen if an animal stays inside a cavity or nest
box during twilight and emerges later, resulting in what would appear to be a late
and atypically fast twilight. Ideally such points should be removed during visual
inspection of the twilights, but in order to make the PF stable to occasional undetected
outliers we have developed two outlier detection approaches. The first approach is
used after tag calibration but before the particle filter run. The main idea of the
approach is to estimate likely longitude of crossing with the equator for each twilight
and then use time series outlier detection software to identify outliers 38], 39]. The other approach is an ‘on the go’ outlier detection technique. This filtering
step is implemented at twilight i when new particle positions were generated but before the resampling for the next
step occurs. At this point, the average distance from i???1???i is compared with i???i?+?1, and if the latter is smaller and mean turning angle at i is less the 100°, then twilight i is considered to be outlier and its likelihood surface is not used in the particle
resampling.

Application

Example 1: simulated stationary tags

We have simulated year-long light data for stationary tags using eqns 9 and 10. The
following parameters for the simulations were estimated from real data collected by
a tag attached to a Tree Swallow while it had been in a known location for ten days,
so they contain a realistic distribution of shading and of slopes.

We began by estimating locations with the threshold method implemented in GeoLight.
For this analysis, we used the month of July as the calibration period. Then, using
the same calibration period, we estimated the locations in FLightR. We did not use
any behaviour or spatial masks and confined our analysis to a reasonably large area
with radius of 1000 km around the tag location for the potential spatial extent required
for the FLightR runs. For all the simulation and real tag runs we used the same priors
of 0.1 for probability of migratory behaviour and 300?±?150 km for distance covered
between twilights. We consider these settings to be suitable for most of the animals
except ones that move longer distances between twilights. In our experience with different
tags, distance priors do not affect results of the model run, if they are generally
correct and wide enough. The probability of migration prior does have an effect, and
if selected too high, the resulting track can become noisy. A value of 0.1 was selected
as a prior on the basis of good results with simulated data and real tracks; we do
not recommend changing it without a specific reason.

To compare the performances of GeoLight and FLightR, we estimated monthly biases (known
minus estimated location) and SD associated with each analysis type for both latitude
and longitude (Table 2). Monthly estimates of latitude were not biased when estimated by FLightR (bias ranged
from ?0.1° to 0.03° due to some rounding error from estimation on a grid) whereas
GeoLight biases in latitude estimation ranged from ?2.41° to 3.63°. Moreover the FLightR
errors on latitude are almost uniform across the year, with SD ranging from 0.23°
to 0.33°, whereas the Geolight SD of monthly errors ranged from 1.02° to 10.96° (Fig. 2). Note that the GeoLight method does not provide estimates for latitude at the time
around equinox without assumption of stationarity. The results in Table 2 show that estimates from the threshold method are biased while FLightR produces unbiased
and consistent estimates.

Table 2. Average bias and SD estimated by GeoLight and FLightR for simulated stationary tags
at 5° N and 55° N. GeoLight calibration was done by July twilights

Fig. 2. Estimated position of simulated stationary tag by classic approach (GeoLight, points)
and FLightR (median estimated latitude with 50 and 95 % credible intervals)

Example 2: tree swallow – 2 min fixes

To demonstrate the application of the FLightR functions to a real-world geolocator
dataset, we analysed the track from a single tree swallow that was fitted with an
Mk12S geolocation logger on 23 June 2011 at the Long Point Bird Observatory in Ontario,
Canada (80.46° W 42.62° N). It was recaptured at the same location and the tag was
removed on 7 June 2012. Geolocator was developed by the British Antarctic Survey (Cambridge,
United Kingdom) and had 15-mm-long stalk positioned at a 30° angle. The tag recorded
maximum light levels at 2 min intervals.

This single track is a small part of a much larger study with many individuals tracked
with geolocators at over ten sites across North America (Bradley et al. in prep.). We provide raw geolocator data and code to estimate positions for both
of the packages online at https://github.com/eldarrak/FLightR/blob/master/examples/tree_swallow_BAS_tag_example/tree_swallow_analysis.Rmd. Note that during the GeoLight run no filters or outlier detection tools were used,
and FLightR detected as outliers and excluded data from ~30 twilights. We present
four different ways of looking at the path to better reconstruct its details and implications.
The reconstructed path from FLightR (coloured symbols in Fig. 3) has the bird staying in the region of Long Point through early July. The bird appears
to have stayed at the same longitude and latitude (Fig. 4), and this is backed up with the highest confidence in the reconstructed path (Fig. 5). On the morning of 13 July, the bird departed Long Point (a diurnal departure indicated
by the orange dot in Fig. 6a) and made a single flight of over 500 km (Fig. 6b, line segment A in Fig. 3) to a stopover site in Virginia (yellow dots in Fig. 3). It stayed in the same general area until 1 or 2 August, when it departed (segment
B Fig. 3) and flew over 200 km to a site near the coastal border of North Carolina and Virginia
(Fig. 3). This movement is less certain in its particulars, having interquartile ranges of
1–2 degrees (Fig. 5). The interquartile ranges of position estimates for the following weeks are near
the equinox (and thus very unreliable for latitude measures) and vary from low (1
degree) to very high (5 degrees; Fig. 5). The bird departed the region on 23 or 24 October, heading south and inland (segment
C), and very likely not stopping until reaching coastal South Carolina, where it remained
until about 2 November. It then departed for Cuba (segment D), spending most of the
rest of the winter there. The one exception to this Cuba winter residency may have
been a departure (segment E) to the vicinity of the Bimini Islands on or about 28
February, returning to Cuba about 5 March, but uncertainties at this time of year,
nearing the equinox, are very high (Fig. 5). Much more certain is that the winter residency ended on 28 March, when the bird
flew from Florida to the coast of North Carolina/Virginia (segment F). It did not
stay there for long, if at all, and continued migrating until it reached Long Point
on 6 April. The interpretation of the locations during this last northward leg are
complicated by the large latitudinal uncertainties associated with the vernal equinox
(Fig. 5). The customary approach using thresholds with GeoLight generally agreed with the
FLightR estimates.

Fig. 3. The track of a tree swallow as estimated from 2-min fixing interval data by classic
approach (GeoLight, grey line and dots) and FLightR. Inset shows Tree Swallow range
in North America. The medians of twilight positions estimated by FLightR are coloured
by the month of a year (colours for each month are illustrated with pie chart). In
June a bird was tagged on the breeding grounds at the Long Point Bird Observatory
in Ontario, Canada. July 19 it left breeding ground and moved to the stopover site
in Virginia (segment A). At 1or 2 August it moved to coastal area on the border of
North Carolina and Virginia (segment B), where remained stationary until end of October.
23–24 October the bird departed towards a stopover site in South Carolina (segment
C), where it stayed for a week and then continued south to Cuba (segment D). The bird
remained in Cuba till the end of March, with one exception for flight to the vicinity
of the Bimini Islands (segment E). On 28 of March bird left wintering grounds and
migrated to the North Carolina/Virginia site (segment F) and after a short stopover
there moved straight to the breeding grounds (segment G). Note that no spatial or
behavioural masks were used for the FLightR run, so positions were allowed to be everywhere.
Raw GeoLight estimates shown on the figure should not be interpreted as a positions
and inference can be made only for most likely location during stationary periods,
new functions in GeoLight 2.0.1 are available now for the estimation of these most
likely stationary locations (Lisovski Liechti, pers. comm). Tree swallow range image
is courtesy of Birds of North America: Cornell Lab of Ornithology

Fig. 4. Longitudes (upper panel) and latitudes (lower panel) of a track of a tree swallow
as estimated from 2-min fixing interval data by classic approach (GeoLight, dots)
and FLightR. The medians of twilight positions estimated by FLightR are shown with
accompanying quartile ranges and 95 % credible intervals. Note absence of the latitudinal
positions from GeoLight during the equinoxes (shown by red vertical lines)

Fig. 5. Precision of the estimates of positions by longitude (orange) and latitude (blue)
of a track of a tree swallow shown by the interquartile ranges

Fig. 6. Medians of probability of migration (upper panel) and migration distance with corresponding
quartiles (lower panel) for a tree swallow estimated by FLightR. Nocturnal migration
is shown in grey and diurnal migration in orange. The circles show transitions which
were characterized by shift of the median latitude and/or median longitude. Vertical
red lines mark equinoxes

Example 3: golden-crowned sparrow – 10 min fixes

Our sample dataset for golden-crowned sparrow was collected from a bird that was tagged
on its wintering grounds in coastal California and tracked to the breeding grounds
on the coast of Alaska 40]. The bird was captured and tagged on 2 February 2010 at the Palomarin Field Station
in coastal California, United States (37.93° N, 122.74° W). It was recaptured at the
same location and the tag was removed on 19 October 2010. The tag was an Mk10S geolocators
developed by the British Antarctic Survey (Cambridge, United Kingdom) with a 15-mm-long
stalk positioned at a 30° angle. The tag recorded maximum light levels at 10 min intervals.
Note that during the GeoLight runs no filters or outliers detection tools were used,
and FLightR detected as outliers and excluded about 40 twilights out of 600.

The reconstructed track for the golden-crowned sparrow indicates a movement from the
wintering grounds in California, north along the west coast of North America to the
breeding grounds in Alaska (Fig. 7). After it was tagged in early February, the bird was sedentary until 18 April, and
clear migration activity was initiated May 13 and ceased in early June (Fig. 8). This is consistent with post-tagging observations of the bird that confirmed that
it was present (and likely sedentary) at the tagging sight until at least March 31.
During the migration both latitude and longitude have high uncertainty (Fig. 9), during the breeding season at Alaska uncertainty in longitude remains high. During
spring migration, daily movements were primarily between 500 and 700 km, with some
as long as 800 km (Fig. 10). Fall migration occurred between 13 September and 14 October (Fig. 10). During both spring and fall migrations, movements were both diurnal and nocturnal.

Fig. 7. The track of a golden-crowned sparrow as estimated from 10-min fixing interval data
by the classic approach (GeoLight, grey line and dots) and by FLightR. Inset shows
golden-crowned sparrow range in North America. The medians of twilight positions estimated
by FLightR are coloured by the month of a year (colours for each month are illustrated
with pie chart). Bird was tagged in 2 February on the wintering grounds in California.
It remained close to the capture site till 13 May and then started northward migration.
After two weeks of migration it arrived on the breeding grounds (~1 June), and remained
there until at least 5 July. It may have moved about 200 km westward after breeding,
though the uncertainty here is high. 10 September bird started migration to the wintering
grounds, where it arrived at in the end of October. Note that no spatial or behavioural
masks were used for the FLightR run, so positions were allowed to be everywhere. Raw
GeoLight estimates shown on the figure should not be interpreted as a positions and
inference can be made only for most likely location during stationary periods, new
functions in GeoLight 2.0.1 are available now for the estimation of these most likely
stationary locations (Lisovski Liechti, pers. comm). Golden-crowned sparrow range
image is courtesy of Birds of North America: Cornell Lab of Ornithology

Fig. 8. Longitudes (upper panel) and latitudes (lower panel) of a track of a golden-crowned
sparrow as estimated from 10-min fixing interval data by classic approach (GeoLight,
dots) and FLightR. The medians of twilight positions estimated by FLightR are shown
with accompanying quartile ranges and 95 % credible intervals. Note absence of the
latitudinal positions from GeoLight during the equinoxes (shown by red vertical lines)

Fig. 9. Precision of the estimates of positions by longitude (orange) and latitude (blue)
of a track of golden-crowned sparrow shown by the interquartile ranges

Fig. 10. Medians of probability of migration (upper panel) and migration distance with corresponding
quartiles (lower panel) for a golden-crowned sparrow) estimated by FLightR. Nocturnal
migration is shown in grey and diurnal migration in orange. The circles show transitions
which were characterized by shift of the median latitude and/or median longitude.
Vertical red lines mark equinoxes

Uncertainty in both longitude and latitude increased during both migrations. Notably,
uncertainty in longitude also increased around summer solstice in the end of June
(Fig. 9). This may be because both twilight lines become almost horizontal and parallel to
each other in summer at high latitudes, thus limiting inference on the longitudinal
position. The extreme case of this happens at the Arctic Circle (66° N) at solstice.
This situation reflects one of the two major limitations in solar geolocation – (1)
determination of latitude near equinoxes at any longitude and (2) determination of
longitude near solstice at high latitudes. The FLightR model was able to cope with
these limitations and provided meaningful positioning results.