A worked example
Below I describe the most common use of the package, downloading recent PhenoCam time series and generating phenological transition dates for a given site and data type. I intentionally disable most automated processing and step through some of the routines to illustrate the workflow which normally happens internally to the main function download_phenocam(). Generated transition date files can be used in later analysis or modelling exercises using for example the phenor R package.
Select a site
A full list of meta-data for all sites can be queried using the list_sites() function.
sites <- list_sites()
head(sites)
#> site lat lon elev
#> 1 acadia 44.37694 -68.26083 158
#> 2 aguatibiaeast 33.62200 -116.86700 1086
#> 3 aguatibianorth 33.60222 -117.34368 1090
#> 4 ahwahnee 37.74670 -119.58160 1199
#> 5 alleypond 40.74284 -73.74304 61
#> 6 alligatorriver 35.78790 -75.90380 1
#> contact1
#> 1 Dee Morse <dee_morse AT nps DOT gov>
#> 2 Ann E Mebane <amebane AT fs DOT fed DOT us>
#> 3
#> 4 Dee Morse <dee_morse AT nps DOT gov>
#> 5 Mary Martin <mary DOT martin AT unh DOT edu>
#> 6 Asko Noormets <noormets AT tamu DOT edu>
#> contact2 date_start date_end nimage
#> 1 John Gross <John_Gross AT nps DOT gov> 2007-03-15 2020-03-20 48623
#> 2 Kristi Savig <KSavig AT air-resource DOT com> 2007-08-16 2019-01-25 36615
#> 3 2003-10-01 2006-10-25 2639
#> 4 John Gross <John_Gross AT nps DOT gov> 2008-08-28 2020-03-20 51302
#> 5 Nicholas Grant <ngrant02 AT fs DOT fed DOT us> 2014-11-05 2015-10-13 1452
#> 6 John King <john_king AT ncsu DOT edu> 2012-05-03 2020-03-20 54856
#> tzoffset active infrared method
#> 1 -5 TRUE N httppull
#> 2 -8 TRUE N httppull
#> 3 -8 FALSE N httppull
#> 4 -8 TRUE N httppull
#> 5 -5 FALSE N ftppush
#> 6 -5 TRUE Y ftppush
#> site_description
#> 1 Acadia National Park, McFarland Hill, near Bar Harbor, Maine
#> 2 Agua Tibia Wilderness, California
#> 3 Agua Tibia Wilderness, California
#> 4 Ahwahnee Meadow, Yosemite National Park, California
#> 5 Alley Pond, Queens, New York
#> 6 Alligator River National Wildlife Refuge, North Carolina
#> group camera_description camera_orientation site_type
#> 1 National Park Service unknown NE III
#> 2 USFS unknown SW III
#> 3 USFS unknown NE III
#> 4 National Park Service unknown E III
#> 5 SmartForests StarDot NetCam SC S II
#> 6 PhenoCam AMERIFLUX StarDot NetCam SC N I
#> flux_data flux_networks flux_sitenames MAT_site MAP_site MAT_daymet
#> 1 FALSE NA NA 7.05
#> 2 FALSE NA NA 15.75
#> 3 FALSE NA NA 16.00
#> 4 FALSE NA NA 12.25
#> 5 FALSE NA NA 11.90
#> 6 TRUE AMERIFLUX US-NC4 16.6 1310 16.75
#> MAP_daymet MAT_worldclim MAP_worldclim
#> 1 1439 6.5 1303
#> 2 483 14.9 504
#> 3 489 13.8 729
#> 4 871 11.8 886
#> 5 1263 11.7 1109
#> 6 1371 16.4 1312
#> dominant_species
#> 1
#> 2
#> 3
#> 4
#> 5
#> 6 Nyssa sylvatica, Taxodium distichum, Nyssa aquatica, Acer rubrum
#> primary_veg_type secondary_veg_type koeppen_geiger ecoregion wwf_biome
#> 1 DB EN Dfb 8 4
#> 2 SH Csa 11 12
#> 3 SH Csa 11 12
#> 4 EN GR Csb 6 5
#> 5 DB Cfa 8 4
#> 6 DB WL Cfa 8 5
#> landcover_igbp
#> 1 5
#> 2 7
#> 3 7
#> 4 8
#> 5 13
#> 6 5
#> site_acknowledgements
#> 1 Camera images from Acadia National Park are provided courtesy of the National Park Service Air Resources Program.
#> 2 Camera images from Agua Tibia Wilderness are provided courtesy of the USDA Forest Service Air Resources Management Program.
#> 3 Camera images from Agua Tibia Wilderness are provided courtesy of the USDA Forest Service Air Resources Management Program.
#> 4 Camera images from Yosemite National Park are provided courtesy of the National Park Service Air Resources Program.
#> 5
#> 6 Research at the Alligator River flux site is supported by DOE NICCR (award 08-SC-NICCR-1072), DOE-TES (awards 11-DE-SC-0006700 and 7090112), USDA Forest Service (award 13-JV-11330110-081) and USDA-NIFA (award 2014-67003-22068).To select a site first download an overview meta-data table of all available sites together with their ROI id’s and vegetation type and a limited set of meta-data parameters.
rois <- list_rois()
head(rois)
#> site lat lon veg_type roi_id_number
#> 1 acadia 44.37694 -68.26083 DB 1000
#> 2 acadia 44.37694 -68.26083 DB 2000
#> 3 ahwahnee 37.74670 -119.58160 GR 1000
#> 4 ahwahnee 37.74670 -119.58160 GR 2000
#> 5 ahwahnee 37.74670 -119.58160 GR 3000
#> 6 alleypond 40.74284 -73.74304 UN 1000
#> description
#> 1 Deciduous trees in foreground center
#> 2 Mixed forest in foreground. Start new timeseries due to camera/FOV change.
#> 3 GR veg type in foreground. Multiple FOV shifts.
#> 4 GR veg type in foreground.
#> 5 GR veg type in foreground.
#> 6 understory plants
#> first_date last_date site_years missing_data_pct
#> 1 2007-03-15 2017-09-20 9.8 7
#> 2 2017-10-11 2020-03-20 2.2 8
#> 3 2008-08-29 2011-10-14 2.8 9
#> 4 2012-05-01 2015-07-01 3.2 0
#> 5 2015-07-28 2020-03-20 4.6 0
#> 6 2014-11-04 2015-10-14 0.7 24Download a PhenoCam time series.
The below code shows you how to download a PhenoCam time series for the “harvard” site, ROI (roi_id) 1 and a time step frequency of 3-days. In this case the default outlier detection and smoothing routines has been disabled and will be run separately in subsequent steps. In normal use these will be enabled by default. The default output directory is tempdir() but any directory can be specified for data management purposes. If default settings are maintained, outlier detection and smoothing will be performed automatically. If so desired phenology dates can be estimated in one pass. In the latter case new data will be written in the same directory as specified for downloading the time series data.
download_phenocam(site = "harvard$",
veg_type = "DB",
roi_id = "1000",
frequency = 3,
outlier_detection = FALSE,
smooth = FALSE,
out_dir = tempdir())
#> Downloading: harvard_DB_1000_3day.csvAfter downloading we read in the data from disk. The data has a header and is comma separated.
df <- read_phenocam(file.path(tempdir(),"harvard_DB_1000_3day.csv"))
print(str(df))
#> List of 10
#> $ site : chr "harvard"
#> $ veg_type : chr "DB"
#> $ roi_id : chr "1000"
#> $ frequency : chr "3day"
#> $ lat : num 42.5
#> $ lon : num -72.2
#> $ elev : num 340
#> $ solar_elev_min: num 10
#> $ header : Named chr [1:24] NA NA NA "harvard" ...
#> ..- attr(*, "names")= chr [1:24] "#" "# 3-day summary product time series for harvard" "#" "# Site" ...
#> $ data :'data.frame': 4549 obs. of 32 variables:
#> ..$ date : chr [1:4549] "2008-01-05" "2008-01-06" "2008-01-07" "2008-01-08" ...
#> ..$ year : int [1:4549] 2008 2008 2008 2008 2008 2008 2008 2008 2008 2008 ...
#> ..$ doy : int [1:4549] 5 6 7 8 9 10 11 12 13 14 ...
#> ..$ image_count : int [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ midday_filename : chr [1:4549] NA NA NA NA ...
#> ..$ midday_r : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ midday_g : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ midday_b : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ midday_gcc : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ midday_rcc : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ r_mean : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ r_std : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ g_mean : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ g_std : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ b_mean : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ b_std : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ gcc_mean : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ gcc_std : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ gcc_50 : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ gcc_75 : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ gcc_90 : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ rcc_mean : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ rcc_std : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ rcc_50 : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ rcc_75 : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ rcc_90 : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ max_solar_elev : num [1:4549] NA NA NA NA NA NA NA NA NA NA ...
#> ..$ snow_flag : logi [1:4549] NA NA NA NA NA NA ...
#> ..$ outlierflag_gcc_mean: logi [1:4549] NA NA NA NA NA NA ...
#> ..$ outlierflag_gcc_50 : logi [1:4549] NA NA NA NA NA NA ...
#> ..$ outlierflag_gcc_75 : logi [1:4549] NA NA NA NA NA NA ...
#> ..$ outlierflag_gcc_90 : logi [1:4549] NA NA NA NA NA NA ...
#> - attr(*, "class")= chr "phenocamr"
#> NULLExpanding 3-day data for maximum phenology resolution
The downloaded time series is of a 3-day resolution. However, to correctly evaluate the phenology on a daily time step the time series needs to be expanded to this one day time step. This can be achieved using the expand_phenocam() function.
Detect data outlier
After reading in the data as a data frame you can apply the outlier detection routine. This routine uses an iterative method to detect outlier values in the Gcc time series. This routine filters out most spurious values due contaminiation by snow, mist, rain or otherwise very bright events. Warnings are suppressed as the routine is iterative and might throw warnings if it does not converge on a solution. This has no implications for the routine and data returned.
Smoothing a time series
After detecting outliers you can smooth the data. This function uses an AIC based methodology to find the opitmal loess smoothing window. Warnings are suppressed as the routine uses an optimization in which certain parameter settings return warnings. This has no implications for the routine and data returned.
Calculating phenophases
Finally, if smoothed data is available you can calculate phenological transition dates. This routine uses a PELT changepoint detection based approach to find meaningful seasonal cycles in the data. By default start of growing season dates are returned. If the reverse parameter is set to TRUE the end of growing season dates are returned. Dates are formatted as unix time and will be provided for three default threshold values (10 / 25 / 50%) of the Gcc amplitude.
start_of_season <- transition_dates(df)
print(head(start_of_season))
#> transition_10 transition_25 transition_50 transition_10_lower_ci
#> 1 14001 14007 14014 13999
#> 2 14361 14366 14374 14357
#> 3 14719 14724 14730 14717
#> 4 15095 15101 15107 15091
#> 5 15453 15460 15467 15444
#> 6 15824 15828 15835 15822
#> transition_25_lower_ci transition_50_lower_ci transition_10_upper_ci
#> 1 14005 14012 14004
#> 2 14365 14372 14364
#> 3 14723 14729 14722
#> 4 15100 15105 15098
#> 5 15459 15465 15456
#> 6 15827 15833 15826
#> transition_25_upper_ci transition_50_upper_ci threshold_10 threshold_25
#> 1 14009 14015 0.37674 0.39191
#> 2 14368 14375 0.37729 0.38890
#> 3 14726 14731 0.37666 0.39202
#> 4 15103 15108 0.37967 0.39426
#> 5 15462 15468 0.38299 0.39800
#> 6 15831 15836 0.38186 0.39409
#> threshold_50 min_gcc max_gcc
#> 1 0.41965 0.36904 0.46329
#> 2 0.41811 0.36808 0.46200
#> 3 0.41718 0.36909 0.46479
#> 4 0.42197 0.37074 0.46826
#> 5 0.42706 0.37250 0.47762
#> 6 0.42516 0.37268 0.47110Alternatively you can use the phenophases() function which is a wrapper of the transition_dates() function. However, as it potentially writes data to disk it needs additional information such as the roi_id, site name etc. The phenophases() function is the function which generated the final data products in the Richardson et al. (2018) paper. If used internally the output will be formatted in unix time, when written to file the dates will be human readable in YYYY-MM-DD format. Both start and end of season estimates will be provided.
With the phenoogy dates calculated we can plot their respective locations on the smoothed time series. In this case the plot will show the 50% amplitude threshold values for both rising and falling parts of the 90th percentile Gcc curve, marked with green and brown vertical lines respectivelly.
plot(as.Date(df$data$date),
df$data$smooth_gcc_90,
type = "l",
xlab = "date",
ylab = "Gcc")
# rising "spring" greenup dates
abline(v = phenology_dates$rising$transition_50,
col = "green")
# falling "autumn" senescence dates
abline(v = phenology_dates$falling$transition_50,
col = "brown")