High-resolution mapping of protected agriculture in Mexico, through remote sensing data cloud geoprocessing

ABSTRACT

In recent years, plasticulture, which is broadly defined as the use of plastics in agriculture, including greenhouse, tunnels and plastic mulch has become a fast-growing production system. While in the 1970s, the system was used on 2.8 M ha worldwide, the area had grown to 38 M ha in 1996/97 and is still increasing at a rapid rate: In Asia and Africa, the use of plastic greenhouses is growing by 15-30% per year. Therefore, plasticulture is expected to rise in many agrarian countries such as Mexico. Nonetheless, there are negative effects of plasticulture, that need to be correctly addressed in order to have responsible rural landscape management. Ergo our study attempts, through cloud geoprocessing of Sentinel-2 imagery, to rapidly provide a reliable and updated high-resolution map of protected agriculture that employs plastics (mostly tunnels and greenhouses). This map can be used as a basis for a protected agriculture monitoring system in Mexico, and eventually help to effectively plan the upcoming rural development, e.g. by providing input for territorial planning, showing the municipal rural practice trends, and monitoring the markets and its economic impacts.

KEYWORDS:

• Protected agriculture
• plasticulture
• remote sensing
• Google Earth Engine
• Sentinel 2

Introduction

Human population is increasing as well as the demand for food and agricultural products, which is expected to grow by 70–100% of the total production of 2010 by 2050 (Godfray et al., 2010). In the context of globalization, many countries, such as Mexico, are reorienting their agriculture towards exports. To achieve these goals, global annual agricultural expansion and improvements in production technology are expected (Tilman, Balzer, Hill, & Befort, 2011).

One of the most common technologies is Protected Agriculture (PA) which consists of infrastructure that modifies the natural environment, to protect the crops from external factors and enhance growth (Jensen & Malter, 1995). These include mulching, controlled environment agriculture, hydroponics, greenhouses, tunnels, etc. In this work, we focused solely on PA that employs plastics (mostly tunnels and greenhouses).

PA has shown to be a reliable option for agriculture, is productive, conservative of water and nutrients, suitable for automatization and easy to implement (Jensen & Malter, 1995). It also increases yields of crops, improves quality, and decreases the incidence of pathogens (Demchak, 2009; Lamont, 2009). These advantages have made PA rise worldwide, even in desertic areas. For instance, it is estimated that plastic-mulched land cover grows to an annual rate of 20% globally (Lu, Di, & Ye, 2014). Therefore, plasticulture can play an important role to satisfy the growing farm products’ demand by increasing food production (Lamont, 2009)

In Mexico, the PA has been incentivized. The subsidies destinated to PA projects reached near 2000 million Mexican pesos (≈ 103,715,000 USD) of the federal budget between 2001 and 2007 (IMAT, 2015). This has resulted in an exponential growth of the total area of PA in Mexico, starting with 100 ha in 1970 and reaching 25,000 ha by 2016 (Bastida Tapia, 2017). Most areas of PA are found in the central and North-western parts of Mexico (states of Sinaloa, Jalisco, Michoacán and Baja California) (SIAP, 2015).

However, PA has negative effects on various environmental aspects and at different scales. It can modify microclimates, alter distributions of pollinators, and create aesthetic pollution of the rural landscape (Jensen & Malter, 1995; Picuno, Tortora, & Capobianco, 2011). Additionally, due to their low production cost, plastics are in high demand in agriculture. They are also largely discarded which is a huge contribution to the plastic waste problem (Jambeck et al., 2015) due to their extremely low biodegradation rate (Meneses, Corrales, & Valencia, 2007; Sanders, Granberry, & Cook, 1996). In addition, PA presents an organic waste disposal problem, because all the roots, leaves and steams that do not represent the final harvest are thrown away (Tolón Becerra & Lastra Bravo, 2010). Also, it requires specialized equipment that means the greater cost to the farmers (Sanders et al., 1996). It has been shown that PA uses around three times more agro-chemicals than non-plastic-protected practices, which increases toxicity, runoff and soil erosion in the surrounding areas (Arnold, Luckenbach, & Unger, 2004; Vázquez, Pardo, Suso, & Quemada, 2006; Wan & El-Swaify, 1999). Similarly, these technologies in many cases need more water because it increases the internal temperature that entails higher evaporation rates, which leads to the use of low-quality water and excessive underground water extraction (Tolón Becerra & Lastra Bravo, 2010). Therefore, as PA has both positive and negative impacts, it is imperative to develop a monitoring system that allows measurement of changes in PA, so the effects of this practice can be properly addressed, leading towards responsible planning of rural landscapes.

