The Temporal-Stability-Based Irrigation MAPping (TSIMAP) Method: A Virtuous Trade-Off between Accuracy, Flexibility, and Facility for End-Users

Abstract

Remote sensing technology is an essential tool for tracking human-induced alterations on the water cycle, among which irrigation prevails. The possibility of obtaining detailed and accurate information on the actual irrigation extent through remote-sensing-based approaches is of paramount importance for water resources management. In this study, an update of the TSIMAP (Temporal-Stability-derived Irrigation MAPping) method, originally developed with satellite soil moisture as an input, is proposed. To demonstrate that the flexibility of the approach does not affect its main strength point (i.e., good accuracy in the face of high simplicity for users), a dual analysis relying on 1 km NDVI (Normalised Difference Vegetation Index) instead of soil moisture is carried out over the Ebro basin (Spain); data delivered by the Copernicus Global Land Service (CGLS) are used. First, results of this work are compared with outcomes from the method’s original implementation obtained over a focus area (denominated “Ebro_CATAR”) through satellite soil moisture. In the proposed configuration relying on NDVI, an overall accuracy (OA) up to 93% is found. Results highlight an increase in OA ranging from +2% to +6% depending on the validation strategy with respect to the TSIMAP implementation relying on soil moisture. Then, a basin-scale application is performed, providing performances still satisfactory (OA = 75%) notwithstanding a higher degree of heterogeneity.

1. Introduction

Irrigation is the prevalent source of anthropogenic water consumption. Its effects have repercussions on multiple dynamics affecting human beings and Earth’s surface [1], such as food production [2,3,4], spread of diseases [5,6], freshwater availability [7,8,9] and climate [10,11,12]. Notwithstanding this, irrigation is still often not considered a hydrological variable in its own right, and data on irrigation dynamics are often missing [13,14,15] or not scientifically organised [16].

The detailed knowledge of actual irrigation extents is crucial for many hydrological and water resource management applications [17,18,19]. Starting from the early 2000s, several irrigation mapping techniques relying on satellite information have been developed. As pointed out by the comprehensive review by Massari et al. [20], various types of satellite retrievals from optical [21,22,23,24], to cite a few, and microwave, see, e.g., [25,26,27,28], sensors can be successfully exploited for mapping irrigated areas at different spatial scales. Methods relying on multi-source input data also exist [29,30,31,32,33,34]. Most of the approaches available in the scientific literature rely on supervised machine learning models requiring training; hence, ancillary ground-truth information is required for their implementation. By contrast, Dari et al. [27] proposed a novel approach for speedy irrigation mapping relying on a combination of statistical features derived from the temporal stability theory [35] and the unsupervised k-means clustering algorithm [36], henceforth defined as TSIMAP (Temporal-Stability-derived Irrigation MAPping) approach. The method does not require any training procedure and, as recently pointed out in a comparative analysis by Bazzi et al. [37], ensures satisfactory performances in face of its simplicity. So far, the method has been developed and implemented with satellite soil moisture at different spatial resolutions as input [27,28]. Nevertheless, it is supposed to work with any kind of input information that is a proxy for irrigation occurrence.

The aim of this work is to actually demonstrate the flexibility of the TSIMAP method, whose applicability is not limited by a single type of input data. In order to achieve this, a twofold analysis carried out over the Ebro basin (Spain) by exploiting 1 km NDVI (Normalised Difference Vegetation Index) from the Copernicus Global Land Service (CGLS) is proposed. First, a comparative analysis with outcomes from the method’s implementation relying on satellite soil moisture by Dari et al. [27] is carried out. The comparison is focused over the same period (2016–2017) and the same portion of the Ebro basin considered in a previous study by the authors. As a second step, a basin-scale application considering the irrigation seasons in the time span of 2015–2020is performed.

2. Materials and Methods

2.1. Study Area

The study is focused over the Ebro basin (Spain), a site that has attracted, in recent years, the interest of many researchers working on the topic of irrigation monitoring through satellite data [26,34,38,39,40,41]. This is because the Ebro catchment is an area rich in irrigation data, a rather rare feature worldwide.

