# Appendix S1: R code, example analysis, tutorial and code for some figures # Prepared by Jason Hill on 3/25/18 in conjunction with a manuscript submitted by Jason Hill & Roz Renfrew to Ecology and Evolution # concerning the migratory movements of Grasshopper Sparrows (Ammodramus savannarum) as estimated from light-level # geolocators. This example will walk through the analysis for a single bird: U300, a male AHY Grasshopper Sparrow # originally captured and fitted with a light-level geolocator at Fort McCoy, WI on 6/9/2018. He was recaptured at the same location, (43.9711, -90.6642), # his territory, on 30 May, 2016. # I choose U300 for this example, because his geolcoator was still functioning when he returned to the breeding grounds in 2016, so we can # estimate the spring migration component as well. # Most of this tutorial was adapted from the many excellent examples from https://github.com/eldarrak/FLightR and from personal communication with # the developer of FLightR: Eldar Rakhimberdiev. Many thanks for this excellent R package. # Also check out the Geolocator Discussion and Support page on Ornithologyexchange.org # You can view several of my questions and answers from package developers there sessionInfo() #R version 3.4.3 (2017-11-30) #Platform: x86_64-w64-mingw32/x64 (64-bit) #Running under: Windows >= 8 x64 (build 9200) #Matrix products: default #locale: # [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 #[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C #[5] LC_TIME=English_United States.1252 #attached base packages: # [1] stats graphics grDevices utils datasets methods base #other attached packages: # [1] rgeos_0.3-26 geosphere_1.5-7 raster_2.6-7 sp_1.2-7 BAStag_0.1.3 FLightR_0.4.7 #loaded via a namespace (and not attached): # [1] Rcpp_0.12.16 pillar_1.2.1 compiler_3.4.3 plyr_1.8.4 tools_3.4.3 #[6] digest_0.6.15 memoise_1.1.0 tibble_1.4.2 gtable_0.2.0 lattice_0.20-35 #[11] png_0.1-7 rlang_0.2.0 mapproj_1.2-5 rgdal_1.2-18 proto_1.0.0 #[16] coda_0.19-1 withr_2.1.2 stringr_1.3.0 RgoogleMaps_1.4.1 maps_3.2.0 #[21] devtools_1.13.5 grid_3.4.3 jpeg_0.1-8 ggmap_2.6.1 ggplot2_2.2.1 #[26] reshape2_1.4.3 magrittr_1.5 scales_0.5.0 MASS_7.3-47 SGAT_0.1.3 #[31] colorspace_1.3-2 stringi_1.1.7 lazyeval_0.2.1 munsell_0.4.3 rjson_0.2.15 # clear your workspace rm(list=ls()) # set working directory to the folder containing this script setwd('c:/Users/Jason Hill/Documents/R/Geolocator Data.2.0/Appendix S3') ################################################################################# ##########Processing the raw data into twilight events########################### ################################################################################# # To work with FlightR, user must identify/define twilights using another package first # here, I use BAStag to define the twilight events and read the .lux file for U300 devtools::install_github("SWotherspoon/BAStag") library(BAStag) df.lux.U300<-readMTlux("U300_19Jul16_153711.lux") df.lux.U300$Light<-log(df.lux.U300$Light) # we will back-transform later, this just makes the light # values temporarily smaller/easier to work with # Define the capture latitude and longitude, FLightR will use this for calibration periods # This male sparrow was captured on his territory. lat.U300<-43.9711; lon.U300<--90.6642 lightImage(tagdata=df.lux.U300, offset = 7, zlim = c(0, 12), dt = 60) # offset simply adjust the Y-axis for ease of viewing # hopefully your data look this nice. Many of the other species' tags I have worked with have looked "not this nice" ############################## # process twilights, FlightR tutorial suggests a twlight threshold value of 1.0, whatever you choose, be consistent twl.U300 <- preprocessLight(tagdata = df.lux.U300, threshold=1.0, offset = 7,lmax=12,map=TRUE) # follow the BAStag documentation for using preprocessLight (it takes some getting used to) save(twl.U300, file="twl.U300.RData") # save this file for re-editing later # You only need to go through that process once. If you come back to this file later, you can reload the results # For re=editing, use the following lines, must write out file path name for some reason #load('c:/Users/...(replace w/ your file path here).../twl.U300.RData') #df.lux.U300<-readMTlux("U300_19Jul16_153711.lux") #df.lux.U300$Light<-log(df.lux.U300$Light) #twl.U300<-preprocessLight(tagdata = df.lux.U300, threshold=1.0, offset = 7,lmax=12,map=FALSE,stage=3,twilights=twl.U300) #save(twl.U300, file="twl.U300.RData") # don't forget to resave the file TAGS.U300.twl<-BAStag2TAGS(df.lux.U300, twl.U300, threshold=1.0) summary(TAGS.U300.twl) # need to backtransform lightlevels TAGS.U300.twl$light<-exp(TAGS.U300.twl$light) save(TAGS.U300.twl, file="TAGS.U300.twl.RData") write.csv(TAGS.U300.twl, file="TAGS.U300.twl.csv") ################################################################################# ##########Data are ready for FLightR############################################# ################################################################################# devtools::install_github('eldarrak/FLightR') library(FLightR) #load the data Proc.data.U300<-get.tags.data("TAGS.U300.twl.csv") summary(Proc.data.U300) str(Proc.data.U300) summary(Proc.data.U300$FLightR.data$Data$gmt) ############################## #User-defined calibration period for this bird/tag # this period represents the time when you belive the bird was stationary at one location # each geolocator records light with varying degrees of sensitivity, so you need to calibrate the light levels for each geolocator to a known location # In our case, that location is this male's territory. In my experience, longer calibration periods perform better in FLightR: # You are essentially giving FLightR more data to work with to calibrate that tag. Eldar Rakhimberdiev and I communicated about the 'appropriate' length # of a calibration period. According to my communications with Eldar, "short" calibration periods may result in non-desired effects. # For example, a week-long calibration period with our data resulted in a lot of due-north movements just prior to fall migration and # due-south movements just prior to spring migration. Eldar said this is indicitive of a calibration period that is "too short" # you will need to use your judgement with your own data plot_slopes_by_location(Proc.data=Proc.data.U300, location=c(lon.U300,lat.U300)) abline(v=as.POSIXct("2015-08-1"),lwd=2,col="DARKORANGE") # here you can see that the red and black + symbols entirely overlap through 1 August. That's good. # For a calibration period, we used day after banding through 1 August for all birds for simplicity # Now we define the calibration periods; you'll have at least one calibration period # But for this particular bird, we can also use a second calibaration period for when the bird was back on the breeding grounds in 2016 with a working geolocator # You won't know that until you run FLightR once for this bird Calibration.periods.U300<-data.frame(calibration.start=as.POSIXct(c(NA,"2016-5-4")),calibration.stop=as.POSIXct(c("2015-8-1",NA)),lon=lon.U300, lat=lat.U300) # 1st period: beginning of file through 1 August 2015 # 2nd period: from 4 May 2016 through end of file Calibration.U300<-make.calibration(Proc.data.U300, Calibration.periods.U300,model.ageing = TRUE) # note the model.ageing value. The geolocator sensor becomes opaque over time due to sun exposure, which may affect the amount of light it gathers # Less light gathering would trick FLightR into thinking the sunrise came later and sunset came earlier than it actually did. # FLightR can estimate the difference in light-gathering capabilities between the first and second calibration periods and linearly account for this undesired effect if you have >1 calibration period. # So the first time you run these data you would have model.ageing turned off (default value) because you'd only have one calibration period ############################## # Spatial extent of grid # here you can define your own spatial grid which restricts movement based on user-defined rules Smart.Grid<-make.grid(left=-110, bottom=0, right=-60, top=49, distance.from.land.allowed.to.use=c(-Inf, 1500),distance.from.land.allowed.to.stay=c(-Inf, 0),probability.of.staying = 0.05) # make sure you set the grid extent large enough to accomindate all likely movement # smaller grids run faster in FLightR, but make sure your results do not have data "stacked" on the edge of your grid; in that case make your grid larger # here we allow a small probability of the bird remainging stationary over water (see manuscript) # it's easy to bring the smart grid into ArcGIS (for example), for further manipulation (see manuscript) #write.csv(Smart.Grid,"Smart.grid.csv") # then edit the attribute table in ArcGIS (add XY data, convert to shapefile), #Then save or copy the entire attribute table into a .csv file for import back into R. # Note, your modified.smart.grid file should have as many rows as your original Smart.Grid file. # this (or something like it) is what Cooper et al. 2016 did with the Kirtland's Warbler paper # bring into R and overwrite original SmartGrid values #dfSG<-read.csv("modified.Smart.Grid.csv",header = T) #length(Smart.Grid[,3]); length(dfSG$Stay2) # they should be the same length, if not, you did something wrong. #Smart.Grid[,3]<-dfSG$Stay2 # overwrite the pixel values in the Smart.Grid with your customized values #plot(Smart.Grid,col=Smart.Grid[,3]) # future extension/idea: you could assign grid values based on habitat quality or availability within that cell ############################## #make pre.run.object ptm_start <- Sys.time() # to measure how fast your system runs all.in.smart.grid.U300<-make.prerun.object(Proc.data.U300, Smart.Grid,start=c(lon.U300,lat.U300), #end = NA, # this is tricky. read the help document carefully, for birds with 1 calibration period you would normally use "end=NA" which indicates that you don't know the location of the bird on the last date Calibration=Calibration.U300,threads=-1,Direction=180, Decision = 0.05) # examine the migration direction settings. str(all.in.smart.grid.U300$Indices$Matrix.Index.Table$Direction) # you see these values are all "180" [i.e., due south], that's a good guess in the fall, but not in the spring. # look at the "days since January 1" value that FLightR automatically creates all.in.smart.grid.U300$Indices$Matrix.Index.Table$yday # the first 411 values correspond with 2015, and after that is 2016 data # you can change all direction values (from 412 onwards to 0 [due north]) all.in.smart.grid.U300$Indices$Matrix.Index.Table$Direction[412:671]<-0 # the priors for movement decision ("decision=") and direction ("direction=") are twilight-specific # you could certainly use external information (e.g., band recoveries) to generate day-specific priors # for example, the movement decision prior could be low for most of the year, and then high during times when your species is known to migrate # Grasshopper Sparrows are stationary for most of the year's twilights, hence the low movement probability ptm_end<-Sys.time() - ptm_start;ptm_end ############################## # run the FLightR particle filter nParticles=1e4 # 1e4 is recommended for test and 1e6 for the analysis # results may change dramatically as you increase from 1e4 to 1e6 particles Result.U300<-run.particle.filter(all.in.smart.grid.U300,nParticles=nParticles,known.last=TRUE, #true b/c we know the location on the last day of data recording check.outliers=T,b=810) # see manuscript for how we choose the 810 km value write.csv(Result.U300$Results$Movement.results,"U300_MovementResults.csv") write.csv(Result.U300$Results$Quantiles,"GIS_Result.U300.csv") # bring the GIS_Result.U300.csv into GIS, for example in ArcGIS, add XY data and export to shapefile save(Result.U300, file="Results.U300.RData") ################################################################################################### ###########################Identifying stationary periods########################################## ################################################################################################### #Stationary.migration.Summary # here we determine the minimum length of what constitutes a stationary period # in our case we defined a stationary period (e.g., a stopover during migration) as two consecutive twilights in the same location # Define prob.cutoff value: the higher the prob.cutoff, the fewer stationary and movement periods you have. # so with 0.05 you might have 10 stationary periods (interestingly, some of these stationary locations might be directly atop each other, more on that later) # with 0.2 you might have 4, for example U300.migration.summary<-stationary.migration.summary(Result.U300, prob.cutoff = 0.05, min.stay = 2) save(U300.migration.summary, file="U300.migration.summary.RData") # note: we were primarily interested in migration stopovers (hence the filename) but this function returns stationary periods throughout the year # I found it helpful to put the $Potential_stat_periods and $Potential_movement_periods data side-by-side in an Excel worksheet while I worked through the rest of this code # note: the values in $Potential_stat_periods and $Potential_movement_periods refer to specific twilights ############################## #for mapping/summarizing purposes I found it easier to work with a reduced version of Result.U300 (which is massive) # create empty dataframe to hold selected results, we just want median lat/lon, quartiles, 95% CRIs, and time # might as well just adapt the original Results$Quantiles since they are the right dimensions (and time is correctly formatted already!) # but we don't need everything, just keep certain columns alt.Result.U300<-Result.U300$Results$Quantiles[c(2:3,5,9:10,12,17:21)] # for FLightR-identified stationary periods, use the lat/lon identified by the stationary.migration.summary function for (i in 1:nrow(U300.migration.summary$Potential_stat_periods)) { # for i in 1 to # of rows (i.e., # of stationary periods) alt.Result.U300$FstQu.lat[U300.migration.summary$Potential_stat_periods$start[i]:U300.migration.summary$Potential_stat_periods$end[i]]<-U300.migration.summary$Stationary.periods$FstQu.lat[i] alt.Result.U300$TrdQu.lat[U300.migration.summary$Potential_stat_periods$start[i]:U300.migration.summary$Potential_stat_periods$end[i]]<-U300.migration.summary$Stationary.periods$TrdQu.lat[i] alt.Result.U300$Medianlat[U300.migration.summary$Potential_stat_periods$start[i]:U300.migration.summary$Potential_stat_periods$end[i]]<-U300.migration.summary$Stationary.periods$Medianlat[i] alt.Result.U300$LCI.lat[U300.migration.summary$Potential_stat_periods$start[i]:U300.migration.summary$Potential_stat_periods$end[i]]<-U300.migration.summary$Stationary.periods$LCI.lat[i] alt.Result.U300$UCI.lat[U300.migration.summary$Potential_stat_periods$start[i]:U300.migration.summary$Potential_stat_periods$end[i]]<-U300.migration.summary$Stationary.periods$UCI.lat[i] alt.Result.U300$FstQu.lon[U300.migration.summary$Potential_stat_periods$start[i]:U300.migration.summary$Potential_stat_periods$end[i]]<-U300.migration.summary$Stationary.periods$FstQu.lon[i] alt.Result.U300$TrdQu.lon[U300.migration.summary$Potential_stat_periods$start[i]:U300.migration.summary$Potential_stat_periods$end[i]]<-U300.migration.summary$Stationary.periods$TrdQu.lon[i] alt.Result.U300$Medianlon[U300.migration.summary$Potential_stat_periods$start[i]:U300.migration.summary$Potential_stat_periods$end[i]]<-U300.migration.summary$Stationary.periods$Medianlon[i] alt.Result.U300$LCI.lon[U300.migration.summary$Potential_stat_periods$start[i]:U300.migration.summary$Potential_stat_periods$end[i]]<-U300.migration.summary$Stationary.periods$LCI.lon[i] alt.Result.U300$UCI.lon[U300.migration.summary$Potential_stat_periods$start[i]:U300.migration.summary$Potential_stat_periods$end[i]]<-U300.migration.summary$Stationary.periods$UCI.lon[i] } # now write out these alt.resuts to a csv file to bring into ArcGIS (add XY data, export to shapefile) write.csv(alt.Result.U300,"alt.GIS.Result.U300.csv") # use the GIS layer to visually ID when the bird started migration and reached the wintering and breeding grounds ################################################################################################### ###########################Define critical periods################################################# ################################################################################################### # Define critical periods for subsequent loops based on the GIS file and the $Potential_stat_periods and $Potential_movement_periods sf<-256# first movement twilight on fall migration ef<-307# first winter twilight (which is also end twilight of last fall movement period) sw<-509# end twilight of first stationary period on wintering grounds ew<-509# last twilight on wintering grounds ss<-510# first movement twilight on spring migration es<-657# end twilight of last spring movement (or stationary) period ################################################################################################### ###########################Fall Migration Summary################################################## ################################################################################################### # nrows = the # of twilights during fall migration, starting with fall movement and ending with first winter period twilight (U300.alt.fall.migration.latlons<-matrix(NA,nrow=(ef-sf)+1,ncol = 2)) # need to add 1 U300.alt.fall.migration.latlons[,1]<-alt.Result.U300$Medianlon[sf:ef] U300.alt.fall.migration.latlons[,2]<-alt.Result.U300$Medianlat[sf:ef] (U300.alt.fall.migration.path <- rbind(c(lon.U300,lat.U300), U300.alt.fall.migration.latlons)) # no need to bind winter location, redundant with last fall migration twilight library(sp) # to measure segments between locations and summarize them sum(spDists(U300.alt.fall.migration.path,segments=1,longlat=TRUE)) # total distance (km) moved between breeding and wintering grounds # we're assuming here that the bird traveled in a straight-line between twilights write.csv(spDists(U300.alt.fall.migration.path,segments=1,longlat=TRUE),"U300.alt.fall.distances.csv") library(raster) # to save fall migration route as a shapefile U300.alt.fall.migration.route<-spLines(U300.alt.fall.migration.path, crs='+proj=longlat +datum=WGS84') shapefile(U300.alt.fall.migration.route,"U300.alt.fall.route") ################ Fall segment summary # for U300, three stationary periods (#s 4-6) correspond with fall migration (U300.alt.fall.segments<-matrix(NA,nrow=(sw-sf)+2,ncol = 3)) U300.alt.fall.segments[1,1:2]<-c(lon.U300,lat.U300); U300.alt.fall.segments[1,3]<-1 # breeding location U300.alt.fall.segments[2:nrow(U300.alt.fall.segments),1]<-alt.Result.U300$Medianlon[sf:sw] U300.alt.fall.segments[2:nrow(U300.alt.fall.segments),2]<-alt.Result.U300$Medianlat[sf:sw] U300.alt.fall.segments[2:nrow(U300.alt.fall.segments),3]<-sf:sw head(U300.alt.fall.segments) # column 3 is the twilight number, which the following code will refer to # bind 1 [breeding grounds lat/lon] with reference to final (i.e., "end") twilights through the 1st winter stationary period fallends<-c(1,U300.migration.summary$Potential_stat_periods$end[4:7]) # in brackets refers to fall stationary periods and first winter stationary period U300.alt.fallsegs <- numeric(length = (length(fallends)-1)) for (i in 1:(length(fallends)-1)){ U300.alt.fallsegs[i]<-sum(spDists(U300.alt.fall.segments[which(U300.alt.fall.segments[,3]==fallends[i]):which(U300.alt.fall.segments[,3]==fallends[(i+1)]),1:2],segments=1,longlat=TRUE)) } # now summarize the movement segments on fall migration # the bird starts at the breeding grounds, and moves from statiory period to stationary period (i.e., stopover to stopver), to the wintering grounds # between stationary periods (i.e., when the bird is moving) U300 passes through all twilight locations that # were not identified as being stationary by FLightR format(round(U300.alt.fallsegs, 2), nsmall = 2) #results in km summary(U300.alt.fallsegs,digits=5) sd(U300.alt.fallsegs); sum(U300.alt.fallsegs) ################################################################################################### ###########################Winter Period Summary################################################### ################################################################################################### # You could easily complete this same process for the winter grounds if your bird moved around # U300 stayed in one stationary period the entire time before leaving for spring migration # so no movement to calculate here # given the uncertainty of winter locations that we reported in the manuscript, I'm not # sure that I'd have a lot of faith in estimates of winter movement ################# # When did bird reach the wintering grounds? # visually examine the GIS data & U300.migration.summary to determine # the first stationary period on wintering grounds, in this case the 7th period winter.lat.U300<-U300.migration.summary$Stationary.periods$Medianlat[7] winter.lon.U300<-U300.migration.summary$Stationary.periods$Medianlon[7] # estimate great-circle distance between wintering grounds and breeding grounds library(sp) # measure segments between locations and summarize them spDists(matrix(c(winter.lon.U300, winter.lat.U300,lon.U300,lat.U300),nrow=2,ncol = 2,byrow = T),segments=T,longlat=T) # this is the shortest possible fall migration flight path for U300 # this is a difference of sum(U300.alt.fallsegs)-spDists(matrix(c(winter.lon.U300, winter.lat.U300,lon.U300,lat.U300),nrow=2,ncol = 2,byrow = T),segments=T,longlat=T) # 299 km from what we actually estimated for U300 # estimate bearing from breeding to wintering grounds library(geosphere) bearing(c(lon.U300,lat.U300),c(winter.lon.U300, winter.lat.U300))+360 #180 is due south from breeding grounds ####################### # Winter space use, 50% and 90% utilization distributions # important to remember that FLightR includes the location uncertainty in these estimates (which is a good thing) # but that also means these space-use estimates are likely much larger that what the bird actually used U300.complete.winter.area<-plot_util_distr(Result.U300,dates=data.frame(as.POSIXct(Result.U300$Results$Quantiles$time[307]), as.POSIXct(Result.U300$Results$Quantiles$time[509])),percentiles=c(0.50,0.90)) writeOGR(U300.complete.winter.area$res_buffers, ".", "U300.complete.winter.area", driver="ESRI Shapefile") library(rgeos) # get centroid of 50 and 90% isopleths with rgeos, do not use the label point [less accurate] gCentroid(U300.complete.winter.area[[1]]) # -101.8832,19.93425 ################################################################################################### ###########################Spring Migration Summary################################################## ################################################################################################### # same process as fall migration summary (U300.alt.spring.migration.latlons<-matrix(NA,nrow=(es-ss)+1,ncol = 2)) # need to add 1 U300.alt.spring.migration.latlons[,1]<-alt.Result.U300$Medianlon[ss:es] U300.alt.spring.migration.latlons[,2]<-alt.Result.U300$Medianlat[ss:es] (U300.alt.spring.migration.path <- rbind(c(alt.Result.U300$Medianlon[ew],alt.Result.U300$Medianlat[ew]), U300.alt.spring.migration.latlons)) # bind last stationary winter location to spring movements sum(spDists(U300.alt.spring.migration.path,segments=1,longlat=TRUE)) # total distance moved between last statinoary period on wintering grounds and end of spring migration (whereever that may be) library(raster) # save spring migration route as shapefile U300.alt.spring.migration.route<-spLines(U300.alt.spring.migration.path, crs='+proj=longlat +datum=WGS84') shapefile(U300.alt.spring.migration.route,"U300.alt.spring.route") (U300.alt.spring.segments<-matrix(NA,nrow=(es-ss)+2,ncol = 3)) U300.alt.spring.segments[1,1:2]<-c(alt.Result.U300$Medianlon[ew],alt.Result.U300$Medianlat[ew]); U300.alt.spring.segments[1,3]<-ew # last wintering location U300.alt.spring.segments[2:nrow(U300.alt.spring.segments),1]<-alt.Result.U300$Medianlon[ss:es] U300.alt.spring.segments[2:nrow(U300.alt.spring.segments),2]<-alt.Result.U300$Medianlat[ss:es] U300.alt.spring.segments[2:nrow(U300.alt.spring.segments),3]<-ss:es head(U300.alt.spring.segments);tail(U300.alt.spring.segments) # column 3 is the twilight number, which the following code will refer to # bind ew [last wintering grounds lat/lon] with reference to final (i.e., "end") twilights through the 1st breeding grounds stationary period # if bird returned to breeding grounds with geolocator functioning, then include the first breeding # grounds stationary period in the next line of code (springends<-c(ew,U300.migration.summary$Potential_stat_periods$end[8:17],es)) # in brackets refers to stationary periods during spring migration U300.alt.springsegs <- numeric(length = (length(springends)-1)) for (i in 1:(length(springends)-1)){ U300.alt.springsegs[i]<-sum(spDists(U300.alt.spring.segments[which(U300.alt.spring.segments[,3]==springends[i]):which(U300.alt.spring.segments[,3]==springends[(i+1)]),1:2],segments=1,longlat=TRUE)) } format(round(U300.alt.springsegs, 2), nsmall = 2) # combine stationary periods <45 km apart (the default minimum movement parameter ["a"] in FLightR) # For U300, that translates to twilights twilights 547-588 & 602-622 alt.Result.U300[547:588,1:10]<-alt.Result.U300[547,1:10] # do not alter the time column alt.Result.U300[602:622,1:10]<-alt.Result.U300[622,1:10] # do not alter the time column # now revisit U300.alt.spring.migration.latlons (U300.alt.spring.migration.latlons<-matrix(NA,nrow=(es-ss)+1,ncol = 2)) # need to add 1 U300.alt.spring.migration.latlons[,1]<-alt.Result.U300$Medianlon[ss:es] U300.alt.spring.migration.latlons[,2]<-alt.Result.U300$Medianlat[ss:es] (U300.alt.spring.migration.path <- rbind(c(alt.Result.U300$Medianlon[ew],alt.Result.U300$Medianlat[ew]), U300.alt.spring.migration.latlons)) # bind last stationary winter location to spring movements sum(spDists(U300.alt.spring.migration.path,segments=1,longlat=TRUE)) # total distance moved between wintering grounds and end of spring migration write.csv(c(spDists(U300.alt.spring.migration.path,segments=1,longlat=TRUE),"test"),"U300.alt.spring.distances.csv") library(raster) # save spring migration route as shapefile U300.alt.spring.migration.route<-spLines(U300.alt.spring.migration.path, crs='+proj=longlat +datum=WGS84') shapefile(U300.alt.spring.migration.route,"U300.alt.spring.route") # revisit spring segment summary U300.alt.spring.segments[2:nrow(U300.alt.spring.segments),1]<-alt.Result.U300$Medianlon[ss:es] U300.alt.spring.segments[2:nrow(U300.alt.spring.segments),2]<-alt.Result.U300$Medianlat[ss:es] U300.alt.spring.segments[2:nrow(U300.alt.spring.segments),3]<-ss:es head(U300.alt.spring.segments);tail(U300.alt.spring.segments) # column 3 is the twilight number, which the following code will refer to # springends stays the same for (i in 1:(length(springends)-1)){ U300.alt.springsegs[i]<-sum(spDists(U300.alt.spring.segments[which(U300.alt.spring.segments[,3]==springends[i]):which(U300.alt.spring.segments[,3]==springends[(i+1)]),1:2],segments=1,longlat=TRUE)) } format(round(U300.alt.springsegs, 2), nsmall = 2) sum(U300.alt.springsegs) summary(U300.alt.springsegs[which(U300.alt.springsegs>0)],digits=4) # summary of segments >0 km sd(U300.alt.springsegs[which(U300.alt.springsegs>0)]) length(U300.alt.springsegs[which(U300.alt.springsegs>0)])-1 # of stopovers # now rewrite these alt.results to a csv file to bring into ArcGIS write.csv(alt.Result.U300,"alt.GIS.Result.U300.csv")