Evaluating Landsat and RapidEye Data for Winter Wheat Mapping and Area Estimation in Punjab, Pakistan

While publicly available, cost-free coarse and medium spatial resolution satellite data such as MODIS and Landsat perform well in characterizing industrial cropping systems, commercial high spatial resolution satellite data are often preferred alternative for fine scale land tenure agricultural systems such as found in Pakistan. In this article, we integrated commercial 5 m spatial resolution RapidEye and free 30 m Landsat imagery in characterizing winter wheat in Punjab province, Pakistan. Specifically, we used 5 m spatial resolution RapidEye imagery from peak of the winter wheat growing season to derive training data for the characterization of time-series Landsat data. After co-registration, each RapidEye image was classified into wheat/no wheat labels at the 5 m resolution and then aggregated as percent cover to 30 m Landsat grid cells. We produced four maps, two using RapidEye derived continuous training data (of percent wheat cover) as input to a regression tree model, and two using RapidEye derived categorical training data as input to a classification tree model. From the RapidEye-derived 30 m continuous training data, we derived Map 1 as percent wheat per pixel, and Map 2 as binary wheat/no wheat classification derived using a 50% threshold applied to Map 1. To create the categorical wheat/no wheat training data, we first converted the continuous training data to a wheat/no wheat classification, and then used these categorical RapidEye training data to produce a categorical wheat map from the Landsat data. Two methods for categorizing the training data were used. The first method used a 50% wheat/no wheat threshold to produce Map 3, and the second method used only pure wheat (≥75% cover) and no wheat (≤25% cover) training pixels to produce Map 4. The approach of Map 4 is analogous to a standard method in which whole, pure, high-confidence training pixels are delineated. We validated the wheat maps with field data collected using a stratified, two-stage cluster design. Accuracy of the maps produced from the percent cover training data (Map 1 and Map 2) was not substantially better than the accuracy of the maps produced from the categorical training data as all methods yielded similar overall accuracies (±standard error): 88% (±4%) for Map 1, 90% (±4%) for Map 2, 90% (±4%) for Map 3, and 87% (±4%) for Map 4. Because the percent cover training data did not produce significantly higher accuracies, sub-pixel training data are not required for winter wheat mapping in Punjab. Given sufficient expertise in supervised classification model calibration, freely available Landsat data are sufficient for crop mapping in the fine-scale land tenure system of Punjab. For winter wheat mapping in Punjab and other like landscapes, training data for supervised classification may be collected directly from Landsat images without the need for high resolution reference imagery.


Introduction
Coarse and medium resolution remote sensing data have advantages over high resolution data due to their spatial coverage, temporal resolution, and availability at near real time [1].For large scale fields (>50 ha) in intensively managed commercial landscapes, Landsat has proven to be sufficient for accurate mapping of crop type [2,3].However, the use of medium spatial resolution remote sensing data for crop characterization in finer-scale land tenure systems, such as found in Punjab in Pakistan, is challenged by a supposed preponderance of mixed pixels [4] that results in higher uncertainty of area estimates and lower map accuracy [5,6].Spatial resolution (pixel size for digital data) determines the area of the smallest separate field that can be identified [7], which is significant for land tenure systems such as those found in Punjab, where the average field size is 5.6 ha [8].One way to improve the performance of medium spatial resolution data is to integrate it with high spatial resolution commercial imagery.High resolution data provide more spatially detailed observations of complicated land tenure systems, thereby potentially improving map accuracies.Areal extent of temporally frequent growing season imagery for high spatial resolution data is still limited, despite recent progress [9].More importantly, the cost of commercial data limits their general use.In this study, we tested the integration of available growing season RapidEye images (5 m) and Landsat 30 m data for Punjab province in Pakistan for winter wheat mapping.
Accurate within-season area estimation of cultivated wheat [6] is crucial to pre-harvest decision making on transportation, storage, and trade [10].The Crop Reporting Service (CRS) of the Agriculture Department of Punjab has the mandate to produce official wheat area and production data in support of policy decisions.CRS uses a labor intensive "area list frame" sampling approach to generate crop statistics months after wheat harvest [11,12], reducing the utility of the data for quick policy and market responses [10,13].An alternative approach to overcome the challenges of timely reporting [8,14,15] and accurate estimation [5,6] of wheat cultivated area is needed.Khan et al. [6] employed 30 m Landsat data to characterize wheat in Punjab using a classification tree algorithm, which yielded an overall accuracy of 76%.With the focus to improve map accuracy in the context of Punjab's small field, multiple cropping agriculture system, we evaluated the use of high spatial resolution commercial data as a training data source for Landsat-scale mapping.The objectives were to: (1) Evaluate which type of training data (continuous versus categorical) and which type of wheat map (percent wheat versus binary wheat/no wheat) produced the best accuracy; and (2) evaluate the benefit of the maps when the maps are used for the purpose of reducing standard errors of wheat area estimates derived from the field sample.The first objective focuses on the common problem of how to produce the most accurate map, whereas the second objective focuses on a practical use of a map which is to provide auxiliary information for increasing precision of area estimates.
Following good practice guidance [16], we implemented a probability sampling design to obtain in situ reference data to assess the accuracy of the wheat maps.Antecedent crop type maps were used to stratify the study area into high, medium and low strata for use in selecting the sample [2,3,6].We also used these field sample data to estimate wheat area for the 2014-2015 growing season.In our area estimation method, the maps are used as auxiliary information incorporated in what is called a "difference estimator".We posited that the more accurate map product would yield a greater reduction in standard error of the wheat area estimate derived from the field sample.This part of our study was designed to assess if sub-pixel percent wheat cover, enabled by high-resolution training data from RapidEye, would serve as better auxiliary information for purposes of reducing the standard error of the area estimate of wheat compared to auxiliary information in the form of a categorical wheat/no wheat map.