Consequently, around the world, there has been a growing interest towards a reliable mapping and detection of plasticulture (Yang et al., 2017). Some of those efforts concentrate in remote sensing using satellite and hyperspectral data, because it offers the advantage that it covers much larger areas, than other techniques along with high accuracy mapping and fast processing. Remote sensing data have been used to map plasticulture in many parts of the world such as: Weifang District (Yang et al., 2017), Xiaoshan District (Wu, Deng, Wang, Ma, & Tahmassebi, 2016) and Xinjiang province (Lu et al., 2014) all in China, Southern Italy (Picuno et al., 2011), Israel (Levin, Lugassi, Ramon, Braun, & Ben-Dor, 2007) and Almeria in Spain (Agüera, Aguilar, & Aguilar, 2006; Aguilar, Nemmaoui, Novelli, Aguilar, & Lorca, 2016; Novelli, Aguilar, Nemmaoui, Aguilar, & Tarantino, 2016), using various remote sensing data as, for example, Landsat 5–8, QuickBird, Sentinel-2 imagery and hyperspectral AISA-ES photography.

In Mexico, there have been very few attempts to map and monitor plasticulture. The official data from the Agroalimentary and Fisheries Information Service (Servicio de Información Agroalimentaria y Pesquera, SIAP) is the only available source for plasticulture at a national scale (SIAP, 2015) (Figure 1). Even so, these data merely show the area of plasticulture per municipality, without mapping the infrastructures. Considering the rapid increase of plasticulture in the country, the updating of these data is far too slow to take effective action concerning landscape management, territorial growth planning, inspection or rapid response (Gorelick et al., 2017; Shimabukuro et al., 2011).

Figure 1. Protected agriculture in Mexico, number of hectares (ha) per municipality for the entire country. Own elaboration based on the data from SIAP (2015)

This study aimed at mapping, in near real time, the current plasticulture infrastructure (august 2017- august 2018) over the entire Mexican territory using Google Earth Engine (GEE), a relatively new software that enables big data analysis (Gorelick et al., 2017).

Materials and methods

We used the available highest resolution free access data, in order to implement a near real-time monitoring system that can be easily replicated and improved.

Material

We obtained Sentinel-2 multispectral imagery, acquired from august 2017 to august 2018, with less than 5% cloud percentage, through the GEE dataset catalogue. These images contain 13 spectral bands representing top of the atmosphere (TOA) reflectance (level 1-C), and an additional band for cloud mask information. Spatial resolution is 10, 20 or 60 m, depending on the spectral band, and revisiting interval is 5 days (European Space Agency, 2015).

We also used ancillary data: [1] a digital elevation model (DEM) at 30 m resolution, provided by the Shuttle Radar Topography Mission (SRTM) (Farr et al., 2007), [2] the night-time data from the Visible Infra-red Imaging Radiometer Suite (VIIRS) Day/Night Band (DNB) for the whole year (Aug. 2017- Aug. 2018) (Mills, Weiss, & Liang, 2013), [3] a digital map (shapefile) of the urban boundaries of localities with more than 15,000 inhabitants from the National Housing Commission (CONAVI, 2017), [4] a digital map of the Mexican physiographic provinces (National Institute of Statistics and Geography (INEGI) 2001), [5] the statistics on per municipality plasticulture area (SIAP 2015) and [6] a land cover/land use map, scale 1/250,000 from INEGI (2012).

Methods

