GEOG 566

         Advanced spatial statistics and GIScience

Archive for Final Project

June 11, 2018

Analysis of flowering phenology as a step towards understanding phenological synchrony

Filed under: 2017,Final Project @ 5:55 pm

BACKGROUND: While the cinnabar moth (Tyria jacobaeae) was intentionally introduced to North America as a biological control insect, its release had unintended consequences as it established on a non-target native perennial herb, Senecio triangularis (Diehl and McEvoy, 1988). Populations of triangularis colonized by the cinnabar moth can experience heavy foliar herbivory (up to 100%) without showing reductions in plant size or seed set (Rodman, 2017). Herbivory of flowers is less common given a mismatch between the flowering time of S. triangularis and peak feeding stages for the moth, but when it does occur, floral/seed herbivory may have direct impacts on population dynamics. Previous work has shown that that larvae experiencing phenological synchrony with S. triangularis flowers decreased seed set by 95% (Rodman, 2017); and that S. triangularis is seed-limited, so that reduction of seed set decreases seedling recruitment to the next generation (Lunde, unpublished data).

Phenological synchrony is defined in this study as the overlap of “virulent” larval stages (late stages: L4 and L5) and “vulnerable” flower stages (early stages: primordia (P), buds (B), young flowers (YF) & flowers (F)). Synchrony varies by site and year, and the environmental drivers causing this variation are unknown. Knowing which populations of S. triangularis would be most likely to experience seed herbivory by cinnabar larvae could help managers track and respond to cinnabar moth presence and the risk to S. triangularis on a site-by-site basis. Because phenological synchrony was a rare occurrence for the five sites represented in 2017 data, most of the analyses described here focus on characterizing the timing of flowering of S. triangularis, while later work will focus on the timing of the L4 and L5 stage “virulent” cinnabar larvae.


RESEARCH QUESTION: How does the timing of “vulnerable” flower stages differ by site? Can I describe and quantify between site-differences in the timing of these stages? And how well do differences in candidate environmental factors account for differences in the occurrence of vulnerable flowers?


DATASET: This project used data collected from July to September, 2017, at five sites in a region of the Willamette National Forest near Oakridge, Oregon. A sixth site, Bristow, was lost due to fire, but is occasionally present in data analyses. Phenology scores (for cinnabar larvae and Senecio triangularis) and environmental variables were taken at each site during surveys at regular 10-day intervals. Some data were taken at the individual plant level, while some variables (such as snowmelt date) could not be measured in a meaningful manner below the site level. Data details:

  1. Population locations: a polygon layer outlining the 5 meadow sites surveyed over the 2017 season.
  2. Phenology: counts of vulnerable and invulnerable flower heads (capitula), and counts of cinnabar larvae in peak feeding stages (4th & 5th instars), collected from random subsets of tagged plants during each survey. Binary response variables created from phenology scores included plant vulnerability (1=vulnerable, 0=invulnerable), larval virulence (1=virulent larvae present, 0=absent), and phenological synchrony (1=co-occurrence of virulent larvae & vulnerable plants, 0 = one or none occurred).
  3. The following environmental variables are included in the phenology dataset: ambient temperature at 36 inches from the ground; soil temperature and soil moisture at 6 inch depth. Soil moisture and temperature readings were taken at each individual plant.
  4. Growing degree days for each site and survey day were modeled using a single triangle method with developmental thresholds set as above 5˚ and below 37˚ C. Growing degree days are defined as accumulated daily heat gain above and below developmental thresholds, here set to 5˚ and below 37˚ C. Developmental thresholds were estimated from previous studies on phenology of alpine flowers (Kudo and Suzuki, 1999); whereas I recognize that 10˚ C may be more commonly used as the minimum temperature for development of insects.


  1. HYPOTHESES: Alpine and subalpine perennial plant species have been shown to vary widely in terms of which environmental factors drive flowering phenology (Dunne et al., 2003). I expect to see differences in flowering phenology between sites. I expected differences between sites to correlate with differences in the abiotic environment: growing degree days, soil temperature, and soil moisture. And in order to test this, I also needed to explore how these three environmental variables were correlated to one another, so that I could understand issues underlying the analysis of these independent explanatory variables.

Because insects are ectotherms whose development is largely constrained by their ability to gain heat externally, the phenology of insects is often predicted on growing degree days (Johnson et al., 2007); and I would also expect to see variation between sites in the timing of occurrence for virulent larvae, perhaps driven more strongly by degree days, and less strongly by soil moisture and temperature. But I was not able to address that hypothesis during the analysis presented here.



  • Lognormal Regression of Flowering Phenology Data – the first approach I used to understand the patterns of flowering phenology present in my 2017 data, I used the fitdistrplus package in r to compare and fit phenomenological curves to counts of vulnerable (F-stage) capitula on given degree days. This required combining degree days of observation date and capitula counts into a single response variable to create a frequency distribution that a density functions could be fit to. I used splitstackexchange and methods described in Emerson and Murtaugh (2012) to create an observation of time (in growing degree days) for each event (F stage capitulum). This method led to the fitting of lognormal curves with the nls() function and parameter estimates given by the fitdist() function within The “goodness of fit” for the lognormal curve fit to all observations was assessed for data from each individual site as an attempt to describe variation in the timing of F-stage flowering by site. Visual comparison of lognormal curves fit to each site gave further insight into left-right shifts of these curves.


  • Assessing correlation in soil temperature and moisture data –in order to examine patterns in variation of soil temperature and moisture, particularly with respect to time, I converted dates and times to POSIXct objects using the lubridate After examining soil temperature against time of day of survey, I used the lm() function to fit linear regressions to this data series, and used the mutate function within the dplyr package to create a new variable of “adjusted soil temperature” normalized for time-of-day (Noon). I fit a linear regression to soil temperature against growing degree days for all sites, and repeated this for each individual site to see if time-of-season had significantly different impacts on soil temperature per site. I also looked at the relationship between soil moisture and time-of-day as well as time-of-season, but found less distinct trends and did not normalize for time-of-day.


  • Logistic Regression of probability of plant “vulnerability” against growing degree day, soil temperature, and soil moistureby simplifying phenology scores per capitulum to a binary response of “vulnerable” or “invulnerable” for an entire plant, I was able to analyze how the probability p of a plant being vulnerable varied with growing degree day, soil temperature, or soil moisture using logistic regressions (Hosmer et al., 2013). I used the glm() function in r to fit logistic regressions (family=binomial) to binary data for growing degree days, adjusted soil temperature (see part 2), and soil moisture. I used the predict() function to test how these fitted models performed in terms of classifying plants as vulnerable or invulnerable with a threshold value at 0.50 (datacamp 2018). I compared predicted response curves for the five sites visually, and made individual pairwise comparisons between two sites by using the coefficient estimates for site contrast parameters (Faustini 2000).


RESULTS: This project was an exploration of the data I had gathered in 2017, and produced a series of comparative figures showing variation in flowering phenology, measured environmental factors, and the relationship between the two. Early in analysis, a frequency distribution of stages counted per site per date showed that phenological synchrony (co-occurrence of F stage or prior with L4 and L5 larvae) was quite rare in the data (Figure 1).

Flowering phenology (as proportion of capitula counted in stage F) was best characterized by lognormal curves using scaled growing degree days (in hundreds). A curve was fit to data from all sites, and then both root mean squared error (RMSE) and mean absolute error (MAE) were calculated for data from each individual site relative to this all-site curve (summarized in Figure 2).

Figure 2. Lognormal curve fit by nls shown against F stage observations by site. Inset shows root mean squared error (RMSE) and mean absolute error (MAE) calculated overall and by site.

Figure 3. Lognormal curves of proportion of capitula found in F stage against growing degree days (hundreds).

This method showed some variation in goodness of fit but didn’t lend insight into the direction or magnitude of “shift” in timing between sites. A visual comparison of individual curves fit to each site showed slight shifts in the peak of each curve by site (Figure 3). This could be partially caused by the fact that a single survey represented different degree days between sites. However, Buck and Waterdog, which were modeled as having the same degree days, had different peaks in this figure, showing some ability to distinguish timing by site using this method.

Analysis of soil temperature found that data had a lurking effect of time-of-day of survey, which was characterized with a linear regression (Figure 4).

Normalizing soil temperature by time-of-day to a common time (Noon) allowed seasonal trends in soil temperature to emerge (Figure 5), where originally the impact of time-of-day obscured these trends (Figure 6).

Figure 5. Trends of adjusted soil temperature across the season (measured in hundreds of growing degree days).

Figure 6. Trends of non-adjusted soil temperature across the season (measured in hundreds of growing degree days)

Soil moisture was not found to be impacted by time of day, but was found to have low variability between most sites except one (Figure 7).

Logistic regression of a per-plant measure of “vulnerability” (as a binary response, 1 = any vulnerable capitula present and 0 = no vulnerable capitula present) against growing degree day showed a clear effect of growing degree day on the probability of plant vulnerability, with increasing degree days leading to a decrease in vulnerability with a clear threshold-like response (Figure 8).

Figure 8. Logistic regression of plant vulnerability (vulnerable = 1) against growing degree days (hundreds) with fitted curve shows a strong switch in probability of vulnerability from high to low with increasing growing degree days.

A similar regression against soil temperature (adjusted for time-of-day) showed that increasing soil temperature also decreased the probability of plant vulnerability, but with a more gradual response rather than a threshold-type response (Figure 9).

Figure 9. Logistic regression of plant vulnerability (vulnerable = 1) against adjusted soil temperature (°F) with fitted curve shows increasing adjusted soil temperature accounts for some decrease in probability of vulnerability.

Per site logistic regression of plant vulnerability against soil temperature showed that models (shown as predicted response curves) fit better for some sites than others, and slopes appeared to differ but not inflection points (Figure 10). This difference could be further resolved by fitting logistic regressions to two sites and analyzing the coefficient of the interaction term (which represents difference in slope of the second site) (Faustini 2000). The poorness of fit for Blair data in this case could be due to one missing survey date of soil temperature measurements.

Figure 10. Logistic regression of probability of plant vulnerability against adjusted soil temperature (degrees F) for all sites (grey dashed curve) and individual sites.

Soil moisture had an opposite relationship with plant vulnerability, with increasing soil moisture associated with increasing probability of plant vulnerability, using the same methods described above (Figure 11). Comparisons of predicted response curves by site for soil moisture showed a variety of slopes depending on the site, with all curves fitting reasonably well (Figure 12).

Figure 11. Logistic regression of probability of plant vulnerability against soil moisture (percent) shows an increase in probability of vulnerability with increased soil moisture.

Figure 12. Logistic regression of probability of plant vulnerability against soil moisture (percent) for all sites (grey dashed curve) and individual sites.


SIGNIFICANCE: Previous work showed that the loss of seeds due to cinnabar larval herbivory could significantly impact population dynamics for Senecio triangularis. The analyses presented here showed the variation in plant phenology, environmental factors, and the relationship between the two for the overall data set and for each site, individually. The biggest finding of this project was that phenological synchrony was too rare in this data set to model directly, as I had originally set out to do. Focusing instead on flowering phenology allowed development of methods and a greater understanding of the data captured in 2017.

For each variable examined, differences could be seen between sites with each site showing distinct behaviors. Further analysis of a full logistic regression model including all terms of interest could help understand the relative contributions of each explanatory variable to the change in probability of plant vulnerability.

For certain variables, it seems I may not have captured a wide enough range of variation in focal explanatory variables to draw clear connections between differences in phenology and these explanatory factors. For instance, soil moisture was largely similar between four sites and then quite distinct for a fifth site, Buck Mountain. Future work will focus on capturing greater variation in the abiotic environment, and hopefully also capturing greater variation in plant phenology.


LEARNING OUTCOMES: Through this project, my understanding and flexibility in R coding has increased by leaps and bounds. This includes restructuring data frames using dplyr pipelines (wide to long data, and collapsing nested data to compare at larger spatial/temporal scales), subsetting at will, and plotting data from every angle (mostly with ggplot). I discovered many useful packages (fitdistrplus, lubridate, splitstackexchange) and delved back into some previously encountered but less-well understood analyses (logistic regressions). I used ArcMap a few times to extract GPS points and other meadow statistics, but will hope to weave together the current analyses with explanatory variables not yet fully explored after the conclusion of this course.


STATISTICAL LEARNING: In terms of statistics used in this project, I learned about fitting phenomenological curves to data and analyzing goodness of fit using RMSE and MAE. I also learned about de-trending data to remove lurking effects such as the time-of-day effect in my soil temperature data. This will be very useful when dealing with other environmental factors that vary with time-of-day, and for which I do not have continuous data. I learned about fitting and analyzing logistic regression models to binary data using visual comparisons as well as numerical comparisons and interpretations for changes in slope and intercept. I was also able to understand and use methods developed by Emerson and Murtaugh (2012) for measuring “time of stage” for phenological data.



