ISPRS Journal of Photogrammetry and Remote Sensing

Mapping high resolution (30-m or better) cropland extent over very large areas such as continents or large countries or regions accurately, precisely, repeatedly, and rapidly is of great importance for addressing the global food and water security challenges. Such cropland extent products capture individual farm ﬁ elds, small or large, and are crucial for developing accurate higher-level cropland products such as cropping intensities, crop types, crop watering methods (irrigated or rainfed), crop productivity


Introduction
Accurate, and precise agricultural cropland extent products over very large areas that map small to large farms are of great importance to assess and monitor global food and water security. They are a critical part of land system studies (Verburg et al., 2013). Such products are also of great importance for assessing global crop water use, crop productivity (productivity per unit of land), water productivity (productivity per unit of water or crop per drop), and food security studies (Foley et al., 2011;Thenkabail et al., 2010;Teluguntla et al, 2015b;Matejicek and Kopackova, 2010). Remote sensing based spatially distributed cropland maps with high spatial resolution provide an efficient way to monitor croplands (Foley et al., 2011;Fritz et al., 2015;Yu et al., 2013). Over the last two decades, several global and regional cropland products have been produced using medium to coarse resolution (250m to 1-km) remote sensing data such as the Advanced Very High-Resolution Radiometer (AVHRR) and the Moderate-Resolution Imaging Spectroradiometer (MODIS) data Kumar et al., 2018;Thenkabail et al., 2009Pittman et al., 2010;Portmann et al., 2010;Siebert and Döll, 2010;Salmon et al., 2015;Waldner et al., 2015Waldner et al., , 2016. These products are very useful for a preliminary understanding of agricultural croplands in terms of their spatial distribution patterns and their characteristics such as crop dominance and cropping intensities. However, the coarse resolution of these products limits their usefulness in assessing small agriculture fields (Teluguntla et al., 2015b;Thenkabail et al., 2010). Further, there are several global to regional land use/land cover (LULC) products produced using multiple remote sensing data in which agricultural croplands is one or more classes. Some examples are: DIScover (Loveland et al., 2000); GLC500m (Friedl et al., 2010); MCD12Q1 (Liang et al., 2015); Globecover (Defourny et al., 2009); FROM-GC ; FROM-GLC ; and Globeland30 (Arsanjani et al., 2016;Chen et al., 2015). However, these products were focused on LULC in which mapping croplands in detail was not the primary objective. Hence, the cropland accuracies suffer (Yang et al., 2017). Further, most of these products are also coarse resolution. Most of these products fail to map individual farm fields, especially when they are small and\or fragmented. Definitions of croplands also vary from product to product, resulting in different results of cropland extent and their characteristics in each of these products. Overall, existing cropland extent products are coarse resolution, lack field level details, and/or are mapped as part of other LULC classes where specific cropland class focus is missing. As a result, uncertainties and errors in cropland locations are very high.
Object-based approaches (Peña-Barragán et al., 2011;Peña et al., 2014) include Hierarchical Image Segmentation software or HSeG (Tilton et al., 2012). A combination of pixel-based and object-based methods have also been recently attempted Chen et al., 2018). However, these methods and approaches were overwhelmingly applied on: (a) multi-temporal moderate resolution (250-m or higher) remotely sensed data, and/or (b) small areas, and/or (c) high-resolution (Landsat 30-m) remotely sensed data with limited multi-temporal images.
Hitherto, availability of cloud-free, high quality images as well as use of multi-temporal, high-resolution data over very large areas for cropland mapping has been challenging and resource intensive. However, these challenges have been overcome through a paradigm shift in remote sensing data collection, management, and processing. First, Landsat-8 Operational Land Imager (OLI) data and Landsat-7 Enhanced Thematic Mapper+ (ETM+) data at 30-m spatial resolution were utilized every 16-days for 3-years (2013-2015) for entire Australia and China. Managing massive volumes of Landsat data for analysis over very large areas is a big challenge when adopting traditional remote sensing approaches that use commercial imaging processing software on workstation PC based systems. No matter how powerful the systems are, the entire process of data analysis including, pre-processing, over very large areas involving 1000's of Landsat images is cumbersome, slow, and tedious. However, in the current era of adopting powerful machine learning algorithms (MLAs) in cloud computing environments such as Google Earth Engine (GEE) these limitations have been overcome allowing planetary scale remote sensing at high spatial resolutions as illustrated by Erickson (2014) and Gorelick et al. (2017). Gorelick et al. (2017) demonstrated that multi-petabyte archive of georeferenced datasets can be combined in the GEE catalog that includes images from Earth observing satellites and airborne sensors (e.g., USGS Landsat, NASA MODIS, USDA NASS CDL), weather and climate datasets, as well as digital elevation models. This system has exceptional data organization, and has enabled geo spatial processing over very large areas. Along with computing and storage resources, GEE also supports major MLAs useful for image enhancement and, image classification, and allows batch processing through JavaScript or Python on Application Program Interfaces (APIs). These capabilities reduce most of preprocessing steps needed in traditional remote sensing approaches. Very recently, several studies have used the GEE platform for largescale (continental, global) mapping .
Thereby, the overarching goal of this study was to map cropland extent in detail (e.g., showing all individual farms whether small or big) using high-resolution (30-m) multi-year (2013-2015) time-series (16day) Landsat-8 OLI data over the entirety of China and Australia. These two countries have very large cropland areas with distinct cropping systems. Australia is a major agriculture producer and exporter with large scale industrial farm fields. Pastoral farming is another major agricultural land use in Australia. Australian farmers produce cereals, legumes and oilseeds on a large scale (ABARES, 2016) for human consumption and livestock feed. In contrast, an overwhelming proportion of Chinese farms are small, fragmented, but often contiguous over large areas due to intensive land use for agriculture. Average crop field size in China is less than a hectare (Samberg et al., 2016); fine resolution satellite data is required to map such crop fields. Landsat data with 30-m spatial resolution is ideal dataset to map cropland extent in China. Whereas the average crop field size in Australia is 100 ha (Samberg et al., 2016) which is much larger than most crop fields in China. Medium resolution sensor such as MODIS at 250-m (1 pixel is approximately 6.25 ha) are highly inadequate to map smaller and/or fragmented crop fields. High-resolution 30-m (1 pixel is approximately 0.09 ha) Landsat-8 OLI 30-m imagery is expected to map small and fragmented farms in addition to large farms. Chinese farmlands are also very diverse, spread across mountains, river banks, and large plains. China's agriculture feeds 1.38 billion people whereas Australia's much smaller population allows it to be a major exporter of food. China is the topmost producer of three major food crops (rice, wheat, and maize) in the world (FAOSTAT). The best cropland products of Australia and China known by the authors are derived from other LULC products. For example, for Australia, croplands were mapped as part of land cover products such as Dynamic Land Cover Data (DLCD) developed by Geo Science Australia (Lymburner et al., 2011). Another cropland product for Australia is from the Australian Collaborative Land Use and Management Program (ACLUMP) by the Department of Agriculture and Water resources (http://www.agriculture.gov.au/abares/aclump). However, there is no unique cropland layer available for Australia and China like United States Department of Agriculture (USDA) National Agricultural Statistical Service (NASS) 30 m Cropland Data Layer (CDL) derived from Landsat for USA (Johnson and Mueller, 2010;Boryan et al., 2011;Han et al., 2012). Due to the above reasons, mapping croplands at high-resolution of these two countries is of great importance.
Thereby, specific objectives of the study were to develop a precise 30-m cropland extent map of Australia and China using 16-day Landsat-8 OLI data for the nominal year 2015 using an RF classifier through the GEE cloud-computing platform. We chose a pixel-based RF MLA approach for this study because MLAs have been particularly useful in classifying large datasets with high spatial and temporal resolutions for large area land cover mapping (Pelletier et al., 2016). Second, this study used a large collection of reference training and validation data obtained from: (a) ground-data through extensive field visits, (b) samples interpreted from sub-meter to 5-m very high-resolution imagery data, and (c) data sourced from reliable secondary sources. These reference datasets were invaluable in training the RF classifier, as well as to assess classification accuracies, and establish uncertainties. Third, computed cropland area statistics from the 30-m cropland product were compared with national and sub-national cropland area statistics obtained from sources such as the Food and Agricultural Organization (FAO) of the United Nations (UN) and National Governmental sources.

Study area
Our study area consisted of two major cropland areas of the world: China and Australia. The two regions were further stratified into four refined agro-ecological zones (RAEZs) (Fig. 1), which helped identify areas of similar farming practices, soil types, and climatic patterns. These RAEZs were modifications of the Food and Agricultural Organization's (FAO's) global agro-ecological zones (GAEZs) which are based on the number of growing degree days, soil, and terrain data, all provided at 10-km spatial resolution (Fischer et al., 2000). However, GAEZs has too many zones for cropland classification purposes as many of these zones only have a very small proportion of crops relative of the overall land area of the zone. Therefore, we refined GAEZs into RAEZs by utilizing Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model Version 2 (GDEM V2) 30-m data (https://www.asterweb.jpl.nasa.gov/gdem.asp), slope Fig. 1. Stratification of the study area into distinct and broad refined agro-ecological zones (RAEZs). The figure also shows the distribution of the reference training data used in the Random Forest machine-learning algorithm. The Random Forest pixel-based supervised machine learning algorithm used in this study was trained using reference training data falling within each of these zones to separate croplands from non-croplands.
derived off 30-m GDEM data, and the proportion of croplands in a zone. A number of smaller RAEZs were combined into larger zones depending on the zone's importance for croplands. This resulted in three broad RAEZs across China, and a single RAEZ for Australia ( Fig. 1). Classification using MLAs was carried out in each of these RAEZs separately and then similar classes across zones were combined.

Datasets
Two major types of data were used in the study. First, we will describe and discuss the satellite sensor data followed by reference training and validation data of Australia and China.

Satellite imagery: Landsat-8 and Landsat-7
In order to cover crop dynamics in different time-periods, Landsat-8 OLI Top of Atmosphere (TOA) 16-day satellite sensor data collections available on GEE cloud for three years (2013)(2014)(2015) were organized for Australia and China. In addition, Landsat-7 was used to fill some data gaps (Irons et al., 2012) over 2013-2015 for China, aiming to provide seamless cloud-free 30-m wall-to-wall data coverage (Fig. 2a). Specific image collections used were LANDSAT/LC8_L1T_TOA and LANDSAT/ LE7_L1T_TOA for 2013-2015 (sometimes referred to as nominal 2015). The Landsat TOA data available on GEE are orthorectified and corrected for solar angle (USGS TOA product guide). The coefficients were derived from USGS TOA product guide and Chander et al. (2009). There is a 16-day revisit time for Landsat-8 OLI and Landsat-7 ETM+ 30-m data. It is not possible to get continuous 16-day cloud free timeseries data for wall-to-wall coverage for any part of the region due to cloud cover. To overcome this limitation and to ensure cloud-free or near-cloud-free wall-to-all coverage, bi/tri-monthly composites (depending on the cloudiness of the countries/regions), were composed (e.g., Figs. 2a and 2b). Finally, 30-m mega-file data-cubes (MFDCs) were created as per the following steps leading to a 48-band MFDC ( . 30-m Data-cube for China composited for 4 time-periods using 16-day data of Landsat-8 and Landsat-7 for the years 2013-2015 (e.g., period 1: Julian days 1 to 90). For each period eight bands (blue, green, red, NIR, SWIR1, SWIR2, TIR1 and NDVI) were composited, taking the median value of a given pixel over the period. From the 4 time-periods, there was a 32-band mega-file data cube in the Google Earth Engine collection. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) a 12-month period. In comparison, tri-monthly time-periods were created for China (Fig. 2b) because a greater number of days were required to achieve cloud-free or near cloud-free wall-to-wall coverage due to greater number of cloudy days over the country. As a result, there were only 4 time-periods (period 1: Julian days 1-90, period 2: 91-180, period 3: 181-270, period 4: 271-365; Fig. 2b) over 12-month period. It should be noted that for each time-period (e.g., 1-60 days) Landsat data from multiple years (2013)(2014)(2015) were used to maximize the chances of pure cloud free pixels over Australia and China. As a result, the process involved gathering all the Landsat-8 16-day images over Australia ( Fig. 2a), and all the Landsat-8 as well as Landsat-7 images over China (Fig. 2b) available for each time-period (e.g., period 1: 1-60 Julian days). The following bands and indices were used for each time-period: blue (0.45-0.51 μm), green (0.53-0.59 μm), red (0.63-0.69 μm), NIR (0.85-0.89 μm), SWIR1 (1.55 1.65 μm), SWIR2 (2.1-2.3 μm), and TIR1 (10.60-11.19 μm); (Table 1, Figs. 2a and 2b) along with Normalized Difference Vegetation Index (NDVI) were time-composited by taking median value of each pixel of each band (Roy et al., 2014). Thereby, for Australia, eight median value bands composed over 6 time-periods resulted in a 48-band MFDC (Fig. 2a). Whereas for China, eight median value bands composed over 4 time-periods resulted in a 32-band MFDC (Fig. 2b). The band stack and time-periods leading to MFDC are shown in Table 1 as well as in Figs. 2a and 2b. All compositions were performed on the GEE cloud-based geospatial platform for planetary-scale data analysis (Gorelick et al., 2017). Landsat TOA products were used instead of surface reflectance (SR) due to the limited temporal availability of Landsat-7 and Landsat-8 SR imagery on GEE for our study regions.

Reference croplands samples for training and validation
Reference training and validation data were collected for Australia and China, from three high quality sources as described below. Part of the training data was used to gather knowledge and train the MLA and the other part was used to test the classes. Validation data were used to assess accuracies, errors, and uncertainties. Data were collected keeping the cropland definition in mind. Croplands were defined in this study as agricultural lands with annual standing crops + cropland fallows + permanent crops (e.g., coffee, tea, rubber) (Teluguntla et al., 2015b). Samples were collected to cover wide range of cropland and non-cropland classes (e.g., Figs. 3a and 3b) with excellent spatial distribution in Australia (Fig. 4a) and China (Fig. 4b). Since the satellite data were collected for the 2013-2015 time-period there are possibilities that croplands in one year or season may remain cropland fallows in another year or season . The cropland areas of Australia and China are well established. Hence croplands changing to non-croplands in Australia and China in such short time period (2013)(2014)(2015) was minimal and was estimated by us to be < 1% of the total cropland areas of these countries. Further, we also observed the multi-year NDVI patters of the sample areas to see whether the sample locations were "pure" (unchanged) or "changed" (showed significant change over time). Reference training and validation data were then appropriately labeled.
3.2.1. Ground reference and validation data through field visits First, we conducted an extensive field survey during September and October of 2014, the peak crop-growing season for crops in Australia. More than 4400 ground samples (e.g., Figs. 3a and 4a) were collected from New South Wales (NSW), Victoria (VIC), South Australia (SA), and Western Australia (WA) states of Australia following the specific guidelines on collecting ground reference data (Congalton, 2015). The sampling sites included various crop fields: such as cereal crops (wheat, barley, and oats), legumes (lupin, lentils, peas, and beans), oilseeds (canola), vegetables, continuous crops (orchard crops), fodder crops (alfalfa and sown pastures), some cropland fallow samples and noncropland samples. Such data collection over a wide range of crop types helped us to capture variability across croplands and helped in the development of robust training datasets. A sample illustration of the ground reference data acquired across Australia is shown in Fig. 3a and the spatial distribution of data across Australia is shown in Fig. 4a. The samples shown in Fig. 4a cover the entire agricultural landscape of Australia except the northern part of the country where we used submeter to 5-m data.
Similarly, we obtained ground reference data for China through collaboration with the Chinese Academy of Agricultural Sciences (CAAS) between 2008 and 2014 time-period. These field data were spatially well distributed throughout agricultural cropland areas of China. A total of 2120 ground samples (e.g., Figs. 3b and 4b) including: (1) location of samples (GPS position, location name, date of collection); crop properties (2) croplands versus non-croplands; (3) irrigated or rainfed; (4) crop intensity (single, double, triple, continuous cropping in 12 months); (5) crop type (major crop types mentioned above, others); and (6) digital photographs of each sample. A sample illustration of the ground reference data distribution across China is shown in Fig. 3b. For both Australia and China, there was a minimum of onekilometer distance between each sample to ensure spatial autocorrelation do not occur (Congalton, 2015). We further verified all ground data samples for croplands versus non-croplands using sub-meter to 5-m very high-resolution data (also see Section 3.2.2).

Very high spatial resolution imagery (VHRI) data
Second, some samples were obtained throughout Australia, and China by interpreting sub-meter to 5-meter very high spatial resolution imagery (VHRI) data of 4 recent years (2012-2015) available to US Government entities through the National Geospatial Agency (NGA) (Marshall and Thenkabail, 2015) and multiple years of imageries obtained from Google Earth Pro. From this, we collected 1210 reference and 1710 validation samples (e.g., Fig. 4a for Australia, Fig. 4b for China).

Tertiary sources
Third, reference data were obtained from national systems and other reliable sources such as Geo Science Australia, Australian Bureau of Agricultural and Resource Economics and Sciences (ABARES).
A summary of sample distribution of the reference training and validation data for Australia and China is provided in Table 2. The   Fig. 4a. Sub-meter to 5-m very high-resolution imagery data for Australia. Illustration of the reference training and validation data of Australia collected using submeter to 5-m very high spatial resolution imagery. Fig. 4b. Sub-meter to 5-m very high-resolution imagery data for China. Illustration of the reference training and validation data of China collected using sub-meter to 5-m very high spatial resolution imagery. P. Teluguntla et al. ISPRS Journal of Photogrammetry and Remote Sensing 144 (2018) 325-340 reference data were divided into two distinct datasets: roughly 50% for training and 50% for validation ( Table 2). The training data were used to develop knowledge base and train the machine-learning algorithms (e.g., Random Forest in this study). The validation data were used to perform accuracy assessments leading to errors of omissions, errors of commissions, and overall accuracies.

Methodology
The goal of this study was to produce accurate Landsat 30-m derived cropland extent product for Australia and China. We used a pixel-based supervised classification approach with the RF classifier on the GEE cloud computing platform to develop cropland extent product for Australia, and China using Landsat 30-m, 16-day time-series data (Figs. 2a and 2b) for the 2013-2015 time-period. An overview of the methodology is shown in Fig. 5.

Random forest machine learning algorithm
We selected the pixel-based supervised Random Forest (RF) machine learning algorithm because it is generally immune to data noise and overfitting and is extremely useful in classifying remote sensing data. Furthermore, RF classifiers can successfully handle high data dimensionality and typically achieve higher accuracies in comparison with other approaches such as maximum likelihood, single decision trees, and single layer neural networks (Belgiu and Drăguţ, 2016;Lawrence et al., 2006;Na et al., 2010). Random Forest algorithm is a machine learning method, in which the RF classifiers construct multiple de-correlated random decision trees that are bootstrapped and aggregated to classify a dataset by using the mode of predictions from all decision trees (Breiman, 2001). The RF classifier is more robust than single decision tree (Chan and Paelinckx, 2008) and easier to be implemented than many other advanced classifiers such as Support Vector Mission (SVM) (Pelletier et al., 2016) as recently established by Xiong et al. (2017a). It uses bootstrap aggregating (bagging) to form an ensemble of decision trees by searching random subspaces from the given data (features) and the best splitting of the nodes by minimizing the correlation between the trees. Random Forest classifiers also provide a quantitative measurement of each variable's contribution to the classification output, which is useful in evaluating the importance of each variable. They provide an internal accuracy assessment by using an 'out-of-bag' (OOB) technique, in which about a third of the data is kept aside as a validation dataset to assess accuracy of the classification. This technique can be used to cross-validate the RF classifier using independent datasets.
Random Forest classifiers available in GEE use six input parameters: (1) number of classification trees, (2) number of variables used in each classification tree, (3) minimum leaf population, (4) bagged fraction of the input variables per decision tree, (5) out-of-bag mode, and (6) random seed variable for decision tree construction. When the number of trees increases, the overall accuracy of classification increases without overfitting (Breiman, 2001). While training sample imbalance can affect the RF classification output by over-fitting the majority class (Breiman, 2001;Chen et al., 2004), various methods such as downsampling the majority class or duplicating the minority samples can provide immunity against over-fitting (Sun et al., 2007). Further, the choice of optimized parameter values using the OOB outputs can eliminate overfitting in the RF classifier. The optimized parameter values were selected by selecting the training samples (Table 2), running the RF algorithm, and testing the classification output for overall, producer's and user's accuracies in error matrix. The goal is to obtain not just the high overall accuracies, but also a good balance of producer's accuracies (or least errors of omissions) and user's accuracies (or least errors of commissions). It is not just the high number of training samples of a class that help attain optimal accuracies, but the purity (e.g., pure cropland samples instead of mixed) of the samples as well.
All supervised pixel-based classifications rely heavily on the input training samples. To discriminate croplands under various environments and conditions, two criteria are very important: 1. RF classifications need to take place in AEZs or RAEZs (e.g., Fig. 1), 2. The sample size (Table 2) of the initial training dataset for the RF classifier needs to be large, especially in complex regions.
All samples were selected to represent a 90-m × 90-m polygon (Table 2 and Section 3 and its sub-sections). First, we made extensive field campaigns in Australia during the 2014-crop growing season when data were collected on precise cropland locations as well as non-cropland locations. This effort led to collection of more than 628 samples spread across Australia (e.g., Table 2). Second, we utilized the ground data from previous efforts for China and other reliable sources. Third, sub-meter to 5-m VHRI, available to us for the entire study region, was used to generate croplands versus non-cropland interpretations by multiple analyses across China, Australia, and a total of ∼1490 data samples were used from these interpretations. To move forward with a larger sample size, an iterative sample selection procedure was introduced with the following steps for training the RF machine-learning algorithm as illustrated in Fig. 5   or abandoned fields), historical Landsat Images and MODIS NDVI time-series are also referenced. All the samples selected to represent a 90-m × 90-m polygon. 5. Loop steps 1-4 with enlarged training dataset until classification becomes stable and high classification accuracy is achieved.
The number of iterations required for the training sample selection is a function of the complexity of the area. China was divided into three zones (Fig. 1) to carry out classification (Fig. 5): the iterative selection was repeated ∼4-5 times, in each of the 3 zones, to improve the initial classified results. Initially, we started with a small number of samples (3 0 0) and slowly increased the sample size till we reached a high level of accuracies. After each iteration, we visually compared the classification result with sub-meter to 5-m VHRI, randomly, at 100 s of locations. If the classification results were not satisfactory, we increased the training samples till we attained satisfactory classification results. Once that was achieved, the accuracy assessment team performed the accuracies using the independent validation data ( Table 2). A product was accepted as final only when the overall, producer's and user's accuracies were adequately high (typically above 80%).