Study Area
Punjab is the most populated province of Pakistan with 56% of the total population and is the second largest in terms of land area with 205,344 km 2 (Figure 1).The province, with 12.4 million hectares (Mha) of cultivated area, is considered the bread basket for Pakistan [17].During Rabi season (December-April) of 2014-2015, official reported area of winter wheat equaled 6.98 million hectares [18].The province is divided into three major agro-ecological zones: (1) the Potohar plateau in the north with rain-fed agriculture, accounting for 10% of the total agricultural area of the province [19]; (2) the arid desert and semi-desert in the south and central region of the province with little agricultural production; and (3) the main irrigated crop growing region of Indus basin [17].Punjab province produces wheat and other staple food [5,6,20] to ensure food security of over 200 million people [21,22] and contributes significantly to region's economic development [23].The livelihoods of approximately 66% of the rural communities in Pakistan are associated with agriculture [24].
Remote Sens. 2018, 10, x FOR PEER REVIEW 3 of 21 (December-April) of 2014-2015, official reported area of winter wheat equaled 6.98 million hectares [18].The province is divided into three major agro-ecological zones: (1) the Potohar plateau in the north with rain-fed agriculture, accounting for 10% of the total agricultural area of the province [19]; (2) the arid desert and semi-desert in the south and central region of the province with little agricultural production; and (3) the main irrigated crop growing region of Indus basin [17].Punjab province produces wheat and other staple food [5,6,20] to ensure food security of over 200 million people [21,22] and contributes significantly to region's economic development [23].The livelihoods of approximately 66% of the rural communities in Pakistan are associated with agriculture [24].

Remotely Sensed Data
In this research, time-series of Landsat 30-m per pixel data and single-date 5-m per pixel RapidEye imagery were used.Given that Landsat data are freely available, all Landsat images within the Rabi winter wheat growing season of 2014-2015 were downloaded and used to generate 750 multi-temporal spectral metrics.Nine RapidEye images from the peak of the growing season were purchased from the data provider.RapidEye images cover 20% of the province area.

Landsat Data
The Landsat imagery for the Rabi growing season starting in December 2014 and ending in March 2015 were used to create a set of multi-temporal spectral metrics for mapping wheat in Punjab.A total of 307 level 1 terrain corrected (L1T) images from 20 WRS2 path/rows (Figure 2) were used including 160 Landsat 8 OLI and 147 Landsat 7 ETM+ scenes.

Remotely Sensed Data
In this research, time-series of Landsat 30-m per pixel data and single-date 5-m per pixel RapidEye imagery were used.Given that Landsat data are freely available, all Landsat images within the Rabi winter wheat growing season of 2014-2015 were downloaded and used to generate 750 multi-temporal spectral metrics.Nine RapidEye images from the peak of the growing season were purchased from the data provider.RapidEye images cover 20% of the province area.