Datacamp. (2018). Logistic Regression in R Tutorial.

Diehl, J.W., and McEvoy, P.B. (1988). Impact of the Cinnabar Moth (Tyria jacobaeae) on Senecio triangularis, a Non-target Native Plant in Oregon. In Proceeding VII International Symposium on Biological Control of Weeds, (Rome, Italy), p.

Dunne, J.A., Harte, J., and Taylor, K.J. (2003). Subalpine Meadow Flowering Phenology Responses to Climate Change: Integrating Experimental and Gradient Methods. Ecol. Monogr. 73, 69–86.

Faustini, J.M. (2000). Stream channel response to peak flows in a fifth-order mountain watershed. PhD Thesis.

Hosmer Jr, D.W., Lemeshow, S., and Sturdivant, R.X. (2013). Applied logistic regression (John Wiley & Sons).

Kudo, G., and Suzuki, S. (1999). Flowering phenology of alpine plant communities along a gradient of snowmelt timing.

Murtaugh, P.A., Emerson, S.C., Mcevoy, P.B., and Higgs, K.M. (2012). The Statistical Analysis of Insect Phenology. Environ. Entomol. 41, 355–361.

Rodman, M. (2017). Non-target Effects of Biological Control: Ecological Risk of Tyria jacobaeae to Senecio triangularis in Western Oregon. Oregon State University.

Exploring historical constraints on fire across the Central Oregon Pumice Plateau

Filed under: 2018,Final Project,Final Project @ 5:45 pm


My original question was how historical fire occurrence varies spatially on the central Oregon pumice plateau. My question morphed into how did constraints on fire including climate, fuel abundance and continuity, and lodgepole pine forest influence the occurrence of fire.


My dataset consists of records of fire occurrence for 52 sample points distributed over an 85,000 ha area.

Study Area & Sample Points

At each sample point the annual occurrence of fire was reconstructed from tree rings during ~1700-1918.  Cross sections were removed from dead trees. Sections were sanded and precisely dated and injuries created by non-lethal low-severity surface fires were dated to their exact year of occurrence.


Fires scars collected at individual sample sites were composited into 1 record of fire occurrence at each sample point.  A range of 8 – 28 fires occurred at each sample point (mean 16).  In the Composite graph horizontal lines represent sample points and vertical tick marks represent fire occurrence.



I hypothesized that climate, time since fire, and lodgepole pine acts as constraints to fire occurrence across the central Oregon pumice plateau.


Earlier investigations demonstrate that fires historically occurred in dry forests in Oregon during drought years.  However fire size has not been related to climate.  It could be that fires of different sizes have additional relationships to previous year climate.  In exercise 2, I checked to see if small fires, large fires, and extensive fires were similarly related to climate.

Recent investigations north and northwest of this study area demonstrated lower fire frequency in lodgepole pine forests and that areas of lodgepole pine forest acted as intermittent barriers to fire spread.

Fuel is an obvious limiting factor to fire spread.  After fire fuel will need to recover sufficiently for fire to spread thus time since fire may be predictive of fire occurrence.


Ex 1 – Mapping fires from Binary Point Data

Prior to exploring constraints on fire occurrence I need to produce maps of fire occurrence.  To do this I used Arc GIS to evaluate using Thiessen Polygons, Kriging, and Inverse Distance Weighting (IDW) to map historical fires. Once I had made maps of fires I visually examined spatial and temporal variation in fire occurrence by creating an animation of fire events over time, mapping fire history statistics, and by identifying fire occurrence groups with cluster analysis and then mapping the groups.

Ex 2 – Determining how fire size and climate were related

I used superposed epoch analysis (SEA) to determine how annual climate and climate in antecedent years was related to fire occurrence.  By breaking fire events into size classes I was able to see if relationships with climate varied with fire extent.

Ex 3 – Using a GLMM to understand the influence of climate, time since fire, and lodgepole pine on annual fire occurrence

I used a generalized linear mixed model (GLMM) to determine the influence of fixed effects (climate, time since fire, and lodgepole) on annual fire occurrence across my study area.  Sample point was included as a random effect to see if the model varied across the study area. I used R2 for GLMM models to determine the relative importance of each fixed effect and the random intercept of sample point in the model.

The Results and the Significance

Ex 1

In exercise 1 I learned the tradeoffs between the 3 difference approaches to mapping fires from binary point data. The Thiessen polygon method provided most efficient, objective, and parsimonious method to map fire perimeters based on the distribution of my sample points.

Fire Extent mapped from Thiessen Polygons

Ex 2

Superposed Epoch analysis demonstrated that extensive fires occurred during years of extreme drought, large fires occurred during average droughts, and small fires were not related to dry or cool wet climate years. No relationships with climate in antecedent years for fires of any size or years without fire were found.


These results demonstrate that climate in antecedent years before a fire is not related to fuel abundance and connectivity. Cool wet years that may have been necessary to produce fine fuels that carry surface fires did not occur before fire years. Only a dry hot year during the fire year was associated with large and extensive fire spread. This suggests that fine fuel production with respect to fire spread is not related to climate. For my investigation this result means that fuel recovery following fire is not moderated by climate after fire. Thus time since fire may be an independent predictor of fire occurrence.

Ex 3

Using the GLMM approach I determined that climate, time since fire (interval), and lodgepole influenced annual fire occurred.  However, lodgepole was weakly significant, and only created a small change in the annual probability of fire. The conditional R2 for the GLMM model demonstrated that the random effect of site accounted for very little of the variance explained in the model. Climate and not interval as I expected explained most of the variance accounted for by the model.

For every 1 unit increase in PDSI the probability of fire occurring decreased by ~25%, for each year without fire (Interval) the probability of fire increased by 7%, and for each percent increase in lodgepole forest around a sample point the probability of fire decreased by 1%.

Significance of Ex 3 – Because lodgepole was significant, but weakly significant I want to explore other metrics that capture how lodgepole varies with respect to each plot.  It could be that lodgepole nearer the plot is more important or it could be its spatial pattern that causes variation in fire occurrence.  I was surprised that interval was less influential than climate. I interpret this as a strong indication that fuel is only limiting for a few short years after fire occurrence in the central Oregon pumice plateau. The small influence of site in the model suggests that bottom-up controls (microclimate, slope, local composition) that are site specific have little influence on fire frequency and fire region.

My next step is to focus on a single explanatory variable and switch to using site as a fixed effect with an interaction with that explanatory variable.  In addition I plan to fit 52 models, one for each site, using a single explanatory variable.  This would switch from a GLMM approach to a GLM approach, and allow me to map coefficients at each site on a map

Software Learning

In R I learned to use the LME4 package, learned how to perform SEA analysis, learned several new scripts for data wrangling, and learned how to make graphical summaries of fire history in ggplot.

In ArcGIS I learned to create animations of fire events and how to map fires with three different interpolations. I also used kriging to summarize variation in fire metrics in exercise 1.

Statistical Learning

I learned the pros and cons of different techniques to map fires, the ins and out of SEA analysis, and got an initial understanding of GLMMs. I plan to spend a lot more time working with GLMMs now that I’ve had this introduction and have found a flexible tiered approach to linear models.


Bolker, B.M.,Brooks,M.E.,Clark,C.J.,Geange,S.W.,Poulsen,J.R.,Stevens,M.H.H. & White, J.S.S. (2009) Generalized linear mixed models: a practicalguide for ecology and evolution. Trends in Ecology & Evolution, 24,127–135.

Hessl A, Miller J, Kernan J, Keenum D, McKenzie D. 2007. Mapping Paleo-Fire Boundaries form Binary Point Data: Comparing Interpolation Methods. The Professional Geographer 59:1, 87-104.

Florian Hartig (2018). DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models. R package version 0.2.0.

Michael H. Prager & John M. Hoenig (2011) Superposed Epoch Analysis: A Randomization Test of Environmental Effects on Recruitment with Application to Chub Mackerel, Transactions of the American Fisheries Society, 118:6, 608-618, DOI: 10.1577/1548-8659(1989)118<0608:SEAART>2.3.CO;2

Nakagawa S. and Schielzeth H (2012) A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods in Ecology and Evolution Vol 4(2):133-142.

Exploring recreational movement behavior through hidden Markov models

Filed under: 2018,Final Project @ 3:55 pm

Research Question Asked

Glacier Bay National Park and Preserve (GLBA), located in southeast Alaska, contains over 2.7 million acres of federally designated terrestrial and marine wilderness (National Park Service, 2015). Recreation users access GLBA Wilderness primarily by watercraft; the park lacks formal trail networks in its wilderness and terrestrial connectivity is fragmented by the park’s water resources. First designated as wilderness in 1980 through the Alaska National Interest Lands Conservation Act, management of the park’s wilderness has been guided by a 1989 Wilderness Management Plan (National Park Service, 1989). Much has changed in Alaska and GLBA since that time, and the park is currently engaged in updating its Wilderness Management Plan to adapt its management practices to modern management contexts. Park managers are particularly interested in developing a better understanding the wilderness experiences of water-based backcountry overnight users in GLBA Wilderness. As such, a dataset of global-positioning system (GPS) tracks of water-based visitor travel patterns were collected during the summer 2017 use season to record the spatial and temporal behaviors of recreationists in GLBA Wilderness.

The movement ecology paradigm provides a useful, organizing framework for understanding and conducting path analysis on GPS tracks of recreationists in GLBA Wilderness. Formally proposed in 2008, the movement ecology paradigm was designed to provide an overarching framework to guide research related to the study of organismal movement, with specific emphasis on guiding questions of why organisms move and how movement occurs through the lens of space and time (Nathan et al., 2008). The framework emphasizes understanding components of the movement, looking for patterns among those components, and understanding meaning behind movement through the underlying patterns (Figure 1). Ultimately, the target understanding is the movement path itself, which can be understood through quantification of the external factors, internal factors, and capacity for movement and/or navigation by the moving organism (Figure 2, Nathan et al., 2008). Through employing a movement ecology approach to the study of overnight kayaker movements in a protected area, individual movement tracks can be broken down into relevant components, the components of the path can be studied for patterns, and ultimately internal and external factors can be explored for influence or explanation of the movement path.

The following two aspects of the movement ecology framework were the focus of this study:

Movement States (Figure 1)– The movement ecology framework operationalizes movement as a series of step lengths and turning angles that can be used to identify underlying behavioral states. This follows the assumption that organisms move for specific reasons, and those reasons are reflected in the distance and direction of travel an organism travels in a set amount of time. The focus of this study was to transform raw GPS data into a path of step lengths and turning angles and to identify underlying behavioral states using step length and turning angle data.

Figure 1. Organizing framework for path analysis approach for this study (figure from Nathan et al., 2008): A) Understanding movement as a series of movement steps and associated turning angles, and B) Understanding values in the series of step lengths and turning angles as characteristic of behavioral states.

External Factors Affecting Movement (Figure 2) – According to Nathan et al. 2008’s framework, external factors are one of the factors influencing the movement path. In this study, external factors were operationalized as landscape-level features that may influence where a recreation travels in Wilderness and the type of behavior in which the recreationist engages. Two environmental variables, distance to shore and bathymetry were explored for influence on the movement paths.

Figure 2. Movement ecology framework (figure from Nathan et al., 2008). Of focus in this study is the exploration to two external factors, bathymetry and distance to shoreline, and their relationship to the movement path and associated behavioral states.

To apply the above-mentioned aspects of the movement ecology framework to my study of recreationist behavior in GLBA Wilderness, the following research questions were asked:

  • What are the mean step length and mean turning angle values for emergent behavioral states observed among the movement patterns of recreational kayakers?
  • How does bathymetry influence the transition probability between emergent behavioral states?
  • How does distance to shoreline influence the transition between emergent behavioral states?
  • Which external factor, distance to shoreline or bathymetry, has more explanatory power for transition probabilities in emergent behavioral states?

Description of Datasets Used

Dependent Variables Dataset:

GPS Dataset: The dataset for analysis was a test group of five GPS tracks taken from a sample of 38 GPS tracks collected during the summer of 2017 (Figure 3). Recreation grade, personal GPS units were administered to a sample of recreationists entering GLBA wilderness via personal kayak, June through August 2017. Study participants were asked to carry the GPS unit for the duration of their trip and return the unit at the end of their trip. GPS units continuously recorded movement throughout the trip.

Temporal Resolution and Extent: Units recorded a GPS point at various intervals, determined as a function of speed of travel. When speed was recorded at 0 miles per hour (MPH), the GPS units recorded an X,Y location point every 60 seconds. When speed was recorded at 1 MPH, the GPS units recorded an X,Y point every 15 seconds. When speed was 2 MPH or greater, the units recorded an X,Y GPS point every 8 seconds. The test group of GPS tracks were collected on trips taken in late June and early July 2017. The recorded tracks had between two and four days of data. For the purposes of this analysis, each of the test batch tracks were down-sampled such that data were aggregated into one-minute time bins, with X and Y data being averaged across the minute. In this way, the temporal resolution of the data during the analysis phase was 1-minute.