We used the information on the cloud band (QA60) to mask cirrus and opaque clouds. The images were stacked, computing the median value of all of them for each band. We choose the median instead of average in order to avoid noise due to outlier values. As a following step, eight indexes, some of them developed to enhance plastic cover, were calculated a 10-m resolution (Table 1).

Table 1. Description of the computed indexes using Sentinel-2 band information. The abbreviations refer to the spectral bands: NIR near-infrared, SWIR short wave infrared, RED, GREEN and BLUE are the bands in the visible parts of the spectrum

Index	Abbreviation	Description	References
Normalized Difference Vegetation Index	(NDVI)	$\frac{NIR-RED}{NIR+RED}$	(Rouse, Haas, & Deering, 1974)
Normalize Difference Build-up Index	(NDBI)	$\frac{SWIR1-NIR}{SWIR1+NIR}$	(Zha, Gao, & Ni, 2003)
Normalized Difference Water Index	(NDWI)	$\frac{GREEN-NIR}{GREEN+NIR}$	(McFeeters, 1996)
Green Normalized Difference Vegetation Index	(GNDVI)	$\frac{NIR-GREEN}{NIR+GREEN}$	(Gitelson, Kaufman, & Merzlyak, 1982)
Plastic Greenhouses Index	(PGI)	$\left\{\begin{array}{cc}0& whenNDVI>0.73\\ 100\ast & \left(\frac{BLUE\ast \left(NIR-RED\right)}{1-mean\left(BLUE+GREEN+NIR\right)}\right)\\ 0& whenNDBI>0.005\end{array}\right.$	(Yang et al., 2017)
Plastic-mulched landcover index	(PMLI)	$\frac{SWIR1-RED}{SWIR1+RED}$	(Lu et al., 2014)
Normalized Difference SWIR2_NIR	(SW2_N)	$\frac{SWIR2-NIR}{SWIR2+NIR}$	(Aguilar et al., 2016)
Normalized Difference SWIR1+ SWIR2_NIR	(SW1+ SW2_N)	$\frac{\left(\frac{SWIR1+SWIR2}{2}\right)-NIR}{\left(\frac{SWIR1+SWIR2}{2}\right)+NIR}$	(Aguilar et al., 2016)

Ancillary data were used to mask areas where plasticulture is unlike to occur. For instance, as extensive plasticulture does not occur on steep areas, a mask was created for areas with a slope above 8% (Ortho Books, 2001). In the same manner, plasticulture is primordially a rural activity (Jensen & Malter, 1995), the Monthly average radiance composite images from the VIIRS nighttime data were summed to produce a one-year “night-light accumulation” layer, which was thresholded at 59 nanoWatts/cm²/sr to create an urban mask. This mask was complemented using the urban boundaries from the CONAVI digital map.

We choose the maximum entropy classifier (Mann, Mcdonald, Mohri, Silberman, & Walker, 2009) over all the other classifiers in GEE, because it enables performing a classification for only two classes (protected agriculture or not) and gives a map of probability ranging, which should be thresholded to obtain a binary final classification. This classifier is widely used in species distribution modelling using georeferenced localities for the presence and pseudo absence of the targeted species (Phillips, Anderson, & Schapire, 2006). The choice of the threshold applied to the probabilities enables users to control the trade-off between omission and commission errors. The classification was carried out using a stratified strategy. Each one of the 15 physiographic provinces was classified independently in order to reduce the diversity of land covers and reduce spectral confusion.

First, we made an attempt to classify a small region between the cities of Morelia and Pátzcuaro, in the State of Michoacán, one of the states where a large amount of PA is reported. We delimited training polygons through visual interpretation of the Sentinel colour composite and Google Earth imagery with the closest temporal match possible (Figure 2). A water mask (NDWI > 0) and a plastic mask with the PGI only considering classifying pixels with similar spectral behavior to plastics films of PA were created. In order to determine the optimum bands combination to accurately map plasticulture, we employed the transformed divergence index (Swain, 1973). This separability index, based on the assumption of normality, uses the mean and variance-covariance matrices of data representing two classes to assess their separability. As GEE does not allow to carry out this analysis, it was necessary to download the Sentinel composite image. Afterwards, we performed a non-supervised classification to produce training data for the non-plasticulture land cover classes able to represent the spectral variability of these classes. The global separability index was focused on the separability of PA versus any other land cover category by weighing higher separability indices involving PA.

Figure 2. Plasticulture in the region between Patzcuaro and Morelia cities, as shown in the bottom left box, both panels display the composite one-year Sentinel-2 image. (a) Bare composite image, the bright white pixels are the plastic films of PA infrastructures. (b) highlighted in red are shown the areas that were detected as PA by our classifier (these vectors are not the ones made by visual interpretation)

After bands selection, we carried out a binary classification (PA/all other non-PA covers) for the Pátzcuaro region using as training data for PA and non-PA, 1000 and 5000 randomly selected pixels, respectively, inside and outside the PA polygons previously visually identified. Finally, we executed the classifier, resulting in an accurate mapping and detection of plasticulture.

After mapping PA in the Pátzcuaro region, we carried out the classification for the entire physiographic Pátzcuaro belongs to. To do so we employ a different training data set by including 5000 random points for non PA covers. A similar procedure was then applied to the other physiographic provinces.

Then, we compared the PA national map to the Sentinel colour composite and Google Earth imagery, evaluating every municipality with a more than 25 ha difference against the SIAP data, and manually cleaned by eliminating commission and drawing omission polygons. The final result (a sample is depicted in Figure 3) was a binary map with PA and non-PA covers.

Figure 3. Protected Agriculture in the State of Sinaloa, according to this study. Location of all plasticulture infrastructures on the national territory, and number of hectares per municipality

Accuracy assessment

The accuracy assessment was carried out following Olofsson, Foody, Stehman, and Woodcock (2013); Olofsson et al. (2014). These authors indicated a strategy to assess accuracy and estimate areas of maps which present rare class(es). Therefore, we used a stratified random sampling, and the sample size (N) was determined using equation (1)(1) $N={\left(\frac{\sum W_i\ast U_i}{S\left(\text{^}O\right)}\right)}^2$(1) .

(1) $N={\left(\frac{\sum W_i\ast U_i}{S\left(\text{^}O\right)}\right)}^2$(1)

where N is the number of samples, S(Ô) the standard error of the estimated overall accuracy that we would like to achieve, W_i is the mapped proportion of area of the category i, and U_i the expected user’s accuracy of class i (proportion of the area mapped as category i which is truly i).

It is recommended between 50 and 100 samples be assigned to the rare classes and then distribute the remaining samples to the other classes proportionally to the class area. However, in the case of our PA map, the PA class is extremely rare (about 50,000 ha for a total territory of 189,991,000 ha). According to equation (1)(1) $N={\left(\frac{\sum W_i\ast U_i}{S\left(\text{^}O\right)}\right)}^2$(1) , N will depend mainly on the user accuracy of the dominant class (non-PA) which is expected to be very high because the class PA is extremely rare. As a consequence, the size of the sample will not be very large. The 50–100 samples assigned to the PA stratum will allow to assess commission error accurately. However, it is doubtful that the remaining samples can enable the detection of PA omission error. For instance, supposing an omission error of 100,000 ha (twice the size of the reported area by SIAP), the probability that a random sample unit outside the mapped as PA area falls into an omitted PA area is (100,000)/(189,991,000–50,000) = 0.0005. Consequently, the probability (p) that a sample of 500 random sampling units does not detect any omission errors is p^500 = 0.77, this probability with a sample of 1000 units, is still p^1000 = 0.59. Therefore, in the most likely case, the only commission will be evaluated, and omission error will not be detected, leading to the underestimation of the rare class area. Paradoxically, in the very unlikely case that one or more sample units detect an omission error, this would lead to an overestimation of the rare class area when correcting the area estimates using the confusion matrix.

As an attempt to deal with this problem, we constructed an additional stratum, in which most of the omission error is likely to occur. For this purpose, we used a 1 km buffer around all the PA polygons of those municipalities where this cover represents more than 2% of the whole municipality area according to the SIAP data. Selected samples were visually interpreted using both the Sentinel colour composite and Google Earth imagery. We also used formulas for estimating accuracy and areas along with the intervals of confidence of these estimators from Card (1982). Finally, we overlaid the map with a land use/cover map dated 2011 in order to assess the impact of PA on land use/cover change.

Results and discussion

The separability analysis showed that combinations using more than five bands do not significantly increase the separability due to the high correlation between bands and indices. Therefore, we selected the five bands combination which enables the highest separability which includes the Blue, NIR, SWIR1, SWIR2 and SW2_N bands. Classifications were carried out for the different physiographic regions and resulting probability images were thresholded with the value of 0.8. We chose this value in order to favour commission error (easy to detect and to correct using the accuracy assessment) and to reduce, as much as possible, omission errors.

Based on expected user accuracies of 0.7 for PA and 0.95 for non-PA categories, we estimated the sampling size in 476 verification plots. Fifty sample units were assigned to the areas mapped as PA, 50 to the buffer area around the area mapped as PA in municipalities where important PA was reported and the remaining units (376) to the rest of the territory. All the plots were selected randomly.

Then, we elaborated the confusion matrix (Table 2) and calculated the omission and commission error after carrying out the adjustments proposed by Card (1982). Taking into account omission and commission errors depicted in the confusion matrix, we estimated the area of PA estimated in 40,311 ha ± 10,231. However, as we show in the Methods section, omission errors outside the buffer area are unlikely to be detected using only 376 samples. So, if we consider that the buffer stratum did not enable us to take into account the totality of the omission errors, then the PA area estimate is likely underestimated. Another aspect is the large confidence interval around the area estimate. In order to reduce this uncertainty, the sample size for the PA and Buffer strata should be increased. However, using 100 samples for these two strata instead of 50 will likely reduce the confidence interval to 7,200 ha only.

Table 2. Raw confusion matrix

	PA	Non-PA (Buffer)	Non-PA (Outside buffer)
Total Area (ha)	47,584	151,273	189,990,788	Total plots	Commission Error	User’s Accuracy
PA	36	0	14	50	28%	0.72 ± 0.12
Non-PA (Buffer)	2	48	0	50	0.46%
Non-PA (outside buffer)	0	0	376	376
Total plots	38	48	390	476
Omission Error	5%	3%
Producer’s Accuracy	0.85 ± 0.17					Overall: 96.64%

Accuracy assessment and area estimates of such maps with both rare and large classes is not an easy task because the omission error of the rare class is difficult to evaluate. As an attempt to overcome this issue, we used a stratified strategy using a stratum where omission error is expected. In the present case, the stratum enables us to detect and evaluate omission error and therefore, after correcting the estimates through the confusion matrix. Without this stratified strategy, the area of PA would have been likely underestimated. However, it is possible that the used stratum was not able to “concentrate” enough the omission error of PA and that some amount of omission error was not taken into account. Other strategies can be envisaged to construct the strata. The first one can be based on the Maximum Entropy classification probability values. As cells with a probability above 0.8 were classified as PA, cells in which probability belongs to an interval between a certain threshold and 0.8 were not classified as PA but present a certain spectral affinity to PA and, therefore, may be prone to be misclassified. Another approach can be based on spatial modelling: maps of suitability for PA can be elaborated using explanatory variables and PA mapped areas (Camacho Olmedo, Paegelow, Mas, & Escobar, 2018; Mas, Kolb, Paegelow, Camacho Olmedo, & Houet, 2014) and areas with high suitability used as strata.

To contrast our results with the data reported by SIAP, we subtracted our municipal PA area estimations from both sources (Figure 4). We found that Michoacán, instead of Sinaloa, is the leading state at PA area; similarly, of the total 2458 municipalities nationwide, 57 had differences (either commission or omission) of 100 ha or more, meaning that both studies give very similar or equal area estimations for most of the municipalities. Figure 5 allows to compare the image and the obtained map in three of the municipalities with most PA area.

Figure 4. Frequency Histogram of the differences between estimates. The X axis is the difference between SIAP and ours municipal PA area estimations in ha while the Y axis is the frequency in the logarithmic scale. Most municipalities have a difference of less than 100 ha, and differences of more than 400 ha are very rare

Figure 5. Plasticulture Infrastructure in different windows of three municipalities with a large PA area. All windows contain two panels: A and B; panel A contains the Sentinel colour composite, where the plasticulture infrastructure is easily discernible, panel B contains the same composite image as background, and an additional layer in bright red which shows the polygons classified as PA. The windows correspond to the following municipalities: 1) Ensenada, Baja California State; 2) Zapotlan, Jalisco; and 3) Zamora, Michoacan

