Earth transformed: detailed mapping of global human modification from 1990 to 2017

Data on the extent, patterns, and trends of human land use are critically important to support global and national priorities for conservation and sustainable development. To inform these issues, we created a series of detailed global datasets for 1990, 2000, 2010, and 2015 to evaluate temporal and spatial trends of land use modification of terrestrial lands (excluding Antarctica). We found that the expansion of and increase in human modification between 1990 and 2015 resulted in 1.6 M km² of natural land lost. The percent change between 1990 and 2015 was 15.2 % or 0.6 % annually – about 178 km² daily or over 12 ha min⁻¹. Worrisomely, we found that the global rate of loss has increased over the past 25 years. The greatest loss of natural lands from 1990 to 2015 occurred in Oceania, Asia, and Europe, and the biomes with the greatest loss were mangroves, tropical and subtropical moist broadleaf forests, and tropical and subtropical dry broadleaf forests. We also created a contemporary (∼2017) estimate of human modification that included additional stressors and found that globally 14.6 % or 18.5 M km² (±0.0013) of lands have been modified – an area greater than Russia. Our novel datasets are detailed (0.09 km² resolution), temporal (1990–2015), recent (∼2017), comprehensive (11 change stressors, 14 current), robust (using an established framework and incorporating classification errors and parameter uncertainty), and strongly validated. We believe these datasets support an improved understanding of the profound transformation wrought by human activities and provide foundational data on the amount, patterns, and rates of landscape change to inform planning and decision-making for environmental mitigation, protection, and restoration.

The datasets generated from this work are available at https://doi.org/10.5281/zenodo.3963013 (Theobald et al., 2020).

1 Introduction

Humans have transformed the Earth in profound ways (Marsh, 1885; Jordan et al., 1990; Vitousek et al., 1997; Hurtt et al., 2006), contributing to global climate change (IPCC, 2019), causing global habitat loss and fragmentation, and contributing to declines in biodiversity and critical ecosystem services (IPBES, 2019). Addressing the consequences of rapid habitat loss and land use change is essential for the implementation of various international initiatives, including the Convention on Biological Diversity 2020 Aichi Biodiversity Targets, the United Nations 2030 Sustainable Development Goals (especially Goal 15; Secretariat of the Convention on Biological Diversity, 2010), the Bonn Challenge (Verdone and Seidl, 2017), and the Global Deal for Nature (Dinerstein et al., 2019). Foundational to addressing these goals is a firm understanding of the rates, trends, and extent of these land use changes. Efforts to date have focused on historical patterns (Klein Goldewijk et al., 2007, 2017; Ramankutty et al., 2008; Ellis, 2018) or have been limited due to the unavailability of contemporary, temporally comparable, and high-resolution (<1 km²) data (Venter et al., 2016; Geldmann et al., 2019a; Kennedy et al., 2019b).

Here we describe a new dataset that maps the degree of human modification of terrestrial ecosystems globally, for recent changes from 1990 to 2015 and for contemporary (circa 2017) conditions. We mapped human activities that directly or indirectly alter natural systems, which we call anthropogenic drivers of ecological stress or “stressors” (following Salafsky et al., 2008; Theobald, 2013). Similarly to other efforts (Sanderson et al., 2002; Theobald, 2010, 2013; Geldmann et al., 2014; Venter et al., 2016; Kennedy et al., 2019b), we augmented remotely sensed data with traditionally mapped cartographic features. This is because remotely sensed imagery has limitations for this application – especially prior to ∼2010 – because it can require human interpretation to classify adequately and can miss development features that are obstructed by vegetation canopy or are small or narrow features (e.g., towers, wind turbines, power lines).

We mapped the degree of human modification based on an established approach that has been applied nationally, internationally, and globally (Theobald, 2010, 2013; González-Abraham et al., 2015; Kennedy et al., 2019b). It uses an existing classification system (Salafsky et al., 2008) to (a) ensure parsimony, (b) distinguish two spatial components (area of use and intensity of use), (c) use a physically based measure that is needed to estimate change (Gardner and Urban, 2007), (d) incorporate spatial and classification uncertainty, and (e) combine multiple stressors into an overall measure that assumes additive relationships among stressors and addresses the correlation among variables (Theobald, 2010). The resulting quantitative estimate of human modification has values ranging from 0 to 1 that support robust landscape assessments (Schultz, 2001; Hajkowicz and Collins, 2007).

To understand temporal landscape change, we calculated the degree of human modification – denoted by H – for the years 1990, 2000, 2010, and 2015 using methods and datasets that minimize noise and bias. Second, we included additional stressors not incorporated previously, including disturbance of natural processes due to reservoirs, effects from air pollution, and human intrusion (Theobald, 2008). Third, we calculated human stressors using data that are of a resolution up to 2 orders of magnitude finer (0.09 vs. 1–86 km²) than past efforts (Ellis and Ramankutty, 2008; Geldmann et al., 2014, 2019a; Haddad et al., 2015; Venter et al., 2016; Kennedy et al., 2019b). This higher resolution reduces the loss of information of the spatial pattern within a pixel; better identifies rare features; facilitates the application of these data for species and ecological processes that often occur at a fine scale; and improves the utility and relevance of these products for policy makers, decision makers, and land use managers.

Calculating H as a real value across the full gradient of landscape change is valuable because it can be applied rigorously to a variety of questions (Theobald, 2010, 2013), including discerning the heterogeneity of human uses that are often lumped within broad classes like “urban”, capturing the extent and pattern of the agricultural lands typically occurring beyond urban centers and protected areas, and delineating areas of low modification, all of which are useful for conservation prioritization and planning efforts (Kennedy et al., 2019b, c). Here, we describe the technical methods and briefly report on results on the temporal trends and current spatial patterns of human modification across all terrestrial lands, continents, biomes, and ecoregions. Because conservation organizations often use these types of data to focus their activities on specific regions (e.g., Jantke et al., 2019), we provide rankings by biome and ecoregion (Dinerstein et al., 2017) and briefly compare our results to other available studies.

2 Methods

2.1 Overview