Spatial Resolution and Extent: At each time interval (described above), the GPS units recorded X and Y coordinates. Coordinates were recorded in decimal degrees. The geographic coordinate system for the data is GCS_WGS_1984. For analysis, the data were projected into the NAD_1983_UTM Coordinate System with a Zone 8N projection. The spatial extent for the dataset is the park boundary for GLBA.

Independent Variables Datasets:

Bathymetry: Bathymetry was incorporated using a 25-meter raster layer of the bathymetry underlying GLBA’s marine Wilderness. The bathymetry data layer was accessed from the publicly available National Park Service data clearinghouse available at the following link The downloaded data layer was added into an existing base map in ArcMap 10.3 that contained the analysis area (Glacier Bay National Park Wilderness) and the aggregated point shapefile of test batch movement data. The Extract Values to Points (Spatial Analyst) tool in ArcGIS was used to join raster cell values underlying the movement data together for analysis.

Distance to Shoreline Dataset: A landcover vector dataset was accessed from the publicly available National Park Service data clearinghouse available at the following link The downloaded data layer was added into an existing base map in ArcMap 10.3 that contained the analysis area (Glacier Bay National Park Wilderness) and the aggregated point shapefile of test batch movement data. Using the Joins and Relates option available in ArcGIS, the point-based GPS data were joined to the polygon (landcover) data. During the joining process, the distance from the landcover layer (the shoreline) was measured from each point, and added to the dataset as a unique field for each point. In this way, the distance from shoreline, in meters, was calculated for each GPS data point for use in analysis.

Figure 3. Map displaying datasets used in analysis. The yellow points on the map represented the aggregated shapefile of movement data used in the analysis. The blue area represents the bathymetry data layer. The distance from each yellow point to the nearest land cover class was calculated through a spatial join to operationalize the distance to shore dataset.


For each research question, the following hypotheses were developed:

  • What are the mean step length and mean turning angle values for emergent behavioral states observed among the movement patterns of recreational kayakers?
    • Two behavioral states will emerge from analysis of the step length and turning angle movement data, one state representing a movement-oriented state in which step lengths are longer and turning angles are narrower and the second state representing a resting-oriented state in which step lengths are shorter and turning angles are wider. Specific values for step length and turning angle for each behavioral state were not hypothesized.
  • How does bathymetry influence the transition probability between emergent behavioral states?
    • As bathymetry increases (i.e., as the depth of the ocean increases), kayakers will be more likely to be in a movement-oriented state rather than a resting-oriented state.
  • How does distance to shoreline influence the transition between emergent behavioral states?
    • As distance to shoreline increases (i.e., as the kayaker is further away from the shoreline), kayakers will be more likely to be in a movement-oriented state rather than a resting-oriented state.
  • Which external factor, distance to shoreline or bathymetry, has more explanatory power for transition probabilities in emergent behavioral states?
    • Distance to shoreline is likely to have greater explanatory power than bathymetry. This hypothesis is sort of an off-the-cuff assumption that because kayakers can see the shoreline and landscape features their behavioral state is more likely to be influenced by distance to shoreline than bathymetry, which would be a harder environmental variable to perceive while kayaking. The hypothesis is therefore based on assumptions about human-perception of environmental characteristics.


The overarching analytic approach for this study was use of hidden Markov models for operationalizing behavioral states within the test batch of movement paths (Langrock, et al., 2012; Michelot, et al., 2016). An R package, moveHMM (Michelot, Langrock, & Patterson, 2016), was the primary analytic tool used. The hidden Markov model approach requires that measured data, in this case time-stamped X and Y coordinates of movement actually represent movement data. The literature on the application of hidden Markov models for movement behavior suggests this requirement can be met by assuming that spatial inaccuracy does not exist within the data (i.e., the measured coordinates represent actual movement behaviors) and by regular time-stamped sampling of the GPS data (i.e., no missing data). Through assuming that the measured data represents a known state (in this case, movement), the hidden Markov model uses patterns in the measured data to reveal the “hidden” underlying states in the data. In this application of hidden Markov models, the hidden states being modeled are two behavioral states defined by combinations of co-occurring step length and turning angle movements. Additionally, covariates can be explored to understand how a co-occurring environmental variable that is changing in space and time with the step length and turning angle movement data may or may not correlate with behavioral shifts between states through hidden Markov models. Data processing and analytic approaches comprised three main phases, described below.

Phase 1 Summary: Relating Environmental Covariates to Step Length and Turning Angle Data (Exercise 2)

ArcGIS tools were used to spatially join environmental covariates (bathymetry and distance to shore) datasets to the point data prior to generating step length and turning angles. ArcGIS was also used as a data visualization tool to generate an overall understanding of the spatial extent of datasets being used.

Phase 2 Summary: Generation of Step Lengths and Turning Angles from GPS Data (Exercise 1)

The “prepData” function in the moveHMM R package was used to convert the series of X and Y coordinates into a series of step length, turning angle, and averaged covariate values. The “prepData” function requires that the GPS data be regularly sampled, and that each observation has a unique numeric ID code associated with the data in the data frame – data processing and summarizing tools in R were used to meet these requirements. Additionally, an R function was used to convert data from a latitude and longitude geographic coordinate system to a projected, UTM coordinate system to generate meaningful step length values of meters through the prepData tool.

Phase 3 Summary: Fitting the moveHMM Models and Evaluating Results (Exercise 3, Parts 1 and 2)

The “fitHMM” function was used from the moveHMM R package to run the hidden Markov models, generate behavioral probabilities, and run an AIC analysis on the fit models.


Step Length and Turning Angle Generation:

The step length and turning angle histograms suggest that across the five test tracks processed using the moveHMM tool, the majority of step lengths across the tracks were between 600-700 meters per minute and kayakers generally traveled in a straight direction turning infrequently. The step lengths were greater than originally anticipated, and after further investigation I learned this is likely because kayakers have the option of taking a day boat back to the visitor center from their backcountry location rather than paddling back to the visitor center. This new information was a surprise to receive so late in the analysis, but likely explains why the step lengths of 600-700 meters per minute occur in the dataset. Future exploration of these techniques on this dataset will need to account for this discovery, likely through elimination of the day boat portion of the GPS track from analysis. The 0-100 meter per minute step length bin likely represents stoppage time. The disparity between the two step length categories suggests that in the case of water-based recreationists, step length may be a good metric for examining changes in behavioral state in future bivariate analyses.

Turning angle histograms suggest that for the most part, kayakers are traveling in a relatively straight direction with little variation away from that direction. When turning movements are made, they tend to turn toward the left. This appears to be the result of making a loop trip in which kayakers initiate their trips following the eastern coast of GLBA inlet and finish their trips following the western coast of GLBA inlet.

Figure 4. Histograms of generating turning angles and step lengths for movement path two.

Best Fitting Models:

Two models were fit to the step length and turning angle data using the moveHMM tool: a two-behavioral state model with distance to shore as a covariate for behavior transition and a two-behavioral state model with bathymetry as a covariate for behavior transition. The decision to model a two-behavioral state model came from review of the step length and turning angle histograms and conversation with my research advisors, instructor, and classmates. Figure 5 reports the model outputs for these model runs.

Figure 5. Results from the model outputs for the best fitting parameters with distance to shore as a covariate (left panel) and bathymetry as a covariate (right panel). The distance to shore model suggests a two-state model, with one state being characterized by shorter step lengths (mean = approximately 1 meter) and wider turning angles (mean = 3.1 degrees) and a second state being characterized by longer step lengths (mean = approximately 224 meters) and narrow turning angles (mean = -0.014 degrees).The bathymetry model suggests a two-state model, with one state being characterized by shorter step lengths (mean = approximately 2.5 meters) and wider turning angles (mean = 3.1 degrees) and a second state being characterized by longer step lengths (mean = approximately 262 meters) and narrow turning angles (mean = -0.007 degrees).

Visual Presentation of Results and AIC Analysis

Results of the AIC analysis to determine the favored model suggest that model 2, the bathymetry model is favored over model 1, the distance to shore model, given AIC values of 72708.23 and 82097.23 respectively. Given these results, the remaining visualizations and interpretation of results is provided for Model 2, with bathymetry as a covariate.

Figure 6. Density histogram of step length states for the Bathymetry Model. The density distribution and model fit curves to not suggest practical significant visual differences in the step length distributions for states 1 and 2. However, the figure suggests that the distribution for state 2, the movement oriented state, has a higher density across step lengths from slightly greater than 0 to approximately 300 meters in length than state 1, which peaks right around a step length of 0.

Figure 7. Density histogram of turning angle states for the Bathymetry Model. Turning angles for State 2 cluster around 0, which is expected given the narrow turning angle movement mean reported for this state. The curve of State 1 is unexpected, given that the mean of the turning angle movement reported for this state is around 3.

Figure 8. Distribution of behavioral states across time for track ID 2. The figures show that the majority of the movement track is spent in state 2, the moving state, and that state 1, the resting state, occurs infrequently throughout the movement track.

Figure 9. Transition probability matrix for influence of bathymetry on transition between or among behavioral states. The four transition probability plots show that as bathymetry increases (water depth increases), the likelihood of staying in state 1, the resting state, decreases dramatically at a small increase in bathymetry and then remains at 0. Similarly, at small increases in bathymetry, the probability of transitions from state 1, a resting state, to state 2, a moving state increases rapidly and then states at a probability of 100%.

Figure 10. Displays three separate outputs for Track 2 that begin to tell a story of what may be going on in the data. A) Example output for location of behavioral states in space for track ID 2. In theory, the figure would show blue to characterize portion of the track during which the individual is in state 2 and orange for portions for the track during which the individual is in state 1. At this scale, state 1 is hard to see, but there are small clusters around the beginning and end of the trip that show some orange track pieces. This figure shows that for individual 2, the majority of the track is in state 2 behavior rather than state 1 behavior. B) A display of the track (the orange color has no relevance) overlayed on a satellite image of the surrounding area. The few instances of state 1 resting behavior in the track co-occur in space with access to a glacier, suggesting that potential landscape features other than bathymetry and distance to shore may be better predictors of behavior state changes. C) Temporal display of track three, with transition in color through the track displaying the passage of time. The color gradient changes somewhat rapidly in a small amount of space where the state 1 behavior occurs.

For each research question, the following hypotheses were developed:

  • What are the mean step length and mean turning angle values for emergent behavioral states observed among the movement patterns of recreational kayakers?
    • Result: A two-behavioral state model was developed, with the following parameter estimates for step length and turning angle by state derived from the best fitting model, the bathymetry model: The bathymetry model suggests a two-state model, with one state being characterized by shorter step lengths (mean = approximately 2.5 meters) and wider turning angles (mean = 3.1 degrees) and a second state being characterized by longer step lengths (mean = approximately 262 meters) and narrow turning angles (mean = -0.007 degrees).


  • How does bathymetry influence the transition probability between emergent behavioral states?
    • Result: As bathymetry increases, the likelihood of staying in a stationary behavioral state, if already in a stationary behavioral state, decreases rapidly. Similarly, the likelihood of transitioning between a stationary state and a movement-oriented state increases rapidly and abruptly as bathymetry increases.


  • How does distance to shoreline influence the transition between emergent behavioral states?
    • Result: As distance to shore increases, the likelihood of staying in a stationary behavioral state, if already in a stationary behavioral state, decreases rapidly. Similarly, the likelihood of transitioning between a stationary state and a movement-oriented state increases rapidly and abruptly as bathymetry increases. (Results not pictured in this exercise)


  • Which external factor, distance to shoreline or bathymetry, has more explanatory power for transition probabilities in emergent behavioral states?
    • Result: Through an AIC analysis comparing the Bathymetry Model and the Distance to Shore Model, the Bathymetry Model was favored over the Distance to Shore Model as the preferred model for modeling two-state movement behavior. This result is counter to my hypothesis that distance to shore would be the favored model.

Despite these findings, I am not overly confident in the model outputs and conclude that additional work is needed to fully apply and understand the use of hidden Markov models for understanding movement behavior. First, the sample of tracks modeled is only representative of 5 movement patterns, and in reviewing these movement patterns they are quite distinct from each other. Additionally, imbedded in the data is the potential for tracks collected while recreationists were taking the dayboat rather than while independently kayaking. This fact not only impacts the step length and turning angle but also the underlying motivations related to external factors that served as the original premise for this work. Additionally, I think the model would perform better on a larger sample of data to help iron-out some of the disparities in the parameter estimates. Finally, the modeling exercise presented assumes that one distribution for each state is appropriate for the population, therefore, any variability in the individual track data is minimized in the model. This is only a helpful approach if the researcher believes that all actors behave in a similar capacity – this is a big assumption for my data and given this modeling exercise I do not necessarily think it holds.