Nonetheless, from those 57 discrepant municipalities, according to our study 32 of them had a less estimated area (omission) and 25 had a greater estimated area (commission) against the SIAP data. Most of these omissions happened in the North of the country, for example, in the municipalities of Elota and Navolato, both in Sinaloa, with “missing” areas of 1676 ha and 1274 ha, respectively, or Mulegue and Ensenada in Baja California (589 ha and 583 ha, respectively). However, an additional examination did not allow to find unclassified PA areas. Probably, due to the arid conditions of these regions, PA is temporal, which leads to a failure in properly detecting these areas using the annual composite image. Although it is possible that PA has decreased, it is unlikely due to increasing global and national trends. Another omission example was Villa Guerrero in the State of Mexico (central region), this municipality is the major producer of flowers in Mexico and presents the singularity of having the PA infrastructures within the urban area. In this case, the missing 762 ha was considered as urban area.

Commissions, on the other hand, occurred mostly in the central region of the country, being Michoacán and Guanajuato the states with more area over-estimated by the present study in comparison with the SIAP data, with 3819 and 1416 additional ha, respectively. No particularly important commission errors were detected in these areas and the observed area differences can be explained by the rapidly increasing PA area in central Mexico, and thus the difference between studies can be attributed to the updating of the present study (2018 imagery) in comparison with SIAP (2015).

