A global land-use data cube 1992–2020 based on the Human Appropriation of Net Primary Production

The LUIcube presents information on area and biomass flows for 32 land-use classes globally at 30 arcsecond resolution and in annual time steps from 1992 to 2020. It was constructed based on the datasets listed in Table 2 and advancements of methods described in^1,51.

Land-use areas

The land-use area dataset was constructed based on a closed-budget approach³³, assigning the total land area to 5 main land-use types, that are further differentiated into 32 land-use classes.

We primarily used the land cover product from the European Space Agency Climate Change Initiative (ESA CCI LC)^9,52 to allocate the land-use types to the grid by defining the correspondence between land cover and land use (Table 3) based on the suitability of certain land-cover classes for a given land use. For cropland and grazing land national area statistics are available from FAOstat³⁴. Here, we used these suitability classes to distribute the national totals sequentially, starting with the most suitable class. Non-productive areas (including snow) were mapped based on the potential productivity and other wilderness areas were derived from^35,36 (see below for details).

Steps 1 through 9 (as shown in Fig. 7) were conducted in sequence, filling up the total land area per grid cell.

Stepwise integration of data inputs for mapping land-use types.

First, we mapped built-up land based on urban areas from ESA CCI LC. Next, the national cropland areas from FAO were allocated to the cropland suitability classes 1 to 3 sequentially. Cropland was then further differentiated into 23 crop-specific LU classes and fallow land, mainly based on the Spatial Production Allocation Model²⁰. Non-productive wilderness was mapped based on the potential productivity, using a threshold of 20gC/m²/yr for NPP_pot³³. Next, wilderness core and periphery were assigned. As a 6^th step, the national totals for permanent meadows and pastures reported in FAOstat were distributed within grazing suitability classes 1 to 3 per country. Since the grazing suitability classes 1 and 2 overlap with the same ESA CCI LC classes as cropland suitability classes 1 and 2, the assigned cropland areas were subtracted before distributing the national FAO total for grazing land. If the remainder of the grazing suitability classes 1 to 3 was smaller than the FAO value, part of peripheral wilderness was reclassified as grazing land. In step 7, forestry was mapped based on the ESA CCI LC classes for closed forests. Remaining open wooded land (step 8), comprising open forests with tree cover below 40% and the remaining land area not accounted for in the preceding steps were categorized as other grazing land (additionally to the FAO grazing land total). These steps are described in detail below and a full list of all LU classes is presented in Table 6.

Land cover product

The compilation of the land-use dataset aimed to reassemble land-cover data available from remote-sensing based on information on how this land is used by humans. While several global land-cover products are available (e.g.^53,54), the ESA CCI land cover product was selected as it provides the longest annual time series starting in 1992. We also chose not to combine multiple land cover products to maintain simplicity and facilitate seamless integration with other ESA products. We defined the correspondence between land cover and land use based on the plant functional types of each land cover class (see Table 3). The ESA CCI LC product was resampled from the original resolution of 300 m to fractional cover at 30 arcseconds, which is approximately 1 km² at the equator, and consists of 9 grid cells from the original resolution. This resolution was chosen (1) to reduce computational requirements, (2) to match the resolution of the change detection module of the ESA CCI LC product⁹, and (3) because it was found to be an appropriate resolution to link HANPP and biodiversity⁵⁵.

Land area and coastlines

The total land area per 30 arcsecond grid cell was calculated by subtracting ESA CCI LC class “210-Water bodies” from the total area of the grid cell. This was the starting point for the distribution of the various land-use classes within each cell.

Built-up land

Built-up land was mapped using the ESA CCI LC class “190-Urban areas”. These areas are not completely sealed, but also include parks or gardens. Additionally, we allocated 5% of the cropland area in each grid cell to built-up land to account for dispersed settlements and infrastructure not captured at the ESA CCI LC resolution of 300 m.

Cropland

The aim of the LUIcube is to track changes in cropland area and patterns across time, warranting consistency with national cropland totals reported by FAO. We thus used year-specific cropland patterns derived from ESA CCI LC adjusted to match the reported national cropland areas from FAOstat.

Data on the total area used as cropland per country⁵⁶ in addition to information on crop-specific area harvested per year⁵⁷ were downloaded from FAOstat on 22^nd May of 2024.

We classified several ESA CCI LC classes as being potentially used as cropland (see Table 3), with decreasing suitability from class1 to class3 and distributed the national total given by FAO (CLarea_FAO) consecutively across these classes in following procedure, which in principle follows the one suggested by²:

If the national total of cropland suitability class1, derived from ESA CCI LC, exceeded CLarea_FAO, class1 cells were assigned to cropland proportionally, with the remainder per cell being available for grazing. If not, all class1 cells were used as cropland. If the remaining area of CLarea_FAO of the country fit into cropland suitability class2 cells, it was distributed in the direct proximity of class1 cells. If not, cropland suitability class3 was utilized in the same manner. The national cropland area was reduced if the sum of suitability class1, class2 and class3 could not cover the value reported by FAOstat.

To differentiate cropland into the production areas for specific crops or crop classes we used data from the Spatial Production Allocation Model (SPAM)^20,28 for 2000⁵⁸, 2005⁵⁹, 2010⁶⁰ and 2020⁶¹, given at a spatial resolution of 5 arcminutes. This dataset is constructed based on subnational statistical data on crop production and uses auxiliary information on crop prices, population density and biophysical suitability to distribute cropping systems across grid cells. Since the thematic resolution was increased over time, we created aggregated crop classes (Table 4) to work with a continuous time series.

For each crop, SPAM provides information on yield, harvested area, and physical area accounting for multicropping within one year. For each gird cell we interpolated linearly between years with data points and kept values constant before the first data point (i.e., 1992–2000). The data obtained this way was integrated with the total cropland area as follows:

The basis cropland map used in SPAM is a synergy map⁶² that combines multiple sources to arrive at a cropland map for circa 2010. While this approach might create a more accurate cropland map for this point in time, tracking changes in cropland area and patterns across time requires using year-specific cropland patterns, which we derived from ESA CCI LC and adjusted to match the reported cropland areas from FAO per country as described above. Hence, we adapted the SPAM layers for harvested area per crop to the LUIcube cropland map of each year. This was done by filling those cropland cells not covered by the SPAM synergy cropland map with the respective average value of harvested area per crop of the closest cells in a radius of 30 cells (~30 km) and omitting the SPAM information of cells not classified as cropland in the LUIcube.

We ensured consistency with the harvested areas reported by FAO through following procedure:

1. 1)

   The minimum fallow area per country was estimated by relating the total harvested area (CLharv_FAO), summed over the 22 crop classes c derived from SPAM, to the total cropland area at the national level (Clarea_FAO), and applying the resulting percentage to the cropland area of each cell to arrive at the area per cell that is available for harvesting (CLharv).

   $$FALLOWmi{n}_{FAO}=CLare{a}_{FAO}-\mathop{\sum }\limits_{c=1}^{22}AREAharveste{d}_{FAO,c}$$

   (1)

2. 2)

   The area available for harvesting (CLharv) was then distributed among the 22 crop classes based on the share of each crop’s physical area in the sum of the physical areas of all crops present per cell (S_area_c) according to the SPAM dataset (adjusted to the cropland pattern and extent of the given year). The resulting crop-specific cropland area per cell (CLharv_c) was then adjusted based on a factor derived on the national level to match the FAO harvested area given for this crop and country (f_{c, FAO}). If a crop in a country was not represented in SPAM, but was reported in FAO we assumed it to be cropped proportionally in all cropland cells of the country.

   $$CLhar{v}_{c,FAO}=(CLharv\ast S\_are{a}_{c})\ast {f}_{c,FAO}$$

   (2)

3. 3)

   Applying the multicropping index given in SPAM to CLharv_c,FAO, we then arrived at the physical area per crop. Due to multicropping, the physical area can be smaller than the harvested area, which would increase the initial fallow areas estimate from step 1.

4. 4)

   While food crops are not cultivated on this fallow area, we assigned the production of fodder crops, for which data is not reported in the FAO database, to 50% of this area. The other 50% were assumed to actually lay fallow.

The multicropping index in cropland cells not identified as such in SPAM was set to the median value of the country. We used the global median value of the respective crop class if no SPAM data for the given country was available.

Wilderness

The land-use type wilderness was defined as areas that are not used by humans and encompasses non-productive and productive areas.

Non-productive areas and snow-covered areas without land use were classified based on the modeling results for the potential NPP (NPP_pot) and allocated to cells where potential productivity was below 20 gC/m²/yr, following the approach by³³.

Information on productive areas not or only rarely used by humans was derived from combining data depicting the human footprint, i.e. human pressures on the environment for 1993⁶³ and annually from 2000 onwards³⁵ with data on intact forest landscapes in 2000, 2013, 2016 and 2020 identified by³⁶, downloaded from⁶⁴. Both datasets were resampled to the target resolution of 30 arcseconds in ArcGIS Pro 3.2.2, linearly interpolated between data points and kept constant before the first, and after the last data point.