Significance of results to science and resource managers

Conceptually, I think the application of hidden Markov models and the use of the moveHMM tool has great potential for exploring movement data outside the realm of animal-based movement, particularly for natural resources managers. The study of human movement in recreation settings through the collection of GPS track data is a relatively new methodological development. As such, analytic methods applied to spatial data in the social sciences to date have focused on providing descriptive summaries of track lengths and travel times and summarizing the data in aggregate. Using a movement ecology approach provides a lens through which to understand individual movement patterns and to look at variation in both space and time within a track. These methods also have the potential to help researcher and managers understand how landscape level features may influence movement patterns and what the characteristics of movement patterns are in certain places and at certain times. As researchers and managers begin to build empirical relationships between behavioral states and landscape level features, managers can work to engage in more directed landscape restoration, visitor use management, and park planning.

Learning: What did you learn about software a) Arc-Info, b) Modelbuilder and/or GIS programming in Python, c) R, d) other?

Through this class, I worked primarily in ArcGIS and R, exploring new tools and developing additional skills in both software packages. Through working in ArcGIS, I explored tools for completing spatial joins with both raster and vector data, data manipulation and formatting challenges and limits, and began to develop a rudimentary process for working between ArcGIS and R. Through the course, I moved beyond my prior knowledge of ArcGIS, using tools through the course exercises that I had not previously used such as the generate random points tool, extract values to points tool, spatial join with distance calculation, and mosaic raster tool.

The majority of my work in this class was performed in R, an analytic programming language in which my own prior experience was through completing homework for STAT 511 taken last term. My proficiency in R grew immensely through this course. Skills I gained through this course include opening, manipulating, and saving excel and csv data files in R, map-based data visualization tools, data summary tools, and several movement based analytic packages including adelhabitat and moveHMM. A central learning opportunity for me in this course was to be exposed to the various ways that R can be used to analyze and manipulate spatial data. Moreover, prior to this class, I had no concept of the depth of the spatial analytical packages available through R. I was introduced to those packages through recommendations by classmates, and the experience has changed how I will approach analysis in the future.

I did not work with Modelbuilder in the course, nor did I work with Python. I have worked with those programs in the past and was happy to develop new skills in R and additional skills in ArcGIS.

Learning: What did you learn about statistics, including a) hotspot, b) spatial autocorrelation, c) regression, and d) multi-variate methods?

Through the presentations of my classmates and through my own project work in this course I expanded my understanding of how to apply statistical concepts to spatial data analysis. Through several presentations given my classmates through Tutorial exercises, I have developed a working knowledge of autocorrelation, including what the analysis seeks to identify, how the term “lag” is operationalized in the analysis, and tools/packages for completing autocorrelation analysis. To date, this concept has been a term that has eluded me, and while I did not work directly with autocorrelation analysis for my own project work, I feel the exposure gained in this class has helped me move forward in understanding what an autocorrelation analysis seeks to accomplish.

For multi-variate statistical methods, I took a deep dive into exploring hidden Markov models as a mechanism for understanding transitions among behavioral states and how the potential for relationships with environmental covariates. I also ran an AIC analysis to identify a favored HMM. I also was exposed to new thinking about how to explore relationships between landscape variables through Exercise 2, beginning to develop a new way of thinking about how landscape level variables are related and why I might be interested in knowing that relationship when looking at the GPS tracks of use. I also was excited to be exposed to geographic weighted regression through Sam’s presentation, and look forward to exploring this analytic tool in the future. I liked the ability of the tool to bring in the spatial component of the data as a sort of third dimension to the analysis.

More than anything, this course has introduced me to new ways of thinking about spatial data generally, and exposed me to the wide range of possibilities available for future analyses. One of the greatest values of the course has been hearing the presentations of my classmates and learning about their research problems, questions, and tools used.


Langrock, R., King, R., Matthiopoulos, J., Thomas, L., Fortin, D., & Morales, J. (2012). Flexible and practical modeling of animal telemetry data: hidden Markov models and extensions. Ecology 93(11): 2336-2342.

Michelot, T., Langrock, R., & Patternson, T. (2017). An R package for the analysis of animal movement data. [online].

Michelot, T., Langrock, R., & Patternson, T. (2016). MoveHMM: An R package for the statistical modelling of animal movement data using hidden Markov models. Methods in Ecology and Evolution 7: 1308-1315.

Nathan, R., Getz, W.M., Revilla, E., Holyoak, R., Saltz, D., & Smouse, P.E. (2008). A movement ecology paradigm for unifying organismal movement research. PNAS 105(49): 19052-19059.

National Park Service. (1989). Wilderness Visitor Use Management Plan: Glacier Bay National Park and Preserve.

National Park Service. (2015). Glacier Bay: Wilderness Character Narrative. Available:

Modelling snow for sheep

Filed under: 2018,Final Project @ 1:21 pm

RESEARCH QUESTION – How do seasonal snow conditions affect Dall Sheep recruitment?

Dall an emblematic species of alpine regions in high latitude North America. Their ranges extend from the mountains of the Yukon Territory, Canada, to the furthest western extent of the Brooks Range in Alaska. Populations of Dall Sheep have declined 21% range-wide since the 1990s with a major mechanism of decline thought to be the increased frequency of extreme spring snow conditions(Alaska Department of Fish and Game, 2014). During the months of April and May mature Dall sheep ewes typically give birth to one lamb. The survival of this lamb is dependent on the mother’s ability to protect it from predators and guide it to accessible forage. If successful, the lamb is ‘recruited’ into the population. Hence, a commonly used metric of animal population growth potential is the mother to child ratio, or in this case the lamb to ewe ratio (hereafter written as lamb:ewe). Extreme spring snow conditions are thought to decrease lamb survival by limiting access to forage, either by deep snow coverage or ice-layers formed in the snow subsequent to rain-on-snow events. The limited forage could cause starvation or increased use of areas where vulnerability to predation is increased. In this project I will examine this question via use of a spatially explicit snow-evolution model, SnowModel, and lamb:ewe ratios from summer sheep surveys.


In this project I used three primary datasets;

Snow / climate dataset; SnowModel (Liston and Elder, 2006)was used to simulate daily snow and climate conditions in 6 different Dall sheep domains where survey data was available. SnowModel was forced with the Modern Era Retrospective Reanalysis for Research and Applications (MERRA2) product (Gelaro et al., 2017). SnowModel effectively downscales temperature, humidity, precipitation, wind speed and direction from a 0.5º by 0.625º to a 30m resolution and physically evolves and distributes a simulated snowpack across digital elevation and landcover layers, derived from the IfSAR DTM distributed by US Geological Survey and the NLCD 2011 product distributed by the Multi-resolution Land Characteristics Consortium respectively (Homer et al., 2015). For 5 of the domains (Brooks Range, Denali, Gates of the Arctic, Lake Clark and Yukon Charley), where available, in-situ data on snow depth and snow water equivalent were used to calibrate and validate the model. The 6thdomain, in the Wrangell St Elias, in-situ data from a March 2017 field campaign was used to calibrate and validate the model and test for model performance (see below). SnowModel is run for the entire period between September 1st1980 and August 31st2017. Daily data for snow and climate above the elevation of shrubline were then aggregated into monthly and seasonal metrics e.g. mean monthly snow depth (m). Seasons in this case were taken as September to November (Fall), December to February (Winter), March to May (Spring) and June to August (Summer).

Figure 1; Map of Alaska with Dall Sheep ranges

Sheep data; Sheep data used here is from annual surveys completed by Alaska Department of Fish and Game, Bureau of Land Management, US Fish and Wildlife Service and National Park Service in the same area as the SnowModel domains. Lamb:ewe ratios were calculated from the number of lambs recorded and the number of ewes and ewe-likes. Survey methods include distance sampling, stratified random sampling, and minimum count methods either from a ground location or fixed wing aircraft.

Climate Indice data; Climate indices of larger scale weather patterns were downloaded from the National Oceanographic and Atmospheric Administration (NOAA) for 7 different indices; the Pacific Decadal Oscillation (PDO), Arctic Oscillation (AO), East Pacific / North Pacific Oscillation (EP/NP), North Pacific Pattern (NP), West Pacific Pattern (WP), Pacific North American Index (PNA), and the North Atlantic Oscillation (NAO) (


The hypotheses for this project fell were delineated by blog post, the three installments of which looked at a) model performance, b) climate indices and snow condition relationships, c) snow condition and lamb:ewe ratio relationships. They are hence;

  1. SnowModel performs best at the elevations and landcover where the greatest amount of in-situ data in available
  2. The influence of different climate indices on snow conditions is not uniform throughout Dall sheep ranges
  3. Spring snow conditions have the greatest impact on lamb:ewe ratios surveyed in summer


For hypothesis A I used a multivariate approach to test the hypothesis, using the FAMD tool of the FactoMineR library in R to test where SnowModel was over or under predicting snow depth by land cover class and other metrics of topography (elevation, slope, aspect and northerness).

For hypothesis B I used autocorrelation and crosscorrelation functions in R to test for patterns within the time series of the snow metrics and in correlation with the climate indices.

For hypothesis C I used a simple approach comparing whether snow conditions were above or below their mean to whether lamb:ewe ratios were above or below their mean to see which conditions and when might best predict levels of recruitment. This then informed a multiple logistic regression model computed in R


Hypothesis A;

The results from my hypothesis A can be best seen in figures 2 to 4;

eigenvalue variance.percent cumulative.variance.percent

Dim.1 3.1793016        24.456166                    24.45617

Dim.2 2.1249171        16.345516                    40.80168

Dim.3 1.4969078        11.514675                    52.31636

Dim.4 1.2530914         9.639165                    61.95552

Dim.5 0.9902287         7.617144                    69.57267

Figure 2; Scree plot and table from the multiple factor analysis

Figure 3; Plot of the importance of the quantitative variables importance to the 1stand 2nddimension of the factor

Figure 4; Plot of the qualititative variables importance to the 1stand 2nddimension of the factor analysis.

The factor analysis revealed that only 40% of the variation in the data could be explained by the first 2 dimensions (figure 2), with elevation being the biggest contributor to the 1stdimension, and the category of SnowModel error (diffCategory) being the biggest contributor to the 2nddimension (figures 3 and 4). However, and perhaps expectedly given its elevational transition and role in in snow accumulation processes, landcover type was a significant contributor in both dimensions. Broadly speaking as the elevation increased, and the landcover went from coniferous forest to prostrate shrub tundra and bare ground, model accuracy increased. Interestingly there is also a pattern of underprediction in bare and coniferous forest landcover and over prediction in prostrate shrub tundra and erect shrub tundra, although the magnitude of the error is greatest in the landcover classes that typically populate lower elevations. As I had greater amounts of in-situ data from higher elevations I could confirm my hypothesis, however the analysis did reveal an under and overpredict pattern between bare ground and prostrate shrub tundra that I didn’t expect.

Hypothesis B;

I conducted analyses for all 6 domains described above, however the autocorrelation of monthly and seasonal snow/climate metrics from 1980 to 2017 and cross-correlations of monthly and seasonal total snowfall from 1980 to 2017 to climate indices, did not reveal any meaningful patterns, though some statistically significant results at certain lags were observed they never rose above ~0.5 for the autocorrelation and ~0.25 for the cross correlation. This suggests that larger term patterns of climate do not explain a large proportion of the inter and intra annual variability in snow conditions in alpine areas. Of note is the stronger influence of different climate indices for different domains, although only weakly significant there appears to be a pattern dependent on latitude and continentality, please refer to blog post 2 to examine this. In the meantime hypothesis 2 can be considered open for further testing. For this post please see the examples from the Wrangell St Elias below as illustrative;

Figure 5; Autocorrelation of monthly variables 1980 to 2017, Wrangell St Elias domain


Figure 6; Autocorrelation of seasonal variables 1980 to 2017, Wrangell St Elias domain


Figure 7; Cross-correlation of monthly variables to climates indices 1980 to 2017, Wrangell St Elias domain

Hypothesis C;

To test hypothesis C a simplistic approach was utilised to confirm the occurrence of the following conditions;

  1. High month/season snowfall and low lamb:ewe ratios
  2. Low month/season air temperature and low lamb:ewe ratios
  3. High month/season snow depth and low lamb:ewe ratios
  4. Low month/season forageable area and low lamb:ewe ratios

These statements are based on what we expect the relationship between snow conditions and sheep recruitment might be – increased snowfall and snow depth, or lower air temperatures and forageable area, produce conditions where greater energy expenditure is required for survival. Dall sheep are in calorific deficit during the snow season so benign conditions mean that ewes reach the lambing period in better condition and are potentially then more able to provide for their lambs, increasing the observed lamb:ewe ratio. An alternative or complementary idea is that conditions during and after lambing are more important as lambs require a narrower range of conditions than adult sheep to survive. From the frequency of occurrence of agreement in these conditions we can select variables for use in a logistic regression model.

