Global maps of cropland extent and change show accelerated cropland expansion in the twenty-first century

Spatiotemporally consistent data on global cropland extent is essential for tracking progress towards sustainable food production. In the present study, we present an analysis of global cropland area change for the first two decades of the twenty-first century derived from satellite data time-series. We estimate that, in 2019, the cropland area was 1,244 Mha with a corresponding total annual net primary production (NPP) of 5.5 Pg C year−1. From 2003 to 2019, cropland area increased by 9% and cropland NPP by 25%, primarily due to agricultural expansion in Africa and South America. Global cropland expansion accelerated over the past two decades, with a near doubling of the annual expansion rate, most notably in Africa. Half of the new cropland area (49%) replaced natural vegetation and tree cover, indicating a conflict with the sustainability goal of protecting terrestrial ecosystems. From 2003 to 2019, global per-capita cropland area decreased by 10% due to population growth. However, the per-capita annual cropland NPP increased by 3.5% as a result of intensified agricultural land use. The presented global, high-resolution, cropland map time-series supports monitoring of natural land appropriation at the local, national and international levels. High-resolution satellite observations provide an accurate and cost-effective solution to monitoring national and global progress towards Sustainable Development Goals. With a new global cropland dataset, this study reports that during the first two decades of the twenty-first century, the global cropland area increased by 9%, whereas the per-capita cropland area decreased by 10% and the per-capita annual cropland net primary production increased by 3.5%.

