Threshold model for detecting transparent plastic-mulched landcover using moderate-resolution imaging spectroradiometer time series data: a case study in southern Xinjiang, China

Abstract. Plastic-mulched landcover (PML) is the land surface covered by thin plastic films. PML has been expanding rapidly worldwide and has formed a significant agriculture landscape in the last two decades. Large-scale PML may impact the regional and global climate, ecosystem, and environment because it changes the energy balance and water cycles of the land surfaces, reduces the biodiversity, deteriorates the soil structure, etc. To study its impact, the spatial and temporal distributions and dynamics of PML have to be obtained. This paper presents a threshold model (TM) for PML detection and mapping with moderate-resolution imaging spectroradiometer (MODIS) time series data. Based on the temporal-spectral features of PML in the early stage of a growing season after planting, a TM was designed with the number of days (d) when the normalized difference vegetation index (NDVI) value is larger than a threshold value (x) as the discriminator. The model has been successfully applied to map PML in southern Xinjiang, China, from the interpolated MODIS NDVI time series (from 90th to 125th day of each year). Results indicate that when TM parameter x is set to 0.2 and d to 8, the overall accuracy and kappa coefficient (κ) are >0.84 and 0.65, respectively. We believe this classification accuracy can meet the PML mapping for large geographic areas. Furthermore, visual comparison between the PML maps from TM classification of MODIS time series and that from the maximum likelihood classification of Landsat ETM+ and OLI images shows they are consistent both in the pattern and location of PML. Therefore, detection and mapping of PML by using MODIS time series with the TM method is feasible. The PML mapping in this study used a cropland mask derived from Landsat images using a maximum likelihood classifier to mask out non-cropland when applying the TM algorithm. The accuracy of such a mask is subject to further study. Because of frequent global coverages of MODIS data, the method presented in this paper could potentially be used for PML detecting and mapping at continental and global scales.