I will present my results only for the Wrangell St Elias domain by month, however seasonal results and other domains have been analysed but not yet interpreted.


Figure 9; Wrangell St Elias results by month

Figure 10; Log regression result

From the simple comparison we can see that spring month (March, April, May) snow conditions that are more hazardous do not always predict low lamb:ewe ratios any more than winter months (December, January, February). From the logistic regression model, where I included the variables by month that had the most frequent occurrence of meeting the conditions, February forageable area and the total amount of November snowfall came out in the final model as being strongest.


Each of my blogposts / hypotheses took me further along towards a better idea about how best to address and answer my research question. Hypothesis A described where SnowModel was performing best, and worst, and gives a quantifiable error to propagate through the analysis (although I have not done that here). Most importantly it gave me reasonable confidence that the model is doing reasonably ok in sheep territory and refined the areas where I can spend energy trying to improve it.

Hypothesis B, and its tests, were murkier, producing know results that really jumped out for their confirmation or rejection of the hypothesis. Further work would be to aggregate the indices and snow seasons further into yearly, or September to August, averages and compare them to lamb:ewe ratios as wells as snow and climate metrics.

Hypothesis C produced a model that rejected the original hypothesis, suggesting that snowfall in the autumn and the area available to forage in February is most important in. predicting lamb:ewe ratios. However, the pseudo R-sq (derived from its AIC score) of this model is not much to shout about and could likely be refined by introducing seasonal values and more variables into the analysis. Logistic regressions weren’t conducted in other domains than the Wrangell St Elias – it would be interesting to compare results to these to examine whether there are range-wide similarities.

The results from hypothesis C can be used to interpret whether years where no sheep survey took place are likely to have below and above average lamb:ewe ratios, giving an indication as to whether the observed population decline is a result of recruitment being affected by increasing occurrence of hazardous snow years.


Using R and Rstudio for the first time was a significant feature of this class for me. I began to appreciate the value of organising data according to the tidyverse philosophy and saw incremental improvement in my abilities to conduct analyses and create figures using the packages available. I’m not yet certain whether it is an overall improvement on Matlab, which I previously used, but I do prefer the ease of ggplot2 over equivalent Matlab tools.

I was pretty unversed statistically at the beginning of this class and have gained an understanding in the principles, use and concepts of Principle Component Analysis, Multiple Factor Analysis, autocorrelation and cross-correlation spatially and temporally, and multiple logistic regression. My future work will be much enhanced as a result.


Alaska Department of Fish and Game, 2014. Trends in Alaska Sheep Populations, Hunting, and Harvests. Division of Wildlife Conservation, Wildlife Management Report ADF&G/DWC/ WMR-2014-3, Juneau.

Gelaro, R., McCarty, W., Suárez, M.J., Todling, R., Molod, A., Takacs, L., Randles, C.A., Darmenov, A., Bosilovich, M.G., Reichle, R., Wargan, K., Coy, L., Cullather, R., Draper, C., Akella, S., Buchard, V., Conaty, A., da Silva, A.M., Gu, W., Kim, G.-K., Koster, R., Lucchesi, R., Merkova, D., Nielsen, J.E., Partyka, G., Pawson, S., Putman, W., Rienecker, M., Schubert, S.D., Sienkiewicz, M., Zhao, B., Gelaro, R., McCarty, W., Suárez, M.J., Todling, R., Molod, A., Takacs, L., Randles, C.A., Darmenov, A., Bosilovich, M.G., Reichle, R., Wargan, K., Coy, L., Cullather, R., Draper, C., Akella, S., Buchard, V., Conaty, A., da Silva, A.M., Gu, W., Kim, G.-K., Koster, R., Lucchesi, R., Merkova, D., Nielsen, J.E., Partyka, G., Pawson, S., Putman, W., Rienecker, M., Schubert, S.D., Sienkiewicz, M., Zhao, B., 2017. The Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2). J. Clim. 30, 5419–5454.

Homer, C., Dewitz, J., Yang, L., Jin, S., Danielson, P., 2015. Completion of the 2011 National Land Cover Database for the Conterminous United States – Representing a Decade of Land Cover Change Information. Photogramm. Eng. 11.

Liston, G.E., Elder, K., 2006. A Distributed Snow-Evolution Modeling System (SnowModel). J. Hydrometeorol. 7, 1259–1276.

Considering Beaver Dam Occurrence based on Stream Habitat and Landscape Characteristics

Filed under: Final Project @ 1:18 pm
  1. Research Question

Q1: How does the Suzuki and McComb Habitat Suitability Index (HSI) relate to observed beaver dams and the West Fork Cow Creek drainage?

 Q2: What other factors explain the selection of suitable habitat?

 At the beginning of this effort, the question I was most interested in was “How well does the Habitat Suitability Index (HSI) developed by Suzuki and McComb (1998) beaver dams and the West Fork Cow Creek drainage located in the South Umpqua River Basin in Southern Oregon.  The HSI suggests that beaver dams are mostly likely to occur where stream gradients are less than 3%, active channel widths are three to six meters, and valley bottom widths are at least 25 meters.

The second question I had was what other variables might explain the selection of some suitable stream reaches for dam building, while others were seemingly ignored.  Other models such as the Beaver Restoration Assessment Tool, developed by McFarland et al (2017) out of Utah weight a number of other variables including flow permanence and vegetation as possible limitations on dam building that I wanted to consider.  However, because of dataset availability I decided, instead to follow another line of inquiry from landscape ecology that considers the influence of habitat size and connectivity.

  1. Datasets:

To answer this question I considered two primary datasets. The first dataset is, Netmap, a proprietary stream network layer generated through a combination of digital elevation models (DEM) and data that was collected through state and federal agency stream surveys. This included estimates of stream gradient, active channel width, and valley bottom width.  The second data set was a collection of survey locations randomly selected from the stream network layer, stratified by locations considered to be suitable and unsuitable for beaver damming based on the HSI criteria .  Given the rarity, however, beaver dam occurrences in a watershed, the sample was weighted toward sites considered suitable for beaver damming. This dataset also included observations from surveys that took place in the fall months of 2017.  In total, 48 beaver dams were observed from the survey locations, principally on two streams in the drainage.


Map 1: West Fork Cow Creek








  1. Hypotheses:

H1: All observed dams will occur in stream reaches classified as suitable by the Suzuki and McComb HSI. 

 H2: Suitable stream reaches that are longer and located in closer proximity to other suitable reaches will are more likely to be occupied than those that are smaller and more isolated.

Figure from Dunning et al. 1992 showing how habitats A and B, while both too small to support a population, may be occupied (A) if in close proximity to other habitats









  1. Approaches:

The most challenging but most informative approaches I used all related to geoprocessing in ArcMap, which I was largely unfamiliar with as a software package prior this class.

Using SQL to identify suitable damming habitat:

The first effort included using SQL to identify reaches in the stream network data that fit the criteria.  The language was a little clunky at first but after a little time was simple enough to select the appropriate thresholds for each of the habitat variables from the HSI .

Converting suitable reaches into habitat patches:
The second approach required that I convert contiguous reaches of suitable habitat into habitat ‘patches’ so that I could eventually calculate the length of each patch and the distance of that patch to its nearest neighboring patch.  While seemingly simple, there was not one tool in Arc to accomplish this task.  All told I used a combination of functions in Arc including Buffer, Dissolve, Join and Spatial Join to accomplish this.

Calculating distance between patches:

To accomplish the last approach, I burrowed a tool designed for traffic engineers that can calculate the distance from one point to another via a specified network.  In most cases the network is a series of roads but for this problem, I used the stream network polylines, since there is literature to support that beavers will have high fidelity to water as they move through the landscape.


  1. Results:

Exercise 1 produced two metrics that I mapped.  The first was the patch length, and the second was distance to the nearest neighboring patch.

Map 2: Suitable dam habitat and survey locations


In Exercise 2, I demonstrated that the habitat criteria from the HSI correctly identified all areas where dams were observed, but that not all reaches meeting the HSI criteria had dams sites.  In short, based on these data, the HSI criteria seemed to be necessary but not ultimately sufficient for damming.




In Exercise 3, I then looked at the relationship between patch length, distance to nearest neighbor, and a final metric that consider the size of the nearest neighbor weighed by the distance to that neighbor. As a result of aggregating continuous stream reaches into single features, the number of patches with observed dams decreases to 4 and was problematic for a logistic regression analysis because the small sample size can lead to Type 1 errors.


Map 3. Habitat Patch Length based on HSI criteria

Map 4. Distance to Nearest Neighboring Patch









  1. Significance

I found two important take-home messages from these efforts. The first is that the HSI developed by Suzuki and McComb (1998) seems to work relatively well for identifying suitable dam habitat in the West Fork Cow Creek in sites, in that observed dams all fell within the criteria thresholds.  However, there is still quite of bit of variation in the suitable habitat locations where beaver dams were observed that cannot be explained by these data.  The second take-home message is that the selection within habitats appears to coincide with the landscape hypothesis that habitats are more likely to be selected based on their size and connectivity.  Due to sample size limitations, however, the strength of evidence for this is not conclusive.

These results are important to managers because there is a growing emphasis on the use of beavers to improve stream characteristics thought to be conducive to salmonids, particularly, Threatened Oregon Coast Coho Salmon (Oncorhynchus kisutch).  To date most management efforts in this regard have focused on relocating beavers from locations where dam building is problematic (e.g. road culverts) to areas where dams will not be a nuisance and accessible to anadromous fish.  Yet, our anecdotal understanding from conversations with managers in the Umpqua Basin is that these efforts have not been guided by an evaluation of suitable habitats for this area.  In other words, this work suggests that, at a minimum, relocation efforts should focus on streams that meet the HSI criteria.

  1. Geo-processing learning

These efforts were helpful to me in two ways.  First, it was useful to consider the HSI criteria relative to the observed dam sites that we found and improves my confidence that these may be necessary criteria for dam building, at least in similar drainages.  Secondly, I had relatively little experience with ArcMap, and what I did was from nearly a decade ago.  These exercises greatly improved my familiarity and confidence using this software.  This is particularly important as we scale this research up to the larger Umpqua Basin where the number of stream reaches increase by a factor 10.  I’ve even begun automating the geo-processesing steps for the patching and OD Cost Matrix in model builder which surpasses any expectations I had in the beginning of the quarter.

  1. Statistical learning

From statistical standpoint I am disappointed that I was not able to consider more sophisticated analyses.  In particular I had hope to apply logistic regression and apply a predictive map to inform my survey selection for the rest of the Umpqua basin but that was not feasible because of my small sample sizes.  As it was, I was limited to very basic analyses such as the permutation test given the small and lopsided data.  Yet, as I turn to defining my protocol for surveys in the rest of the Umpqua I have a much better appreciation for the relative rarity of dam site occurrences and that I will need to generate my sample accordingly so that I can develop a larger dataset to allow more robust statistical procedures.  Despite this, it was helpful to see how other students approached their problems and consider what of those tools I may apply to future datasets.


Dunning, J. B., Danielson, B. J., & Pulliam, H. R. (1992). Ecological Processes That Affect Populations in Complex Landscapes. Oikos, 65(1), 169–175.

Suzuki, N., & McComb, W. C. (1998). Habitat classification models for beaver (Castor canadensis) in the streams of the central Oregon Coast Range. Retrieved from





Using correlation-based techniques to investigate population trends in Bull Kelp (Nereocystis luetkeana) in southern Oregon

Filed under: Final Project @ 9:52 am

The Question: I was initially looking to explore correlation between a kelp canopy coverage data set and a suite of environmental variables. My question morphed into examining the correlation between canopy cover and two temperature data sets.

The Data: I used three datasets to investigate this question. The first was a 35-year time series of kelp canopy cover in southern Oregon derived from Landsat satellite imagery. This dataset was shared with me by my colleague Tom Bell and reported the percent cover of each 30mx30m pixel whenever cloud free imagery was available. The other two datasets were two time series of temperature on the Oregon coast, one derived from satellite imagery and the other from direct measurements of intertidal temperature near Port Orford. The satellite-derived data set was over 35 years long, but came at a resolution of 0.1 degree raster cells, meaning that it represented offshore sea temperatures more so than local, nearshore temperatures. The local, intertidal dataset was measuring local, nearshore water temperatures that were within 20 miles of most of my kelp beds. That time series only covered about half of my time series however, from late 1999 to the present.


