Estimating Carrying Capacity with Big Data (ECCwBD)

SUMMARY

Estimating Carrying Capacity with Big Data (ECCwBD) summarizes remotely sensed MODIS 500m datasets of gross primary productivity (GPP) and net photosynthesis (PsnNet) from Greenup through senescence to estimate available forage within a user specified area of interest (AOI) and timeframe (Running et al. 2004, Running et al. 2015, Running and Zhao 2021, Friedl et al. 2022, LP DACC 2019). ECCwBD offers worldwide insights into understanding large scale ecological carrying capacity dynamics via user designated environmental nutrient constraints and ungulate nutrient intake requirements. Environmental thresholds are based on user input options for average percent digestibility of preferred forage types and associated crude protein content (Stringham et al., 2003). Input options for ungulate daily intake rates (kg) and crude protein composition by weight within AOI for a specified number of ungulates, up to five ungulate species and three groups per species, can lead to strong inference via readily approximated growth and offtake estimates based on growing season MODIS total PsnNet (EPA 2019, Leeuw J., et al 2019, RMEF 2020, Running et al. 2004, Sierra Club 2010, Steenweg et al. 2016). 

Estimated offtake is spatially evenly distributed as average offtake per 500m pixel. More specifically, mean ungulate offtake is evenly dispersed within AOI on grassland or shrubland terrain slope < 15% (Steenweg et al. 2016). Forage abundance and quality is represented by percent cover of AOI rangeland per floral species multiplied by species' respective percent crude protein content. Resultant product is summed for all designated rangeland forage species. Ungulate carrying capacity is calculated by first multiplying the number of ungulates per age group by daily intake rate (kg) necessary to maintain fitness. Resultant product (kg/day) is multiplied by ungulate crude protein intake requirement (%) and then divided by rangeland crude protein content (%) of designated dominant, preferred forage types. Quotient is nutritional based daily intake (kg/AOI m^2/day), which is in turn multiplied by the number of days per year ungulates are present; result is total annual intake (kg/AOI m^2/time). Total available forage is then growing season PsnNet (kg*C/AOI m^2/time) minus total annual intake (kg/AOI m^2/time). Proportion of available forage to bison offtake, times the number of animals equals relative adjustment estimate. Growing season available forage less ungulate offtake needed to build fat reserves leaves remaining forage to sustain a given number of ungulates through winter conditions into the following spring Greenup. If the proportion of remaining forage to ungulate offtake is positive and within fitness thresholds, ungulates may be added and vice versa if negative. ECCwBD can further assist conservation initiatives via herd fitness and resource stability optimization by estimating ungulate stocking rates per forage utilization levels, e.g., 25%, 65%, 87% and a theoretical 100% utilization (Leeuw J., et al 2019, Steenweg et al. 2016). Ungulate interspecific and sex ratios required to facilitate species abundance, richness and genetic diversity are beyond this scope. 

In order to answer AOI specific concerns surrounding impaired ecological processes and possible ecological state shifts, we offers services to inform the interpretation of ECCwBD results. Our standard requires set and repeatable analytical approaches related to analyses outside of Earth Engine, e.g., measurable rangeland AOI specific ecosystem structure and function impacts under variable ungulate utilization levels.  Please see more about our services here. 

JUSTIFICATION

ECCwBD summarizes photosynthesis within any user specified area of interest (AOI) worldwide thanks to the Earth Engine Code Editor, Client Libraries & Catalog of MODIS remotely sensed data (Gorelick et al. 2017, NASA 2020, Running et al. 2004). Significant correlations and low standard deviations exist between growing season net photosynthesis (PsnNet) and net primary productivity (NPP) on Northern Hemisphere grasslands and shrublands above and within the Tropic of Cancer (Running et al. 2015, Running and Zhao 2021). This strong relationship also exists for Southern Hemisphere grasslands within the Tropic of Capricorn. Southern Hemisphere growing seasons above the Tropic of Capricorn overlap years where Greenup and senescence are in separate years, e.g.  Dec. 2021 and  Jan.-Mar. 2022 (Friedl et al. 2022, LP DACC). Survivingenious is working to solve this known issue for the Southern Hemisphere.