With a total drainage area of about 86,000 km², the Ebro basin is the largest Mediterranean catchment of Spain; it is prevalently characterised by a cold semi-arid climate with an oceanic climate area in its Northern portion [42], where several reservoirs store water from the Pyrenees to be delivered to the crops located downstream in the central valley through a complex network of irrigation canals [39]. The analysis presented in this study is carried out on the basin’s flat agricultural portions only, which are shown in Figure 1. Such a domain includes both irrigated and non-irrigated areas; it is outlined through the 100 m resolution Corine Land Cover data set referring to the year 2018 (CLC 2018), distributed by the Copernicus Land Monitoring Service (https://doi.org/10.2909/960998c1-1870-4e82-8051-6485205ebbac (accessed on 13 November 2019). The land cover types that are considered are listed in the figure’s legend. It is important to specify that the information from the CLC 2018 is only used to delineate the spatial domain of agricultural areas (rainfed and irrigated) to be considered; as a benchmark for the actual extent of irrigated areas, we use more site-specific and accurate data sets as described in Section 2.3. Moreover, a double-scale analysis is carried out. First, the area enclosed by the red box of Figure 1 is considered in order to compare results with outcomes of Dari et al. [27]; such portion, whose extension is 10,000 km², encloses a wide irrigated area shared between the Catalonia and Aragon regions and is henceforth denominated “Ebro_CATAR”. Then, a basin-scale application is performed (see Section 2.3 for further details).

2.2. The TSIMAP Method: Theory and Input Data

The TSIMAP method was first proposed by Dari et al. [27]. It relies on the temporal stability concept originally formulated by Vachaud et al. [35] for optimising soil moisture measurement networks [43,44,45,46]. Morbidelli et al. [47] also explored the possibility of adapting the concept to raingauge network design. The theory relies on statistical features able to describe the spatiotemporal dynamics of the input data set. Hence, if applied to an irrigation proxy variable, it is possible to use the temporal stability theory to detect where irrigation actually occurs. With this intuition, the approach was originally proposed by exploiting satellite soil moisture as an input [27]. In this study, the method is implemented using NDVI data from the CGLS (https://land.copernicus.eu/global/products/ndvi (accessed on 29 November 2023)), retrieved through the VEGETATION sensor aboard the SPOT (Satellite Pour l’Observation de la Terre) platform, subsequently substituted, from 2013 onwards, by PROBA-V (PRoject for On-Board Autonomy-VEGETATION) [48]. Specifically, product version 2.2 characterised by a spatial resolution of 1 km and sampled at a temporal frequency of 10 days is used [49]. Data from 2015 to 2020 are considered.

Considering the NDVI value observed at the pixel $x$ and time step $t$, $NDV{I}_{xt}$, and the time-varying spatial average, $\overline{{NDVI}_t}$, the relative differences, $R{D}_{xt}$, are defined as

$${RD}_{xt}=(NDV{I}_{xt}-\overline{{NDVI}_t})/\overline{{NDVI}_t}$$

(1)

While the relative differences reveal features of spatial patterns of the input data set, the temporal anomalies, $T{A}_{xt}$, provide information on temporal dynamics. $T{A}_{xt}$ values indicate, for each pixel, the magnitude of the discrepancy with its temporal mean, $\overline{{NDVI}_x}$; they can be calculated as

$${TA}_{xt}=(NDV{I}_{xt}-\overline{{NDVI}_x})/\overline{{NDVI}_x}$$

(2)

If averaged in time during the highest intensity irrigation season, which for the study area coincides with the period spanning from May to September [27], the relative differences and the temporal anomalies can be used as diagnostic indices for detecting where irrigation occurs. In fact, during the irrigation seasons, irrigated crops are expected to show well-developed vegetation conditions, which can be translated into higher values of relative differences and of temporal anomalies with respect to rainfed areas. On the basis of such an assumption, the temporal-stability-derived indices can be used as input features for the k-means clustering method [36]. It allows to group $n$ data points into $k$ clusters on the basis of their similarity, expressed by the minimum Euclidian distance with each cluster’s centroid. The k-means is an unsupervised clustering algorithm; hence, it does not require any training, but the user is asked to interpret the clusters’ meaning. Further details can be found in Dari et al. [27,28]. In this work, the mean relative differences (RD) and the mean temporal anomalies (TA) during the highest intensity irrigation seasons are used as input data points for k-means clustering. The final output consists in binary maps of rainfed and irrigated areas at 1 km spatial resolution.

2.3. Work Logic and Ground Truth Data Sets

A twofold analysis is carried out. First, a comparison between results obtained by forcing the TSIMAP approach with NDVI data and outcomes of the method’s first implementation relying on soil moisture [27] is performed. The analysis iscarried out over the same basin’s portion considered in Dari et al. [27], the Ebro_CATAR domain (10,000 km²), enclosing the largest irrigated portion of the entire basin. In this first analysis, the same validation strategy proposed in the previous work by the authors is implemented to ensure that the comparison is as fair as possible. The crop maps produced over the domain of the Geographic Information System for Agricultural Parcels (SIGPAC) and distributed by the Department d’Agricultura, Ramaderia, Pesca i Alimentació are used. This data set, henceforth referred to as SIGPAC data, is available at the plot scale and for each year at http://agricultura.gencat.cat/ca/serveis/cartografia-sig/aplicatius-tematics-geoinformacio/sigpac/mapa-cultius/ (accessed on 9 October 2019). It is noteworthy that the data set is available for the Catalonia region only. The pre-processing chain to aggregate SIGPAC data from the original resolution to 1 km is fully described in Dari et al. [27]; the only difference adopted in this work is that natural or forested areas are not considered, so, this kind of pixels are masked out. Two types of information contained in the SIGPAC data set are exploited, resulting in an equal number of validation assessments. First, the irrigated/non-irrigated attribute is exploited (first validation). Even though the irrigated areas of the Ebro basin are quite static in time, as thoroughly explained in Dari et al. [27], the irrigation attribute in the SIGPAC data set means that the irrigation network can reach a specific field, but it does not necessarily imply that irrigation actually occurs there. Hence, the information referring to summer crops only, which are very likely irrigated during the irrigation season, is considered in the second validation, in which all the other crops are considered as rainfed. Similarly to what was performed in Dari et al. [27], this first analysis covers the 2-year period of 2016–2017. It is important to mention that in the previous work by the authors, experiments considering three clusters (rainfed, natural/forested and irrigated areas) were carried out. To ensure a consistent comparison with the authors’ previous study, the best performing classification in reproducing irrigated areas resulting in that work is re-implemented here by considering only two clusters, i.e., rainfed and irrigated areas. In this way, results obtained in Dari et al. [27] and in this work are fully comparable.

As a second step, a basin-scale application is proposed. Hence, the method is applied over the whole domain of agricultural areas (see Figure 1). In this case, as a ground truth information, the map of the basin’s irrigated areas according to the Hydrological Plan of the Ebro basin is used. The map, produced by the Ebro Hydrological Confederation (CHEBRO, https://www.chebro.es/ (accessed on 15 April 2020)), is dated 2004; hence, it could miss some irrigation portions recently developed. Nevertheless, the main irrigation districts were already developed at that time. Such information, henceforth defined as the CHEBRO map, is used to validate a static map obtained by applying the TSIMAP method to all the irrigation seasons of the considered period, i.e., from 2015 to 2020. This choice is led by the reasonable assumption of irrigation extents quite static in time at the basin scale and also by the nature of the available ground truth information. For validating the basin-scale irrigation map, all the pixels of the CHEBRO map are considered as irrigated (since the map only provides the edges of irrigated areas), while the other pixels of the agricultural areas within the Ebro basin are considered as rainfed. Figure 2 provides a scheme summarising the main features of the analyses carried out.

3. Results

3.1. Comparison with Original Implementation Relying on Soil Moisture

In this section, results obtained over the smaller domain considered in Dari et al. [27] and validated with detailed ground truth information from the SIGPAC data set are presented. Outcomes referring to the irrigation seasons of 2016 and 2017 are shown in Figure 3 and Figure 4, respectively. Both figures are organised in six panels: Figure 3a and Figure 4a show the location of the Ebro_CATAR area with respect to the basin extension, Figure 3b, Figure 4b and Figure 3c, Figure 4c show the mean RD and the mean TA, respectively, during the highest-intensity irrigation season (May–September), Figure 3d, Figure 4d and Figure 3e, Figure 4e provide the SIGPAC ground truth maps according to the two validation strategies, and Figure 3f and Figure 4f show the rainfed/irrigated map resulting from the TSIMAP method, respectively. In Figures from Figure 3d and Figure 4d to Figure 3f and Figure 4f, light brown and blue green are used to indicate rainfed and irrigated areas, respectively. We note that for Figures from Figure 3b and Figure 4b to Figure 3f and Figure 4f, data are shown over the portion falling into the Catalonia region only, as this is the domain where benchmark information from the SIGPAC data set is available, respectively. Even with different magnitudes, for both irrigation seasons of 2016 and 2017, the mean RD and TA of NDVI prove themselves to be promising diagnostic indices for detecting irrigated areas. In fact, the maps of Figure 3b, Figure 3c and Figure 4b and Figure 4c show consistent spatial patterns with higher (mostly positive) values over areas where irrigation is expected to occur and lower (mostly negative) values over the surrounding dryland area extending to the east, respectively.

For both irrigation seasons of 2016 and 2017, the mean RD and TA are used as input data points to perform the k-means classification and produce the 1 km maps of rainfed and irrigated areas shown in Figure 3f and Figure 4f, respectively. Figure 5a provides the confusion matrices synthesising the accuracy of the resulting maps according to the two validation approaches, i.e., comparing maps in Figure 3f, Figure 4f with those shown in Figure 3d, Figure 4d and Figure 3e, Figure 4e. Overall, a satisfying accuracy in properly reproducing rainfed and irrigated portions is obtained in both irrigation seasons of 2016 and 2017. According to Validation 1, 88% and 87% of irrigated areas are properly reproduced by the TSIMAP method in 2016 and 2017, respectively. Almost all the rainfed areas are properly classified, since in both irrigation seasons of 2016 and 2017 an accuracy of 97% is reached. According to Validation 2, the same results are obtained for both irrigation seasons considered. The accuracy in properly classifying irrigated areas rises to 98%, but this happens at the expense of the accuracy in reproducing rainfed areas, which drops to 84%.

Dari et al. [27] explored different combinations of input parameters derived from satellite soil moisture to map irrigation. The authors found that the standard deviation of RD and mean TA calculated for the 1 km DISPATCH (Disaggregation based on Physical and Theoretical scale Change) downscaled SMAP (Soil Moisture Active Passive) product [50] were the best diagnostic indices for reproducing irrigated areas, reporting a maximum accuracy of 78% reached according to Validation 2 considering the irrigation season of 2017. However, as explained in Section 2.3, such a result refers to the capability in reproducing three different clusters, and not in discerning between rainfed and irrigated areas only as in this case. For this reason, the experiment offering the best accuracy according to Dari et al. [27] is repeated here by considering two clusters only and results are summarised by the confusion matrices of Figure 5b. The differences with the implementation relying on NDVI proposed here thin out, especially in the capability of the method in reproducing actual irrigated areas; accuracies differ by 1–2%, still in favour of the new implementation. Higher differences are found with respect to the capability of recognising rainfed areas. The classification based on satellite soil moisture provides accuracies lower than those relying on NDVI by 10% and 8% (and by 3%) for both the validation experiments in 2016 (2017).

3.2. Basin-Scale Application

The TSIMAP method is then applied at the basin scale to produce a static map of rainfed and irrigated agricultural areas across the years of 2015–2020. The k-means classification is carried out by using as input data points the mean RD and TA during the irrigation seasons (May–September) of the period of 2015–2020. The resulting maps are shown in Figure 6; consistent spatial patterns of higher average RD and TA values are found over areas known to be irrigated, such as agricultural lands within the Ebro_CATAR domain, the narrow area in the basin central valley running along the main reach of the Ebro and the river delta [39].

The accuracy of the map of rainfed/irrigated areas derived from the input indices shown in Figure 6 is assessed through comparison with the CHEBRO map. Results are summarised in Figure 7. Figure 7a shows the map produced by the TSIMAP method, Figure 7b provides the ground truth map and Figure 7c contains the confusion matrix. As can be depicted by visual inspection of maps in Figure 7a and Figure 7b and by accuracies reported in Figure 7c, irrigated areas are satisfactorily reproduced at the basin scale also, which is characterised by higher heterogeneity. The percentage of irrigated areas properly classified is equal to 81% of the total. Higher uncertainty exists in properly classifying rainfed areas, for which the accuracy drops to 72%.

4. Discussion

The aim of this work is to showcase the flexibility of the TSIMAP method originally proposed in a previous work by the authors with satellite soil moisture data as an input. However, the method is expected to provide satisfactory results with any input observation that can be an irrigation proxy and for which the temporal stability concept applies, without losing its characteristic of ease of use that distinguishes it from other methods available in the literature [37].

For this reason, the approach is implemented with 1 km NDVI data from the CGLS, as several studies linking irrigation dynamics with NDVI are available in the literature [23,51,52,53]. The comparison against outcomes from the method’s first implementation reveal a slight increase in the capability of properly classifying irrigated areas and more substantial improvements in detecting rainfed areas when forcing the approach with NDVI observations. In terms of overall accuracy (OA), a value of 93% and 92% is reached according to Validation 1 for the irrigation seasons of 2016 and 2017, respectively. For Validation 2, the OA is 88% for both years. Over the same area, Elwan et al. [34] reported an accuracy on the order of 87% through a supporting vector machine classification fed with statistical features from Sentinel-1 and Sentinel-2 as an input. Gao et al. [26] reached an accuracy of 81% using Sentinel-1-derived features over the same domain, but three classes (irrigated crops, irrigated trees, and non-irrigated fields) were considered in that work. In both works by Gao et al. [26] and Elwan et al. [34], methods requiring training were adopted, while the TSIMAP approach did not require this step. It is important to point out that results from Gao et al. [26] and Elwan et al. [34] refer to a finer spatial scale (at the field level) with respect to the 1 km scale adopted in this work.

In the basin-scale application, an OA equal to 75% is reached. Even though a decreasing accuracy can be encountered when enlarging the considered domain and thus facing higher heterogeneity, it is important to point out that results of the two analyses carried out are not directly comparable since different ground truth data sets are used. A closer look at NDVI dynamics can help in assessing results. Figure 8 shows the rainfed and irrigated areas properly classified (true rainfed and true irrigated, respectively), the irrigated areas erroneously classified as rainfed (false rainfed) and the rainfed areas wrongly classified as irrigated (false irrigated). The 14-day aggregated time series of NDVI averaged over each one of the domains above are also shown in Figure 8b, where the main watering seasons (May–September) are highlighted in light grey. False rainfed areas are scattered and often located adjacent to the edges of irrigated areas; hence, it is possible that the 1 km sampling adopted here is not fine enough to properly detect complex patterns. However, the NDVI signal over those areas is similar to the one detected over true rainfed portions (see light brown and pink lines in Figure 8b). The green blue line referring to true irrigated areas clearly shows a second NDVI peak during the period May–September attributable to irrigation occurrence. False irrigated areas (mostly located in the northwestern side) represent the highest source of uncertainties. For those areas, the NDVI time series (black line in Figure 8b) do not show a clear peak during the irrigation season, but values are comparable to those detected over true irrigated areas. This can be an index of areas where NDVI is systematically higher than the basin average (see also Figure 6), thus leading to wrong classification. Nevertheless, uncertainties could come from the ground truth data set as well. As explained in Section 2.3, the CHEBRO map dates back to 2004; hence, it is possible that recently developed irrigated areas are detected by satellite NDVI but are still classified as rainfed in the benchmark data set.

5. Conclusions

The analyses carried out in this work demonstrate the flexibility of the TSIMAP method, able to provide optimal performances in face of its simplicity for end-users [37]. The approach’s applicability is not limited to the use of satellite soil moisture as an input, as the method can be fed with any irrigation-proxy data (e.g., vegetation indices, evapotranspiration estimates, land surface temperature records). To showcase this, an implementation relying on 1 km NDVI is proposed. Results outperform the method’s first implementation by Dari et al. [27], with increases of OA equal to +6% in the irrigation season of 2016 according to both validation procedures implemented and to +2% and +3% according to Validations 1 and 2, respectively, considering 2017. A basin-scale application is performed as well, producing a static map across irrigation seasons within the period of 2015–2020 providing satisfactory results despite a higher heterogeneity. The method’s flexibility allows potential operational applications at different spatial resolutions [28]. For instance, the TSIMAP approach represents a useful tool to refine irrigation water use estimates from satellite data, as hinted by Kragh et al. [54] and Dari et al. [39]. Further assessments of the performances of the proposed approach with different types of input (e.g., merging soil moisture and vegetation indices) data with respect to varying climatic contexts are the foreseen next steps of this research.