The Hypotheses: Worldwide, warming temperatures have been implicated in global declines of kelp populations (Wernberg et al.,2016; Filbee-Dexter et al., 2016). Perhaps even more important than the direct effects of temperature is the fact that ocean temperature is closely correlated with nutrient availability, which is of crucial importance to these fast-growing primary producers. However, in many local or regional studies, temperature does not necessarily emerge as one of the top drivers of kelp growth and biomass. For example, in southern California giant kelp biomass was not significantly impacted by the extreme 2013-2015 warm water anomalies that impacted the Pacific coast of North America (Reed et al., 2016). Furthermore, several studies on giant kelp (Macrocystis pyrifera) in southern California have found that wave intensity as one of the primary factors controlling giant kelp biomass (Bell et al., 2015, Parnell et al., 2010). Biotic interactions, such as competition with understory kelp and herbivory by urchins, can also tightly control kelp populations (Dayton and Tegner, 1984). Overall, untangling the relative importance of abiotic drivers of kelp populations, is a highly context dependent business. In southern Oregon, I expect temperature/nutrient availability to have an effect on kelp biomass but for that effect to be moderated by the presence of urchins, wave events, and climate oscillations.


The Methods: To investigate the relationship between kelp and temperature, I employed three methods: 1) Autocorrelation and cross correlation to investigate temporal correlation within and between the two biggest kelp patches in southern Oregon.
2) Graphical representations and cross correlation to examine patterns between satellite-derived sea surface temperatures and local, directly-measured nearshore temperatures.

3) Interpolation of my kelp time series to monthly temperatures using polynomial splines in order to utilize wavelet and cross-wavelet analysis to look for areas of shared power within and between the kelp time series and the two temperature time series.

The Results and the Significance

Exercise 1: Autocorrelation and cross correlation were quite limited in the kelp patches I looked at, even at the yearly scale (e.g. Fig 1). This suggests that a) the species responds quickly to changing environmental conditions rather than to holding on to momentum from previous population sizes and 2) that local factors (less than 20km) are more important in driving kelp population size than local factors

Figure 1: Cross Correlation between maximum annual kelp cover at Rogue Reef and Orford Reef over a 35 year time series. The lag is in years.

While this finding was one of the simplest, it was also one of the most significant for me. This finding suggests that focusing on local, patch-specific dynamics will be important in untangling environmental drivers of kelp in Oregon. I’ve already utilized this finding when responding to feedback from managers. An ODFW employee suggested that my satellite-derived kelp time series was incorrect because it showed moderately sized populations since 2014. He said that Orford Reef, one of the largest kelp beds in the state, had been practically non-existent in past years, so there was no way regional populations were at anything other than historical lows. I looked into patch-specific dynamics over the past 4 years and found that, while Orford Reef had indeed gone essentially to zero since 2014, other patches in Oregon had recovered in that same time period and were bolstering regional kelp cover numbers.

Exercise 2: Overall, the relationship between SST and nearshore temperatures was fairly consistent. The temperatures matched very well in the winter, but then nearshore temps warmed up less and more slowly in the summer months. The average difference between mean annual temperature for the two datasets was consistent as well, usually staying between 1.1-1.5 degrees (Fig 2). The rankings of the warmest versus coolest years were also very similar between the datasets.

Figure 2: Difference in annual mean temperature between satellite derived sea surface temperature (black) and measured, nearshore temperatures (red) from 2000-2016. Note how consistent the difference between the two is.

While I can use satellite SST to infer average nearshore temperatures, it will not give me good insight into the variability of nearshore temperature, which could be important to kelp population regulation. If I want to try to incorporate this variability into my analyses, it may be useful to look further into the effect of upwelling and terrestrial water cycles on nearshore temperatures. However, considering a) that much of the interannual variability in temperatures (hottest to coolest years) is similar between satellite and nearshore temperatures and b) how much time it might take to more finely predict local temperatures from satellite temperatures, I think that the satellite temperatures may be good enough to use for most of my analyses.

Exercise 3: According to wavelet analysis, strong annual summer peaks in canopy were not a consistent pattern in my kelp time series, but rather only in the 1985-1991, 1999-2000, and 2014-2017 periods (Fig 3). This inconsistency extends to the cross wavelet analysis of kelp and temperature, where only certain years have high-power, annual cycles. Another important finding from this exercise was that sometimes there appear to be high-power interannual oscillations between kelp and temperature, suggesting climate oscillations may also be influencing kelp cover.

For me, there are two important takeaways from this analysis. First, it is particularly interesting to see a high degree of power in the 2014-2017 period because in these years a) we have had a boom in urchin populations, b) we have had anomalously high water temperatures and c) kelp in northern California has collapsed during this period. All of these other factors suggest that annual kelp population oscillations would have a smaller amplitude, not a larger. And second, the high power periods and years for temperature (both satellite and intertidal) were not necessarily the same as those for kelp. To me, this indicates something other than a 1:1 relationship between temperature and kelp. This suggests that I need to explore other environmental variables and that I might need to explore them in multivariate analysis rather than multiple univariate analyses.

Figure 3: A) Wavelet analysis of the kelp canopy time series and B) Cross-Wavelet analysis of kelp canopy and satellite temperature. Areas of high-power are in red. The black line in A represents the power ridge and the black arrows in B represent whether the two are in phase (right), out of phase (left), x leading (down) or x lagging (up).

I think my next steps with this project will be to use PCA to try and identify particular environmental factors that are most strongly correlated with kelp canopy cover. Incorporating my findings from Exercise 1, I will do it utilize PCA not only on the regional Oregon kelp time series, but also using the local population time series from individual patches. In this way, I will be able to look at regional drivers of kelp canopy cover as well as how the extent to which local drivers may differ a) between patches and b) from the regional drivers.


The Learning: I almost exclusively utilized R for this class, and gained familiarity with interpolation techniques as well as a number of very user-friendly, useful packages for performing correlation analyses. Most of the techniques I utilized in this class dealt with autocorrelation, in one form or another, but I learned about a number of other spatio-temporal statistics via the student presentations. In particular, I’m excited to utilize Geographically Weighted Regression to better understand local kelp patch correlations with various environmental variables and PCA to identify environmental variables that may be mostly closely correlated with kelp cover.

The References:

Bell, T. W., Cavanaugh, K. C., Reed, D. C., & Siegel, D. A. (2015). Geographical variability in the controls of giant kelp biomass dynamics. Journal of Biogeography, 42(10), 2010-2021.

Dayton, Paul K., and Mia J. Tegner. “Catastrophic storms, El Niño, and patch stability in a southern California kelp community.” Science 224.4646 (1984): 283-285.

Filbee-Dexter, K., Feehan, C. J., & Scheibling, R. E. (2016). Large-scale degradation of a kelp ecosystem in an ocean warming hotspot. Marine Ecology Progress Series, 543, 141-152.

Reed, D., Washburn, L., Rassweiler, A., Miller, R., Bell, T., & Harrer, S. (2016). Extreme warming challenges sentinel status of kelp forests as indicators of climate change. Nature communications, 7, 13757.

Wernberg, T., Bennett, S., Babcock, R. C., de Bettignies, T., Cure, K., Depczynski, M., … & Harvey, E. S. (2016). Climate-driven regime shift of a temperate marine ecosystem. Science, 353(6295), 169-172.





June 9, 2018

Relationships between environmental features and recreationist behavior in Grand Teton National Park

Filed under: 2018,Final Project @ 7:12 pm

The research question that you asked.

Broadly, my research question at the beginning of the course was: “what spatial and temporal patterns emerge from day-use hikers in Grand Teton National Park?”

I appreciated starting out with a broad, exploratory question as it allowed me to think creatively and learn a variety of approaches for analyzing and conceptualizing human behavior using GPS data. As the course progressed, and as I examined patterns within the data, my research question became more specific:

“What are the relationships among the spatial and temporal behavior of recreationists and the environmental features within the recreation area?”

A description of the dataset you examined, with spatial and temporal resolution and extent.

Throughout the course, I used a sub-set of five GPS tracks (i.e. five hikers) collected on July 19, 2017. This subset came from a collection of 652 GPS tracks of day-use visitors at String and Leigh Lakes in Grand Teton National Park. The GPS units were distributed to a random sample of visitors between July 15 – September 8, 2017. Each intercepted visitor was asked to carry the GPS unit with them throughout the duration of their visit at String and Leigh Lakes. When deploying the units, study technicians also recorded the total number of people in the group, and the intended destination for their day visit. To maintain independence between samples, only one GPS unit was given to each group.

The GPS units used in this study were Garmin eTrex 10 units. These units collected point data every 5 seconds. The GPS tracks were saved as point features for analysis in ArcGIS so that each visitor’s hiking path can be represented by a series of points. The positional accuracy of these units can vary up to 15 meters. However, the Garmin units were calibrated with a high accuracy Trimble GPS unit which indicated a low average positional error of 1.18 meters.

Hypotheses: predictions of patterns and processes you looked for.

I hypothesize that people will spend more time, have shorter step lengths, and more acute turning angles the closer they are to a water feature.

This hypothesis is grounded in an assumption that summertime recreationists are drawn to open viewscapes, particularly cooling water features like lakes, waterfalls, and streams. The recreation site the data were collected in contains stunning lakes that are couched right on the edge of the Teton mountain range. Perhaps hikers will feel compelled to stop and sightsee the closer they are to these water features, and will have less stopping behavior the further away they are from these features.

Approaches: analysis approaches you used.

I sourced data layers representing vegetation cover and elevation from I used tools in ArcMap to build and calculate attributes for analysis. All analyses were conducted in R.

Exercise 1: (a) I used R to plot spatially explicit graphs of human movement through space and time. This allowed me to visually examine how people were using the recreation site and provide context for subsequent analysis.

(b) I created histograms representing the proportions of the step lengths and turning angles of the five recreationists. This approach provided me with a better understanding of the distributions and characteristics of the response variable.

Exercise 2: I generated box plots and scatter plots that represented the relationships among various environmental features — my independent variables — in the study area. The variables I examined were: vegetation, elevation, and distance to water.

Exercise 3: I conducted OE analysis (observed vs. expected) to examine if hikers were spending significantly more or less time in certain vegetation types despite being near water. I followed this with a chi-square test to see if the differences in distributions of observed vs expected were statistically significant.

Final analysis: (a) I ran the data through a hidden Markov Model package in R “moveHMM” to identify the probability of people changing from state 1 (short steps, acute turning angles) to state 2 (long steps, obtuse turning angles) as a function of their distance to water.

(b) I also attempted doing a multiple linear regression on the data. Unfortunately I couldn’t normalize the distribution of step length, realizing this approach may not be an adequate choice for this dataset.

Results: what did you produce — maps? Statistical relationships? other?

Throughout the course I produced numerous graphs representing recreationists movements through space and time. Figure 1 demonstrates a simple plot of the five tracks I worked with. This plot indicates most people in this sample stayed close to shore and relatively close to the parking lot. This makes sense as the trail hugs the lake shore.

Figure 1. A visual representation of the five tracks used throughout the analysis. All tracks collected on July 19, 2017.

Figure 2 summarizes the distribution of step lengths and turning angles for all five hikers. Step length is calculated as the distance between one GPS point to the subsequent GPS point (all GPS points have a temporal resolution of 1 minute). This figure suggests that these hikers typically walked straight, and nearly 40% of the time their step lengths, or distances between points, were 41-60 meters.


Figure 2. Distributions of turning angle (left) and step length (right).

Once I had a better understanding of the features of my response variable, I was curious to learn more about the environmental variables that I thought could be influencing human behavior. I also wanted to explore these variables to check for any confounding relationships between human movement and distance to water (my overarching hypothesis).

Examining the environmental features suggested that elevation may not be an influential variable in relation to the behavior of the recreationists (Figure 3). The range in elevation within the hikers’ movement paths was only 10 meters. However, the results did indicate that the predominant vegetation type for the hikers was coniferous woodland, and that conifer was primarily located near the lake shore. Perhaps the presence of conifer is playing a role in the amount of time the hiker spends near water. In other words, is the conifer deterring a person from stopping, even though they are near water? These questions encouraged me to examine the relationships among the recreationists time spent near water and the vegetation type they were in.

Figure 3. Box plots representing vegetation type relative to elevation (top left) and distance to water (top right). Elevation plotted against distance to water.

Before diving into analyzing the relationship between time, distance to water, and vegetation, I first calculated the overall proportions of the length of the hikers’ tracks grouped by distance class, and the time spent in each distance class. Figure 4 indicates that the five hikers spent 75% of their time 0 -20 meters from the lake. Additionally, of the entire track length, over 60% of the track is 0 – 20 meters from the lakeshore. From an ocular check, it looks like my hypothesis is being supported: people are spending more time near the water. But, is this a function of the lake feature or some other variable?

Figure 4. Proportion of time spent in the area compared to the length of the track (left). Frequency of vegetation types along every distance class (right).

The next set of figures illustrate the expected amount a time a person would spend in specific vegetation types within each distance class if they were traveling at a constant rate compared to the observed amount of time that person actually spent in each vegetation type within each distance class.  Coniferous woodland was the first vegetation type I examined.