Regarding latitudinal gradients north of the Tropic of Capricorn, MODIS 500m net photosynthesis (PsnNet) summed between Greenup and senescence day is highly correlated with predicted PsnNet (see scatterplot below). When compared with net primary production (NPP) growing season PsnNet is slightly less, meaning it slightly underestimates carrying capacity when compared with the NPP based estimate. Though, differences in primary productivity between NPP and season PsnNet are in part attributable to late or early detection of actual Greenup and/or senescence day by the MODIS 8-day intervals. Primary productivity outside the growing season can also reflect photosynthesis by woody species or other non-primary forage types low in nutrients. This means the NPP based carrying capacity estimate is likely biased high in terms of available nutritious forage types. Growing season PsnNet does account for most of annual NPP and is made available shortly after acquisition making it more up to date than NPP. Actual PsnNet is verified here with GPP and with NPP where available (Running et al. 2004, Running et al. 2015, Running and Zhao 2021, LP DACC 2019). Here, prediction of the prior year's total growing season PsnNet per 500m pixel provides a basis for calculating change in productivity per variation in environmental predictor variables. A deprecated Generalized Linear Mixed Model (GLMM) is applied here to predict growing season PsnNet (Bolker et al. 2009). Predicted PsnNet is then verifiable with actual growing season PsnNet as soon as it becomes available, post senescence.

Selected model type for ECCwBD PsnNet prediction is the GLMM with a gaussian family distribution and identity link function.  A normal distribution in the response and, in this case, for most predictor variables means the GLMM is deprecated to a linear mixed model, though here we continue referencing GLMM to lesson confusion. A selected GLMM is fitted with a random effect and fixed predictor variables with Shiny interactive options to select from 23 models fitted to combinations of up to 5 variables. There are 13 variable predictors and 13 static predictors from which to select. In the example below, a fitted GLMM with a random effect from average growing season PsnNet 2001-2020 and additive fixed effects from mean MODIS snow cover Jan.-Greenup day, mean MODIS leaf area index Jan.-Greenup day, and a direct interaction term between latitude and longitude. In a GLMM an intercept is calculated per random effect group, which accounts for and is distinguished by within group variation. Where a local intercept per group is applied with generalized slope coefficients for predictor variables (Bolker et al. 2009). GLMM comparison is accomplished by referencing Akaike information criterion (AIC) as a measure of model residual deviance. This estimate of error is important when comparing models because an AIC closer to zero means a lower residual deviance from zero, i.e., less model error (Bolker et al. 2009). If we plan to select for the most significant model by referencing AIC, then we will statistically see greater correlation between actual PsnNet and predicted PsnNet. Keep in mind, spatial temporal variation leads to varying degrees of AOI predictive accuracy. Moreover, results may vary where different combinations of variables per selected model capture greater or lesser degrees of variation. Theoretically, there is less variation in localized vegetative community productivity and increased variation within larger extents, a phenomenon known as spatial auto-correlation (Cliff and Ord 1973). Where there is vegetative heterogeneity or extents of homogeneity across the landscape Shiny interactive GLMM fitting may help capture residual deviance. Therein, where prior year's GLMM predicted PsnNet is highly correlated with actual season PsnNet we can expect a statistically significant (> 0.5) correlation for the following, current year of interest.

