Accelerated global cropland expansion and primary production increase in the 21st century

Spatiotemporally consistent data on global cropland extent is a key to tracking progress 10 toward hunger eradication and sustainable food production 1,2 . Here, we present an analysis 11 of global cropland area and change for the first two decades of the 21 st century derived 12 from satellite data time-series. We estimate 2019 cropland area to be 1,244 Mha with a 13 corresponding total annual net primary production (NPP) of 5.5 Pg C yr -1 . From 2003 to 14 2019, cropland area increased by 9% and crop NPP by 25%, primarily due to agricultural 15 expansion in Africa and South America. Global cropland expansion accelerated over the 16 past two decades, with a near doubling of the annual expansion rate, most notably in 17 Africa. Half of the new cropland area (49%) replaced natural vegetation and tree cover, 18 indicating a conflict with the sustainability goal of protecting terrestrial ecosystems. From 19 2003 to 2019 global population growth outpaced cropland area expansion, and per capita 20 cropland area decreased by 10%. However, the per capita annual crop NPP increased by 21 3.5% as a result of intensified agricultural land use. The presented global high-resolution 22 cropland map time-series supports monitoring of sustainable food production at the local, 23 national, and international levels. Global and

cropland map time-series supports monitoring of sustainable food production at the local, 23 national, and international levels. 24 Global population growth and increasing standards of living inevitably cause the expansion and 25 intensification of global agricultural land use to fulfill growing demands for food, biofuel, and 26 other commodities 3-5 . In turn, agriculture expansion and intensification threaten ecosystem 27 functioning and lead to species extinction through habitat loss and fragmentation 4,6,7 . The United 28 Nations 2030 Sustainable Development Goals (SDGs) call for balancing increasing agricultural 29 production with maintenance of ecosystem services 8 . Implementation of SDGs to improve food 30 security, protect freshwater and terrestrial ecosystems, and mitigate climate change require 31 national policies and international cooperation that are based on consistent, independent, and 32 timely data on agriculture extent and productivity 2 . Spatiotemporally consistent satellite 33 observations provide the most accurate and cost-effective solution for global agricultural land 34 use mapping and monitoring 9 . Satellite data have been shown to enable national and global 35 agriculture mapping [10][11][12][13][14] . However, no globally consistent cropland time-series data at locally 36 relevant spatial resolutions exist to date. 37 Here, we present a global cropland extent and change dataset that can serve as a tool for 38 monitoring national and global progress towards SDGs. We define cropland as land used for 39 annual and perennial herbaceous crops for human consumption, forage (including hay), and 40 biofuel. Perennial woody crops, pastures, and shifting cultivation are excluded from the 41 definition. The fallow length is limited to four years for the cropland class. Our definition is 42 consistent with the arable land category reported by the Food and Agriculture Organization of 43 (with no crop cultivation for more than 4 years). In North and South America cropland expansion 85 through the conversion of pastures and long fallows was more common (75% and 61%, 86 respectively) than through clearing of natural vegetation. 87 Abandonment or conversion to other land uses affected 10% of the 2003 cropland area (115.5 ± 88 24.1 Mha). Of that area, 52% was converted into pastures or abandoned (Extended Data Table   89 2); such conversions may be temporary and followed by crop recultivation years later. Industrial 90 and residential construction and infrastructure development was the second largest driver of 91 gross cropland loss, responsible for 16% of the total cropland area reduction. In Southeast Asia, 92 35% of cropland reduction was due to urban sprawl. A portion (13%) of 2003 cropland was 93 converted to permanent woody crops or aquaculture, with the highest proportion of such 94 transitions in Southeast Asia (28%). Flooding caused by surface water increase, water erosion, 95 and reservoir construction affected cropland area on all continents (3% total reduction). The

125
The global MODIS-derived annual NPP within cropland area (Extended Data Fig. 4) increased 126 by 25% between 2003 and 2019 (from 4.4 Pg C yr -1 to 5.5 Pg C yr -1 , Fig. 1). South America had 127 the highest NPP increase (by 0.38 Pg C yr -1 , or 88%) followed by Africa (by 0.29 Pg C yr -1 , or 128 50%) (Extended Data Table 5). The per capita annual crop NPP also increased globally by 3.5%, 129 balancing the per capita cropland area reduction. Two processes contributed to the global crop 130 NPP increase, namely the increase of cropland area and the increase in crop productivity per unit 131 area. We found that the mean NPP per unit area within stable crops increased by 10%, from 402    above the global average include the largest grain and soybean exporting countries (Australia,194 Argentina, Brazil, USA, Russia).  The global boundaries for the cropland mapping were informed by the USGS global cropland 253 map 11 . The cropland mapping extent was defined using the geographic 1°×1° degree grid. We 254 included every 1°×1° grid cell that contains cropland area according to the USGS map. Small 255 islands were excluded due to the absence of Landsat geometrically corrected data 256 ( Supplementary Information Fig. 1).

257
The cropland mapping was performed at four-year intervals (2000-2003, 2004-2007, 2008-2011, 258 2012-2015, and 2016-2019). Using a long interval (rather than a single year) increased the  Table 1.  At the second stage, we used the 924 tiles that had been classified as crop/no-crop and the 2016-299 2019 metric set to train a series of regional cropland mapping models. The classification was 300 iterated by adding training tiles and assessing the results until the resulting map was satisfactory.