Introduction
Plastic covered farmland, or plasticulture, in the broad sense, is defined as the use of plastics in agriculture. 1,2,3,4The purposes of using plastic cover in agriculture are to mitigate the threats of insects, crop diseases, coldness, heat, drought, and strong rainfall and wind, and to improve the productivity. 4,5Plastic-mulched landcover (PML), greenhouse or walk-in tunnels, and small tunnels are the three main types of plasticulture.This paper only deals with PML, the dominate type of plasticulture in term of area coverage.
China is one of the major nations in the world that applies plastic mulch on a large scale.According to recent statistics, China used 1760 kiloton (kt) of thin plastic film as mulch in 2005, 2200 kt in 2010, and expects it to reach 2600 kt in 2015, and the PML area in China has grown from 4200 ha in 1981, 11,966,873 ha in 2003, to 28,000,000 ha in 2010, ∼17% of the total farmland. 6,7The Ministry of Agriculture of China predicted that the PML area will reach 33,000,000 ha in 2015. 7The PML area is also growing significantly in Italy, 8 Israel, 9 Spain, 10 and North America. 11Plasticulture areas have been expanding rapidly worldwide and now represent an important cultivated landscape. 5uch a large-scale landcover change must impact the climate, ecosystem, and environment because plasticulture changes the energy balance and water cycles of the land surfaces, reduces the biodiversity, deteriorates the soil structure, etc. 5 Yet studies on the impact of such a landcover change are rare because the information about the spatial and temporal distributions of the plasticulture area around the world is unavailable.The first step toward the understanding of the impact will be the mapping of the spatial and temporal distribution of plasticulture.Remote sensing is the only available approach for mapping plasticulture over a large geographic region (e.g., all of China, or all of North America).
Recent studies of mapping plasticulture mainly concentrate on high spatial resolution (at the meter level) remote sensing imagery to extract greenhouse information.For example, Levin et al. 9 claimed that there were three absorption bands (centered at 1218, 1732, and 2313 nm, respectively) for clear plastic, which used 1 m resolution AISA-ES hyperspectral image data to achieve a >90% detection accuracy for clear plastic cover and an ∼70% for black plastic cover.Agüera et al. 10,12 used the maximum likelihood (ML) classification method to extract greenhouse locations from both Quickbird and IKONOS images and achieved satisfactory detection accuracy.Tarantino and Fugorito 13 successfully mapped rural areas with widespread plastic covered (tunnel/greenhouse) vineyards using both spectral information and spatial texture from very high spatial resolution true color aerial images.Carvajal et al. 14 developed an artificial intelligence neural network to identify greenhouses using 2.44 m resolution QuickBird imagery.
A few studies of mapping plasticulture based on medium-resolution (at a tens of meters level) imagery have also been conducted.Colby and Keating 15 attempted to classify greenhouse areas using the parallelepiped method and ML from Landsat TM images, but failed to get the expected results.Thunnissen and Wit 16 also failed to extract the expected accuracy of mulched landcover information using Landsat TM and SPOT imagery.Lu et al. 5 developed a decision-tree classifier with rules obtained from analyzing the spectral characteristics of transparent PML on Landsat-5 TM images in the study area in Xinjing and then applied the classifier for extracting the PML information in 1998, 2007, and 2011, respectively.The classifier requires using Landsat TM images in two temporal phases, one in the early growing session when the transparent plastic film is just applied and another during the peak of the growing season.The results show that the decision-tree classifier is effective and temporally stable for extracting PML from Landsat-5 TM images.
PML is by far the largest type of plasticulture in term of the area it covers, and in China, 95% of plasticulture exists in mulch form. 7The time for applying plastic film as mulch is normally concentrated in a short period (e.g., within one week) during the planting stage of a growing season.Due to the long revisiting time of Landsat, it is hard to collect cloud-free images for mapping PML over a large geographic area using the method described in Ref. 5.
The moderate-resolution imaging spectroradiometer (MODIS), onboard the NASA EOS satellites, was designed to provide high-quality observations of land surfaces, atmosphere, and oceans with frequent temporal coverage. 17Because of frequent global coverage, MODIS allows obtaining cloud-free images for a large geographic area during the short planting stage, although such images usually have low spatial resolution.If an algorithm can be developed to successfully extract PML from MODIS images, annual PML mapping for a large geographic area (e.g., national, continental, or even global) will become possible.
In this paper, we present a threshold model (TM) method for mapping PML using MODIS images.The model has been tested in the Xinjiang Uyghur Autonomous Region of China, where 100% of the 1,484,100 ha of cotton fields were mulched by transparent plastic film (TPF). 18,19Methodology Supervised classifiers, such as ML, K-nearest neighbor, support vector machines, spectral angle mapper, artificial neural networks (ANN), and decision tree, have proven to be efficient classifiers for landcover mapping. 20,21,22,23However, most of these supervised classification methods require all classes that occur in a training image to be exhaustively labeled. 24This will not only increase the classification cost since the process of gathering training samples or otherwise labeling training samples is very expensive in terms of time and manpower, 25 but also be unnecessary in those studies in which only a specific class needs to be extracted.Therefore, one-class classification methods, 26 which try to detect a specific class and reject the others, have been employed in landcover mapping 27 and have proved to be effective. 28,29ne example of a one-class classification method is the TMs, which can be expressed as follows: If (the discriminative features of a pixel meet the threshold conditions) Then Assign the pixel to the special class; Else Assign the pixel to the other class.
To some extent, TMs are simple, efficient one-class classifiers, which classify the whole image into the specific class and the other class via threshold conditions.Therefore, the key for successfully using TMs for landcover mapping is to correctly select the discriminative features (or simply the discriminators) and set the correct threshold values for the conditions.In this study, we applied TMs to detect PML from MODIS images since we just focus on the PML class and ignore the other classes.
Because of its manageable data volume and high temporal resolution (covering the entire Earth surface every one to two days), MODIS images have been widely used in monitoring regional land surface processes. 30,31In this study, both the spectral and temporal features of MODIS images are explored for determining the discriminators and threshold values for mapping the PML with the threshold method.
Spectrally, plastic film has very distinct features that can be used to set the threshold conditions.Figure 1 shows the spectral curves of typical plastics, such as clear polyethylene (Clear_PE), Black_PE, and white polyvinyl chloride (White_PVC), where Clear_PE is a type of TPF, in comparison with water and vegetation by using the U.S. Geological Survey (USGS) digital spectral library splib06a.From Fig. 1, we can find that the TPF has very distinct, high, and almost constant reflectance from visible to near-infrared wavelengths, which correspond to MODIS bands 1 and 2. Therefore, high reflectance values for these two bands during the planting stage signify the possible PML with TPF.
The spectral features above are infused with the temporal features to form two temporalspectral features for PML detection with MODIS.The first is the transition of the land surface from nonplastic cover to plastic cover during the planting season.Famers put the new clear plastic mulch on the field quickly in the first several days of the planting season.This results in the significant increase of the surface reflectance in the visible/near-infrared bands within a short period of time.The MODIS time series should be able to show such a transition for the PML pixels.The second is the high reflectance in the visible/near-infrared bands during the first several weeks after planting.In the first several weeks of a growing season, the satellite sensors view mostly the plastic films because of the low canopy coverage of the ground.With the development of crops, more and more land surface is covered by crop canopy and the spectral features of PML will gradually transit from TPF spectral features (high reflectance values in both visible and near-infrared bands) to vegetation spectral features (high reflectance in the infrared band and low in the visible band).With these two temporal-spectral features in hand, we can select the discriminators and threshold values.
Based on the analysis above, we can conclude that the best time to detect PML is from the time of planting to sometime before the ground is completely covered by canopy (we call this time period the detectable period), and PML information can be extracted from band 1 and band 2 and the normalized difference vegetation index (NDVI) time series of MODIS data from the detectable period.The threshold condition is If (the number of the accumulated days >d when the discriminator's value of a pixel in MODIS time series <x during the detectable period) Then Assign the pixel to PML class; Else Assign the pixel to the other class, where x is the threshold value of the discriminator, and d is the number of the accumulated days that the pixel in the MODIS time series meets the threshold condition during the detectable period.In this study, the discriminator, d, and x are obtained by analyzing the MODIS time series for the detectable period discussed in Sec. 4 of this paper.
3 Study Areas and Data

Study Areas
In order to test the threshold model methods for mapping PML, we chose southern Xinjiang, China, (Fig. 2) as our research area for three reasons.First, >11% of PML in China is located in Xinjiang, which accounts for >40% of the total annual cotton production in China where plastic mulch is used in 100% of the cotton fields.Second, southern Xinjiang is the key cotton plantation area in Xinjiang, accounting for >70% of cotton acreage in the province. 32Third, the size of individual cotton fields or their aggregation in southern Xinjiang is big enough to be classified in MODIS images.To simplify the problem, we ignored the subpixel PML, which commonly occurs at the 250 m spatial resolution MODIS images, and defined PML as pixels or areas covered by >50% plastic mulch and non-PML as <50% covered by plastic mulch when we aggregated Landsat PML data to MODIS resolution for ground truth.
The study area, which includes Shule, Shufu, and Aketao counties, is centered at 39°1 8'8.15"N, 76°3 '33.85"E with a total of 3230 sq.km and has an arid and continental climate with a warm temperature and low annual precipitation.The average temperature is 11.3°C; the annual rainfall is ∼60 mm; 33 and the altitude ranges between 1244 and 1353 m.Cotton, winter wheat, and spring corn are the major crop types.Cotton is usually planted between early and middle April, and winter wheat is planted before the middle October, while spring corn is planted between late April and early May (Fig. 3).

MODIS time series data
MODIS has 36 bands.Among them, two have a spatial resolution of 250 m, five 500 m, and the rest have a spatial resolution of 1000 m.For PML detection, a higher spatial resolution is preferred.Therefore, the red band (620 to 670 nm, band 1) and near-infrared band (841 to 876 nm, band 2) are selected to explore temporal-spectral features for PML detection.By using the GeoBrain system developed by the Center for Spatial Information Science and System, George Mason University, USA, 35 MODIS surface reflectance daily L2G global 250 m Sin Grid V005 products (MOD09GQ MODIS data products) for the study area were subset and downloaded for covering the detectable period from the 85th day to the 150th day in 2009, 2013, and 2014.A batch tool, which we developed with ArcGIS ModelBuilder, was used to reproject MODIS data to the Universal Transverse Mercator coordinate system using nearest-neighbor resampling so that they can be coregistered with Landsat images used in this study.
Three time series of 250 m spatial resolution images were produced: red band (band 1), nearinfrared band (band 2), and the NDVI, which is defined as the following: 36  Fig. 3 Phenological calendar of major plantation in southern Xinjiang. 34here R is the value in the red band of an image pixel and NIR is the value in the near-infrared band of an image pixel.In the study, we use MODIS band 1 as the R band and band 2 as the NIR band.

Landsat images
In order to obtain the training and testing data as well as the cropland mask for this study, we manually interpreted the Landsat images for the same years as the MODIS data.Landsat images were freely downloaded from the website of Geospatial Data Cloud, China, 37 and the USGS official website. 38Table 1 lists the images from Landsat 7 and Landsat 8 used in this study.
Because of the malfunction of the Landsat 7 ETM+, image LE71490332009119SGS01 lost some scanning lines.In this study, we repaired the lost scanning lines by using the gap-fill algorithm provided by Geospatial Data Cloud website. 37

Ground Truths from Landsat Remote Sensing Images
One of the most economical methods for obtaining ground truths is to interpret the available relatively high-resolution satellite images (e.g., Landsat 4-5 TM, Landsat 7 ETM, or Landsat 8 OLI at 30 m spatial resolutions) since PML can be clearly visible in those images (as shown in Fig. 4).In this study, the ground truths for training and testing the TM methods were manually interpreted from the Landsat images and then aggregated to the MODIS resolution.
Only part of the study area with a large concentration of PML and non-PML were interpreted.
The overlay of the visually interpreted PML and non-PML from Landsat images with the MODIS time series generated the ground truths at MODIS resolution by assigning a MODIS pixel to PML if more than half of the Landsat pixels within the MODIS pixel are  PML pixels.For each year, ∼25;000 MODIS pixels in the study area were identified as the ground truths in this way.Among them, ∼10;000 pixels were PML pixels and the rest were non-PML pixels.
In addition to the agricultural land, the study area also contains other landuse/landcover types.Since the objective of this research is to test the effectiveness of MODIS time series data with TM methods for extraction of PML over agriculture land, a cropland mask was needed to mask out the non-agriculture land in the study area.A classification algorithm, the ML, was used to obtain the cropland mask because ML is one of the widely used classifiers with relatively high classification accuracy. 39In this study, training datasets for the six landcover classes in the study area, including PML, non-PML planted land, fallow land, built-up area, water body, and bare land (such as Gobi and sand), were manually interpreted from the Landsat images and used to train the ML classifier.The trained ML classifier was used to automatically classify the Landsat images for the entire study area.The classification results were further processed with majority/minority analysis and visual correction to filter out the small spatial aggregation of classes.Table 2 shows the final classification accuracy of Landsat images.The classes of PML, non-PML planted land, and fallow land were merged into the cropland class and other three classes into the non-cropland class.Then the result was aggregated to the MODIS resolution by using the majority rule to form the cropland mask at MODIS resolution.

Determination of Threshold Condition and Value
There are multiple temporal-spectral features of MODIS time series useable as the discriminator for PML mapping.In order to determine the best one as the discriminator as well as its associated d, and x, we need to analyze the performance of different features.By using PML and non-PML training sets, we calculated the means of band 1 reflectance, band 2 reflectance, and NDVI for PML and non-PML for each day.Therefore, three pairs of daily mean spectral reflectance (DMSR) curves, which are band 1, band 2, and NDVI, are obtained from the MODIS time series for each year.
Figure 5 shows the results of the comparison of the DMSR time series between PML and non-PML using original MODIS data in 2009, 2013, and 2014.From Fig. 5, we can find that in all three years, band 1 reflectance of PML is generally higher than non-PML, while both band 2 and NDVI are just the opposite.However, there were 35 days in 2009, 43 days in 2013, and 35 days in 2014 when values of band 1, band 2, and NDVI for both PML and non-PML were very close.By checking the quality control (QC) layer of MODIS data, it was found that the dates when the values for both PML and non-PML were very close were covered by cloud either fully or partly.
Although the study area was located in a dry region, there was frequent large cloud cover during the planting seasons.In order to remove the impact of cloud in this study, visual inspection of MODIS images was performed to find out the MODIS images with large cloud cover, which were not used in the calculation of the DMSR curve.Instead, a rebuilt image was constructed for the date of a discarded MODIS image by linearly interpolating the reflectance of individual pixels from reflectances of the corresponding pixels in good images of the nearest earlier and later good days indicted by the MODIS QC layer.For the period from the 90th to 125th day of a year, 57, 71, and 54% of MODIS images in 2009, 2013, and 2014, respectively, were discarded due to large cloud cover.After the interpolation, the means were recalculated to reconstruct the DMSR for both band 1 and the NDVI time series (Fig. 6).Band 2 time series was Fig. 5 Comparison of daily mean spectral reflectance (DMSR) time series between PML and non-PML using the original data.
Fig. 6 Comparison of DMSR time series between PML and non-PML using the interpolated data.
not reconstructed because there was no significant difference in band 2 between PML and non-PML.The reason for no significant difference is that the non-PML area in this study was mainly winter wheat, which has a high near-infrared reflectance in May.Because a large percentage of rebuilt images was used in the DMSR calculation in 2013, the curve was smoother than in other years.Figures 6(a), 6(b), and 6(c) show that, from the early sowing season to the end of the early growing season (the detectable period), the DMSR of band 1 for PML are obviously greater than that of non-PML (winter wheat and spring corn fields).The reason is that, from the time the farmers put the new transparent plastic mulch on the fields in the sowing season, the plastic mulch, which has a higher spectral reflectance in band 1, is visible in the early growing season.By further analyzing Figs.6(a), 6(b), and 6(c), we find that the PML DMSR curves of band 1 fluctuate around 0.2, while non-PML 0.15, and this pattern of difference was the same for all the three years.Figures 6(d), 6(e), and 6(f) illustrate that the NDVI values of non-PML are above that of the PML in three years, since cotton phenology differs from that of winter wheat, namely winter wheat is in the green-up stage, while cotton is in sowing and emerging stages.Figures 6(d), 6(e), and 6(f) also show that NDVI values for non-PML in the three years of study fluctuated between 0.2 and 0.5, but generally went up with the progress of the growing season.However, the NDVI values for PML can be divided into two stages.The first stage is from the 90th day to the 125th day when the NDVI fluctuated between 0.1 and 0.25, and the second stage is from the 126th day to the 150th day when the NDVI went up with the progress of the season.Therefore, it is clear that we can use NDVI values from the 90th to 125th day (a total of 31 days), which includes the sowing and seedling emergence periods of cotton, to discriminate PML and non-PML.
To determine the best d value, there were 31 training sessions with d from 1 to 31.For each training session, the accuracy of PML detection was calculated by using the testing sets of ground truths.Then d value with the highest detection accuracy was selected as the d value of the model.
To sum up, the PML curves in NDVI of three years are all under non-PML, despite their different ranges of DMSR time series, and there is a significant difference (∼0.1) between the PML and non-PML of DMSR from the 90th to the 125th day in all three years; therefore, NDVI DMSR from the 90th to the 125th day is the best discriminator to detect the PML among the three time series.

Detecting and Mapping PML
The statement that NDVI DMSR curves for PML fluctuate around 0.2 in the three years does not mean that NDVI threshold value of detecting PML should be set to x < 0.2.The reason is that DMSR values are the average spectral reflectance values of all pixels in the training sets, which In other words, this threshold value will lead to the misclassification of PML pixels as non-PML pixels.
In order to avoid the misclassification, we used the threshold model to identify non-PML instead of PML in the study area and then subtracted non-PML from the croplands' landcover images to get the PML area.According to above analysis and referring to Figs. 6(d), 6(e), and 6(f), among TM, x can be 0.2.Therefore, the TM for PML detection can be written as follows: If (the number of the accumulated days >d when the NDVI value of a pixel in MODIS time series <0.2 during the 95th to 125th day of a year) Then Assign the pixel to non-PML class; Else Assign the pixel to PML class We developed a java program of TM to detect PML from MODIS NDVI time series from the 90th to the 125th day.Using PML and non-PML ground truth, the classification results were evaluated by overall accuracy (OA) and κ, the kappa coefficient 40 in ENVI 5.0 software, and the results of the TM experiments are given in Table 3.
Table 3 shows that, in 2009, when threshold parameter x is 0.2 and the accumulated days parameter d is 7, the classification with the TM method has achieved the highest accuracy with OA 0.848 and κ 0.660, while in 2013 and 2014, when the threshold parameter x is 0.2 and the accumulated days parameter d is 9 and 8, respectively, classification accuracy reaches their peak with OA 0.898 and 0.864 and κ 0.796 and 0.679, respectively.It also shows that we can use the same TM to detect PML with a high accuracy in three diverse years with x ¼ 0.2 and d ¼ 8. OA values in all the three years are >0.84 and κ > 0.65.Therefore, we can conclude that TM is an effective method for mapping PML with MODIS time series data.
In order to further examine the effectiveness of the TM method, we visually compared the PML maps from MODIS with the TM method and those from Landsat with the ML classification method.To reduce the paper length, this paper only shows the MODIS and Landsat 8 OLI results for 2014 (Fig. 7).By comparing Figs.resulting from MODIS is a little bit larger than that from Landsat 8 OLI.In addition, a mismatch map between Fig. 8(b) and Fig. 8(c) was created by overlaying the two PML maps with the method described in Refs.41 and 42 and is shown in Fig. 8(d).We attribute the difference to classification errors in both the ML classifier and TM classifier.

Conclusion and Discussion
From this study, we found that the interpolated NDVI DMSR time series curve in the early growing season is a good discriminator for mapping PML in Xinjiang since NDVI values of PML are all lower than that of non-PML in 2009, 2013, and 2014, respectively, and the difference between PML and non-PML curves is more significant than that in band 1 and band 2.
We designed TM experiments and successfully applied this TM to map PML from MODIS interpolated NDVI time series (from the 90th to the 125th day).Results of TM experiments show that, when the TM parameter x is 0.2 and d is 8, the OA and κ of PML images detected by TM are >0.84 and 0.65, respectively.Therefore, we believe that this accuracy can meet the needs for mapping PML over a large geographic area.Furthermore, visual comparison shows that detected locations of PML and non-PML from MODIS are consistent with the results from ML classification of Landsat images.Therefore, TM detection and mapping of PML is feasible.The PML mapping in this study is based on a croplands' map derived from Landsat images using the ML classifier, and the accuracy of ML classifier products needs further study.

Fig. 1
Fig. 1 Spectral characteristics of typical plastics in comparison with water and vegetation (U.S. Geological Survey Digital Spectral Library splib06a).

Fig. 7
Fig. 7 Comparison of classification results in the study area in 2014 (PML: white, non-PML: gray, non-cropland: black): (a) maximum likelihood (ML) classification map with Landsat 8 OLI and (b) threshold model classification map with moderate-resolution imaging spectroradiometer (MODIS) time serials data.
7(a) and 7(b), we find the distribution of PML and non-PML in the two figures is very consistent.Meanwhile, we also provided Landsat 8 OLI false color composite, ML classification results with Landsat 8 OLI, and TM classification results with MODIS time serials data for part of the study area (see Fig. 8).Comparing Fig. 8(b) with Fig. 8(c), we found that the spatial distribution pattern of PML in both classification results is almost identical, but the acreages of PML

Fig. 8
Fig. 8 Comparison of classification results in local area in 2014: (a) Landsat 8 OLI (RGB: bands 6, 5, 4); (b) ML classification results with Landsat 8 OLI; (c) TM classification results with MODIS time serials data; (d) error of classification by comparing OLI results with MODIS results.

Table 2
Classification accuracy assessment of Landsat remote sensing images.Hang, and Di: Threshold model for detecting transparent plastic-mulched landcover. . .

Table 3
Design and results of threshold model experiments using moderate resolution imaging spectroradiometer normalized difference vegetation index time series from 90th to 125th day.OA, overall Accuracy; italics values, the highest classification accuracy for a given year; bold values, the highest classification accuracy for a given d value and classification accuracy for the best d value.Lu, Hang, and Di: Threshold model for detecting transparent plastic-mulched landcover. . .
means some PML pixels may have NDVI values >0.2.