Actual available forage is then season PsnNet multiplied by forage crude protein content and digestibility, minus estimated total ungulate intake per number of animals within AOI/all species/all age group/pixels/time on grassland & shrubland terrain slope < 15% (Steenweg et al. 2016, RMEF 2020, EPA 2019, Sierra Club 2010, Leeuw J., et al 2019). Carrying Capacity adjustment is net actual available forage divided by total ungulate intake multiplied by number of animals (see Instructions). Predicted season PsnNet can be used to estimate carrying capacity at the beginning of the growing season, whereas actual season PsnNet can be used at the end of the growing season with greater confidence. In this manner, accuracy and validity of predicted PsnNet is verifiable with actual season PsnNet for any rangeland north of the Tropic of Capricorn where MODIS productivity estimates and predictor variables are available (Running et al. 2004). Here ungulate offtake is used interchangeably with ungulate intake. In United States estimates, grassland and shrubland classes are from the USGS National Land Cover Database; classes 52 & 71 respectively (Yang et al. 2018). Terrain slope & aspect are derived from USGS 10m NED (USGS). 

Alternatively, grassland and shrubland NPP may be used to estimate the landscape ungulate carrying capacity. In years where NPP is unavailable, predicted NPP can be used to estimate net annual photosynthesis; where predicted NPP results from 2020 PsnNet to NPP regression coefficients. In order to establish a point of reference for NPP verification, the 2020 NPP (updated) is regressed against total season PsnNet to estimate slope and y-intercept regression coefficients and to verify one another. When season PsnNet is adjusted with regression coefficients from 2020 NPP to season PsnNet and regressed once more against NPP their correlation approaches 1.00 and covariance remains high for 2020. However, inter-annual variability in productivity poses an issue for years with different growing conditions and/or utilization levels. Inter-annual variability control is accomplished by first checking if regression of MODIS net primary productivity (NPP: MOD17A3HGF.061) over MODIS season net photosynthesis (PsnNet: MOD17A2H.006) shows a significant linear relationship (Running et al. 2015, Running and Zhao 2021, LP DACC 2019). If collinearity is high, then multiplying growing season PsnNet by regression coefficients from NPP over season PsnNet gives bases for annual NPP carrying capacity estimates for current year of interest. For example, if input panel date range is set to 2020 for NPP & PsnNet, ECCwBD can verify the 2020 coefficient adjusted PsnNet (predicted NPP). Therefore, given little change in cover-type, 2021 carrying capacity can be estimated with a strong degree of confidence using 2020 NPP to 2020 season PsnNet regression coefficients. Predicted NPP is auto generated during ECCwBD run-time.

To improve carrying capacity estimates under actual utilization levels, ECCwBD application works for a single ungulate species or combinations of bison and elk and/or cattle and/or sheep and/or other animals within AOI. Optional settings for 'Other animals' daily intake rate per age group along with the number of animals per age group make it possible to estimate carrying capacity for at risk species, such as the critically endangered black rhinoceros (Diceros bicornis). Both seasonal PsnNet and coefficient adjusted PsnNet less estimated total ungulate offtake can offer an ecologically informed decision support tool for adaptable management strategies (Leeuw J., et al 2019). 

STATISTICAL CONSIDERATIONS

ECCwBD geospatial summarization and statistical inference application can assist the scientific toolset for Big Data analyses with complete sampling of all 500m pixel centers, 5000 pixels or less, of all environmental layers used in these analyses (Gorelick et al. 2017, LP DACC 2019, NASA 2020, Running et al. 2004). For rangelands with greater than 5000 pixels, where the fraction of pixels sampled represents population values the overall standard deviation, standard error and margin of error are low for all layers (Barde and Barde 2012).  For instance, Yellowstone Park's northern range is approximately 999.58 km^2 of which 562.87 km^2 is on grassland and shrubland represented by approximately 3577 MODIS 500m pixels. All grassland and shrubland net photosynthesis (PsnNet) pixels for the northern range within Yellowstone are illustrated by the histogram below. Further, the scatterplot below shows actual 2022 growing season PsnNet over the predicted 2022 growing season PsnNet for Yellowstone's northern range. Where extents have more than 5000 pixels, using a fraction sampled underestimates carrying capacity in terms of actual AOI total growing season PsnNet. In such cases, inference is drawn from statistical measures of population sample data, inductive reasoning, and reference to ECCwBD Earth Engine application results panels where ungulate capacity estimates are based on all AOI pixels. 