The human footprint values range from zero to 50, based on information on the extent of infrastructure, population density and agriculture. To identify unused areas, we only considered areas with a score of zero, i.e. areas with no discernable human presence. Since two different data sources were used for 1993 (data from⁶³) and for 2000 onwards (data from^35,65) we established continuity from 1993 to 2000 and beyond by adapting the 1993 layer. We first created the average of 1993 and 2000 and only counted cells with a value below 1 as wilderness and further set cells to zero that also have a value of zero in 2000.

Through the combination of both datasets, we differentiated between core wilderness and peripheral wilderness. Core wilderness was assumed to be without land use, while peripheral wilderness might be subject to infrequent land uses such as extensive, sporadic wood fuel gathering, mushroom picking, wild-honey-collecting, etc. Within the forest zone extent provided by⁶⁴ core wilderness grid cells were defined as those, where both the human footprint was zero and an intact forest landscape was present. Outside the forest zone a human footprint score of zero was the only criterium for assigning core wilderness. We defined peripheral wilderness as those areas within the forest zone extent where only one of the two datasets indicated wilderness and also all areas remaining in cells where core wilderness covers only a fraction of the cell.

Grazing land, open wooded land and forestry

We classified several ESA CCI LC classes as being potentially used as permanent pastures and meadows (see Table 3), with decreasing suitability from class1 to class3 and distributed the national total given by FAO consecutively across these classes in the same procedure as for cropland. Since we assumed that some of these classes can potentially be used as both cropland and/or grazing land, the areas already assigned to cropland were subtracted from the areas available for grazing land in each step. Additionally, we re-classified peripheral wilderness areas to grazing land, in the rare cases where the FAO target value could not be met within suitability classes 1 to 3.

ESA CCI LC classes with tree cover between 15% and 40% were aggregated to “open wooded lands”. Following the approach that areas which are not unambiguously assigned to other land uses might also be used for grazing³³, we included these open wooded lands in the grazing land-use type, but also allowed wood extraction through wood fuel collection.

We hence differentiated grazing land into two classes: grazing class “notrees” consists of grazing areas with grazing suitability 1, which were identified as agricultural areas or grassland by ESA CCI LC and have (almost) no trees. The remaining grazing land is categorized as grazing class “owl” (open wooded land), where woody vegetation is present and can also be harvested and the human influence on the landscape is less than in grazing class “notrees”.

Land predominantly used for forestry was mapped based on ESA CCI classes of closed forests (>40% tree cover based on the plant functional types given in the ESA CCI LC product description). We differentiate coniferous and non-coniferous forests.

Land-use intensity: biomass harvest, harvest residues, and land-use induced changes in NPP

Based on the land-use area dataset all components of the framework of Human Appropriation of Net Primary Production (HANPP) were calculated using adapted methods from^1,37,51,66. This framework discerns three different types of NPP: potential NPP (NPP_pot), i.e., the NPP that would prevail in the absence of land use but with current climate; the NPP of the actually prevailing vegetation (NPP_act), and the NPP that remains intact in ecosystems after harvest (NPP_eco). The difference between NPP_pot and NPP_act is the human-induced change in NPP due to land conversions (HANPP_luc); the difference between NPP_act and NPP_eco is HANPP_harv, i.e., the biomass extracted or killed during harvest. HANPP is the sum of HANPP_luc and HANPP_harv, equals the difference between NPP_pot and NPP_eco, and integrates output intensity and system-level intensity. HANPP can also be expressed as fraction of NPP_pot, then we denote it as HANPP%.

NPP_pot of all land-use classes

NPP_pot, the potential productivity in the absence of human land use, is the reference point of the HANPP framework and was derived from the Dynamic Global Vegetation Model (DGVM) LPJ-GUESS⁶⁷ in a run that assumed no land use or land management but used historical climate forcing data from 1900–2015 as in the preceding global HANPP study by¹. We set negative NPP values to zero and calculated 5-year moving averages to reduce the influence of annual fluctuations.

