Our analysis had two major components: 1) compiling individual and cumulative development threats globally, and 2) locating and prioritizing where cumulative threats pose a risk to terrestrial natural habitats. We projected development threats for nine sectors on terrestrial lands: urban and agricultural expansion, fossil fuels (conventional oil and gas, unconventional oil and gas, and coal), renewable energy (solar, wind, and biofuels), and mining. Future resource development potentials for each sector were created from publicly available global datasets (see S1 Table) and relatively ranked based on either the amount of unexploited resources (i.e. for fossil fuels, renewables, and mining) or estimated future area expansion derived from past trends (i.e. for urban and agriculture).

Calculating individual and cumulative development threat

Sector development threat rankings were based on the locations of unexploited or potential resources necessary to support development and/or estimates of land predicted to be modified (S1 Table). More specifically, these relative threat assessments were derived from synthesizing fifty global datasets: urban and agricultural expansion (n = 3; refs [13,35,36]), fossil fuels (n = 30; refs [37–66]), renewables (n = 12; refs [67–79]), and mining (n = 5; refs [80–84]) and then for each sector aggregating values to a 50 km2 grid cell. This cell resolution was selected due to the varying scales of source data (S1 Table) and the flexibility this resolution provided for aggregation. Future threats from resource development potentials for each sector were relatively ranked from 1 to 100 across the globe, excluding Antarctica and any 50-km analysis grid cells with greater than 50% overlap with marine environments. We then summed the individual sectors scores to produce a cumulative global threat map (see section below on combining individual sector threats and Fig 2). We projected all spatial data to a Mollweide projection to minimize area distortion except for between-feature distance calculations in which we used the Two-Point Equidistance projection. Unless otherwise specified, we used ArcGIS v.10.2 with the Spatial Analyst Extension [85] to perform all spatial data development, procedures and analyses.

Area-ranked threat scores. We relatively ranked threat scores using an equal-area rank method [86,87]: (1) where C r is the ranking of the target cell value, C i is the count of all grid cells with values less than the target cell value, f i is the number of cells with the target cell value, and N is the total number of cells in the study extent. To ensure ranking consistency (i.e. top ranked cells all equal to 100) across sectors, we rescaled all equal-area rankings from 1 to 100. The area-ranking approach created uniform distributions of scores per sector, such that equal bin ranges represented an equal area on the globe, therefore allowing for similar distributions and equal weighting across threat sectors. Prior to our selecting the area-ranking approach we tried several normalization approaches (e. g., log, log-log, square-root, cubic and min-max scaling) [88], but these transformations failed to create normally distributed values from the generally right-skewed individual sector threat scores and caused some to have more weighting than others.

Urban expansion. We used published maps of urban expansion probabilities by 2030 [13] (S1 Table). Maps were developed by first forecasting an aggregate amount of urban expansion per defined global regions from probability density functions of projected GDP and urban population. Then the aggregate amount of expansion was spatially distributed using a spatially-explicit land-change model with covariates slope, distance to roads, population density, and land cover. Given that our analysis focused on future development expansion into existing natural areas, we excluded areas already classified as urban and calculated the mean urban expansion probabilities for each 50-km grid cell. We then area-ranked mean probabilities of urban expansion for an urban development threat score (Fig 6). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 6. Projected future development threat of urban expansion. Area-ranked threat scores based on mean probabilities of global urban expansion by 2030, after excluding current urban areas. https://doi.org/10.1371/journal.pone.0138334.g006