Additionally, I was introduced to a neat way to conceptualize differences in distributions, particularly when the distributions represent different units of measurement (i.e. time vs. length). Creating an odds ratio essentially divides the values of time vs. length to give me a ratio. A value less than 1 indicates that people spent less time in the distance class than they would have if they had been traveling at a constant pace. A value greater than one indicates the people spent more time in distance class than they would have if traveling at a constant pace.

Figure 5  indicates that people were spending less time in conifer despite being closer to water! The further away from water, the more time they spent in conifer. The chi-square test indicates these differences in distributions are statistically significant (X2  = 90.359 df = 7, p-value < .001).

Figure 5. Observed/Expected for coniferous vegetation types grouped by distance class.

Figure 6 represents the observed versus expected distributions for herbaceous vegetation types (ex: open meadows). Interestingly, it appears that the hikers spent more time in meadows than was expected. It also appears that meadows are prevalent further away from shore. Back in Figure 5 we learned that people are spending, on the whole, more time in areas 0-20 meters from shore and 61-80 meters from shore. Perhaps for the 61-80 meters distance class, this can be explained by the presence of open meadows. However, the chi-square test indicates that these distributions are not significantly different from one another (X2= 7.4523 df = 3, p-value = .06).

Figure 6. Observed/Expected for meadow vegetation types grouped by distance class.

Figure 7 is very revealing as it indicates that people are spending more time in residential facilities that are close to shore. This includes areas with picnic tables, paved sidewalks, and some areas of the parking lot. This makes sense as I imagine people are drawn to the amenities that are offered. Also, the chi-square test indicates that the differences in these distributions are significant (X2 = 197.86, df = 3, p-value < .01) Seeing these results, I wonder if next time I should modify my hypothesis and examine if/ how distance from vehicle explains behavior.

Figure 7. Observed/Expected for residential facilities grouped by distance class.

These results suggested that environmental variables may play a larger role in human movement and behavior than originally hypothesized. As a final step in this exploratory journey I experimented with the moveHMM package to identify what the probabilities are of a person changing from state 1 (short step lengths, acute turning angles) to state 2 (long step lengths, obtuse turning angles) as a function of their distance to water. It is important to note that I still don’t have a comprehensive understanding of hidden Markov Models. Therefore, the following figures represent only a fraction of the output generated from this package. Further, the following figures aim to demonstrate how outdoor recreation scientists could potentially use these methods to better understand human behavior.

Figure 8 represents the two state model distributions of step length and turning angle. Figure 9 shows us that, indeed, as the distance to water increases, the likelihood a person will transition from state 1 (slow) to state 2 (fast) increases. According to this output my hypothesis is being supported!

While HMMs provide a statistically rigorous framework for incorporating covariates, allow for the autocorrelation commonly experienced with GPS data, and enable researchers to make inferences about changes in behavioral states, my own knowledge of these processes remains limited. Although it was relatively simple for me to push the data through an HMM model, this final step in my analysis still leaves me with a lot of questions. Further, I was unable to build an HMM model that controls for the vegetation type; thus I have a feeling that only including distance to water as a covariate may not be telling the whole story.

Figure 8. State-dependent distributions in the 2 state model.

Figure 9. Effect of the covariate ‘distance to water’ on the transition probabilities

As a final exercise, I was curious to see if it was possible to do a multiple linear regression on my own as another way to address the hypothesis. First, I checked for autocorrelation in the data. The output below represents one track.

Figure 10. Autocorrelation on distance to water (above) and step length (below).

The acf informed me that the spacing at which observations are no longer autocorrelated is about 2 minutes. I subsequently took a subsample of track to only select points that are greater than two minutes apart. After re-sampling I re-tested for autocorrelation – both variables are no longer autocorrelated (Figure 11).

Figure 11. Autocorrelation on re-sampled dats on distance to water (above) and step length (below).

I then tested for normality. The distributions were skewed, so I did a square root transformation on the distance to water. This normalized the distribution and a Shapiro-Wilk normality test provided a p-value > .05, suggesting the distribution is not statistically different from the normal curve.

Figure 12. Density plot of distance to water before (above) and after (below) transformation.

Unfortunately, I wasn’t able to normalize step length. Doing a log transform or square root just made the skewness even more severe.

Figure 13. Density plot of step length. The responses are not normally distributed.

At this point, I realized that doing a regression was perhaps not the best approach – especially since I wasn’t able to normalize the distribution of the response variable (step length). However, despite the cliff hanger ending, I intentionally included this final exercise in the blog post to demonstrate the use and consideration of autocorrelation on GPS data.

Significance. What did you learn from your results?  How are these results important to science? to resource managers?

These results revealed to me that environmental features such as vegetation and open viewsheds have an influence on the behavior of recreationists. Further, the results indicate that vegetation type  also plays more of a role in the behavior of recreationists than originally anticipated. I was surprised to see how much vegetation type influenced the amount of time the hikers spend in the area.

Overall, these results also demonstrate how recreationist activity can be measured and analyzed to develop deeper understandings of behavior.

These results are important to both science and resource managers. Parks and protected area land managers strive to provide a quality user experience while also protecting natural and cultural resources. Accurately understanding how people move and behave in a recreation system allow for more informed management decision making. For example, understanding what environmental conditions influence behavior could indicate a need for additional infrastructure, signage, or educational initiatives depending on the management objectives for the area. Additionally, by using statistical tools to analyze these relationships the results can have predictive power for managers.

In the scientific and academic communities, applying spatial methods to outdoor recreation science allows for a more accurate understanding of how people move, experience, and interact in outdoor spaces. By integrating GIScience with other common social science techniques in outdoor recreation — such as surveys, observations, and interviews — scientists glean richer results that can support and contribute to existing theory, generate deeper understandings about human behavior, and inspire additional studies.

Your learning: what did you learn about software (a) Arc-Info, (b) Modelbuilder and/or GIS programming in Python, (c) R, (d) other?

I worked exclusively in ArcMap and R. I went into this class with very limited knowledge and minimal confidence with the software; now, I am coming out of the class with more confidence in my ability to learn and figure it out. I appreciated the exploratory, self-guided structure of the class which provided space to play around with the technology and make mistakes. Even more so, I appreciated the insight of Julia and Laura, as well as the tips, tricks, and support from fellow classmates.

I learned how to work with GPS data in ArcMap and source data layers that were relevant to my research question. I learned various tools in ArcMap that allowed me to make calculations (i.e. TrackAnalyst tool and spatial join tools).

I learned how to work with spatial data in R and became more efficient at aggregating and manipulating dataframes for analysis. I also learned how to represent the data in a way that is meaningful to other audiences.

What did you learn about statistics, including (a) hotspot, (b) spatial autocorrelation (including correlogram, wavelet, Fourier transform/spectral analysis), (c) regression (OLS, GWR, regression trees, boosted regression trees), and (d) multivariate methods (e.g., PCA)?

I learned how to do expected vs. observed analysis and chi-square test for significant differences in the distributions. I learned more about the application of hidden Markov Models in developing probability models representing changes in human behaviors. I also was able to work a little with autocorrelation.

Beyond my project I learned quite a bit from listening to the tutorials of other classmates. I am now (ever so slightly) more familiar with concepts like spatial and temporal autocorrelation, kriging, and geographically weighted regression. Additionally, it was extremely helpful working with Susie who had a similar dataset and was also using a hidden Markov Model approach.


Michelot, T., Langrock, R., Patterson, T. (2017). An R package for the analysis of animal movement data. Available:

June 8, 2018

Interannual variation in phenology and productivity at a C3 and a C4 grassland

Filed under: 2018,Final Project @ 12:48 pm

1. Research question


Climate change is altering the production and distribution of plant species (Kelly & Goulden, 2008); however, the response of grasses to climate change is relatively understudied compared to woody plants, especially in tropical regions (e.g., Schimel et al. 2015). Globally, grasslands and savannas are estimated to comprise 30 percent of non-glacial land cover, and many important crop species are grasses (Still et al. 2003). Grasses are predicted to respond more quickly to climate change because of their short life span and dispersibility, and are aggressive invaders that can exacerbate fire cycles, cause land cover transitions, and alter carbon cycling (D’Antonio and Vitousek 1992, Angelo and Daehler 2013).

The community composition of a grassland mediates its response to its environment, and is critical to consider in forecasting the climate change impacts (Knapp et al., 2015). Grasses with the C4 photosynthetic pathway, in contrast to the ancestral C3 pathway, have comparatively higher resource-use efficienciesight, water, and nitrogen, especially under high temperatures. As a result, species using the C3 or C4 pathway will have distinct responses to a warming climate and rising CO2 (Collatz et al., 1992, 1998; Lloyd & Farquhar, 1994; Suits et al., 2005). Understanding the differential environmental constraints on the phenology, productivity, and distribution of grasses of different functional types will be critical in predicting how they will respond to future climate change, as well as identifying areas that may be prone to invasion from non-native grasses.


My initial research question, was, broadly: How are phenology and productivity of natural grasslands related to local climate variables, and how do these relationships differ between C3 and C4 grasslands?

Specifically, the research questions I answered were:

  1. Is interannual variation in production the result of multi-year patterns in production (autocorrelation)?
  2. How does the timing of the relationship between production and local climate variables differ between a C3 and a C4 grassland, and between years?

2. Data

My study sites are two eddy covariance (EC) flux tower locations in natural grassland areas located ~90 miles apart in eastern Kansas. The sites experience nearly identical climates, but the first is a natural tallgrass prairie composed of 99% C4 grass at Konza Prairie Biological Station outside Manhattan, KS, while the second is a replanted agricultural field composed of 75% C3 grass at the University of Kansas Field Station, outside Lawrence, KS. Because the two sites experience a very similar climate, despite being in distinct ecoregions, I hypothesize that photosynthetic type strongly controls differences in production and phenology at each site.

At these sites, production is measured first from satellite data using the normalized differential vegetation index (NDVI), which essentially measures the greenness of vegetation. In this analysis, I use NDVI derived from the MODIS satellites, which has 16-day temporal resolution, and is available from 2001 to present.

Production is also measured at the EC flux towers, and is calculated from ecophysiological equations that use ground-based measurements of atmospheric gas concentrations and meteorological data. The eddy covariance (EC) flux approach uses tower-mounted instruments to measure atmospheric concentrations of water and CO2, as well as air temperature, solar radiation, and other environmental data. All measurements are taken continuously every 30 minutes. EC flux data reflect the “footprint,” or area upwind of the tower where the instruments are mounted. The footprint varies with wind speed and direction, but averages about ~250m2. EC flux data span from 2008-2015.

The main metrics I am interested in are: daily total gross primary production (GPP) and the mean daily light use efficiency of production (LUE). GPP is reported in units of  gC•day^-1•m^-2. LUE is calculated as the daily sum of gross primary productivity (GPP) divided by the daily sum of the amount of incoming radiation, or photosynthetic photon flux density (PPFD). Daily LUE is converted to units of gC•MJ^-1•day^-1•m^-2 from units of µmol CO2• µmol photon^-1•day^-1• m^-2 using the molecular weight of carbon and Planck’s equation.

3. Hypotheses

My hypotheses were:

  1. Interannual variation in the phenology of production will not be strongly autocorrelated. Interannual variation in production is more strongly controlled by local climate variables
    1. a. Phenology will differ in the C3 vs C4 sites: the C4 site will have a longer growing season and a higher mean and maximum production.
  2. The C3 and the C4 site will have distinct cross-correlative relationships between production and local climate. They will differ in the timing and the strength of this relationship, and will respond distinctly to interannual variation in climate.

4. Approaches