In Mexico, PA is mostly concentrated on several areas which need to be monitored to properly manage plasticulture down effects. The overlaying procedure between our maps of PA and a land use/cover map dated 2011, showed that 93% of the PA was established on croplands, mainly irrigated areas. Therefore, almost all PA establishes on areas already with an agricultural use and occur as a shift of agricultural activities. A deep analysis of PA distribution would allow to develop a projection under prospective scenarios.

Concluding remarks

Knowing the scope that can be achieved with GEE and the goals behind its development (Kumar & Mutanga, 2018), we hope that this work can set an antecedent towards a national accurate PA monitoring system and, ultimately, contribute to understanding rural changes and the promotion of sustainable rural management.

This study is, as far as we know, the first attempt to map plasticulture infrastructures and PA at a national level using cloud geoprocessing. The maps we obtained present an acceptable accuracy and show a consistency, with the approximations of the tendencies for Mexico, and the SIAP data. Our approach can be improved in the future, for instance by using Bottom Of the Atmosphere (BOA) imagery, employing more accurate training data for the classifier, calculating the separability for each physiographic province, experimenting with other classifiers, developing a more suited accuracy assessment, including eventually a phase of ground-truthing.

Geolocation information

The study area involves the entire territory of Mexico.

Acknowledgments

This research was possible thanks to a residency, within the framework of an international cooperation agreement recently settled, between the Environmental Geography Research Centre (CIGA) at the National Autonomous University of Mexico (UNAM), and the Rural and Environmental Studies Faculty (FEAR) at the Pontifical Xavierian University in Bogota, Colombia (PUJ-Bogota) that aims at promoting scientific collaboration between both countries. This study was supported by the project PAPIME (PE301919) “Herramientas para la enseñanza de la Geomática con programas de código abierto” (DGAPA - UNAM). We also thank the careful revision of our manuscript. The final corrections were done during a sabbatical stay of the second author with the support of PASPA-UNAM.

Data availability statement

The corresponding shapefile and metadata produced as a result of this study can be found at http://lae.ciga.unam.mx/recursos/Agricultura_Protegida_2018.rar

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the PAPIME-DGAPA-UNAM [PE301919].