Landsat Data
The Landsat imagery for the Rabi growing season starting in December 2014 and ending in March 2015 were used to create a set of multi-temporal spectral metrics for mapping wheat in Punjab.A total of 307 level 1 terrain corrected (L1T) images from 20 WRS2 path/rows (Figure 2) were used including 160 Landsat 8 OLI and 147 Landsat 7 ETM+ scenes.Four spectral bands, including red (0.626-0.693 µm), near infrared (NIR; 0.776-0.904µm), and short-wave infrared (SWIR1; 1.567-1.784µm and SWIR2; 2.097-2.349µm) were used as input data from the Landsat 7 ETM+ images.For the Landsat 8 OLI images, the corresponding bands were used: red (0.630-0.680 µm), NIR (0.845-0.885 µm), SWIR1 (1.560-1.660µm), and SWIR2 (2.100-2.300µm).The shorter wavelength visible blue and green spectral bands were not employed due to their greater sensitivity to atmospheric effects [25].Thermal infrared bands (ETM+ 10.40-12.50µm and TIRS 10.60-11.19µm) were used for multi-temporal spectral metrics production (see below), but were not included as variables for mapping.Normalized Difference Vegetation Index (NDVI) [26] and Normalized Difference Water Index (NDWI) [27] values were calculated for all observations as inputs to generate metrics [6], which were created using all cloud-free Landsat observations within the wheat growing season.These multi-temporal spectral metrics were designed to capture spatial and reflectance variations within the wheat growing season and facilitate large area mapping [28].
Multi-spectral data for each Landsat image were initially converted to top-of-atmosphere reflectance [29].Reflectance data were then normalized using Moderate Resolution Imaging Four spectral bands, including red (0.626-0.693 µm), near infrared (NIR; 0.776-0.904µm), and short-wave infrared (SWIR1; 1.567-1.784µm and SWIR2; 2.097-2.349µm) were used as input data from the Landsat 7 ETM+ images.For the Landsat 8 OLI images, the corresponding bands were used: red (0.630-0.680 µm), NIR (0.845-0.885 µm), SWIR1 (1.560-1.660µm), and SWIR2 (2.100-2.300µm).The shorter wavelength visible blue and green spectral bands were not employed due to their greater sensitivity to atmospheric effects [25].Thermal infrared bands (ETM+ 10.40-12.50µm and TIRS 10.60-11.19µm) were used for multi-temporal spectral metrics production (see below), but were not included as variables for mapping.Normalized Difference Vegetation Index (NDVI) [26] and Normalized Difference Water Index (NDWI) [27] values were calculated for all observations as inputs to generate metrics [6], which were created using all cloud-free Landsat observations within the wheat growing season.These multi-temporal spectral metrics were designed to capture spatial and reflectance variations within the wheat growing season and facilitate large area mapping [28].
Multi-spectral data for each Landsat image were initially converted to top-of-atmosphere reflectance [29].Reflectance data were then normalized using Moderate Resolution Imaging Spectroradiometer (MODIS) top-of-canopy reflectance data composite as a normalization target [28].To create the cloud-free growing season MODIS reference data, all 16-day MODIS composites from 2000 through 2011 were ranked by NDVI value.Using 16-day composites corresponding to NDVI ranks between 50th and 90th percentile, we calculated mean surface reflectance for each spectral band.For individual Landsat image the per-pixel difference between top-of-atmosphere spectral reflectance and MODIS surface reflectance composite was used to (i) apply reflectance bias adjustment, which largely accounted for atmospheric scattering; and (ii) apply a cross-track adjustment to account for effects of surface anisotropy.To perform the cross-track adjustment, the reflectance bias was modelled as a function of sensor view angle per band; the derived relationship was then to apply bias adjustment to all pixels of the Landsat image.The radiometric normalization algorithm reduced between and within image reflectance differences caused by variation in atmospheric conditions, dates and surface reflectance anisotrophy [30].Using the method from [28], quality assessment codes were assigned to each pixel to reflect its probability to be a cloud-free land or water observation.
All viable observations within the defined study period were aggregated and used to derive a set of spectral-statistical derivations called multi-temporal spectral metrics.Metrics represent a generic feature space that facilitate large area mapping and have been used extensively with Advanced Very High Resolution Radiometer (AVHRR) and MODIS data [31,32] and more recently with Landsat data [6,28].Landsat-based metrics are calculated between a start and end date without direct relation to the day of the year, in this case the Rabi growing season of December, January and February (Figure 3).Multi-temporal spectral metrics are statistical derivations of all good quality assessed pixels and have been shown to capture salient phenological information for mapping land cover [33,34].Landsat-based metrics enable mapping the same crop type within a broad geographic context and over large regions [32,35].Using metrics instead of single-date images allowed us to create a complete, wall-to-wall dataset without observation gaps, which is well-suited for regional classification.To create a set of multi-temporal spectral metrics, we ranked all cloud-free Landsat observations per pixel corresponding to (i) band reflectance value; (ii) NDVI value; (iii) NDWI value; and (iv) brightness temperature value.For each spectral band and ranking method, minimal, 10th, 25th, 50th, 75th, and 90th percentile, and maximal values were recorded as a set of metrics.In addition, we calculated and recorded mean reflectance for all values between minimum and 10th percentile, minimal and 25th, 10th and 25th, 25th and 50th, 50th and 75th, 25th and 75th, 10th and 50th, 75th and maximal, 90th percentile and maximal, 10th and 90th percentile, and minimal and maximal values for each band and ranking.Multi-temporal spectral metrics allowed us to collect spectral information corresponding to major phenological stages (e.g., peak of growing season, warmest time of year, etc.).In addition to metrics, time-sequential monthly composites for December to March were created based on a median value taken from all cloud/shadow-free observations within each calendar month to facilitate image interpretation and assignment of training data.

RapidEye Data
RapidEye, a constellation of five satellites, delivers high spatial resolution data [36] that enable fine scale mapping [37].RapidEye 5 m spatial resolution 1B images acquired between January and March 2015 were procured.From available RapidEye data collections during the Rabi season of 2015, nine images within the peak and end of the growing season were used in this study (Table 1/Figure 4).Images collected during the peak of growing season were most useful for crop type interpretation.Late growing season images provided enough information for crop identification even though wheat crops were ripe in the southern part of the province.Images outside of the mid-January to late March interval were not useful in separating wheat from other land cover [38].The RapidEye data were resampled to a 5 m raster grid nested to Landsat 30 m pixels.All spectral bands (blue, 0.44-0.51µm; green, 0.52−0.59µm; red, 0.63−0.69µm; red edge, 0.69−0.73µm; and NIR, 0.76−0.85µm) were used as inputs to classification.Three spectral bands (red, red edge and NIR) were used for visual interpretation of wheat crops.The red edge band sensitivity to chlorophyll The RapidEye data were resampled to a 5 m raster grid nested to Landsat 30 m pixels.All spectral bands (blue, 0.44-0.51µm; green, 0.52−0.59µm; red, 0.63−0.69µm; red edge, 0.69−0.73µm; and NIR, 0.76−0.85µm) were used as inputs to classification.Three spectral bands (red, red edge and NIR) were used for visual interpretation of wheat crops.The red edge band sensitivity to chlorophyll content [39] was a useful input in facilitating wheat cover identification [40].The RapidEye data were used to derive per pixel percent wheat area and per pixel categorical wheat/no wheat training data for the respective regression and classification tree analyses.
content [39] was a useful input in facilitating wheat cover identification [40].The RapidEye data were used to derive per pixel percent wheat area and per pixel categorical wheat/no wheat training data for the respective regression and classification tree analyses.