We calculated the degree of human modification using the Direct Threats Classification v2 (Salafsky et al., 2008; http://cmp-openstandards.org, last access: 28 December 2019), which defines a stressor as the proximate human activities or processes that have caused, are causing, or may cause impacts on biodiversity and ecosystems. Table 1 lists the specific stressors and data sources we included in our maps: urban and built-up, crop and pasture lands, livestock grazing, oil and gas production, mining and quarrying, power generation (renewable and nonrenewable), roads, railways, power lines and towers, logging and wood harvesting, human intrusion, reservoirs, and air pollution.

Table 1 Overview of stressors, datasets, spatial resolution, and years data were available and used in the maps of human modification. Stressor classification levels in parentheses correspond to those within the Direct Threats Classification v2 (Salafsky et al., 2008). Acronyms of source data are in bold in the “Source” column for reference throughout the paper. For each stressor, the years 1990–2015 are used for change analysis, and ∼2017 is a compilation of all stressors that represents “current” conditions with the median year of 2017.

^a Based on interpolation. ^b Used major roads (i.e.,  highways) for 2019.

To estimate temporal change in H from 1990 to 2015, we followed established criteria (Geldmann et al., 2014) and included 11 stressors for which we could obtain global data with a fine-grained resolution (≤1 km²) and that provided consistent and comparable repeated measurements, especially in regards to the data source, methods used, and appropriate time frame (Table 1). We included current major roads and railways as a static layer in the temporal maps because in most cases some form of road existed prior to our baseline year of 1990 (except for the relatively rare, though important, new highways constructed).

To estimate the current amount of H circa 2017 (median = 2017, min = 2012, max = 2019), we included three additional stressors, comprising grazing, oil and gas wells, and power lines. We note that we did not map stressors for invasive species or pathogens and genes, geologic events, or climate change. This was because suitable temporal global data were not available to capture stressors due to invasive species or pathogens and genes, the majority of geological events is not directly caused by humans, and climate change is better modeled as a separate process distinct from the effects of direct human activities and there is a plethora of research on this topic (Geldmann et al., 2014; Titeux et al., 2016).

Table 2 Estimates of the intensity value for each stressor. “Best” estimates were determined from Brown and Vivas (2005)¹, Theobald (2013)², Kennedy et al. (2019b)³, or expert judgment⁴ and are bracketed by a minimum and maximum range, following the lowest–highest–best-estimate elicitation procedure to reduce bias (McBride et al., 2012). Results presented here use the best estimate, while minimum and maximum estimates are used to specify the range of possible randomized intensity values in the uncertainty analysis.

^∗ Assumed width of roads and railways (meters) provided in parentheses. Use of roads is incorporated into estimates of human intrusion.
${}^{\ast \ast }$ Forest loss due to wildfire was not included because of the challenges in understanding human causation and suppression, especially over a global extent. Also, loss due to urbanization was not included in this stressor because it is incorporated directly into the built-up stressor.
${}^{\ast \ast \ast }$ Minimum value is half of best; maximum is twice the best.

For each stressor s we quantified the degree of human modification as

$$\begin{array}{}\text{(1)}& H_s=F_s\cdot p\left(C_s\right)\cdot I_s,\end{array}$$

where F_s is the proportion of a pixel occupied (i.e., the footprint) by stressor s, p(C_s) is the probability that a stressor occurs at a location to account for spatial and classification uncertainty, and I_s is the intensity. Importantly, F and I have a direct physical interpretation (Gardner and Urban, 2007), are well bounded, and range from 0 to 1 and values are a “real” data type. Consequently, H provides the basis for unambiguous interpretation to assess landscape change (Hajkowicz and Collins, 2007; Riitters et al., 2009). Specific formulas used to map raw stressor data as indicator layers are provided below. Table 2 details our estimates of intensity values for each stressor (modified from Theobald, 2013; Kennedy et al., 2019b), which is used to differentiate land uses that have varying impacts on terrestrial systems (e.g., grazing is less intensive than mining). Our intensity values were informed by standardized measures of the amount of nonrenewable energy required to maintain human activities (Brown and Vivas, 2005) and found to generally correlate with species' responses to land use where examined (Kennedy et al., 2019b).

We generated datasets that represent temporal changes between 1990 and 2015 and for current (∼2017) conditions by combining stressor layers using the fuzzy algebraic sum (Bonham-Carter, 1994; Malczewski, 1999; Theobald, 2013), which is calculated as

$$\begin{array}{}\text{(2)}& H=\mathrm{1}-{\prod }_{s=\mathrm{1}}^n\left(\mathrm{1}-H_s\right),\end{array}$$

where n is the number of stressors (s) included. Of critical importance, the fuzzy sum formula is an increasive function that calculates the cumulative effects of multiple stressors in a way that minimizes the bias associated with nonindependent stressors and assumes that multiple stressors accumulate (Theobald, 2010, 2013; Kennedy et al., 2019b). This differs substantially from simple additive calculations that are commonly used (Halpern et al., 2008; Halpern and Fujita, 2013; Venter et al., 2016) but assume that stressors are independent and results in a metric that is sensitive to the number of stressors included in the model (Malczewski, 1999).

We mapped human modification of all terrestrial lands (excluding Antarctica) and included lands inundated by reservoirs but excluded other rivers and lakes. An often overlooked but critical aspect to understanding human modification is how water is mapped, especially for the interface between land and coastlines, lakes, reservoirs, and large rivers. We mapped nonreservoir areas dominated by water (i.e., oceans, lakes, reservoirs, and rivers) by processing data on ocean from the European Space Agency's Climate Change Initiative program (ESA CCI; 0.15 km, circa 2000) and on surface waters using the Global Surface Water dataset (GSW; 30 m; Pekel et al., 2016). We identified inland water bodies (i.e., lakes, reservoirs, rivers, etc.) using ESA CCI nonocean pixels that were at least 1 km from the interior of the land–ocean interface. We identified interior water pixels using GSW with at least 75 % water occurrence from 1984 to 2019 and that were at least 0.0225 km² in area (to remove small lakes, ponds, and narrow streams). As a result, inland water bodies and the ocean–land interface are more distinct, more consistent, and better aligned.

We summarized our estimates of human modification across all terrestrial lands, biomes, and ecoregions (defined by Dinerstein et al., 2017) and here report median (H_median) and mean (H_mean) statistics. We summarized results of temporal trends using the mean annualized difference (H_mad), calculated as the mean values across each analytical unit (e.g., biomes, ecoregions) of the annualized difference assuming a linear trend (H_ad):

$$\begin{array}{}\text{(3)}& H_{\mathrm{ad}}=(H_u-H_t)/(u-t),\end{array}$$

where u and t are the years of the datasets (e.g., u=2015, t=1990) and u>t. When discussing trends between 1990 and 2015, we emphasize the mean statistic because it better captures locations where H values have changed (mostly increasing over time), partly due to land uses with high values (e.g., urbanization ∼0.8) that are not well represented in the median statistic. We calculated the increase in H, or conversely the amount of natural habitat loss, as the per-pixel value times the pixel area, summed across a given unit of analysis. This assumes that any increase in the level of human modification causes natural land loss regardless of the original H level. We also report the median statistic because, as is typical of spatial landscape data, the distribution of H values is skewed to the right. Finally, we compared our results of the H_mad to those calculated on the human footprint (HF; for 1993–2009; Venter et al., 2016) and the temporal human pressure index (THPI; for 1995–2010; Geldmann et al., 2019a).

2.2 Stressors mapped

2.2.1 Urban and built-up

To map built-up areas that are typically found in urban areas and dominated by residential, commercial, and industrial land uses, we used the most recent version of the Built-up Grid from the Global Human Settlements Layer dataset (GHSL R2018A; Pesaresi et al., 2015; Corbane et al., 2018). The degree of human modification that is contributed by built-up areas, H_bu, is

$$\begin{array}{}\text{(4)}& H_{\mathrm{bu}}=F_{\mathrm{bu}}\cdot p\left(C_{\mathrm{bu}}\right)\cdot I_{\mathrm{bu}},\end{array}$$

where F_bu measures the proportion of the area of a pixel classified as built-up, p(C_bu) applies the GHSL-reported confidence mask (for 2014) for locations of the built-up areas (for the target year; Pesaresi et al., 2015), and I_bu is the intensity factor specified in Table 2.

Table 3 Probability of a land cover type being classified as cropland or pasture, calculated using the producer's accuracy, which is how often features on the ground are classified or the probability that a certain pixel is classified as a given land cover class. Probabilities of being cropland or pasture cover type (C_cp) are adjusted based on patch size (A) for patches with A<1 km², where $p\left(C_{\mathrm{cp}}\right)=C_{\mathrm{cp}}\cdot A_{\mathrm{cp}}^{\mathrm{2}}$.

2.2.2 Agriculture

We mapped agriculture stressors by identifying land cover classes associated with crop and pastureland from ESA CCI land cover datasets (ESA CCI, 2017; Pérez-Hoyos et al., 2017; Li et al., 2018) available at 0.09 km² for 1992, 2000, 2010, and 2015. We merged the cropland and pastureland stressors because these two classes are combined in the ESA land cover data and they are challenging to distinguish even at higher resolution (∼30 m; Wickham et al., 2017). To incorporate classification errors associated with all cover classes, we multiplied the footprint F_cp=1.0 with the probability that a pixel with cover class C was found to be cropland or pasture, p(C_cp), by interpreting reported accuracy assessment results (ESA CCI, 2017, in Table 3). To reduce the effects of scattered pixels that have some probability of being mapped as cropland–pastureland (e.g., misclassified pixels of high-elevation tundra or alpine areas), we multiplied p(C_cp) by the proportion of lands estimated to be cropland from the Unified Cropland Layer (Waldner et al., 2016), v, so that

$$\begin{array}{}\text{(5)}& p(C_{\mathrm{cp}}{)}^{\prime }=p\left(C_{\mathrm{cp}}\right)\times v\end{array}$$

and also reduced the value of p(C_cp)^′ based on patch size A, assuming that accuracy declines rapidly with cropland–pastureland small patches (A<1 km²) using

$$\begin{array}{}\text{(6)}& p{\left(C_{\mathrm{cp}}\right)}^{\prime \prime }={\left(p{\left(C_{\mathrm{cp}}\right)}^{\prime }\right)}^{\mathrm{2}},A<\mathrm{1}.\end{array}$$

We then calculated H_cp as

$$\begin{array}{}\text{(7)}& H_{\mathrm{cp}}=F_{\mathrm{cp}}\cdot p{\left(C_{\mathrm{cp}}\right)}^{\prime \prime }\cdot I_{\mathrm{cp}}.\end{array}$$

We developed spatially explicit estimates of agricultural intensity based on land management, such as cropping and number of rotations, tilling, and cutting operations, because these activities typically vary geographically (van Asselen and Verburg, 2012; Kehoe et al., 2017). We followed existing methods (Chaudhary and Brooks, 2018) to estimate three intensities of agricultural land use – minimal, light, and intense – and then mapped them using cover types from the Global Land Systems v2 dataset (GLS; Kehoe et al., 2017) by estimating intensity values (I) for each of the agricultural intensity types (Table 2). Although GLS v2 represents conditions circa 2005, we incorporated temporal changes by weighting the proportions of agricultural lands from the time-varying ESA CCI land cover datasets.

To estimate the modification associated with the grazing of domestic livestock (H_au), we used Gridded Livestock of the World v3 (Robinson et al., 2014; Gilbert et al., 2018a–i) which maps the density of animals per square kilometer (G) for eight types of livestock (j): buffaloes, cattle, chickens, ducks, goats, horses, pigs, and sheep. To calculate the overall footprint of grazing (F_au), we summed the weighted densities by global averages of livestock unit (LU) coefficients (w_i=0.84, 0.67, 0.01, 0.01, 0.10, 0.84, 0.23, 0.10, listed, respectively, for each livestock species stated above). We used a lower threshold found at 10 % to remove values <1.0 LU km⁻² (similar to Jacobson et al., 2019) and 1000 LU km⁻² as an upper threshold because it is a common breakpoint between grazing and industrial feedlots (Gerber et al., 2010). We assumed (here and below unless otherwise provided) no uncertainty (p(C_au)=1.0), because we lacked explicit data to do so. We then log₁₀ transformed and max-normalized (Kennedy et al., 2019b) to obtain 0–1 values and calculated the mean H_au using a 10 km radius moving window to reduce the effects of the coarser-resolution pixels:

$$\begin{array}{}\text{(8)}& {\displaystyle }{F}_{\mathrm{au}}={\sum }_{j=\mathrm{1}}^{\mathrm{8}}{G}_j{w}_j,max\left(\mathrm{1000}\right),min\left(\mathrm{1}\right),\text{(9)}& {\displaystyle }{H}_{\mathrm{au}}=\left(\right(\log \left(F_{\mathrm{au}}+\mathrm{1}\right))/\log (\mathrm{1000})\cdot p(C_{\mathrm{au}})\cdot I_{\mathrm{au}}.\end{array}$$

2.2.3 Energy and extractive resources

To estimate stressors associated with extractive energy production, we mapped gas flares derived from nighttime lights using data from the Visible Infrared Imaging Radiometer Suite from the Suomi National Polar-orbiting Partnership (VIIRS; Elvidge et al., 2013). Roughly 90 % of gas flares occur at locations where oil and gas are extracted (Elvidge et al., 2016). We used point data processed specifically to identify gas flares in VIIRS for 2012 and 2013 (Elvidge et al., 2016). For each flare, we approximated a footprint of 0.057 km² per wellhead (Allred et al., 2015). It is common to approximate the footprint of points (and lines) using a simple “buffer”, which implicitly assumes no location error and no distance decay from the point of origin. Such a buffer approach essentially centers a cylinder on each data point, where volume (V) equals the approximate footprint and height (h) and a perfect certainty of 1.0. Here, however, we assumed some uncertainty in the location of the point and that the effects associated with a feature such as an oil or gas wellhead diminish with distance. That is, rather than use a cylinder with volume V (or similarly a simple uniform buffer away from linear features, e.g., power lines or roads), we used a conic-shaped kernel centered on the point to calculate the uncertainty p(C_og), where the height of the cone h=0.5 represents a conservative estimate of spatial accuracy (Theobald, 2013). We derived the cone radius D=0.329 km by setting V to the footprint of 0.057 km²:

$$\begin{array}{}\text{(10)}& D=\sqrt{\left(\mathrm{3}/h\right)V/\mathit{\pi }}.\end{array}$$

Thereby the uncertainty parameter for each point is calculated using

$$\begin{array}{}\text{(11)}& p\left(C_{\mathrm{og}}\right)=\mathrm{3}h/\mathit{\pi }{D}^{\mathrm{2}}.\end{array}$$

We assigned the value of p(C_og) that overlapped the center of each pixel, with max p(C_og)=1.0. Human modification was then calculated as

$$\begin{array}{}\text{(12)}& H_{\mathrm{og}}=F_{\mathrm{og}}\cdot p\left(C_{\mathrm{og}}\right)\cdot I_{\mathrm{og}}.\end{array}$$

2.2.4 Mines and quarries

To estimate modification due to mines and quarries, we derived locations represented as points from a global mining dataset (n=34 565; S&P, 2018; Valenta et al., 2019). We retained surface mines that were constructed, where construction had started, were in operation, were in the process of being commissioned, or had residual production (n=22 705). For the temporal change analysis, we removed locations that did not have a specified year of construction (n=3634). We calculated the mean disturbed area and associated infrastructure of a mine by intersecting mine point locations with 441 623 polygons that represent footprints of quarries and mines (OpenStreetMap, 2019). For four types of mines – coal, hard rock (bauxite, cobalt, copper, gold, iron ore, lead, manganese, molybdenum, nickel, phosphate, platinum, silver, tin, uranium oxide, and zinc), diamonds, and other (antimony, chromite, graphite, ilmenite, lanthanides, lithium, niobium, palladium, tantalum, and tungsten) – we estimated the mean area (a) to be 12.95 km² (n=647) for coal, 8.54 km² (n=860) for hard rock, 5.21 km² (n=39) for diamonds, and 3.40 km² (n=27) for other. Finally, following Eqs. (8) and (9), we calculated p(C_m) for each of the four mining types using D values of 4.973, 4.038, 2.548, and 3.154 km, respectively, and calculated H_m as

$$\begin{array}{}\text{(13)}& H_{\mathrm{m}}=F_{\mathrm{m}}\cdot p\left(C_{\mathrm{m}}\right)\cdot I_{\mathrm{m}}.\end{array}$$

2.2.5 Power plants

To estimate the effects of where energy is produced, we mapped the location of power plants represented as points (n=29 903; WRI, 2019). For the temporal change analysis, we removed locations that did not have a specified year of construction (n=16 288). We estimated p(C_pp) using a conic-shaped kernel (Eqs. 8 and 9) and h=0.5. We mapped both nonrenewable energy forms (H_ppn; coal, oil, natural gas) and renewable energy forms (H_ppr; geothermal, hydro, solar, wind), where we assumed F_pp=1 and calculated a single p(C_pp) for both nonrenewable and renewable energy sectors with D_pp=1224 m (following Theobald, 2013):

$$\begin{array}{}\text{(14)}& {\displaystyle }{H}_{\mathrm{ppn}}=F_{\mathrm{pp}}\cdot p\left(C_{\mathrm{ppn}}\right)\cdot I_{\mathrm{ppn}},\text{(15)}& {\displaystyle }{H}_{\mathrm{ppr}}=F_{\mathrm{pp}}\cdot p\left(C_{\mathrm{ppr}}\right)\cdot I_{\mathrm{ppr}}.\end{array}$$

2.2.6 Transportation and service corridors

For transportation, we mapped roads and railways using OpenStreetMap highway linear features (OpenStreetMap, 2019). We calculated the footprint for the following transportation types: major (motorway, primary, secondary, trunk, link), minor (residential, tertiary, tertiary link), and two-track roads and railways as

$$\begin{array}{}\text{(16)}& {\displaystyle }{F}_{\mathrm{rr}}={\sum }_{i=\mathrm{0}}^c\left(w/\mathit{\alpha }\right)\cdot \mathit{\mu },\text{(17)}& {\displaystyle }{H}_{\mathrm{rr}}=F_{\mathrm{rr}}\cdot p\left(C_{\mathrm{rr}}\right)\cdot I_{\mathrm{rr}},\end{array}$$

where w is the estimated width of a road of type i from Table 2, α is the pixel width (i.e., 300 m), and μ=0.79 to adjust for the fractal dimension of road lines crossing cells (Theobald, 2000) because road lines often cross pixels at random angles. If a divided highway is represented as two separate lines, then each is represented independently. Also, if a cell has two or more roadway types cross it (e.g., where a secondary road joins a highway), the fuzzy sum of H_rr for both roads is calculated. Note that use of roads is incorporated into the “human intrusion” stressor (described below).

To map the modification associated with aboveground power lines (H_pl), we used

$$\begin{array}{}\text{(18)}& H_{\mathrm{pl}}=F_{\mathrm{pl}}\cdot p\left(C_{\mathrm{pl}}\right)\cdot I_{\mathrm{pl}},\end{array}$$

where F_pl is calculated using a 500 m buffer (Theobald, 2013), p(C_pl) is calculated using h=0.5, and I_pl is the estimate of intensity.

To estimate a stressor associated with electrical infrastructure and energy use (H_nl), we mapped “nighttime lights” using the Defense Meteorological Satellite Program Operational Linescan System (DMSP OLS; Elvidge et al., 2001) “stable lights” dataset. We included this as a distinct stressor from the energy extraction stressor (oil and gas flares, discussed above) because gas flares are derived by finding anomalies (high values) in the images rather than from the stable lights product and the footprints associated with the flares are an extremely small fraction of the overall extent of energy infrastructure.

To maximize temporal consistency, we used the intercalibrated DMSP OLS dataset and extended their (Zhang et al., 2016; Li and Zhou et al., 2017) approach for 2013 (using a=1.01, b=0.00882, $c=-\mathrm{0.965}$; Zhang et al., 2016). DMSP OLS values, L, are expected to range from 0 to 63, but because max values differed yearly (ranging from 57.87 to 66.16), we normalized all images (1992–2013) to range from 0 to 1.0 using the max-adjusted value for each year (L^′). To reduce the effects of noise in the images in areas with low light and in high northerly latitudes, we removed nighttime light values when $L^{\prime }<\mathrm{0.077}$ – that is, we set values to null when they were below the 25th percentile of the global terrestrial distribution compared to the often-used noise threshold of L=5 (following Elvidge et al., 2001).

To adjust for interannual spatial-misalignment errors (Elvidge et al., 2013), we adjusted the normalized DMSP image for 2013 to align with the 2013 VIIRS product by identifying sharply contrasting and consistent signals at 10 locations (n=10) distributed across the continents. We then visually compared each of the images from 1992 to 2012 to the DMSP image for 2013 and shifted the images to align them (averaged shift in meters – x=359.5, y=476.2). To further reduce interannual variability, we averaged image values at each pixel using a 3-year “tail” and used a rank-ordered-centroid weighting (Roszkowska, 2013) such that the spatially aligned and temporally smoothed nightlight value Y for year t is

$$\begin{array}{}\text{(19)}& Y_t=\left(L_t^{\prime }\cdot \mathrm{0.62}\right)+\left(L_{t-\mathrm{1}}^{\prime }\cdot \mathrm{0.26}\right)+\left(L_{t-\mathrm{2}}^{\prime }\cdot \mathrm{0.12}\right).\end{array}$$

Finally, to reduce the blooming effects and to take advantage of the higher-quality VIIRS-based nightlights (i.e., higher spatial resolution, reduction in saturated pixels), we sharpened DMSP nightlight values y_t using the VIIRS brightness value y to be proportional to the ratio of the DMSP values:

$$\begin{array}{}\text{(20)}& Y_t^{\prime }=Y_t\cdot \left(L_t^{\prime }/L_{\mathrm{2013}}\right).\end{array}$$

We then transformed $Y_t^{\prime }$ following Kennedy et al. (2019b), capping values above 126.0 (the 99.5 percentile of global values):

$$\begin{array}{}\text{(21)}& H_{\mathrm{nl}}=\left({\log }_{\mathrm{10}}\left(\mathrm{1}+Y_t^{\prime }\right)/\mathrm{2.104})\right)\cdot p\left(C_{\mathrm{nl}}\right)\cdot I_{\mathrm{nl}}.\end{array}$$

2.2.7 Logging

To estimate stressors on forested lands, we used maps of forest loss (Curtis et al., 2018) associated with commodity-driven deforestation, shifting agriculture, and forestry. (Note that we excluded as stressors wildfire because of the challenges of attributing wildfires to human causation – especially over global extent – and urbanization because it is measured directly by the built-up stressor.) We then identified locations where forest was lost due to one of the three mapped stressors (using v1.6, updated to 2018; Hansen et al., 2013) prior to the year of our estimated human modification map and applied the intensity value associated with that stressor. Thus,

$$\begin{array}{}\text{(22)}& H_{\mathrm{fr}}=F_{\mathrm{fr}}\cdot p\left(C_{\mathrm{fr}}\right)\cdot I_{\mathrm{fr}},\end{array}$$

where F_fr is pixels of forested loss in a given year and I_fr is an estimate of intensity associated with the cause of forest loss.

2.2.8 Human intrusion

We estimated human intrusion (H_i) using a method that builds on and extends accessibility modeling (Nelson, 2008; Theobald, 2008, 2013; Theobald et al., 2010; Weiss et al., 2018; Nelson et al., 2019). Human intrusion (a.k.a. “use”; Theobald, 2008) uses central place theory (Alonso, 1960) and integrates human accessibility throughout a landscape from defined locations, typically along roads and rails, as well as to off-road areas from urban areas (Theobald et al., 2010; Esteves et al., 2011; Theobald, 2013; Larson et al., 2018).

Accessibility measured in travel time in minutes is calculated from each mapped settlement point j (e.g., cities, towns, villages) from GRUMP v1.01 and GPW v4 (CIESIN, 2017, 2018). This approach is much less sensitive to arbitrary thresholds of city or town size (e.g., 50 000 residents), often used due to computational constraints (e.g., Nelson, 2008; Weiss et al., 2018). Second, to estimate intrusion of people into adjacent areas from a given settlement, we estimated the number of people (using population estimates at settlement j) at a given location (X; population density is people per square kilometer) following the assumption that the human density halved with every 60 min traveled (Theobald, 2008, 2013). The resulting intrusion map for each settlement was then summed to account for typical overlaps of intrusion from nearby settlements. We assumed that there is a limit at very high population densities, and so we capped the maximum value of intrusion, X, at 1 000 000 and then max-normalized using a square-root transform:

$$\begin{array}{}\text{(23)}& {\displaystyle }{F}_{\mathrm{i}}=X^{\mathrm{0.5}}\cdot \mathrm{0.001},\text{(24)}& {\displaystyle }{H}_{\mathrm{i}}=F_{\mathrm{i}}\cdot p\left(C_{\mathrm{i}}\right)\cdot I_{\mathrm{i}}.\end{array}$$

Note that accessibility was calculated using estimates of travel time along roads and rails, as well as off-road through different features of the landscape, using established travel time factors (Tobler, 1991) and presuming walking off-trail or via boats on freshwater or along ocean shoreline (Nelson, 2008; Theobald et al., 2010; Weiss et al., 2018; Nelson et al., 2019). This included effects of international borders following Weiss et al. (2018), and accessibility to lands was calculated across oceans.

2.2.9 Natural systems modification

Dams and their associated reservoirs flood natural habitat and strongly impact the natural flow regimes of the adjacent rivers (Grill et al., 2019). We mapped the footprint of reservoirs F_r created from 6849 dams from the Global Reservoir and Dam database (GRanD v1.3; Lehner et al., 2011; http://globaldamwatch.org/grand/, last access: 28 December 2019).

$$\begin{array}{}\text{(25)}& H_{\mathrm{r}}=F_{\mathrm{r}}\cdot p\left(C_{\mathrm{r}}\right)\cdot I_{\mathrm{r}}\end{array}$$

Because there are some potential analyses that would benefit from treating all water bodies consistently, we provided an additional version with all water bodies masked out in the dataset.

2.2.10 Pollution

We estimated the stress of air pollution by using data on nitrogen oxides (NO_x) through time from the Emission Database for Global Atmospheric Research (EDGAR v4.3.2; Crippa et al., 2018). We selected NO_x because it is a strong contributor to acid rain and fog and tropospheric ozone and because atmospheric levels are predominantly from human sources (Delmas et al., 1997). We used the 99th percentile (46 750 Mt) as the maximum value and then max-normalized (F_nox) and adjusted using the intensity value I_nox:

$$\begin{array}{}\text{(26)}& H_{\mathrm{nox}}=F_{\mathrm{nox}}\cdot p\left(C_{\mathrm{nox}}\right)\cdot I_{\mathrm{nox}}.\end{array}$$

2.3 Uncertainty and validation analyses

To understand the uncertainty in our results associated with our estimated intensity values (Table 2), following Kennedy et al. (2019c), we recalculated H, where I_s was randomized (n=50) between the minimum and maximum intensity values for each stressor. We then calculated the per-pixel mean and standard deviation for the 50 randomizations at 1 km² resolution for computational efficiency and provided corresponding maps.

We also assessed the accuracy of our maps following validation procedures described in Kennedy et al. (2019a, b, c). Because historical ground truth human modification data in a comparable form are not widely available, we restricted our analysis to test the contemporary conditions of human modification (∼2017 map) that included all stressor layers. We used an independent validation dataset from Kennedy et al. (2019b) that quantified the degree of human modification from visual interpretation of high-resolution aerial or satellite imagery across the world. We selected plots using the Global Grid sampling design (Theobald, 2016), a spatially balanced and probability-based random sampling that was stratified on a five-class rural-to-urban gradient using stable nighttime lights 2013 imagery (Elvidge et al., 2001). Within each of the 1000 ∼1 km² plots, we selected 10 simple random locations to capture rare features and heterogeneity in land use and land cover (for a total of 10 000 subplots), which were separated by a minimum distance of 100 m. The spatially balanced nature of the design maximizes statistical information extracted from each plot, because it increases the number of samples in relatively rare areas that are likely of interest (in contrast to simple random sampling) – especially for urbanized and growing cities (Theobald, 2016).

2.4 Processing platform

We processed, modeled, and analyzed the spatial data using the Google Earth Engine platform (Gorelick et al., 2017). We calculated all distances and areas using geodesic algorithms in decimal degrees (EPSG 4326). We summarized areas and percentages after projecting the data to the Mollweide equal-area projection (WGS84) to simplify calculations. All datasets and maps conform to the Google Earth Engine terms of service. We used program R 3.6.1 (R Core Team, 2019) to generate Fig. 2.

Table 4 Summary of estimates of the degree of human modification (H) and the mean annualized difference between 5- or 10-year increments for which change over time can be calculated (1990, 2000, 2010, and 2015) and H values for the contemporary dataset (∼2017, all stressors). Mean annualized mean difference is calculated as the mean value across the continents of the difference in H values divided by the number of years (e.g., $H_{\mathrm{mad}}=(H_{\mathrm{2015}}-H_{\mathrm{1990}})/\mathrm{25}$).

Table 5 A comparison of the mean annualized difference in human modification values for changes from 1990 to 2015 (H; 1990–2015), human footprint (HF; 1993–2009; Venter et al., 2016), and the temporal human pressure index (THPI; 1995–2010; Geldmann et al., 2019b). Mean annualized mean difference is calculated as the mean value of the difference in H values divided by the number of years (e.g., $H_{\mathrm{mad}}=(H_{\mathrm{2015}}-H_{\mathrm{1990}})/\mathrm{25}$), for each continent. Note that Oceania extends below Papua New Guinea (excluding the country of Australia).

3 Results

Below we describe the temporal and spatial trends of human modification by continent (Table 4), biome (Table 5), and ecoregion (Fig. 2).

Figure 1 A comparison of the recent trends in human activities by ecoregion using the mean annualized difference estimated by (a) human modification (H; from 1990 to 2015), (b) human footprint (for 1993–2009; Venter et al., 2016), and (c) temporal human pressure index (for ∼1995–2010; Geldmann et al., 2019a). Note that interactive maps are available at https://davidtheobald8.users.earthengine.app/view/global-human-modification-change (last access: 29 August 2020).

Figure 2 Graphs of the ratio of natural land loss (2015:1990) and contemporary (∼2017) degree of human modification (denoted as HM) for each of the 14 biomes and their ecoregions, globally. Note that ecoregions with change ratios ≥3.0 are placed on the maximum x-axis value (3.0).

3.1 Changes from 1990 to 2015

The mean value of H for global terrestrial lands increased from 0.0822 in 1990 to 0.0946 in 2015, a percentage change of 15.04 % overall and 0.60 % annually (Table 4). This equates to 1.6 M km² of natural lands lost – 178 km² daily. Increases in human modification occurred globally and in both urban and rural locations. We found that the largest increases in the H_mad occurred in Oceania, followed by Asia and Europe. Australia had the lowest increase followed by North and South America (Table 4). The biomes that exhibited the greatest increases in modification were mangroves, tropical and subtropical moist broadleaf forests, and tropical and subtropical dry broadleaf forests, while the biomes with the smallest increases were tundra, boreal forests and taiga, and deserts and xeric shrublands. Maps of changes in the H_mad between 2015 and 1990 for each ecoregion are shown in Fig. 1a, relative to the HF (Fig. 1b) and THPI (Fig. 1c). Figure 2 shows the ratio of natural land loss between 1990 and 2015, for each ecoregion and grouped by biome, in the context of the contemporary extent of human modification. We found most ecoregions (n=814) had increased in human modification, while the few (n=32) that had decreased were concentrated in higher latitudes and in more remote areas. We also found that changes in the H_mad have increased over time, from 0.0004 to 0.0005 to 0.0006, during 1990–2000 to 2000–2010 to 2010–2015. The percent change has also increased over time from 0.51 % to 0.59 % to 0.68 %.

3.2 Contemporary extent

We found that about 19.1 M km² (±0.0013) of natural lands was lost by ∼2017 – about 14.6 % of land globally (Table 4). South America was the most transformed (28.7 %), followed by North America (16.8 %), while Australia (5.0 %) and Africa (10.7 %) were the least transformed. Broad-scale patterns of the extent of human modification in ∼2017 are shown in Fig. 3. Note that “natural lands lost” was calculated using the continuous value of H, rather than approximations based on classifying the distribution.

Figure 3 The degree of human modification for circa ∼2017: (a) globally; (b) in central America; (c) in Europe, and (d) in Oceania. Note that interactive maps are available at https://davidtheobald8.users.earthengine.app/view/global-human-modification-change (last access: 29 August 2020).

Terrestrial lands with very low levels of human modification (H<0.01; Kennedy et al., 2019a; Riggio et al., 2020) are concentrated in less productive and more remote areas in high latitudes and are dominated by inaccessible permanent rock and ice or are within tundra, boreal forests, desert regions, and to a lesser extent montane grasslands. Table 5 shows that the biomes with the highest levels of H in ∼2017 were temperate broadleaf and mixed forests (H=0.3744); tropical and subtropical dry broadleaf forests (H=0.3317); and Mediterranean forests, woodlands and scrub (H=0.2903). The five least modified biomes were tundra (mean H=0.0023), boreal forests and taiga (H=0.0213), deserts and xeric shrublands (H=0.0572), and montane grasslands and shrublands (H=0.0894).

Following thresholds from Kennedy et al. (2019b), we found that in ∼2017, 51.0 % of global lands had a mean value of H≤0.01 (i.e., very low human modification), 13.3 % had a mean of $\mathrm{0.01}<H\le \mathrm{0.1}$ (low), 21.0 % had a mean of $\mathrm{0.1}<H\le \mathrm{0.4}$ (i.e., moderate), 12.3 % had a mean value of $\mathrm{0.4}<H\le \mathrm{0.7}$ (high), and 2.4 % (3.2 M km², ±0.0003) had a mean of $\mathrm{0.7}<H\le \mathrm{1.0}$ (very high). By area, these results by class amount to the following: very low = 66.8 M km² (±0.0067), low = 17.4 M km² (±0.0017), moderate = 27.6 M km² (±0.0028), high = 16.1 M km² (±0.0016), and very high = 3.18 M km² (±0.0003). Also, we found that 17.6 % had a mean value of H<0.0001 (23.0 M km², ±0.0023) and 4.2 % had H<0.00001 (5.5 M km², ±0.0006).

Table 6 A summary of the data, methods, and results comparing the degree of human modification (HM; this paper), degree of human modification 1 km (HM1k; Kennedy et al., 2019a, b, c), human footprint (HF; Sanderson et al., 2002; Venter et al., 2016), and temporal human pressure index (THPI; Geldmann et al., 2019b). Also see discussion of comparison in Kennedy et al. (2019a, c), Venter et al. (2019), and Riggio et al. (2020). Data source acronyms are provided in Table 1.

NA: not available

3.3 Comparisons

We compared our work to earlier efforts (summarized in Table 6) to determine if overall trends and extents were generally consistent with similar priorities of biomes and ecoregions. Globally, the H_mad from 1990 to 2015 (t=1990, u=2015) was 0.0005, while for the HF and THPI it was higher (HF_mad=0.0006, THPI_mad=0.0008). Perhaps more important is that the variability of the mean annualized difference values in the HF and THPI was 2.3 and 3.2 times that of H. By continent, we found that the H_mad increased the most in Oceania, followed by Asia, Europe, Africa, South America, North America, and Australia. Continental ranks by the THPI followed H roughly, though the HF differed more substantially (Table 5). The H_mad increased for all continents, but the HF_mad showed declines in modification for Europe and South America, while the THPI_mad showed a decline for North America.

Table 7 Summary of results by biome, comparing trends using the mean annualized difference for the human modification (H_mad), human footprint (HF_mad; Venter et al., 2016), and the mean temporal human pressure index (THPI_mad; Geldmann et al., 2019) score. Also provided are estimates of the proportion of terrestrial lands modified as estimated from Kennedy et al. (2018, 2019b; H1k) and the HF (score was max-normalized to rescale to 0–1). The THPI dataset characterizes only change, and so estimates of the proportion of lands modified in 2010 could not be provided. Mean annualized mean difference is calculated as the mean value across the continents and globally of the difference in H values divided by the number of years.

We also found the ranking of biomes by mean annualized difference for the HF and THPI was fairly different from ranks developed from H values (Table 7). Of the three biomes with the largest increase for the H_mad, two were also identified by the HF (tropical and subtropical dry broadleaf forests and tropical and subtropical moist broadleaf forests) and none by the THPI. Of the five biomes with the largest increase for the H_mad, three were also identified by the HF and THPI. The biomes that had the greatest disagreement among the ranking of H, the HF, and the THPI were mangroves, tropical and subtropical coniferous forests, and tropical and subtropical dry broadleaf forests.

The biggest differences in rankings between the H and the HF were for temperate and broadleaf mixed forests (and see comparisons of H1k and the HF in Kennedy et al., 2019b, c; Riggio et al., 2020). The HF was estimated to result in a 12.3 % modification for an earlier date (∼2009; Venter et al., 2016) and is lower likely because fewer stressors were included, because of its additive combination method, and because of its strongly right-skewed distribution caused by max-value normalization. The ranks of the extent of modification by biomes, however, were very similar in H, H1k, and the HF. In general, H had intermediate modification levels compared to H1k and the HF, with H1k levels being slightly higher (difference between 0.00 min and 0.09 max and average difference of 0.05 by biome) and the HF being slightly lower (difference between 0.00 min and 0.13 max and average difference of 0.04 by biome; Table 7). Results for ecoregions shown in Fig. 1 are even more striking, as the mean annualized difference values for the HF and THPI were inconsistent with our results. Of the 814 ecoregions that had increases in the H_mad, a decrease in modification was found for 201 ecoregions for the HF_mad and 202 for the THPI_mad; and for the 32 ecoregions that were found to have decreases in the H_mad, an increase in modification was found for 20 in the HF and 22 in the THPI.

Finally, the global estimate for H1k was likely higher than H because H1k did not limit the livestock stressor at LU km⁻² <1.0, used a slightly higher value for the low-threshold on the electrical infrastructure and energy use stressor (i.e., nightlights), and reported results that incorporate uncertainty into estimates of intensity. Furthermore, global modification from farming was estimated at 37 % for 2000 (Ramankutty et al., 2008) compared to 14.6 % with H. The difference with our results is largely due to their (Ramankutty et al., 2008) mapping of the area land cover types without differentiating the intensity of the impact of those cover types (crop and pasture).

3.4 Uncertainty and validation analyses

We addressed uncertainty in our results by incorporating the parameter p(C_s) for every sector s to best quantify uncertainty in its spatial location and classification as detailed in Sect. 2.2; for example, we adjusted p(C_cp) by directly incorporating measured confusion among land cover types using the results from the accuracy assessment of the land cover dataset (from Eq. 4). Additionally, we incorporated uncertainty by calculating the global mean for each of the 50 randomizations, which across the 50 iterations was 0.1434 ($\text{SD}=\pm \mathrm{0.0076}$) and ranged from 0.1243 to 0.1612. Thus, the global mean of 0.1461 obtained using our best-estimate intensity values was in line with our uncertainty results. We also mapped the per-pixel variance (standard deviation) to examine the spatial pattern of uncertainty (Fig. 4). The locations of the highest levels of uncertainty tend to be in more highly developed landscapes.

Figure 4 A map of the uncertainty as a result of randomizing the intensity factors when calculating the degree of human modification for 2017, showing the per-pixel standard deviation of 50 randomized maps. The highest levels of uncertainty tend to be in more highly developed landscapes (minimum = 0.0, median = 0.0, mean = 0.009, standard deviation = 0.014, maximum = 0.186).

We found strong agreement between H for ∼2017 and our validation data (r=0.783), with an average root mean square error of 0.22 and a mean absolute error of 0.04, for the 926 ∼1 km² plots (9260 subplots). There were 726 plots within ±20 % agreement, while for 161 plots H was estimated to be higher (and for 39 plots lower) than our visual estimate from the validation data. Estimates of H were biased high, likely because the stressors for human intrusion and electrical infrastructure (based on nighttime lights) are not readily observable from the aerial imagery used to generate the validation data. Our results here are consistent with our earlier findings (Kennedy et al., 2019a, b, c).

4 Data availability

The datasets generated from this work are available at https://doi.org/10.5281/zenodo.3963013 (Theobald et al., 2020), which includes the land–water mask used to support subsequent analyses. Extracts of specific geographic areas can be obtained by contacting the authors. All other datasets used in our work are open-source data cited herein.

5 Discussion

5.1 Summary

We found rapid and increasing human modification of terrestrial systems, resulting in the loss of natural lands globally. Our findings foreshadow trends and patterns of increased human modification, assuming future trends in the next 25–30 years continue as they have recently. Thus, our study reinforces calls for stronger commitments to help reduce habitat loss and fragmentation (Kennedy et al., 2019b; Jacobson et al., 2019) – which should be considered in conjunction with current commitments to reduce CO₂ emissions through the Paris climate accord (Baruch-Mordo et al., 2019; Kiesecker et al., 2019). We believe that the comparisons of ecoregions and biomes offer valuable contextual information that provides initial guidance on conservation strategies that may be most appropriate (Kennedy et al., 2019b). Also, it is important to consider the feasibility of achieving ecoregional (Dinerstein et al., 2017) or ecosystem (Jantke et al., 2019) representation goals as well as additional stresses caused by climate change (Costanza and Terando, 2019). We emphasize that although global, continental, biome, and ecoregional summaries provide a general understanding of trends and patterns, the high resolution of H and its gradient nature support robust estimates of change in human modification within a country and within an ecoregion, which are essential for tracking progress toward international and national conservation commitments (Mace et al., 2018), especially when placed within a broader structured decision-making framework (Tulloch et al., 2015).

Our datasets of human modification provide the most granular, contemporary, comprehensive, high-quality, and robust data currently available to assess temporal and spatial trends of global human impacts on landscapes. Our work is grounded in a structured classification of stressors, uses an internally consistent model, evaluates uncertainty, and incorporates refinements to minimize the effects of scaling and classification errors. Our validation approach uses an independent and spatially balanced random-sample design to provide strong support for the quality of our findings (Kennedy et al., 2019a).

Our overarching goal in producing and publishing these datasets is to support detailed quantification of the rates and trends, as well as the current extent and pattern, of human modification and to understand the degree of human modification across the continuum from low (e.g., wilderness) to high (e.g., intense urban). Beyond the basic findings presented here, we believe there are many potential applications of these datasets, including examining temporal rates and trends of land modification in and around protected areas (e.g., Geldmann et al., 2019b), estimating fragmentation for ecoregions and biomes (Kennedy et al., 2019b; Jacobson et al., 2019), and evaluating conservation opportunities and risks (e.g., the conservation risk index; Hoekstra et al., 2005). We also note that the human modification approach allows, in a straightforward and logically consistent way, for the inclusion of additional stressors and higher-resolution datasets that may become available over time or may be available for specific, local areas.

5.2 Caveats

As with any model, we recognize there are limitations to our work. We did not include data for all human stressors, largely because of incomplete global coverage or coarse mapping units (Klein Goldewijk et al., 2007; Geldmann et al., 2014), an inability to discern human-induced versus natural disturbances, or uncertainty in the location and directionality of their impact (e.g., the impact of climate change on terrestrial systems; Geldmann et al., 2014). In particular and discussed in Kennedy et al. (2019b, c), changes to land cover due to ecological disturbance events, such as wildfires or flooding, are not included in our analysis because of the difficulty in separating natural from human-caused disturbances – yet, we recognize that the broad extent of wildfire in particular could have strong implications. We did not include climate data as a stressor in this product to keep our analysis manageable and tractable. Although we attempted to map each stressor comprehensively, we recognize that some datasets may have missing features, particularly for mine and oil/gas wells – though large mines and concentrated oil/gas fields have been mapped quite well. For more integrated analyses, our data product should be used in combination with datasets of impacts due to climate change (e.g., Parks et al., 2020).

Stressors that are particularly important to improve include effects of grazing (currently coarse data and very broad expanse), pasture land, invasive species, and climate change (especially wildfire and effects of sea-level rise), and we encourage future work to focus on developing appropriate datasets and approaches to include or better capture these stressors. Key datasets we believe should be improved include transportation networks, including logging roads (e.g., van Etten, 2019), that are comparable through time; livestock grazing, rangelands, croplands, timber plantations, and pasturelands and their intensity of use; resource extraction (especially mining footprints); and temporal trends in gas flares, utility-scale solar and wind installations (Dunnett et al., 2020), and electrical substations.