Croplands versus non croplands knowledge capture
Knowledge is very important to train a RF algorithm. First knowledge (e.g., Figs. 6a and 6b) was captured (i.e., extracted values from all bands used in this study) using training-samples (  (Figs. 2a and 2b). For illustration, we have shown knowledge distinguishing croplands, non-croplands, and pastures for all training samples for Australia plotted in a box plot (e.g., Fig. 6a). This shows which bands and periods are able to separate croplands from non-croplands. For example, Fig. 6a shows which period is more appropriate to separate croplands from non-croplands. The "Zoom-in" part of NDVI clearly shows that separability of croplands from non-croplands is maximum during period 4 and period 5, which are during the peak growing season of crops in Australia. Similarly, separability of croplands from pastures is maximum during period 3 and period 6.

Cloud computing on Google Earth Engine/processing platforms
We used GEE cloud computing for the pixel-based RF machine learning algorithm for cropland classification. The Landsat archive available on GEE is already pre-processed for atmospheric and topographic effects, which saved us much effort in data download and preprocessing time. We used JavaScript API in the GEE code editor. The zonal boundaries and all training samples were imported into GEE via Google fusion tables.

Results and discussions
There are four distinct components under this section. It starts with a presentation of the 30-m Landsat derived cropland extent product for Australia and China. This is followed by an assessment of accuracies, errors, and uncertainties of the cropland extent products. Then cropland areas, at the national and sub-national level are presented and discussed. Finally, the cropland extent product is spatially compared with other cropland products generated with coarser spatial resolution data.

30-m cropland extent product for Australia and China
The study produced a 30-m cropland extent map for the nominal year 2015 derived using Landsat 16-day time-series data for Australia (Fig. 7a) and China (Fig. 7b) by using the pixel-based supervised Random Forest (RF) machine learning algorithm for cropland classification. Initially, the RF machine learning algorithm, discussed in Section 4, were trained to separate croplands versus non-croplands for each of the zones (Fig. 1) based on knowledge generated using reference data (Section 3 and its sub-sections). The machine learning algorithm described in Section 4 and its sub-sections were then run on the GEE cloud-computing environment using Landsat data (Figs. 2a and 2b) for each of the zones to separate croplands versus non-croplands. The process was iterated and the training samples (e.g., Figs. 6a and 6b) in the algorithms were adjusted several times until we achieved optimal Reflectance/ BT/ NDVI Bands Fig. 6b. Croplands versus non-croplands knowledge capture to train random forest machine learning algorithm in Landsat bands for the Zone 2 (China) of the study area (Band 1-4 (blue), band 5-8 (green), band 9-12 (NIR), band 13-16 (red), band 17-20 (SWIR1), 21-24 (SWIR2). 25-28 (TIR1), and 29-32 (NDVI). Note: TIR1 (thermal infrared band 1) values are scaled for convenience. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 7a. Cropland Extent Product at 30-m for Australia (left panel) with illustrative zoom in view for a location (right upper) and corresponding high-resolution imagery (right lower). This product is made available for visualization at www.croplands.org. The data are downloadable from LP DAAC: https://www.lpdaac.usgs. gov/node/1282. results of croplands versus non-croplands (Figs. 7a and 7b below). The products are called the Global Food Security-support Analysis Data @ 30-m Cropland Extent for Australia (Fig. 7a), and China (Fig. 7b). This product is available through Land Processes Distributed Active Archive Center (LP DAAC https://www.lpdaac.usgs.gov/node/1282) for download  and the same dataset is also available for visualization at www.croplands.org.

Accuracy assessment
Accuracy assessment is a key component of map production using remote sensing data. The cropland extent maps for Australia, and China were assessed for their accuracies using error matrices with independent validation dataset Yadav and Congalton, 2017) that was not available to the producer of this dataset. Accuracy error matrices were produced separately for each of the zones as well as for the entire country.
Australia: A total of 900 stratified, randomly distributed validation samples were used to determine the accuracy of the final cropland extent map of Australia. An error matrix (Table 3) was generated for Australia providing producer's, user's, and overall accuracies (Congalton, 1991(Congalton, , 2009Congalton and Green, 2008;Nhamo et al., 2018).
For the whole of Australia (zone# 1), the overall accuracy was 97.6% with producer's accuracy of 98.8% (errors of omissions = 1.2%) and user's accuracy of 79.0% (errors of commissions = 21%) (Table 3) for cropland class. These results clearly imply the high level of confidence in differentiating croplands from non-croplands for Australia. High producer's accuracy from the statistical assessment using stratified, randomly distributed samples across Australia (98.8%) suggest that almost all croplands of Australia were mapped accurately (with only 1.2% missing). However, as a trade-off, there is a commission error of 21% (meaning non-cropland mapped as croplands). When using any classification this trade-off is inevitable. If we reduce errors of commissions, errors of omissions will increase and vice versa. In this trade-off, it is better to capture almost all croplands even when we have some non-croplands also mapped as croplands. This provides a more complete cropland mask for higher level products.
China: A total of 1972 validation samples, combined for all three zones in China, were used to determine the accuracy of the final cropland extent map of China. Error matrices (Table 4) were generated for each of the zones separately and for all of China providing producer's, user's and overall accuracies. For whole of China, the overall accuracy was 94% with producer's accuracy of 80% (errors of omissions = 20%) and user's accuracy of 84.2% (errors of commissions = 15.8%) (Table 4) for the cropland class. When considering all three zones (Table 4), 95.1% of total cropland area (TCA) of China is in zone #2 as shown in Table 4, 3.9% of total cropland area of China is in zone #3, and only 1% of total cropland area of China is in zone #4. Thereby, the accuracy in Zone #2 is the most important. This zone had very balanced and high producer's and user's accuracies with a producer's accuracy of 83.3% (errors of omission = 16.7%) and user's accuracy of 84.2% (errors of commission of 15.8%). The overall accuracies ranged between 91% and 98% for all zones. These results clearly imply the high level of confidence in differentiating croplands from non-croplands for China.
Within this project, the accuracy team conducted an assessment to determine the optimal number of samples required to perform accuracies. They concluded a total of 250 samples per zone as optimal, beyond which accuracies asymptote. Indeed, the total samples for training the classifiers and for accuracy assessments were significantly higher than the minimum of 250 samples. For example, accuracy assessments had: 900 samples for Australia, 1972 samples for China with the three zones in China having: 1184 for zone 2, 439 for zone 3, and 349 for zone 4. So, all zones have many more validation samples than the required optimal 250 samples.
It is possible to further increase accuracies of classification-especially the user's and producer's accuracies (Table 3 and 4) through number of measures. First, by breaking the zones (Fig. 1) to smaller ones; particularly zone 2 (eastern part of China) which contains approx.  95% of the total cropland area in China. This was not done in this study because of a requirement from the project independent accuracy team to minimize zones. Increasing the number of zones requires additional reference data to be acquired. This was not feasible within the scope of this project; however, is recommended in future work of this magnitude. Second, increased collection of reference training and validation data from agricultural fields of greater diversity (e.g., croplands along the mountains, flooded rice, fallows of various types, permanent crops of various types) is required to ensure that the classification algorithms are robust. Third, in addition to pixel-based RF MLA other pixel-based classifiers such as the support vector machines (SVMs) as well as the object-based algorithms such as the recursive hierarchical segmentation (RHSEG) will help resolve some of the uncertainties from classification results obtained from single algorithm . Fourth, the Google Earth Engine (GEE) is a powerful platform for image collection, organization, and classification over very large areas such as large countries, continents, and the entire planet through cloud computing platform. Whereas, GEE helps in powerful big-data management, and fast computing, there are several limitations on how the GEE processes the big-data using the MLAs. For example, for some classifiers, instead of classifying data of every pixel, it does "sampling" of the entire population of the big-data resulting in uncertainties in classification output and reduction in accuracies. There is also certain positional accuracy issue that sometimes can cause uncertainties on precise geographic location of a sample, especially when the pixel is 30-m resolution. To overcome this, we used a 90-m by 90-m sampling, yet, some uncertainties still need to be better understood. There are other ways of further improving accuracies and reducing uncertainties such as improved temporal resolution of imagery such as integration of global forest maps (Hansen et al., 2013) and global water mask (Carroll et al., 2009). Overall, Australia had nearly 2000 and China nearly 4000 training plus validation samples ( Table 2). Given that the 30-m Landsat data over such large areas have billions of pixels, this is not a very large sample, but probably the best achievable given the complexity of the large areas and resources required. However, larger sample for both training and validation, especially from diverse croplands ranging from high mountains to lowland deltas that take into account various subclasses within standing crops, cropland fallows, and permanent crops will help produce products with even greater accuracies.
Similar methods to those used here were also applied to other continents of the world. For example, in a recent publication on Africa, Xiong et al. (2017a) reported a weighted overall accuracy of 94%, with a producer's accuracy of 85.9% (or omission error of 14.1%), and user's accuracy of 68.5% (commission error of 31.5%) for the cropland class across Africa. Xiong et al. (2017a) used two pixel-based machine learning algorithms (RF and SVM) along with an object-based classification (RHSEG). However, mapping croplands in Africa is much more complex than mapping croplands of Australia or China due to fragmented and mixed cropping in savannas and rainforests. As a result, the user's accuracies of Africa were relatively low compared to our study in Australia and China. Further, uncertainties in cropland training and validation samples is also somewhat higher in Africa than for Australia or China.

Cropland areas and comparsion with statistics from other sources
Apart from producing a map, calculation of cropland area statistics is an important component of the 30-m cropland extent product. Australia has two distinct cropland classes (croplands and managed pastures). Total cropland areas of the two classes (croplands + pasture) were estimated as 65.7 million hectares (Mha) of which croplands alone was 35.1 Mha for the year 2015. In contrast, cropland extent area of China was 165.2 Mha. Table 5 shows country-wise cropland area statistics generated from this study for Australia and China that were compared with other sources such as the national census data based MIRCA2005 (Stefan Siebert and Portmann, personal communication;Portmann et al., 2010) which was also updated in the year 2015 (Table 5) and national level official cropland statistics from Bureau of Statistics for China and Australia.  However, 30-m cropland extent of Australia and China (www. croplands.org) produced by this study is at high resolution. Every 30m pixel is 0.09 ha. Thereby, it is feasible to capture areas at the subnational level, indeed even at individual farm level, small or big. This is a huge advantage over other existing cropland products. As a result, we calculated cropland areas of China at sub-national level for the 31 provinces of China and compared them with cropland areas from national statistics of China (http://www.stats.gov.cn/tjsj/ndsj/2014/ indexeh.htm) as shown in Fig. 8. The relation between province wise areas from both sources show very high correlation with a slope value of 1.008 and R 2 value of 0.85. Among 31 provinces, two provinces Henan and Gansu shows lower-estimation relative to national statistics. This was due to some of the mountainous agriculture that was not completely captured by the RF algorithm in our study. Two provinces in the northeastern region of Heilongjiang, and Inner Mongolia shows higher-estimation in our study relative to national statistics (Fig. 8). This was due to some of the pasture which was classified as cropland in our cropland product.

Comparison of croplands mapped at different spatial scales/resolutions
Spatially, we have compared cropland products generated at three spatial resolutions with sub-meter to 5-m VHRI. The three-different spatial resolution cropland extent products are: 1. 30-m cropland extent product (this study), 2. 250-m cropland extent product , and 3. 1-km cropland extent product (GFSAD1-km, Teluguntla et al., 2015b) Comparison of the above three products with VHRI ( Fig. 9) showed the advantage of mapping croplands at 30-m finer spatial resolution (Fig. 9). We have presented 3 examples in Fig. 9. It is obvious that the field boundaries so well captured in 30-m cropland extent product (this study) are missing in 250-m and 1000-m. Further, omissions and commissions are high in coarser 250-m and 1-km products whereas the 30-m provides near precise estimates.

Conclusions
The study demonstrated a "paradigm-shift" in deriving high resolution 30-m cropland extent products over very large areas (e.g., Australia, China) using multi-year (2013-2015), 16-day Landsat data, random forest machine learning algorithms, large volume of reference training and validation datasets from multiple sources by utilizing the big-data management and processing power of Google Earth Engine cloud-computing platform. These methods and approaches demonstrated the ability to map croplands rapidly and accurately at con-tinental\big-country scale involving petabyte volume big-data. The resulting 30-m Landsat-derived cropland extent products of Australia and China for the nominal year 2015 had overall accuracies of over 94%. The cropland class was mapped with producer's accuracy of 98.8% (errors of omissions = 1.2%) for Australia and 80% (errors of omissions = 20%) for China. The user's accuracies were 79% (errors of commissions = 21%) for Australia, and 84.2% (errors of commissions = 15.8%) for China. The study also demonstrated the ability to compute sub-national province level cropland area statistics that explained 85% variability when compared with national statistics of China. The 30-m cropland extent product mapped complex series of    Fig. 9. Spatial comparison of cropland products for precision. Comparison of cropland extent products across three different spatial resolutions (30-m, 250-m, and 1000-m) with < 5-m very high-resolution imagery.