Agriculture expansion. We calculated agriculture expansion rates using a 2000–2011 time series of global cropland and pasture maps following the methods of Ramankutty et al. [35] (S1 Table). We then: 1) summed for each year and 5 arc-minute (approx. 10x10-km) grid cell the total agricultural area in cropland and pasture (hereafter ag), 2) calculated the yearly fraction of area in ag within each grid cell, and 3) linearly regressed the 12-year time series and used the slope parameter as cell-specific rate of ag expansion. To focus on areas of potential development, we limited our analyses to only those cells with positive rates (slopes), and then averaged the rates of expansion within a resampled 50-km rectangle (corresponding to our threat analysis scale), which in effect accounted for higher likelihood of expansion into neighboring cells. To estimate the fractional area of agriculture expansion by 2030 for each 10-km grid cell, we resampled averaged rates back to a 10-km resolution and multiplied the averaged expansion rates by 19 (representing 19 years from 2012–2030). In cases where the fractional area of ag expansion for 2030, current ag land (i.e. 2011), and urban areas [36] summed to be greater than one (i.e., greater than the entire grid cell), we adjusted the fractional area of ag expansion as the maximum potential land conversion in the cell by subtracting the fractional areas of current ag and urban areas from one. Finally, we calculated the mean fractional area at a resampled 50-km grid resolution and area-ranked this mean value (Fig 7). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 7. Projected future development threat of agricultural expansion. Area-ranked threat scores based on estimates of fractional amount of agricultural expansion by 2030 extrapolated from 2000–2011 cropland and pasture time series maps. https://doi.org/10.1371/journal.pone.0138334.g007

Conventional oil and gas. For conventional oil and gas, our analysis used undiscovered volumes produced by the U.S. Geological Survey (USGS) for those global, geologic provinces which either currently contribute or are estimated to contribute in the future to the world’s reserves [37,39] (S1 Table). We augmented these global USGS assessments with more detailed national-level assessments available for the U.S. [38] and Australia [40], which resulted in a total of 305 geologic provinces worldwide with undiscovered oil and gas volume estimates. From this total, we excluded provinces with zero undiscovered volume (n = 8), and those provinces identified as having only offshore development [41,42] or not having at least 50% of the province overlapping land (n = 55). For each of the remaining 242 provinces we calculated the million barrels of oil equivalent (MBOE) of undiscovered oil, liquid natural gas, and natural gas with a conversion factor of 6 MBOE per each billion cubic feet of natural gas. We summed these values to quantify development potential per province and assigned this total MBOE value to overlapping 50-km grid cells with 50% or more of the cell intersecting a province. We then area-ranked cells based on this total MBOE value (Fig 8). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 8. Projected future development threat of conventional oil and gas. Area-ranked threat scores based on province-level estimates of undiscovered million barrels of oil equivalent for oil, natural gas, and liquid natural gas resources. https://doi.org/10.1371/journal.pone.0138334.g008

Unconventional oil and gas. For unconventional oil and gas, we focused on resources found in shale and other sedimentary formations but did not include any coal-bed methane resource as this latter form of development was estimated separately (see below). We used global [43] and U.S. [38] assessments of technically recoverable unconventional oil and natural gas (S1 Table). For non-U.S. regions, we geo-referenced and digitized basin maps from the global assessment [43] and linked resource volumes listed in the assessment to each basin (n = 98). For the U.S., we relied on spatially available data on basin location and resource estimate [38] (n = 17). We combined both datasets and for each basin, and converted all technically recoverable oil, natural gas and liquid natural gas volume estimates to billion barrels of oil equivalent (BBOEs) and summed these values for a total basin-specific resource estimate. Finally, we converted these basins to a raster with a 50-km resolution grid and area-ranked cells based on the summed BBOE value (Fig 9). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 9. Projected future development threat of unconventional oil and gas. Area-ranked threat scores based on basin-level estimates of technically recoverable billion barrels of oil equivalent for unconventional oil, natural gas, and liquid natural gas resources. https://doi.org/10.1371/journal.pone.0138334.g009