Accuracy improves for AOIs with less than 5000 pixels where covariance is positive and correlation approaches 1.00; season PsnNet and predicted season PsnNet are increasing together with higher degrees of approximation near 1.00 (see Instructions, Schober et al. 2018). Included are 41 environmental predictor variables; mean MODIS cycle 1 Greenup and senescence 2001-2020, winter mean MODIS leaf area index, fraction of photosynthetically active radiation, snow cover and albedo from January 1st through Greenup, OpenLandMap 250m soil organic carbon content, CHIRPS 5566m  winter and growing season total precipitation, NLCD 30m national land cover data, USGS 10m national elevation, ALOS 30m global elevation data, and ERA5 11.1km level 2 winter mean soil temperature and snow depth, and ERA5 11.1km mean volumetric soil water content from January 1st through Greenup and for Greenup through senescence. (Hall et al. 2016, Funk et al. 2015,  LP DAAC, Muñoz 2019, Tomislav and Wheeler 2018, USGS, please see NOTICE). Also available are EE reduced (bicubic resampling) 500m elevation along with elevation derived slope and cosine aspect to account for topographical variation at the MODIS scale (Google 2023). Currently also available is the option to  sample up to 5000 pixels with user specified sample points. 

Run ECCwBD and click "View data" to generate a scatterplot of 500m point data with coordinates. Please see instructions for more on sampling and launching ECCwBD analyses application in R Shiny Server. Be aware there are edge effects where a 500m MODIS pixel overlaps at least two land cover classes, e.g., open water and/or evergreen forest cover types. Forage estimates are less influenced by edge effects for larger grassland/shrubland ranges and for homogenous regions where ranges are within grassland/shrubland cover types. These areas are less influenced by edge effects for the same reasons a random sample size of up to 5000 pixels is sufficient for large rangelands. 

Moreover, high- and low-end pixel values outside two standard deviations are offset by each other and by the estimated ungulate daily intake rate. For generalization purposes, ungulate daily intake is set to be consistent throughout the year since seasonal intake rates may differ due to variable energy expenditures and metabolic rates where terrain and snow conditions limit mobility and forage access (Steenweg et al. 2016, Belile Thesis). Belile is working to resolve this known issue concerning winter biotic and abiotic multi-covariate interactions and influences on actual forage access and availability (see Works in progress). ECCwBD input options for proportion of year ungulates are present and the option to set 'other animal' daily intake rates for all age groups can help lesson inter-annual and intra-annual residual error due to forage availability, access and utilization levels. By generalizing AOI productivity and utilization per 500m pixel residual error influences of possible small-scale heterogeneity in available forage caused by inconsistent cover types and utilization patterns is lessoned. Normal probability density distributions of natural phenomena across landscapes is theoretically adequately represented (Ross 2017). 

Please reference R Shiny exploratory data analyses tabs for verification of MODIS growing season PsnNet and other layers (LP DAAC 2019, NASA 2020). Currently, all analysis takes place with 500m scaling for all pixels for all layers, aside from OpenLandMap 250m soil organic carbon content, CHIRPS 5566m precipitation, NLCD 30m national land cover data, USGS 10m national elevation, ALOS 30m global elevation data, and ERA5 11.1km level 2 soil temperature, volumetric soil water and snow depth, (Funk et al. 2015, Muñoz 2019, Tomislav and Wheeler 2018, USGS, please see NOTICE). Also available are EE reduced (bicubic resampling) 500m elevation, slope and aspect to account for landscape variation at the MODIS scale (Google 2023).

ECCwBD development and dial-in was performed using Yellowstone National Park's northern range. First image below shows AOI sum productivity over time and distribution of growing season net photosynthesis for the northern range within Yellowstone National Park (LP DAAC). Second image shows 2022 carrying capacity estimates for Yellowstone's northern range. The eccwbd application algorithm was ran with default settings to generate these results. Right click and select 'Open image in new tab'.