The computational requirements of such a global model run currently limit the feasible spatial resolution to 30 arcminutes, which is a lot coarser than the 30 arcsecond resolution aimed for in the LUIcube. While previous studies have resolved this mismatch in resolution with simple bilinear interpolation¹, we applied an alternative approach that takes fine-scale local conditions and their influence on NPP into account. This approach downscales the coarse-scale DGVM results based on NPP patterns derived from the empirical MIAMI model⁶⁸ utilizing data on temperature and precipitation available at 30 arcseconds⁶⁹ while retaining the original NPP totals on the level of ecoregions. The approach is described in detail in⁷⁰.

HANPP on built-up land

On built-up land, which is dominated by unproductive, sealed areas, but also includes urban green spaces such as parks and gardens, HANPP_harv was estimated at 1/6^th of NPP_pot and HANPP_luc was estimated at 2/3^rd of NPP_pot¹.

HANPP on cropland

The national level of HANPP_harv per crop is calculated based on the crop production reported in FAOstat³⁴ (downloaded on 10^th June of 2024) and factors to account for losses and byproducts accruing during harvest^51,66. The spatial distribution of HANPP_harv within a country follows the pattern of production derived from the SPAM database. For cells not considered as cropland in the SPAM database, but classified as cropland in the LUIcube, we assumed the yield to be at the lower end of the range within a country as to not overemphasize the production in these grid cells, and set it to 80% of the 5^th percentile to ensure also a reduced yield in cases where one yield level is given for the whole country. If the SPAM database did not provide information on the production of a given crop in a country the global median yield of this crop was applied. A moving average was calculated to smooth out the edges of the 5-arcminute resolution of the SPAM layers.

NPP_act was then calculated by accounting for pre-harvest losses (as in¹) and compared to the NPP_pot on the respective cropland area per cell to arrive at HANPP_luc.

HANPP of industrial roundwood and wood fuel harvest

As described in^1,51 factors were used to derive HANPP_harv (the sum of all biomass felled, i.e. stem wood, bark, branches, foliage, understory, roots etc.) from the national level wood harvest statistics from FAOstat³⁴ (downloaded on 10^th June of 2024), which reports roundwood volumes (excluding bark). NPP_act of woodlands was assumed to be equal to NPP_pot, i.e. HANPP_luc is assumed to be zero^32,43.

The LUIcube offers a finer thematic resolution than previous publications on HANPP, differentiating between industrial roundwood (IR) and wood fuel (WF) harvest as well as between coniferous and non-coniferous harvest and introducing open wooded lands as parts of grazing land where also wood extraction takes place.

The procedure for allocating national data to the grid worked as follows: IR harvest was assumed to predominantly occur in closed forests designated for forestry, while WF harvest was assigned to closed forests and open wooded lands proportionally to NPP_act. As in previous publications, the pattern of harvest was assumed to follow NPP_act patterns and the limit for wood harvest was set to 70% of NPP_act^1,37.

First, harvest occurring in closed forests was allocated, starting with the allocation of IR harvest to the grid. In countries where the IR harvest exceeded the available NPP_act of the respective tree type (i.e., coniferous or non-coniferous, hereafter termed tree type1 for illustrative purposes), the WF harvest was moved to open wooded lands of the same tree type1. Any still exceeding IR harvest was allocated to the closed forests of the other tree type2 (Fig. 8a, step 1a), as the tree type classification is not assumed to be 100% accurate. If the closed forest of the other tree type2 could not provide enough NPP_act for IR and WF harvest, the assigned WF harvest was moved to open wooded lands of tree type2 first (step 2a). Then, the remaining exceeding IR harvest was added to the harvest allocated to open wooded lands of tree type1 (step 3a).

Allocation procedure for wood harvest, exemplified for coniferous tree type. The procedure starts with closed forests (a) and then proceeds to open wooded lands (b).

Second, the harvest in open wooded lands was allocated. If the assigned WF harvest exceeded the NPP_act available in open wooded lands of the respective tree type1, it was re-allocated to be harvested in closed forests of the same tree type1, to the extent NPP was still available for harvest (step 1b). Any remaining WF harvest was then allocated to closed forests of the other tree type2 (step 2b). If the available NPP_act in these classes could still not cover the harvest, as would be the case if steps 1-3a had been triggered, it was moved to open wooded lands of the other tree type2 (step 3b), and was reduced if this was not sufficient.

HANPP of livestock grazing