My approach was to:

  1. Calculate annual phenology metrics from NDVI at each site, using the R package [greenbrown](
    • Phenological metrics include: the timing of the start of season and end of season, the length of the season, and mean and peak NDVI over the season. The phenological metrics are calculated by fitting curves to the annual growth cycle at each site.
    • After calculating the annual phenology metrics, I calculated autocorrelation using the acf() function n R to assess whether annual differences in phenology were a product of change over time, or cyclical trends.
  2. Calculate cross-correlation between GPP, LUE, and soil water content at each site, using the R function ccf(). Assess how cross-correlation differs between water-years  at each site.

5. Results


Fig. 1: Start, end, and length of growing season over time at each site.

Fig. 2: Mean (MGS) and Peak growing season value of NDVI at each site.

The C3 site has a consistently longer growing season than the C4 site, with an earlier start of season and an later end of season. However, based on the NDVI data, the mean and peak growing season values look like they may be the results of cyclic patterns. I fit an autocorrelation function to assess whether there were interannual patterns in the phenological indices.

Fig. 3: Autocorrelation function for start, end and length of season, at each site. Sites are plotted with distinct line types (KFS, C3 =solid; Konza, C4 = dashed).

Fig. 4: Autocorrelation function for mean and peak growing season value, at each site. Sites are plotted with distinct line types.

In the autocorrelograms above, the dashed lines represent confidence intervals for determining whether autocorrelation is statistically significant. Because the autocorrelation function rarely, and marginally, rises above the confidence intervals for any of the phenological indices, this indicates that the phenology we observe is more likely to result from interannual climatic variation– rather than cycles in production.

Production and cross-correlation with climate

After inferring that interannual variation in phenology is likely not the result of temporal autocorrelation of the indices, I investigated the relationship between production and climate. Specifically, I examined the the timing and magnitude of the relationship between GPP and soil water content. Importantly, I am comparing the timing of cross-correlation between water years, which are measured from October to September. Water years, in contrast to calendar years, are a more biologically meaningful division to use when examining the relationship between plant growth and water variables.

For this project, I will just show the relationships between GPP and soil water content. Because LUE is derived from GPP, the cross-correlation between LUE and soil water content looks similar to the relationship between SWC and GPP.

Fig. 5: Annual time series of GPP (gC / m^2/ day) and soil water content (% * 10), for 2010-2013, plotted by water year. The thick lines represent the monthly average of total daily GPP for each site, in each year. The dashed lines represent the average daily soil water content. The daily total of GPP for each site and for each year is potted in thin, light-colored lines in the background.

Inspecting the time series of GPP and soil water content, by water year, we can perceive a clear hump in SWC that precedes a hump in GPP. Calculating cross-correlation between SWC and GPP, using the ccf() function in R, allows us to quantify this relationship.

In my phenological questions, I acknowledge the caveats of NDVI data– which capture distinctions in the timing, but not the magnitude of difference between production of C3 and C4 grasses (Fig. 7). I was curious whether cross-correlation between NDVI and soil water content would reveal similar or different trends from the cross-correlation between GPP and soil water content.

Fig. 6: Time series of GPP and NDVI at the C3 and C4 sites, from 2010-2013. NDVI is multiplied by 10 in order visualize on a similar scale to GPP. Mean monthly GPP is plotted in the thick, solid lines; the daily average of GPP for each site is plotted in the background using thin, lighter-colored lines.

Importantly, NDVI does not accurately capture the magnitude of difference between C3 and C4 production. While C3 and C4 grasses have distinct relationships between NDVI and GPP, NDVI alone is not always a great indicator of production differences between the two photosynthetic types.

Fig. 7: Cross-correlation strength vs. lag for GPP and soil water content, for 2010 -2013. The cross-correlation between NDVI and soil water content for each site in each year is plotted in lighter-colored lines in the background of each plot.

The vertical lines above are plotted at x = 0, and x = 100, to simplify visual comparison. Dashed lines represent confidence intervals.

When we examine how cross-correlation between GPP and soil water content varies between sites and between years, we see a clear shift at Konza Praire, the C4 grassland in 2012. This shift indicates that GPP and LUE are peaking earlier in relation to soil water content: with a lag of about 30 days, as opposed to 100 days at the Kansas Field Station, and at Konza in previous years.

Because 2013 followed the severe drought year, 2012, this shift suggests that the C4 grasses at Konza may be more flexible in shifting the timing of their water use in response to drought or other irregular moisture conditions.

Importantly, NDVI reveals the same patterns of cross-correlation between production and soil water content.

6. Significance.

First, my results confirm that C3 and C4 grasses respond distinctly to similar climatic variation. The shift in the timing of cross-correlation of GPP with soil water content in 2013 (Fig 7) is evidence that the same climatic change will affect these photosynthetic pathways in different ways. We need to continue using both experimental and observational research to quantify the interacting effects of rising temperatures, increased aridity, and elevated CO2 on the production of C3 and C4 grasses.

Second, my results enlighten me that though NDVI does not capture differences in the magnitude of production at C3 vs C4 sites, it does capture differences in phenology / timing of production. This means that NDVI, though a coarse metric of production, is valid to use for assessing variation in the timing of production. Particularly because NDVI is widely used and easy to obtain and calculate, this distinction will be helpful in future phenology research. However, importantly, phenological metrics that calculate mean, max, or total growing season productivity– i.e., metrics that describe the magnitude of production, rather than the timing of production– should be interpreted from NDVI with caution. Ideally, other production indices should be used to assess variation across functional types in the magnitude of production.

7. My learning, w/r/t software

Over the course of this class, I have become much more confident using R Markdown to edit and generate reports for publication. This blog post was drafted and formatted using R Markdown, with only minimal modifications to fit the blog format. I have also learned:

  • how to gapfill time series data using the function zoo::na.approx, which uses linear interpolation to replace NAs in a time series
  • how to convert outputs from acf() and ccf() functions into data frames that are easily plottable using ggplot

8. My learning, w/r/t statistics

I became much more comfortable generating and interpreting correlograms. I used autocorrelation and cross-correlation to determine which variables to proceed with in my analysis.

Future directions will involve using the information I now have about cross-correlationto select variables for future regressions. E.g., I learned that GPP is cross-correlated with soil water content, and as a result, I am more interested in investigating the shifts in the timing of that relationship, than the magnitude. I learned that GPP is more strongly correlated, than cross-correlated, with air temperature and soil temperature. As a result, will use linear mixed models to assess the mean annual relationships between daily total GPP and means of these variables.

Ultimately, this work prepares me to quantify and discuss the similarites and differences in the timing of key phenological and production events for a C3 vs a C4 grassland.



Angelo, Courtney L., and Curtis C. Daehler (2013). “Upward expansion of fire‐adapted grasses along a warming tropical elevation gradient.” Ecography 36.5 , 551-559.

Collatz G, Ribas-Carbo M, Berry J (1992) Coupled Photosynthesis-Stomatal Conductance Model for Leaves of C4 Plants. Australian Journal of Plant Physiology, 19, 519.

Collatz GJ, Berry JA, Clark JS (1998) Effects of climate and atmospheric CO2 partial pressure on the global distribution of C4 grasses: Present, past, and future. Oecologia, 114, 441–454.

D’Antonio, Carla M., and Peter M. Vitousek (1992). “Biological invasions by exotic grasses, the grass/fire cycle, and global change.” Annual review of ecology and systematics 23.1 : 63-87.

Kelly AE, Goulden ML (2008) Rapid shifts in plant distribution with recent climate change. Proceedings of the National Academy of Sciences, 105, 11823–11826.

Knapp AK, Carroll CJW, Denton EM, La Pierre KJ, Collins SL, Smith MD (2015) Differential sensitivity to regional-scale drought in six central US grasslands. Oecologia, 177, 949–957.

Lloyd J, Farquhar GD (1994) 13C Discrimination during CO₂ Assimilation by the Terrestrial Biosphere. Source: Oecologia, 994, 201–215.

Schimel, David, et al. (2015) “Observing terrestrial ecosystems and the carbon cycle from space.” Global change biology 21.5 : 1762-1776.

Still CJ, Berry JA, Collatz GJ, DeFries RS (2003) Global distribution of C3 and C4 vegetation: Carbon cycle implications. Global Biogeochemical Cycles, 17, 6-1-6–14.

Suits NS, Denning AS, Berry JA, Still CJ, Kaduk J, Miller JB, Baker IT (2005) Simulation of carbon isotope discrimination of the terrestrial biosphere. Global Biogeochemical Cycles, 19, 1–15.

June 6, 2018

Exploring the relationship between floating guidance structures, hydraulics, and the location of behavior changes in fish

Filed under: Final Project @ 7:43 pm

Question, dataset, and approach:

This research investigated the hydraulic and behavioral impacts of a floating guidance structure in an experimental channel on juvenile Chinook salmon. Three exercises were conducted to determine if any relationships exist between channel hydraulics (which were stationary, measured at discrete locations and interpolated to become spatially-continuous) and the locations of behavior changes (as determined by a behavioral change point analysis tool, “smoove”). Exercise 1 examined differences in the spatial distribution of behavior changes at 20, 30, and 40 degree guide wall angles. Exercise 2 determined which hydraulic variables were most predictive of the location of behavior changes angles using multivariate methods. Exercise 3 used geographically weighted regression (GWR) to examine spatial variation in the relationship between each of the above variables and the location of behavior changes.


1) Because the hydraulics created by guide wall angles of 20, 30, and 40° differ from one another, we expect the location of behavior changes to vary in space with angle.

2)  Because similar flume experiments have found hydraulic thresholds in velocity gradient for fish behavior, we hypothesize that velocity gradient (or another hydraulic variable) will predict the downstream location of behavior changes at all 3 guide wall angles.


Although the spatial distribution of the location of behavior changes appears to vary, no pattern of statistical significance can be concluded from these data (Figure 1). Because potential hydraulic thresholds are created at increasingly downstream positions along the guide wall with increasing angle, it was hypothesized that the average location of behavior changes similarly appear farther downstream with increasing angle. Although such a pattern appears to exist, limited observations withhold the ability to draw statistically significant conclusions from purely the spatial distribution of behavior changes (Figure 2).

Figure 1. Two displays of behavior changes and the 95% confidence intervals surrounding their spatial distributions. The data in the graphic on the left, presented in Exercise 2, was displayed backwards.
Figure 2. Although average distance downstream of behavior changes varies with guide wall angle, no statistical significance was found in this study.

A comprehensive analysis of channel hydraulics and the location of behavior changes found turbulent kinetic energy (TKE) gradient, water speed, and velocity gradient to be potential predictors of the spatial distribution of behavior changes. A multivariate regression analysis relating the X and Y coordinates of a behavior change with water speed, TKE, TKE gradient, velocity gradient, and acceleration found only TKE to not be significantly correlated with the location of a behavior change (p-value > 0.05). A Principal Component Analysis (PCA) further determined that water speed, as the largest component of the first principal component, may be an important predictor of behavior change (Figure 3). However, PCA’s ability to only model one dependent variable (distance downstream) warranted analysis using Partial Least Squares regression (PLS). PLS found that velocity gradient and TKE gradient best predict the spatial distribution of behavior changes, while water speed contributed the least to each component axis (Figure 4).

Figure 3. Principal Component Analysis (PCA) indicated that water speed is an important hydraulic variable in predicting the downstream location of behavior changes.

Figure 4. Partial Least Squares (PLS) regression indicated that velocity gradient and TKE gradient are perhaps good predictors of the spatial distribution of behavior changes in both the X and Y directions.

Water speed, TKE gradient, and velocity gradient were found to vary spatially in their relationship with the location of behavior changes in a geographically weighted regression (GWR). GWR differs from other regressions by estimating local relationship coefficients for every point in the dataset, rather than assigning a global coefficient across the entire dataset. In this way, local variation in the relationship between independent and dependent variables can be seen. The relationship of TKE gradient and water speed with the location of behavior changes was both strongly negative and strongly positive at fine local scales (Figures 5 and 6). However, velocity gradient demonstrated a consistent positive relationship with the downstream location of a behavior change (Figure 7). That is, if a behavior change occurred at high (low) velocity gradients, we can expect its location to be relatively far downstream (upstream). Although variability in these data are a result of a fish’s past experience, physiology, and perception of stimuli, velocity gradient may most consistently predict the location of behavior changes in this experiment.

Figure 5. Geographically weighted regression depicted local variation in coefficients between water speed and the portion of the boom passed that span positive and negative values.

Figure 6. Spatial variability exists in the relationship between TKE gradient and the downstream location of behavior changes.

Figure 7. Velocity gradient is the only hydraulic variable that showed a consistent relationship with the distance downstream of a behavior change. This implied that if a behavior change occurred at high velocity gradients, we could expect its location to be relatively far downstream, and visa versa.


These results may verify previous experiments, which found velocity gradient to be an important predictor of the location of behavior changes in fish. If true, designs for improving fish passage at dams should aim to minimize velocity gradients created by floating guidance structures to promote safe bypass routes. However, substantial variation exists in these results (e.g. the spatial distribution of behavior changes does not vary significantly with guide wall angle, as hypothesized; no one hydraulic variable consistently predicted of the location of behavior changes across all 3 guide wall angles). Perhaps scientists should focus future studies on hydraulic stimuli as experiential (e.g. signal-to-background noise, cumulative to a threshold) rather than absolute (e.g. one magnitude of a hydraulic variable as a single threshold).

My learning:

For this course, I was forced to revisit R in order to run complicated statistical analyses. Although I was familiar with multiple linear regression analyses already, more sophisticated approaches (multivariate linear regressions, PCA, PLS, and GWR) were new to me. Furthermore, my ability to troubleshoot in Python in order to visualize the results of my analyses (creating multi-dimensional confidence intervals for Exercise 1, for example) was tested and grew.

April 2, 2018

Final Project Sample

Filed under: Final Project @ 2:40 pm

© 2019 GEOG 566   Powered by WordPress MU    Hosted by