G lobal population growth and increasing standards of living inevitably cause the expansion and intensification of global agricultural land use to fulfil growing demands for food, biofuel and other commodities [1][2][3] . In turn, agriculture expansion and intensification threaten ecosystem functioning and lead to species extinction through habitat loss and fragmentation [3][4][5][6] . The United Nations' 2030 Sustainable Development Goals (SDGs) call for balancing increasing agricultural production with maintenance of ecosystem services 7 . Implementation of SDGs to improve food security, protect freshwater and terrestrial ecosystems, and mitigate climate change requires national policies and international cooperation that are based on consistent, independent and timely data on agriculture extent and productivity 8,9 . Spatiotemporally consistent satellite observations provide the most accurate and cost-effective solution for global agricultural, land-use mapping and monitoring 10 . Satellite data have been shown to enable national and global agriculture mapping [11][12][13][14][15][16][17] . However, no globally consistent, multidecadal, cropland time-series data at locally relevant spatial resolutions (30 m per pixel) exist to date.
In the present study, we present a global cropland extent and change dataset that can contribute to monitoring national and global progress towards SDGs. We define cropland as land used for annual and perennial herbaceous crops for human consumption, forage (including hay) and biofuel. Perennial woody crops, permanent pastures and shifting cultivation are excluded from the definition. The fallow length is limited to 4 years for the cropland class. Our definition is largely consistent with the arable land category reported by the Food and Agriculture Organization (FAO) of the UN 18 . To create the cropland dataset, we utilized the consistently processed 30 m spatial resolution Landsat satellite data archive 19 from 2000 to 2019. The Landsat time-series data were transformed into multitemporal metrics that characterize land surface phenology. These metrics were used as independent variables for a machine-learning classification to map global cropland extent. The classification models were locally calibrated using extensive training data collected by visual interpretation of freely available, high-spatial-resolution satellite images. We used a probability sample, stratified based on the Landsat-based global cropland maps, to estimate cropland area and its associated uncertainty, and to analyse pathways of land-use conversion. Sample reference data were collected through visual interpretation of Landsat time-series data and higher-spatial-resolution satellite images. Cropland maps were integrated with the Moderate Resolution Imaging Spectroradiometer (MODIS)-derived annual net primary production (NPP) 20 as a proxy variable for analysing crop productivity. The analysis was performed in 4-year epochs (2000-2003, 2004-2007, 2008-2011, 2012-2015 and 2016-2019). We created one cropland map per epoch (five maps in total), hereafter referred to by the last year of the interval (for example, the 2019 map represents the 2016-2019 epoch).
America had the largest relative cropland gain (by 37.1 ± 8.7 Mha, or 49%). Australia and New Zealand, as well as south-west Asia, displayed moderate cropland expansion (<0% of the 2003 area). North America, Europe and north and south-east Asia featured small net cropland area change but pronounced gross cropland gain and loss, which balanced each other at the continental scale.  Table 2). Of that total, 11% represents dryland conversion through irrigation, mostly found in south-west and south-east Asia and North America. The largest proportions of natural vegetation conversion to croplands (excluding dryland irrigation) were found in Africa (79% of all gross cropland gain area), south-east Asia (61%) and South America (39%). Cropland expansion is a major factor of forest loss and wildland fragmentation 4,23,24 , which illustrates a conflict with SDG 15, specifically, the SDG's targets to halt deforestation and degradation of natural habitats 9 . The other half of cropland expansion (51%) was due to pasture conversion and recultivation of abandoned arable land. Nearly all cropland expansion in Australia, New Zealand, Europe and northern Asia was found within pastures and long fallows (with no crop cultivation for >4 years). In North and South America, cropland expansion through the conversion of pastures and long fallows was more common (75% and 61%, respectively) than through clearing of natural vegetation 24,25 .
Abandonment or conversion to other land uses affected 10% of the 2003 cropland area (115.5 ± 24.1 Mha). Of that area, 52% was either converted into pastures or abandoned ( Table 2); such conversions may be temporary and followed by crop recultivation years later. Industrial and residential construction and infrastructure development were the second largest driver of gross cropland loss, responsible for 16% of the total cropland area reduction. In south-east Asia, 35% of cropland reduction was due to urban sprawl. A portion (13%) of 2003 cropland was converted to permanent woody crops or aquaculture, with the highest proportion of such transitions in south-east Asia (28%). Flooding caused by surface water increase, water erosion and reservoir construction affected the cropland area on all continents (3% total reduction). The remaining 16% of cropland reduction represented tree plantations or restoration of natural vegetation after cropland abandonment.
Cropland dynamics on the continental and national scales. The global Landsat-based, cropland map time-series is complementary to the sample analysis in characterizing global area dynamics (Fig.  3). The sample analysis showed high accuracy of the global cropland maps with variability between regions and lower accuracies for change dynamics ( Table 3). The cropland map time-series allowed us to disaggregate change over time and conduct national-scale analyses. Global cropland expansion accelerated over the past two decades, with a near doubling of the annual expansion rate from 5.1 MHa per year to 9.0 Mha per year ( Table 4). The change in annual cropland expansion rates highlights differences between cropland establishment in Africa and South America. In Africa, cropland expansion accelerated from 2004-2007 to 2016-2019, with a more than twofold increase in annual expansion rates. In contrast, cropland expansion in South America decelerated by 2019, with an annual expansion rate reduced by almost half compared with the 2004-2007 interval.
At the national level, the USA had the largest cropland area by 2019, closely followed by India and China (Supplementary Table  4). The largest net cropland increases were found in Brazil (by 23 Fig. 3b). The differences between national cropland estimates for selected countries may be attributed to different factors. We suggest that, in Russia, where crop abandonment is widespread and not fully documented, the arable land is overestimated by the FAO. In Brazil and Paraguay, the Copernicus cropland fraction dataset shows almost twice the size of cropland area compared with our estimate. This overestimation is partly due to misclassification of pastures as croplands by the Copernicus dataset.  South America had the highest NPP increase (by 0.38 Pg C year −1 , or 88%) followed by Africa (by 0.29 Pg C year −1 , or 50%) ( Table 5). The per-capita annual cropland NPP also increased globally by 3.5%, balancing the per-capita cropland area reduction. Two processes contributed to the global cropland NPP increase, namely the increase in cropland area and the increase in crop primary production per unit area. We found that the mean NPP per unit area within stable croplands (croplands presented over the entire 2000-   through fallow land recultivation. The cross-boundary distribution of major cropland areas and synchronous cropland dynamics illustrate the importance of international cooperation to ensure global progress towards SDGs. Global cropland maps provide spatial context on national, cross-boundary and local cropland dynamics, reflecting the history of land tenure, national policies and abrupt events such as natural and man-made disasters (Extended Data Fig. 6). In eastern Europe, the Baltic states and Russia's Kaliningrad region featured cropland expansion through recultivation of long fallows abandoned after the breakdown of the USSR, whereas the cropland area in neighbouring Poland and Belarus was relatively stable. Cereal, forage and hay production land of the northern Great Plains has different dynamics within Canada, where we observed land abandonment or conversion to permanent pastures, and the USA, where land management has been intensified. The irrigated croplands in Saudi Arabia declined after the depletion of groundwater resources and the implementation of state policies to discourage water-intensive crop production 27 . The 30-m spatial resolution of the cropland maps supports the analysis of local dynamic factors, for example, cropland abandonment after radioactive contamination following the 2011 nuclear disaster on the Fukushima Daiichi nuclear power plant in Japan.
Changes in total and per-capita mapped cropland area from 2003 to 2019 demonstrate the variability of national responses to the need for increased food production to feed a growing population (Extended Data Fig. 7). For most countries with moderate cropland area gains, we observed small decreases in per-capita cropland area. In many African nations (for example, Cameroon, Chad, Tanzania and Uganda, among others) the relatively large cropland area increases compensated for population growth and resulted in small changes in per-capita cropland area. In other countries, cropland increase was not adequate to follow population growth, causing a substantial decrease of cropland per capita (for example, in Ethiopia, Nigeria, Pakistan, Senegal and Tajikistan). Per-capita cropland area decreased almost twofold in Niger, which experienced high population growth and slow cropland expansion. Per-capita cropland area reduction can be an indicator of food insecurity in poor countries that rely on subsistence agriculture, whereas rich countries like Saudi Arabia can compensate for cropland area decline with food imports 28,29 . Several African countries with rapid cropland increase (Angola, Cote d'Ivoire, Democratic Republic of the Congo, Mozambique and Zambia) and South American countries with industrial export-oriented agricultural expansion (Brazil, Bolivia, Paraguay and Uruguay) increased per-capita cropland area. The Baltic states of Lithuania and Latvia had the largest increase of cropland per capita due to cropland gain through recultivation of agricultural lands abandoned in the 1990s, coupled with a sharp population decline (>20% reduction since 2000). Despite their small size, these countries are among the top 15 global wheat exporters.
More than three-quarters (77%) of the global population live in regions with per-capita cropland area and cropland NPP below the year 2019 global average (Extended Data Fig. 8). The lowest per-capita 2019 cropland NPP was in south-west Asia (40% of the global average), which has decreased by 7% since 2003. This finding is aligned with the decrease of per-capita cereal production in western Asia by 14% from 2000 to 2019, reported by the FAO 22 . The per-capita cropland area in south-east Asia in 2019 was half the global average. In contrast, the per-capita cropland areas and cropland NPP in North America, Europe and north Asia in 2019 were twice the global average. South America had nearly threefold higher per-capita cropland NPP than the global average, and it has increased by 59% since 2003. Although per-capita cropland area and  Regions with per-capita cropland area and cropland NPP above the global average include the largest grain-and soybean-exporting countries (Australia, Argentina, Brazil, USA and Russia). Regional accuracies (Table 3) highlight the limitations of the Landsat-based cropland maps. North and South America, which are dominated by large-scale industrial farming, have the highest accuracies. In Europe, Asia and Africa, the global map underestimates cropland area due to spatial resolution limitations in mapping heterogeneous landscapes. Cropland maps in Australia and New Zealand overestimate cropland area due to the inclusion of intensively managed permanent pastures, which are not always separable from crops using Landsat data. In addition, mapping change was shown to be more difficult, with accuracies generally lower across all regions. A probability-based sample analysis is the recommended good practice approach 30 to estimating the extent of and change in land cover and land use, including croplands. The global cropland map time-series enables a higher sampling efficiency through stratification at the subnational, national and global scales. The difference between our sample-based and map-based cropland area estimates and the arable land area reported by the FAO partly related to the definitional inconsistency. The FAO country reports may include unused arable land and other agricultural land uses 10,18 , whereas our estimates represent the actively cultivated cropland area.
The annual MODIS NPP is the only publicly available, globally consistent data that reflect recent changes in primary production within croplands. NPP correlates with crop yield, biomass production and carbon sequestration, although variation of crop types and management practices precludes direct estimation of the crop yield from cropland NPP 31 . The MODIS dataset has been shown to underestimate NPP compared with process-based model estimations, especially for irrigated crops 32 . The difference in spatial resolution between Landsat-based cropland maps and MODIS-based NPP data may impede the analysis of cropland primary production within heterogeneous landscapes. Our Landsat-based, cropland extent time-series data can provide a useful input for improved NPP and crop yield modelling at higher spatial resolution and with better precision.
High-resolution, satellite-based synoptic data on cropland extent and change provide the basis for tracking progress towards sustainable food production and reduction of the environmental impact of agriculture expansion, as well as for applying crop condition monitoring to support decision-making 33 . Cropland extent is a key variable required to estimate emissions from agriculture and is, therefore, a part of the essential climate variables required for monitoring and modelling the Earth's climate 34 . Locally relevant cropland map time-series enable the monitoring of land-use conversion within high-conservation-value ecosystems and protected areas 35 . The cropland extent map, integrated with other high spatial and temporal resolution data, such as forest change 36 and surface water extent 37 , can provide a comprehensive overview of human-induced environmental change, which supports assessing the progress of individual countries towards SDGs.

Cropland-mapping extent and time intervals.
The global boundaries for the cropland mapping were informed by the US Geological Survey (USGS) Global Food Security-Support Analysis Data at 30 m (GFSAD) 11 . The cropland mapping extent was defined using the geographic 1° × 1° grid. We included every 1° × 1° grid cell that contains cropland area according to the GFSAD. Small islands were excluded due to the absence of Landsat geometrically corrected data ( Supplementary Fig. 1).
The we implemented the criterion of the maximum fallow length: if an area was not used as cropland for >4 years, it was not included in the cropland map for the corresponding time interval.
Landsat data. We employed the global 16-day normalized surface reflectance Landsat Analysis Ready Data (Landsat ARD 19 ) as input data for cropland mapping. The Landsat ARD were generated from the entire Landsat archive from 1997 to 2019. The Landsat top-of-atmosphere reflectance was normalized using globally consistent MODIS surface reflectance as a normalization target. Individual Landsat images were aggregated into 16-day composites by prioritizing clear-sky observations.
For each 4-year interval, we created a single annualized gap-free 16-day observation time-series. For each 16-day interval, we selected the observation with the highest near-infrared reflectance value (to prioritize observations with the highest vegetation cover) from 4 years of Landsat data. Observations contaminated by haze, clouds and cloud shadows, as indicated by the Landsat ARD quality layer, were removed from the analysis. If no clear-sky data were available for a 16-day interval, we filled the missing reflectance values using linear interpolation.
The annualized, 16-day time-series within each 4-year interval were transformed into a set of multitemporal metrics that provide consistent land-surface phenology inputs for global cropland mapping. Metrics include selected ranks, inter-rank averages and amplitudes of surface reflectance and vegetation index values, and surface reflectance averages for selected land-surface phenology stages defined by vegetation indices (that is, surface reflectance for the maximum and minimum greenness periods). The multitemporal metrics methodology is provided in detail 19,38 . The Landsat metrics were augmented with elevation data 39 . In this way, we created spatially consistent inputs for each of the 4-year intervals. The complete list of input metrics is presented in Supplementary Table 1.
Global cropland mapping. Global cropland mapping included three stages that enabled extrapolation of visually delineated cropland training data to a temporally consistent, global cropland map time-series using machine learning. At all three stages, we employed bagged decision tree ensembles 40 as a supervised classification algorithm that used class presence and absence data as the dependent variables, and a set of multitemporal metrics as independent variables at a Landsat ARD pixel scale. The bagged decision tree results in a per-pixel cropland probability layer, which has a threshold of 0.5 to obtain a cropland map.
The first stage consisted of performing individual cropland classifications for a set of 924 Landsat ARD 1° × 1° tiles for the 2016-2019 interval ( Supplementary  Fig. 1). The tiles were chosen to represent diverse global agriculture landscapes. Classification training data (cropland class presence and absence) were manually selected through visual interpretation of Landsat metric composites and high-resolution data from Google Earth. An individual supervised classification model (bagged decision trees) was calibrated and applied to each tile.
At the second stage, we used the 924 tiles that had been classified as cropland/ other land and the 2016-2019 metric set to train a series of regional cropland mapping models. The classification was iterated by adding training tiles and assessing the results until the resulting map was satisfactory. We then applied the regional models to each of the preceding 4-year intervals, thus creating a preliminary time-series of global cropland maps.
At the third stage, we used the preliminary global cropland maps as training data to generate temporally consistent global cropland data. As the regional models applied at the second stage were calibrated using 2016-2019 data alone, classification errors may arise due to Landsat data inconsistencies before 2016. The goal of this third stage was to create a robust spatiotemporally consistent set of locally calibrated cropland detection models. For each 1° × 1° Landsat ARD tile (13,451 tiles total), we collected training data for each 4-year interval from the preliminary cropland extent maps within a 3° radius of the target tile, with preference to select stable cropland and non-cropland pixels as training. Training data from all intervals were used to calibrate a single decision tree ensemble for each ARD tile. The per-tile models were then applied to each time interval, and the results were post-processed to remove single cropland class detections and omissions within time-series and eliminate cropland patches <0.5 ha. Manual masks to remove map artefacts (for example, cropland overestimation over temperate wetlands and flooded grasslands) were applied in some regions to improve the map quality. The final global cropland map time-series are available at https://glad.umd.edu/dataset/croplands. Sample analysis. The sample analysis had two objectives: to estimate cropland area and its associated uncertainty and to assess cropland map accuracy. Sample interpretation and sample-based analysis were done only for the start (2003) and the end (2019) of the cropland-mapping interval. Accuracies of intermediate cropland maps (2007, 2011 and 2015) were not assessed, but were considered to be similar to those of the 2003 and 2019 maps due to implementation of the same classification model and consistently processed Landsat data 41 . The analysis was performed separately for each of the seven regions outlined in Extended Data Fig. 1, as well as globally. The regional boundaries were aligned with national boundaries to enable comparison with national data. Only land pixels were considered; pixels labelled as permanent water and snow/ice in the Landsat ARD data quality layer were excluded. In each region, we selected five strata based on the map time-series corresponding to stable croplands, cropland gain and loss, possible cropland omission area and other lands ( Supplementary Tables 2 and 3). The possible cropland omission stratum (stratum 4) includes areas where omission errors are probable, specifically pixels that were not mapped as cropland and either (1) were identified as crops by the GFSAD 11 or (2) had the decision tree-based cropland probability between 0.1 and 0.5. We randomly selected 100 sample units (Landsat data pixels) from each stratum (500 samples pixels per region, 3,500 in total).
Sample interpretation was performed visually using available remotely sensed data time-series, including Landsat ARD 16-day data, composites of selected multitemporal metrics and high-resolution images provided by Google Earth (Supplementary Fig. 2). Each sample pixel was interpreted by two experts independently and the disagreements were discussed and resolved by the research team. The interpretation legend includes the 2003-2019 cropland dynamics categories and land-use transition types. The sample reference data and interpretation results are available at https://glad.umd.edu/dataset/croplands. Area estimation. The sample-based area estimation was performed following previously published methods 42,43 . The 2003 and 2019 total cropland area, stable crops, gross cropland loss and gain, and net change were estimated within each region separately, and for the entire world using equation (1). The area and the total number of Landsat pixels for each region and each stratum are provided in Supplementary Table 3. For each of the 100 sample pixels sampled in each stratum, Proportion of land-use trajectories. We analysed the land-use trajectories of cropland loss and gain using reference sample data within cropland gain and loss strata only. Inclusion of sample pixels from other strata where cropland change was detected would have inflated the area of land-use trajectories that these pixels represent (that is, if a sample pixel from a stable cropland stratum was interpreted as cropland gain due to forest clearing, including the proportion of forest clearing from this large stratum, it will dominate the total regional estimate). The proportion of each land-use trajectory (within cropland gain and loss separately) was estimated from the sample and reported as the percentage of the total gain or loss along with its s.e.m. (Table 2). A combined ratio estimator for stratified random sampling 43 was employed to estimate the percentages (equation (3)).
where: R is the estimated class proportion expressed as a percentage; H the number of sampling strata; where: s.e.
(R) is the s.e.m. of the estimated proportion expressed as a percentage; N h the total number of pixels in stratum h; n h number of sample pixels in stratum h; A hxh the estimated total area of cropland loss/gain expressed in area units; and s 2 yh and s 2 xh the sample variances in stratum h; and s xyh the sample covariance in stratum h estimated as follows: Map accuracy. The map accuracy metrics include overall accuracy (the proportion of correctly mapped sample pixels), user's accuracy of the cropland class (which reflects the cropland class commission) and producer's accuracy of the cropland class (which reflects the cropland class omission) 42 . All accuracy metrics and respective s.e.m.s are presented as percentages (Table 3).
To estimate overall accuracy, we defined y u = 1 if pixel u is classified correctly and y u = 0 if pixel u is classified incorrectly. The estimator for overall accuracy is then expressed by equation (5), and s.e.m. for overall accuracy is computed using equation (6).
where: Ô is the estimated overall accuracy, expressed as a percentage; H the number of sampling strata; N h the total number of pixels in stratum h; N the total number of pixels in the reporting region; and ȳ h = ∑ u∈h yu/n h the sample mean of the y u values in stratum h. s.e.
(Ô) is the s.e.m. of the overall accuracy, expressed as percentage; n h the number of sample pixels in stratum h; and s 2 yh the sample variance: . For estimating user's accuracy of the croplands class, we defined y u = 1 if sample pixel u is correctly mapped as cropland, otherwise y u = 0, and x u = 1 if sample pixel u is mapped cropland, otherwise x u = 0. For the producer's accuracy, we defined y u = 1 if sample pixel u is correctly mapped as cropland, otherwise y u = 0, and x u = 1 if sample pixel u is interpreted as cropland, otherwise x u = 0. The estimator of the user's accuracy and producer's accuracy was then expressed as a ratio estimator (equation (7)) and their s.e.m. calculated using equation (8), which are similar to equations (3) and (4), except that the strata were weighted by their total number of pixels (N h ) rather than the areas (A h ) for the purposes of map accuracy assessment (with pixel being the primary mapping unit): where R is the estimated user's/producer's accuracy, expressed as a percentage.
(R) is the s.e.m. of the estimated user's/producer's accuracy, expressed as a percentage.
Cropland NPP. The cropland NPP was evaluated using the globally consistent Collection 6 MODIS-based, annual year-end gap-filled NPP product (MOD17A3HGF 20 ). The product provides the sum of total daily NPP through the year at a 500-m spatial resolution (kg C m −2 year −1 ). The annual NPP data were resampled to our Landsat ARD data grid and were overlaid with the corresponding 4-year cropland maps to calculate total and per-unit area NPP for each region and each year. We used average annual NPP for each 4-year interval, except for the 2000-2003 interval, where a 3-year average was used instead to avoid using the year 2000 when MODIS data were incomplete. The s.d. of the annual estimates is provided as an uncertainty metric.

March 2021
Corresponding author(s): Peter Potapov Last updated by author(s): Aug 19, 2021 Reporting Summary Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Portfolio policies, see our Editorial Policies and the Editorial Policy Checklist.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted

Software and code
Policy information about availability of computer code Data collection The satellite data processing, classification, and sample analysis were done using freely available GLAD ARD Tools V1.1 (glad.umd.edu/ard).
Source Landsat data obtained from public sources (USGS data depository) Data analysis Data analysis performed using GLAD ARD Tools V1.1 (glad.umd.edu/ard) and other open-source packages (GDAL, R, QGIS) For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Portfolio guidelines for submitting code & software for further information.

Data
Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A description of any restrictions on data availability -For clinical datasets or third party data, please ensure that the statement adheres to our policy All global cropland maps, cropland dynamic maps, and sample data available at https://glad.umd.edu/dataset/croplands/. The global Landsat ARD data, metric generation and image classification codes are available from https://glad.umd.edu/ard.

March 2021
Field-specific reporting Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection. All studies must disclose on these points even when the disclosure is negative.

Study description
Global satellite-based wall-to-wall cropland mapping and sample-based area estimation for 2000-2019 time interval. Cropland maps were integrated with the Moderate Resolution Imaging Spectroradiometer-derived annual net primary production (NPP) as a proxy variable for analyzing crop primary production.

Research sample
We utilized the consistently processed global Landsat satellite data archive from 2000 to 2019 for both wall-to-wall mapping and sample analysis. Global MODIS-based annual NPP within cropland extent used as proxy variable to analyze crop productivity.

Sampling strategy
Stratified random sampling design, with Landsat data pixels serving as a sampling unit. We randomly selected 100 sample units (Landsat data pixels) from each stratum (500 samples pixels per region, 3,500 in total). The strata areas and pixel numbers are provided with the manuscript. The entire database of sample data are provided on-line (https://glad.umd.edu/dataset/croplands/)

Data collection
Mapping and sample analysis were performed by the authors. Mapping was done through extrapolation of manually collected training data at the global extent using machine learning tools. Sample interpretation was performed visually using available remotely sensed data time-series, including Landsat ARD 16-day data, composites of selected multi-temporal metrics, and high-resolution images from Google Earth.
Timing and spatial scale The analysis was performed in four-year intervals (2000-2003, 2004-2007, 2008-2011, 2012-2015, and 2016-2019). The analysis was performed at the global extent. The spatial resolution of the cropland maps is 0.00025 degree per pixel, or approximately 30 m at the Equator.