Coal. For coal, we combined tabular data of 2008 coal reserve estimates (million short tons) at the country-level [59] (S1 Table) with spatial data identifying coal-bearing areas for 65 countries (S2 Table). Spatially explicit data for 39 countries were available [44–48] while for the remaining 26 countries we geo-referenced existing digital maps and digitized coal-bearing areas [49,51–58,61–66] (S2 Table). We intersected all coal-bearing areas with the 65 country boundaries, calculated for each area its proportion that contributed to the countries’ overall total, and assigned individual coal reserve values per area by multiplying this proportion times the total country coal reserves. For four of the five top coal-producing countries (U.S., China, Australia, and India), we were able to further refine reserve estimates with published local government estimates [49,50,57,58]; we used the same attribution procedures based on proportion of overlap of coal-bearing areas with state or province boundaries and multiplied by the reserve estimates. Additionally 29 countries had coal-bearing areas but did not have any country or local reserve estimate. Due to these coal-bearing areas having some potential development threat, we assigned each of these coal-bearing areas with the lowest reserve value for all calculated areas of one thousand short-tons. Finally, we converted these coal-bearing areas to a raster with a 50-km resolution grid and area-ranked cells based on the reserve estimates (Fig 10). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 10. Projected future development threat of coal. Area-ranked threat scores based on coal basin reserve estimates in million short tons attributed form country- and state-level coal reserve data. https://doi.org/10.1371/journal.pone.0138334.g010

Wind. We used three main characteristics to estimate wind power development: 1) wind resources, 2) land suitability based on accessibility and physical restriction, and 3) economic feasibility based on electricity demand and distribution. Following other utility-scaled wind siting analyses [89–91], we synthesized and scaled each characteristic separately on its likelihood to support wind development. For wind resources we used annual averaged wind speed measured as m/s at 80m above Earth’s surface [67] (S1 Table), and restricted the analyses to wind speeds ≥ 6.4 m/s identified as most feasible for utility-scaled development. We then min-max normalized these wind speeds from 0.01–1. To create an overall binary land suitability map, we excluded land cover categories of rock and ice, artificial areas, water and wetlands [69], urban areas [36], and slopes > 20 degrees [70] and restricted all remaining lands to be within 80 km of an existing roads [71]. This produced an initial 300-m resolution raster identifying suitability that we resampled to 900-m resolution (3x3 cells) summing the binary results. Using a conservative approach where only resampled cells with a value of 9 (i.e. fully developed) were classified as suitable for wind power development, we resampled this result to 1 km through a simple bilinear process setting only those suitable cells to a value of 1. To account for economic feasibility, we used proximity to demand centers and existing power plants based on the inverse Euclidean distances (i.e., smaller, straight-line distances result in higher feasibility) from large urban areas [72], defined as greater than 10,000 people, and current utility-scaled power producing locations [73] identified by power plants producing ≥ 5 MW (n = 15,782) and hydropower plants [74] (n = 1541). We created each distance raster using a Two-Point Equidistance projection, projected them to a Mollweide projection with a bilinear sampling technique, and then rescaled values from 0.001–1 with 1 being those cells nearest to the feature (i.e. large urban areas or power plants). These two distance raster datasets were combined by calculating the average for each cell. To account for possible government incentives and the proven ability to develop wind power, this average was then doubled for those cells falling within a country that already produces wind power [75,76]. For consistency with the other three factors, this final feasibility raster (with values ranging from 0.001 to 2) was then re-scaled back to 0.001 to 1. We combined the outputs layers of wind resources (5-km resolution), land suitability (1-km resolution), and economic feasibility (5-km resolution) by multiplying the three metrics into the final wind development threat score, maintaining a grid cell size of 1 km. We then resampled this product to a 50-km resolution grid summing all cell values from the smaller 1 km raster and area-ranked this summed value (Fig 11). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 11. Projected future development threat of utility-scale wind power. Area-ranked threat scores based on combined metric of wind resources (m/s), land suitability, and economic feasibility for wind power development. https://doi.org/10.1371/journal.pone.0138334.g011