301
We then applied the regional models to each of the preceding four-year intervals, thus creating a 302 preliminary time-series of global cropland maps.

303
At the third stage, we used the preliminary global cropland maps as training data to generate 304 temporally consistent global cropland data. Because the regional models applied at the second 305 stage were calibrated using 2016-2019 data alone, classification errors may arise due to Landsat 306 data inconsistencies before 2016. The goal of this third stage was to create a robust 307 spatiotemporally consistent set of locally calibrated cropland detection models. For each 1°×1° 308 Landsat ARD tile (13,451 tiles total) we collected training data for each four-year interval from 309 the preliminary cropland extent maps within a 3° radius of the target tile, with preference to 310 select stable crop and no-crop pixels as training. Training data from all intervals were used to 311 calibrate a single decision tree ensemble for each ARD tile. The per-tile models were then 312 applied to each time interval, and the results were post-processed to remove single cropland class 313 detections and omissions within time-series and eliminate cropland patches below 0.5 ha.

314
Manual masks to remove map artifacts (e.g., crop overestimation over temperate wetlands and 315 flooded grasslands) were applied in some regions to improve the map quality. The final global 316 cropland map time-series are available at https://glad.umd.edu/dataset/croplands/.

318
The sample analysis had two objectives: to estimate cropland area and its associated uncertainty 319 and to assess cropland map accuracy. The analysis was performed separately for each of the 320 seven regions outlined in Extended Data Fig. 1, as well as globally. The regional boundaries 321 were aligned with national boundaries to enable comparison with national data. Only land pixels 322 were considered; pixels labeled as permanent water and snow/ice in the Landsat ARD data 323 quality layer were excluded. In each region, we selected five strata based on the map time-series 324 corresponding to stable croplands, cropland gain and loss, possible cropland omission area, and 325 other lands ( Supplementary Information Tables 2 and 3). The possible cropland omission stratum 326 (stratum 4) includes areas where omission errors are probable, specifically, pixels that were not 327 mapped as cropland and either (a) were identified as crops by the USGS cropland map 11 or (b) 328 had the decision tree-based cropland probability between 0.1 and 0.5. We randomly selected 100 329 sample units (Landsat data pixels) from each stratum (500 samples pixels per region, 3,500 in 330 total).

331
Sample interpretation was performed visually using available remotely sensed data time-series, 332 including Landsat ARD 16-day data, composites of selected multitemporal metrics, and high-333 resolution images provided by Google Earth ( Supplementary Information Fig. 2). Each sample   354 We analyzed the land-use trajectories of cropland loss and gain using reference sample data 355 within cropland gain and loss strata only. Inclusion of sample pixels from other strata where 356 cropland change was detected would have inflated the area of land use trajectories that these 357 pixels represent (i.e., if a sample pixel from a stable cropland stratum was interpreted as cropland 358 gain due to forest clearing, including the proportion of forest clearing from this large stratum will 359 dominate the total regional estimate). The proportion of each land use trajectory (within cropland 360 gain and loss separately) was estimated from the sample and reported as the percent of the total 361 gain or loss along with its standard error (Extended Data  where -standard error of the estimated proportion expressed as percentage; -total number of pixels in stratum h; -number of sample pixels in stratum h;

Proportion of land use trajectories
-the estimated total area of cropland loss/gain expressed in area units; and are the sample variances in stratum h and is the sample covariance in stratum h estimated as follows: 367 The map accuracy metrics include overall accuracy (the proportion of correctly mapped sample 368 pixels), user's accuracy of the cropland class (which reflects the cropland class commission) and 369 producer's accuracy of the cropland class (which reflects the cropland class omission) 33 . All 370 accuracy metrics and respective standard errors are presented as percentages (Extended Data 371 Table 3).

372
To estimate overall accuracy, we defined if pixel u is classified correctly and if 373 pixel u is classified incorrectly. The estimator for overall accuracy is then expressed by Eq. 5, 374 and standard error for overall accuracy is computed using Eq. 6.

375
Eq. 5 where -estimated overall accuracy, expressed as percentage; H -number of sampling strata; -total number of pixels in stratum h; N -total number of pixels in the reporting region; is the sample mean of the values in stratum h.

376
Eq. 6 where -standard error of the overall accuracy, expressed as percentage; -number of sample pixels in stratum h; -the sample variance:

Eq. 8
where -standard error of the estimated user's/producer's accuracy, expressed as percentage. year -1 ). The annual NPP data were resampled to our Landsat ARD data grid and were overlaid 393 with the corresponding four-year cropland maps to calculate total and per unit area NPP for each 394 region and each year. We used average annual NPP for each four-year interval, except for the  Fig. 1) are aligned with country boundaries, we were able to summarize the 402 regional population totals from national data. The population data were related to our sample-403 based (for global and regional estimates) and map-based (for national estimates) cropland area to 404 estimate per capita cropland area and change. Similarly, we related regional crop NPP to 405 population data to estimate per capita crop NPP for 2003 and 2019.

550
OA stands for overall accuracy, UA for user's accuracy, and PA for producer's accuracy.

551
Standard errors of accuracy metrics are shown in parenthesis.