Topographic Data
Topographic data were also used as independent variables, including elevation, slope and aspect.These layers were computed from 30 m spatial resolution void-filled seamless Digital Elevation Model (DEM) derived from the National Aeronautics and Space Administration (NASA) Shuttle Radar Topography Mission (SRTM) (downloaded from CGIAR-CSI: http://srtm.csi.cgiar.org).

Wheat Mapping: Regression and Classification Trees
We related the 5 m RapidEye multispectral data to visually derived categorical training data via a bagged decision tree algorithm to obtain wheat/no wheat maps at the 5 m spatial resolution within each image footprint (Figure 5).The individual 5 m RapidEye classifications were aggregated and resampled to 30 m spatial resolution by calculating the percent of 5 m wheat pixels within each 30 m

Topographic Data
Topographic data were also used as independent variables, including elevation, slope and aspect.These layers were computed from 30 m spatial resolution void-filled seamless Digital Elevation Model (DEM) derived from the National Aeronautics and Space Administration (NASA) Shuttle Radar Topography Mission (SRTM) (downloaded from CGIAR-CSI: http://srtm.csi.cgiar.org).

Wheat Mapping: Regression and Classification Trees
We related the 5 m RapidEye multispectral data to visually derived categorical training data via a bagged decision tree algorithm to obtain wheat/no wheat maps at the 5 m spatial resolution within each image footprint (Figure 5).The individual 5 m RapidEye classifications were aggregated and resampled to 30 m spatial resolution by calculating the percent of 5 m wheat pixels within each 30 m output cell.The 30 m Landsat multi-temporal spectral metrics were related to 30 m RapidEye-derived categorical and percent wheat training data via respective bagged classification and regression tree models [6,28,[41][42][43].The classification and regression tree software implemented by our research group follows the algorithm described by [44].Classification trees employ an entropy measure, referred to as deviance, to split a multi-dimensional space of dependent variables into successively more homogeneous subsets [45].In regression trees, a sum of squares criterion is used to split the independent variables into successively less varying subsets of a continuous variable [45].The best univariate split that yields the maximum reduction in entropy or variability is selected to build the classification or regression tree model, respectively [45].We terminated each classification or regression tree when additional splits decreased model deviance or variability by less than 0.001 of the deviance of the total training set population [43,[46][47][48].To further avoid overfitting, we used a set of seven bagged tree models, each derived from a 20% random sample of training pixels.Each classification tree yielded a per pixel probability of wheat cover class membership, and each regression tree yielded a per pixel percent wheat cover estimate (Figure 5).For both the classification and regression tree models, the per-pixel median of the seven model outputs was selected as the final result [28].To derive the categorical wheat classification (from the regression tree output), the pixel was labeled as wheat if the median value from the seven regression tree outputs was equal to or greater than 50%.
Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 21 output cell.The 30 m Landsat multi-temporal spectral metrics were related to 30 m RapidEye-derived categorical and percent wheat training data via respective bagged classification and regression tree models [6,28,[41][42][43].The classification and regression tree software implemented by our research group follows the algorithm described by [44].Classification trees employ an entropy measure, referred to as deviance, to split a multi-dimensional space of dependent variables into successively more homogeneous subsets [45].In regression trees, a sum of squares criterion is used to split the independent variables into successively less varying subsets of a continuous variable [45].The best univariate split that yields the maximum reduction in entropy or variability is selected to build the classification or regression tree model, respectively [45].We terminated each classification or regression tree when additional splits decreased model deviance or variability by less than 0.001 of the deviance of the total training set population [43,[46][47][48].To further avoid overfitting, we used a set of seven bagged tree models, each derived from a 20% random sample of training pixels.Each classification tree yielded a per pixel probability of wheat cover class membership, and each regression tree yielded a per pixel percent wheat cover estimate (Figure 5).For both the classification and regression tree models, the per-pixel median of the seven model outputs was selected as the final result [28].To derive the categorical wheat classification (from the regression tree output), the pixel was labeled as wheat if the median value from the seven regression tree outputs was equal to or greater than 50%.

RapidEye-Derived Training Data
The nine RapidEye images procured for Punjab Pakistan from January, February and March were visualized as false color composites of NIR, red edge and red spectral bands.The training data for the RapidEye images consisted of manually labelled wheat and non-wheat pixels with visual image-interpretation of each image based on spectral properties and landscape context [28,30].The

Landsat-Derived Map Products
The 30 m continuous (Figure 6) and categorical training data derived from the RapidEye imagery were used to create four wheat maps for Punjab (Figure 5).Map 1, the 30 m percent wheat per pixel map produced by relating the RapidEye derived 30 m continuous training data to 30 m Landsat multi-temporal spectral metrics.Map 2 was a 30 m binary wheat/no wheat classification derived from

Landsat-Derived Map Products
The 30 m continuous (Figure 6) and categorical training data derived from the RapidEye imagery were used to create four wheat maps for Punjab (Figure 5).Map 1, the 30 m percent wheat per pixel map produced by relating the RapidEye derived 30 m continuous training data to 30 m Landsat multi-temporal spectral metrics.Map 2 was a 30 m binary wheat/no wheat classification derived from the percent wheat map (Map 1) with a per pixel threshold labelled as wheat if the percent wheat value exceeded 50%.From categorical training data, we produced two wheat/no wheat classifications.Map 3 was a 30 m binary classification derived using training data based on a 50% wheat/no wheat threshold defined for all RapidEye-derived training pixels at 30 m scale.Map 4, was a second binary classification map produced using only 30 m percent wheat RapidEye-derived training pixels that were either ≥75% (wheat) or ≤25% (no wheat) (Figure 7).That is, for Map 4 we eliminated those pixels from the training data that were ambiguous in regard to whether wheat or not wheat.Map 4 represents a map derived from "pure pixel" training data.The wheat area for Punjab for each map was calculated simply by summing the area proportion of each pixel mapped as wheat in the case of the percent wheat map or summing the area of all pixels mapped as wheat in the case of the categorical wheat/no wheat maps.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 21 the percent wheat map (Map 1) with a per pixel threshold labelled as wheat if the percent wheat value exceeded 50%.From categorical training data, we produced two wheat/no wheat classifications.Map 3 was a 30 m binary classification derived using training data based on a 50% wheat/no wheat threshold defined for all RapidEye-derived training pixels at 30 m scale.Map 4, was a second binary classification map produced using only 30 m percent wheat RapidEye-derived training pixels that were either ≥75% (wheat) or ≤25% (no wheat) (Figure 7).That is, for Map 4 we eliminated those pixels from the training data that were ambiguous in regard to whether wheat or not wheat.Map 4 represents a map derived from "pure pixel" training data.The wheat area for Punjab for each map was calculated simply by summing the area proportion of each pixel mapped as wheat in the case of the percent wheat map or summing the area of all pixels mapped as wheat in the case of the categorical wheat/no wheat maps.

In Situ Data for Area Estimation and Map Validation
To estimate map accuracy and wheat area from the reference classification, an independent in situ data set was collected using a stratified two-stage cluster sample.An initial population of 2261 blocks, (each block 10 km × 10 km) covering all of Punjab was created.Blocks with less than 25% of their area inside the provincial boundary of the Punjab were excluded resulting in 105 blocks being removed.Each block was assigned to one of three strata based on wheat area as determined from a previous Landsat-derived wheat classification for the Punjab produced by [6].These three strata were employed: low wheat (blocks with 0-22% wheat), medium wheat (22-47%) and high wheat cover (>47%).These strata boundaries were selected so that the low wheat stratum represented half of the population (i.e., the low wheat boundary was the median wheat percent cover of the population of blocks), the high wheat stratum included the blocks with the greatest wheat cover so that the total area of wheat in the high stratum represented 50% of all wheat area of the map used to create the stratification, and the medium wheat stratum had all remaining blocks not assigned to the low or high wheat strata.The purpose of the stratification was to form strata that were internally relatively homogeneous but with blocks in different strata having different wheat percent cover and to capture the diversity in wheat crop across the province.From the population of 2156 blocks, 1035 were assigned to the low stratum, 680 to the medium stratum, and 441 to the high stratum.
Considering logistical constraints in field data collection and recognizing spatial and phenological variation in wheat crop across Punjab, seven sample 10 km × 10 km blocks from each stratum were randomly selected.The rationale for allocating equal effort to the three strata was to quantify the full range of different wheat densities, particularly given the small total sample size due to the high cost of field data collection.If proportional allocation had been used, the low wheat stratum would have received the largest sample size.We did not believe the effort spent in the low stratum would have been beneficial.Within each selected sample block, we used simple random sampling to select ten 30 m pixels for field visit (Figure 8).The sample pixels that fell outside of the provincial boundary (5 sample pixels total) were not visited resulting in a total of 205 sample pixels visited of the 210 sample pixels selected from the 21 sample blocks.
For each sampled pixel, we quantified the proportion of wheat area in the field.For non-cropland land cover, we incorporated information from Google Earth™ in interpreting cover such as trees and built-up area.The percentage of wheat in each of the sampled pixels, as observed in the field, was recorded.Error matrices for the four maps were calculated using wheat percent for the sampled pixels from the field (i.e., the reference data) compared to the continuous and binary wheat characterizations (i.e., the wheat maps).Because the reference value was the proportion of wheat area within each sample pixel, the error matrix cell entries were quantified differently from the more typical case in which the reference data represent a hard classification (i.e., proportion of wheat is 0 or 1).The formulas for estimating accuracy of the different maps and for estimating area of wheat based on the field interpretations are those presented as appendix in Potapov et al. [50].
Olofsson et al. [16] noted the value of using the reference sample data as the basis for estimating area.We examined two approaches for estimating wheat area from the field sample data.In the first approach, we estimated area using only the sample interpretations (i.e., we did not use information from the four wheat maps in the estimator of area).This approach is sometimes called a "direct estimator" because the estimate is directly from only the field sample data.However, it is generally possible to reduce the standard error of the area estimate relative to the direct estimator by incorporating information from the wheat maps.In this study, we used a "difference estimator" [51] for this purpose.The difference estimator is applied to the stratified two-state cluster sampling design and (Appendix in Potapov et al. [50]) provides the specific formulas for the difference estimator and its estimated standard error.The starting point for the difference estimator is the total area mapped as wheat from one of the four maps produced.Because the wheat map has classification error, the total area based on the map is expected to be biased.We can remove that bias using the reference sample data, specifically the difference between the proportion of area of wheat (for each sampled pixel) from the map and the proportion of area according to the reference data.A sample estimate of the bias adjustment of mapped wheat area is then obtained.Because the complete coverage wheat map is used in the difference estimator of wheat area, we expect the map information to translate into smaller uncertainty of the wheat area estimate.The standard error of the difference estimator or wheat area provides quantitative evidence for the reduction in uncertainty achieved.We produced four area estimates using the difference estimator, one for each of the four wheat maps (the reference sample data are the same for all four estimates).The four area estimates produced using the difference estimator will depend on the wheat map incorporated in the difference estimator.The standard errors of the estimates indicate how much benefit was gained by using each map in the difference estimator of area, and we can also assess the value of the maps by comparing the standard error of the difference estimates to the standard error of the direct estimate.Smaller standard errors are preferable.
Remote Sens. 2018, 10, x FOR PEER REVIEW 12 of 21 map is used in the difference estimator of wheat area, we expect the map information to translate into smaller uncertainty of the wheat area estimate.The standard error of the difference estimator or wheat area provides quantitative evidence for the reduction in uncertainty achieved.We produced four area estimates using the difference estimator, one for each of the four wheat maps (the reference sample data are the same for all four estimates).The four area estimates produced using the difference estimator will depend on the wheat map incorporated in the difference estimator.The standard errors of the estimates indicate how much benefit was gained by using each map in the difference estimator of area, and we can also assess the value of the maps by comparing the standard error of the difference estimates to the standard error of the direct estimate.Smaller standard errors are preferable.In addition to the sample-based validation of the wheat maps based on the in situ reference data, we performed an inter-comparison of wheat area estimates with the wheat crop forecasts of Punjab Crop Reporting Service (CRS) using their crop reporting system [18] and wheat area estimates from Space and Upper Atmosphere Research Commission (SUPARCO) [52].Provincial and district level wheat area data from the CRS were compared to our wheat maps for district-level area comparison.These estimates were used as auxiliary information available for validation of our mapped wheat area estimates for the province.

Intercomparison of Wheat Area Estimates
The direct sample estimate of the total area of wheat in the Punjab based on the field sample data was 7.17 Mha for the 2014-2015 Rabi growing season.By comparison, the total area of cultivated wheat in Punjab estimated from the four difference estimators (i.e., a difference estimator was produced using each of the four wheat maps as the auxiliary data) varied from 7.25 Mha to 7.76 Mha (Figure 9), indicating that the general effect of the difference estimator was to increase the estimated area relative to the direct estimator that incorporated no map information.Although we have a small sample size, both the Landsat maps and CRS wheat estimates were within the 95% confidence interval of the field-based estimate (Figure 9).Of all map-based area estimates, Landsat trained with first and last quartile RapidEye classification data was closest to the CRS estimate (Table 2).
The difference estimator increased the total wheat area relative to the direct estimator, and with one exception, the difference estimator reduced the standard error relative to the direct estimator.Incorporating the wheat map information via the difference estimator reduced the standard error of the area estimate to approximately 80% of the standard error of 1.08 Mha obtained for the direct estimate that does not incorporate the wheat map information.Consequently, a substantial reduction in standard error can be achieved by incorporating information from the wheat map into the sample-based estimate of wheat area.
Remote Sens. 2018, 10, x FOR PEER REVIEW 13 of 21 In addition to the sample-based validation of the wheat maps based on the in situ reference data, we performed an inter-comparison of wheat area estimates with the wheat crop forecasts of Punjab Crop Reporting Service (CRS) using their crop reporting system [18] and wheat area estimates from Space and Upper Atmosphere Research Commission (SUPARCO) [52].Provincial and district level wheat area data from the CRS were compared to our wheat maps for district-level area comparison.These estimates were used as auxiliary information available for validation of our mapped wheat area estimates for the province.

Intercomparison of Wheat Area Estimates
The direct sample estimate of the total area of wheat in the Punjab based on the field sample data was 7.17 Mha for the 2014-2015 Rabi growing season.By comparison, the total area of cultivated wheat in Punjab estimated from the four difference estimators (i.e., a difference estimator was produced using each of the four wheat maps as the auxiliary data) varied from 7.25 Mha to 7.76 Mha (Figure 9), indicating that the general effect of the difference estimator was to increase the estimated area relative to the direct estimator that incorporated no map information.Although we have a small sample size, both the Landsat maps and CRS wheat estimates were within the 95% confidence interval of the field-based estimate (Figure 9).Of all map-based area estimates, Landsat trained with first and last quartile RapidEye classification data was closest to the CRS estimate (Table 2).
The difference estimator increased the total wheat area relative to the direct estimator, and with one exception, the difference estimator reduced the standard error relative to the direct estimator.Incorporating the wheat map information via the difference estimator reduced the standard error of the area estimate to approximately 80% of the standard error of 1.08 Mha obtained for the direct estimate that does not incorporate the wheat map information.Consequently, a substantial reduction in standard error can be achieved by incorporating information from the wheat map into the samplebased estimate of wheat area.

Accuracy of the Wheat Maps
The four wheat maps produced from the RapidEye training data are shown in Figure 10.The overall accuracies of the four wheat maps were relatively similar ranging from 87% to 90% (Table 2).

Accuracy of the Wheat Maps
The four wheat maps produced from the RapidEye training data are shown in Figure 10.The overall accuracies of the four wheat maps were relatively similar ranging from 87% to 90% (Table 2).User's accuracies of wheat were generally higher than producer's accuracies for wheat for the four map products indicating that omission error of wheat was generally more problematic than commission error of wheat.
District-level wheat area from the CRS final wheat report for 2014-2015 was compared to our map products in Figure 11.The R 2 values for our four map products do not reveal a performance difference, whether using sub-pixel training and percent wheat mapping (Map 1) or homogeneous wheat/no wheat training for wheat classification (Map 4).Results suggest that commercially acquired high resolution data as a training input have no advantage over classical methods of training pure pixels into categorical classes.For most districts, the map estimates agreed with the CRS estimates except for the Rahim Yar Khan district, which is a desert zone.The map products showed higher wheat acreage than the CRS estimate for this district.

Discussion
The four wheat maps we produced for Punjab, demonstrated relatively similar overall, user's and producer's accuracies supporting the idea that 30 m wheat/no wheat mapping was a viable approach, and a key advantage of such an approach is that it can be implemented without the expense or processing of high resolution commercial data.We employed the same training data pool in all cases, from percent wheat per pixel training data to a subset of only pure yes/no training data (≤25% and ≥75% test).Our other key objective was to determine the benefit of the different wheat maps when used for the purpose of reducing the standard error of sample-based wheat area estimate.Specifically, we compared the standard errors when each map was incorporated in a difference estimator of wheat area.The pure pixel map (i.e., Map 4) combined with the difference estimator yielded a standard error of 0.85 Mha, which was nearly the same as the standard error of 0.82 Mha achieved by the percent wheat map used in the difference estimator.For the purpose of reducing the standard error of the area estimate, the percent wheat map was not substantially better than the wheat/no wheat map, so based on this case study, the evidence does not indicate additional value is gained by producing the percent wheat map.Moreover, incorporating the wheat/no wheat map in the difference estimator of wheat area substantially reduced the standard error relative to the standard error (1.08 Mha) of the direct estimator (the direct estimator is what would be used if we had no wheat map).Therefore, we can conclude that an image analyst with experience in identifying winter wheat cover could create a wheat/no wheat map that would yield a substantial reduction (nearly 20%) in standard error for the area estimate of wheat based in the field sample data.
An analysis of specific map errors based on comparison to field data revealed some of the complexity of the Punjab agricultural landscape.Many of the omission errors were for wheat under orchards, which is a fundamental limitation, not easily overcome using satellite data.Other omission errors were in rain-fed wheat production areas of semi-arid environments due to the lack of dense wheat cover and a higher contribution of soil background reflectance.Commission errors were associated with other winter crops, such as clover [6].While landscape heterogeneity impacted map accuracies less, factors such as crop type, intercropping, and intensification also reduced overall accuracies.It is these kinds of mapping limitations that require the combined use of the map with field data to achieve precise area estimates.
Figure 12 shows the three metrics that contributed most to the decision tree models estimating wheat cover, illustrating the importance of longer wavelengths in characterizing wheat.Figure 13 shows full-resolution subsets of RapidEye and Landsat input data and wheat maps.Differences in spatial resolution between RapidEye and Landsat imagery are clearly evident.However, the percent wheat map products derived from the aggregated 5 m wheat/no wheat training data exhibited a discrete appearance made up of largely 100% wheat and 100% no wheat at 30 m training cells.Results from Landsat (Figure 13e-h) were visually congruent with the training data and each other.Given that Punjab is characterized by relatively fine-scale land tenure, there is an expectation that mixed pixels would challenge crop characterization using medium spatial resolution data such as Landsat (30 m).For example, SUPARCO, the Space and Upper Atmosphere Research Commission of Pakistan, employs high spatial resolution SPOT 5 m data for mapping crops due to this fact.However, our results show that winter wheat is a fairly homogeneous land cover target.When examining spatial contiguity of 100% wheat pixels from our percent wheat map, we found an effective field size of 42 ha.The spatial contiguity of cultivated winter wheat was also evident in the 5 m RapidEye data.Figure 14 illustrates the histogram of 5 m wheat/no wheat training data aggregated to 30 m spatial resolution, compared to predicted 30 m percent wheat.Despite the fine-scale land tenure, winter wheat was such a dominant crop that 30 m data captured it as a largely wheat/no wheat categorical (binary) theme.Crop diversity in the Kharif (summer monsoon) season is much greater and may require sub-pixel or finer-scale categorical characterizations.However, for winter wheat, Landsat was a sufficient input.
The demonstrated approach could easily complement the conventional "Village List Frame" method used by the Crop Reporting Service of Punjab by providing an interim assessment in lieu of the more field-intensive and delayed CRS reporting.Less intensive field sampling, guided by an in-season wheat/no wheat map could provide timely information on wheat production and support key decisions regarding wheat management, transportation and storage.The delayed results from the standard CRS method preclude such timely decision making.Importantly, such a capability can be implemented using freely available Landsat data, making the approach a cost-effective improvement for pre-harvest policy decision making.
Remote Sens. 2018, 10, x FOR PEER REVIEW 17 of 21 the more field-intensive and delayed CRS reporting.Less intensive field sampling, guided by an inseason wheat/no wheat map could provide timely information on wheat production and support key decisions regarding wheat management, transportation and storage.The delayed results from the standard CRS method preclude such timely decision making.Importantly, such a capability can be implemented using freely available Landsat data, making the approach a cost-effective improvement for pre-harvest policy decision making.

Conclusions
Wheat maps derived from Landsat data for the 2014-2015 Rabi growing season in the Punjab province of Pakistan correspond closely to official statistics and field validation data.The high accuracies of the wheat maps supports the utility of the maps for depicting the spatial distribution of wheat and potential benefits to the maps for reducing standard error of area estimate when the wheat maps are combined with stratified sampling of field data.For many landscapes characterized by small field sizes and fine-scale land tenure, an evaluation of freely available medium spatial resolution data should be made before presuming that commercial high resolution data are required.Calibration of mapping algorithms (training data) can likely be delineated directly without the need for high resolution reference imagery.Other landscapes where a single crop dominates, such as soybean in northeast China, or paddy rice in parts of Southeast Asia, may be similarly characterized.The value of our approach is in the comparatively low level of resources needed to implement the method and the ability to derive pre-harvest crop area estimates.For Pakistan, such results are important for improving wheat crop management, storage and transportation, as well as decisions on exports and imports of wheat grain in addressing food security issues.

Figure 1 .
Figure 1.Map of the study area, Punjab Province (Pakistan).

Figure 1 .
Figure 1.Map of the study area, Punjab Province (Pakistan).

Figure 2 .
Figure 2. Landsat data volumes used in this study.

Figure 2 .
Figure 2. Landsat data volumes used in this study.

Figure 5 .
Figure 5. Flow chart of the research.

Figure 5 .
Figure 5. Flow chart of the research.

Figure 7 .
Figure 7. Percent wheat training data of Figure 6 thresholded at ≤25% and ≥75% to derive categorical wheat/no wheat training data set.

Figure 7 .
Figure 7. Percent wheat training data of Figure 6 thresholded at ≤25% and ≥75% to derive categorical wheat/no wheat training data set.

Figure 8 .
Figure 8. Stratified random sample blocks and sample pixels from Punjab.Figure 8. Stratified random sample blocks and sample pixels from Punjab.

Figure 8 .
Figure 8. Stratified random sample blocks and sample pixels from Punjab.Figure 8. Stratified random sample blocks and sample pixels from Punjab.

Figure 9 .
Figure 9. Wheat area estimates from field data, the four RapidEye-trained Landsat-derived map products, and official data from CRS and SUPARCO.For the Landsat-derived map products, both the map pixel counts and difference estimator results are shown.Uncertainties of ±1 standard error are shown for the area estimates derived from the probability sample of field (in situ) data.

Figure 9 .
Figure 9. Wheat area estimates from field data, the four RapidEye-trained Landsat-derived map products, and official data from CRS and SUPARCO.For the Landsat-derived map products, both the map pixel counts and difference estimator results are shown.Uncertainties of ±1 standard error are shown for the area estimates derived from the probability sample of field (in situ) data.

Figure 10 .
Figure 10.Results of the four wheat maps classified from Landsat time-series data and using 30 m training data derived from RapidEye.(Map 1) Percent wheat map; (Map 2) 50% thresholded percent wheat map; (Map 3) Binay wheat map derived with RapidEye categorical training; (Map 4) binary wheat map derived with training ≤25% and ≥75%.

Figure 11 .
Figure 11.district-level wheat area from Landsat-basedwheat maps and CRS wheat area estimates for Punjab.Figure 11.District-level wheat area from Landsat-basedwheat maps and CRS wheat area estimates for Punjab.

Figure 11 .
Figure 11.district-level wheat area from Landsat-basedwheat maps and CRS wheat area estimates for Punjab.Figure 11.District-level wheat area from Landsat-basedwheat maps and CRS wheat area estimates for Punjab.

Figure 12 .
Figure 12.Composite of the top three Landsat metrics used to estimate wheat cover: mean of SWIR (1.6 µm) reflectance values between the 10th and 50th percentile ranked by thermal reflectance values (in Red), mean of NIR reflectance values between the 10th and 50th percentile ranked by thermal reflectance values (in Green), and the mean of the last three NDWI reflectance values (in Blue).

Figure 12 .
Figure 12.Composite of the top three Landsat metrics used to estimate wheat cover: mean of SWIR (1.6 µm) reflectance values between the 10th and 50th percentile ranked by thermal reflectance values (in Red), mean of NIR reflectance values between the 10th and 50th percentile ranked by thermal reflectance values (in Green), and the mean of the last three NDWI reflectance values (in Blue).

Figure 13 .
Figure 13.Top row: (a) RapidEye false color composite of red-NIR-red edge data; (b) wheat/no wheat 5 m map (green = wheat, yellow = no wheat); and (c) percent wheat from 5 m map aggregated to 30 m spatial resolution.Middle row: (d) Landsat false color composite of red-NIR-SWIR (1.6 µm);(e) percent wheat from Landsat; and (f) wheat/no wheat map of percent wheat thresholded at 50%.Bottom row: (g) wheat/no wheat map from yes/no RapidEye training thresholded at 50%; and (h) wheat/no wheat map from yes/no RapidEye training thresholded at ≤25% and ≥75%.

Figure 14 .
Figure 14.Histogram of RapidEye-derived wheat training data aggregated to percent cover at a 30 m spatial resolution and Landsat-derived percent wheat cover for the same grid cells used as training data.

Figure 14 .
Figure 14.Histogram of RapidEye-derived wheat training data aggregated to percent cover at a 30 m spatial resolution and Landsat-derived percent wheat cover for the same grid cells used as training data.Figure 14.Histogram of RapidEye-derived wheat training data aggregated to percent cover at a 30 m spatial resolution and Landsat-derived percent wheat cover for the same grid cells used as training data.

Figure 14 .
Figure 14.Histogram of RapidEye-derived wheat training data aggregated to percent cover at a 30 m spatial resolution and Landsat-derived percent wheat cover for the same grid cells used as training data.Figure 14.Histogram of RapidEye-derived wheat training data aggregated to percent cover at a 30 m spatial resolution and Landsat-derived percent wheat cover for the same grid cells used as training data.

Table 1 .
RapidEye data used in development of winter wheat training data.

Table 1 .
RapidEye data used in development of winter wheat training data.

Table 2 .
Comparison of accuracies of four map products, two derived using continuous training data and two derived using wheat/no-wheat training data, along with final area estimates.

Table 2 .
Comparison of accuracies of four map products, two derived using continuous training data and two derived using wheat/no-wheat training data, along with final area estimates.