Solar. For solar, we followed a similar approach to wind resources where we considered three main characteristics to estimate solar development: solar resources, land suitability, and economic feasibility. Utility-scaled solar power produces electricity using two main types of technologies: concentrating solar power (CSP) and photovoltaic (PV). Each technology is optimally implemented at different solar radiation levels measured as Global Horizontal Irradiance (GHI), where PV development is best implemented at GHI ≥ 182 Watts/m2 and CSP development is best implemented at GHI ≥ 217 W/m2. Therefore for solar resources, we used GHI data [68] (S1 Table) to create two solar resource grids, one for CSP and one for PV, where we included only cells with GHI values ≥ 217 and ≥182 W/m2, respectively. We then normalized values from 0.01–1 for each CSP and PV grid, summed the two resulting grids into one solar resource availability output, and normalized results again to a scale of 0.01–1. For the remaining procedures we followed the steps as described above for wind resources with two exceptions: slopes were classified as ≤ 3 degrees [91] and we doubled feasibility scores based on countries producing solar power [77,78]. Similarly to our wind threat, we multiplied the three development factors (i.e. solar resource, suitability and feasibility) to produce one solar development threat value, resampled this threat value to a 50-km resolution grid via summation, and area-ranked these summed values (Fig 12). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 12. Projected future development threat of utility-scale solar power. Area-ranked threat scores based on combined metric of solar resources (W/m2), land suitability, and economic feasibility for solar power development. https://doi.org/10.1371/journal.pone.0138334.g012

Biofuels. Using crop-specific data for yield and harvested area [79] we focused biofuel production analysis on six first-generation biofuel crops (maize, soybean, sugarcane, rapeseed, sunflower, and oil palm), which make up the vast majority of commercial biofuel production [92] and have mature commercial markets, well-understood technologies, and therefore the potential to accelerate indirect land use change [93]. We assessed development threat by combining the maps of fractional area of cropland expansion by 2030 as described above (see Agriculture expansion) with maps of potential biofuel production measured in gallons of gasoline equivalents (GGE). To derive the latter, we first defined 100 crop-specific climate bins based on temperature and precipitation. To capture a range of yields within each climate bin, each climate bin had 1% of the total harvested area for each crop. Within each bin, the maximum potential yield (tons/ha) was defined as the area-weighted 95th percentile yield (i.e. 95% of harvested area within that bin had a lower yield). This methodology is described in more detail in Licker et al. [94] and Mueller et al. [95]. Yields were converted to GGEs using defined values (S3 Table). We then estimated potential expansion of each biofuel crop by mapping the full extent of the 100 crop-specific climate bins in the previous step. The driest climate bin at each temperature range was removed from the analysis as these bins represent extreme growing conditions requiring very intensive irrigation (e.g., Sahara Desert, or interior Australia) and are less likely to be developed. We then generated a maximum potential GGE map (10-km resolution) by combining all six biofuel crops while maintaining the highest GGE value for grid cells where crops overlapped. For final biofuel development threat, we multiplied the potential GGE map by the fractional area of cropland expansion by 2030, resampled the result to a 50-km resolution grid summing potential GGE values, and area-ranked the summed values (Fig 13). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 13. Projected future development threat of first generation biofuels. Area-ranked threat scores based on values of maximum potential gallons of gasoline equivalent multiplied by fraction of agriculture expansion by 2030. https://doi.org/10.1371/journal.pone.0138334.g013

Mining. We combined three main sources (S1 Table) to identify unexploited mineral deposits and quantify mining development threat: 1) Global Mineral Resources Data System [80], 2) Global Minerals Deposits update of 2011 [81–83] and 3) World Geoscience Database [84]. To discern patterns of future potential development, we removed current or past mining locations and any duplicate locations of the same mineral, resulting in a global dataset of occurrence or prospect deposits (n = 116,594). We created a global map that summed the number of unexploited deposits within a 50-km2 cell grid. Due to sampling bias towards the U.S. where 74% of deposits occurred, we area-ranked the number of mining occurrences within the U.S. separately from non-U.S. regions and merged the resulting grids into a final mining development threat map (Fig 14). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 14. Projected future development threat of mining. Area-ranked threat scores based on number of minerals and geologic materials deposit occurrences and prospects. https://doi.org/10.1371/journal.pone.0138334.g014