The national level livestock feed demand was calculated as a “grazing gap”, i.e. the difference between livestock feed demand and the sum of feed crops and crop residues used as feed, as described in^51,66. The feed demand of roughage consuming livestock was calculated based on the factors and linear regressions between average daily feed intake per head and average national milk yield and carcass weight from⁷¹. However, we here defined the feed demand of cows and cattle as the minimum between the estimate for milk cows and beef cattle, as opposed to previous estimates which used the maximum^1,37,51,72, to increase the consistency with other estimates described in⁷³. National amounts of feed supply from market feed crops were sourced from Supply Utilization Accounts provided by FAO³⁴. Supply Utilization Accounts of by-products from oil production (oil cakes, e.g., soybean cake) are currently not reported by FAO and hence were estimated: the production of oil was extrapolated from vegetable oil production via technical conversion factors. Total imports and exports of oilcakes for each country were estimated by aggregating the bilateral trade data of the FAO (prioritizing importer reported before exporter reported values, see⁷⁴). The supply of oilcakes was calculated as production plus imports minus exports and was entirely attributed to animal feed. The use of crop residues as feed was included as in⁵¹ but adjusted by introducing thresholds to ensure that minimum 30% of the total feed demand is covered by roughage (crop residues and grazing) and minimum 15% by grazing alone.

NPP_act on grazing lands was estimated by taking NPP_pot as a starting point and reducing it in the case of grazing-induced degradation⁷⁵, and in the case of productivity changes due to ecosystem conversions, i.e. the replacement of forests by pastures and meadows. The portion of NPP_act available for grazing was calculated based on data on the proportion of non-palatable woody and belowground biomass⁷⁶. Boosting of NPP due to fertilization was allowed in countries which applied more than 5% of their overall fertilizer use to grasslands each year⁷⁷ and only assumed to occur on the grazing land class “notrees”.

The grazing mapping algorithm from the national level to the grid was adapted in this study. While previously grazing HANPP_harv was assigned to the grid based on NPP_pot patterns and the assumption that grazing intensity is disproportionally higher on high-productivity grazing lands, this study integrates human population density and the grazing land classification as additional factors to achieve a clearer separation between intensively used pastures and extensively used, remote rangelands.

The influence of population density on livestock distribution was derived from the Gridded Livestock of the World dataset^78,79 and the Global Human Settlement Layer⁸⁰, by plotting ruminant density against population density and calculating coefficients of a 2^nd order polynomial function (Fig. 9).

Plot of grid-cell values of population density VS. ruminant density with derived best-fitting polynomial function of 2^nd order.

The multiplication of NPP available for grazing and the influence of population density (based on the derived polynomial function) was used to assign a potential for grazing to each grid cell. Within the national distribution of “grazing potentials” we assigned a minimum grazing intensity of 40% of NPP to the 10% of cells with the lowest grazing potential, a maximum grazing intensity of 80% of NPP (due to seasonal constraints of productivity⁷⁶) to the 15% of cells with the highest grazing potential, and linearly interpolated the grazing intensity in-between, creating the baseline grazing intensity gradient (see Fig. 10). If the biomass grazable under these intensity constraints was sufficient to cover the calculated national grazing demand, the factor derived from the ratio between biomass demand and availability (below 1) was applied to all grazing land cells to arrive at the actually grazed biomass per cell. If the grazing demand could not be covered with the baseline grazing intensity gradient, it was adapted: the minimum grazing intensity and the percentage of cells to which the maximum grazing intensity was allocated were increased at the same time until the grazing demand could be met or a minimum grazing intensity of 75% was reached and 85% of cells were assigned the maximum grazing intensity (see Fig. 10). This adaptation flattens the gradient along which grazing harvest is distributed among the grazing land cells within a country and increases the share of cells with a high grazing intensity. If the grazing demand could still not be met with this adaptation and no fertilization of grazing lands was indicated grazing harvest was reduced to the maximum level determined by these constraints.

Grazing intensity (GI) gradient used to distribute grazing harvest per country. The baseline GI gradient (in black) was adapted by step-wise increasing the minimum GI and the percentage of cells with maximum GI (in orange). Boosting of NPP available for grazing on grazing land class “notrees” was allowed if fertilization was indicated (in blue).

For countries where increases of NPP_act on grazing land due to fertilization were indicated based on⁷⁷, the grazing intensity gradient was adapted only to the extent to restrict the boosting of NPP required to meet the grazing demand to a factor below 2. NPP increases were only applied on grazing land class “notrees”. In extreme cases (e.g., India) when even the highest possible grazing intensity could not meet the calculated grazing demand grazable NPP on class1 grazing lands was boosted by a maximum factor of 3.

NPP_act on grazing lands of class “notrees” was subsequently increased by the additional NPP created through fertilization.