Satellite-Based Mapping of Cultivated Area in Gash Delta Spate Irrigation System, Sudan

Abstract

In this study, a simple methodology for mapping the seasonal cultivated area of the Gash Delta Spate Irrigation System based on satellite images was developed. The methodology combined information from multiple bands to characterize the land surface in terms of spectral indices (e.g., Normalized Difference Vegetation Index (NDVI), and surface temperature (Ts)). Visual interpretations of a conveniently selected image were undertaken to identify and select sample points of interest. The NDVI and Ts values (computed from multi-date images that represented the crop growing period) of the sample points were used to developed typical NDVI and Ts plots. By analyzing these plots and the cropping calendar, an NDVI and Ts threshold-based algorithm was developed to extract the cultivated area of a given season. Analysis of the developed algorithm showed that it was simple, easily modifiable, and had interpretable rules and threshold values. Comparing the extracted cultivated area with the field report area showed a promising application of the methodology to map and estimate the cultivated area from only remote sensing data.

1. Introduction

As competition for water increases and commitments for sustainable ecosystems grow, an increasing awareness has been observed in recent times for making the best use of water, a scarce and valuable resource for all economic activities [1,2]. Irrigated agriculture, which produces approximately 40% of the world’s production, consumes about 70% of the world’s total available freshwater [3]. Irrigation practices can play a significant role in agricultural production to meet the projected food demand in several parts of the world by maintaining or increasing crop yields under changing climatic conditions. As a result, accurate information on irrigated agricultural water use, and its spatial extent and variation, is essential for water resource and crop productivity assessments [4]. One type of information that requires accurate estimation is the acreage of an irrigation system. It is an important input parameter for many irrigation performance indicators such as water use, water productivity, water rights, impact assessment, and performance diagnosis, among others [5]. Despite their importance in estimating irrigation performance indicators, irrigation system data including irrigated/cultivated area are seldom measured in a regular and reliable manner to cover the entire irrigation system. According to [6], the lack of proper and effective management of irrigation systems is related to the conventional method of data collection through field observation, which is difficult, time-consuming, and inadequate in temporal and spatial coverage.

The advancement of satellite technology in terms of spatial, temporal, spectral, and radiometric resolution offers an opportunity to capture the agricultural information of a large area at frequent intervals [7,8,9]. Satellite observations provide reliable, economical, and synoptic data of the Earth’s surface. It has been able to provide information with varying degrees of success and accuracy on: crop types [10,11], yield [12,13], and evapotranspiration [14,15]. Remote sensing techniques have also become a viable option and are extensively used for mapping irrigated/cultivated areas [16,17,18,19,20,21]. Ozdogan et al. [22] provided a thorough review of the sensors (e.g., Moderate Resolution Imaging Spectroradiometer (MODIS), Satellite Pour l’Observation de la Terre (SPOT), Land remote-sensing Satellite (LANDSAT), etc.) and methods used for mapping irrigated/cultivated areas at local, regional, continental, and global scales.

At a local scale, the most common mapping methods are visual interpretation and the digital classification of single or multi dated images [22]. The most commonly used methods of the latter method include multi-stage classification, unsupervised clustering, density slicing, and decision tree classification. Most of the digital image classification approaches depend on statistical similarities to a predefined set of conditions based on observations or existing maps. The requirement of enormous amounts of training and ancillary data makes these methods of classification impossible from remote sensing data alone [23]. Moreover, these methods depend on the magnitude and temporal evolution of spectral vegetation indices derived from red and near infrared reflectance [8,19]. However, some studies [24,25,26,27] have shown the potential of the combined analysis of spectral vegetation indices with surface temperature observations for studying surface energy exchange processes. A simple and completely remote sensing based land cover classification scheme was proposed by [23]. They employed the seasonal dynamics of the Normalized Difference Vegetation Index (NDVI) and surface temperature to characterize different land cover types such as barren, shrub, forests, water, crops, and others. We adopted their classification logic, with necessary modifications, to develop a methodology to extract the seasonal cultivated area of the Gash Delta Spate Irrigation System from remotely sensed data.

The methodology deals with visual interpretation and spectral and thermal analysis of carefully selected Landsat-7 Enhanced Thematic Mapper Plus (ETM+) images that fall within the given season. This analysis, plus information on the cropping calendar, was focused on developing a simple and easily modifiable and interpretable threshold-based algorithm to extract the cultivated area. That is, although there are voluminous remote sensing-based indices and classification approaches, in this study, however, we sought to develop a methodology that was simple and intuitive for the target users. In doing so, the main objective of the study was to investigate the capability of visual interpretations of a conveniently selected satellite image in identifying and selecting sample points and the importance of combining the spectral and thermal information in uniquely characterizing the land cover of the study area.

2. Study Area and Data

2.1. Study Area

The Gash Delta Spate Irrigation System (GDSIS), with a net command area of 100,000 ha, is located in Kassala State, Sudan, close to the border with Eritrea (Figure 1). It is an oasis in a surrounding desert with relative humidity ranging from 20–50%, annual rainfall from 180–280 mm (from Hadaliya to Kassala), and an average temperature from 26 °C in winter to 42 °C in summer [28].

The total command area of GDSIS is divided into six irrigation blocks: Kassala, Mekali, Degain, Tendelai, Metateib, and Hadaliya. Each block is further divided into a number of irrigation units known locally as “Misga”. Annually, about one-third of the total Misgas are planned for irrigation and subsequent cultivation, while the remaining two-thirds are fallow for that particular season. Every season, land preparation, which involves clearing of mesquite (Accacia proposes, a thorny and evergreen tree), land leveling, and routine maintenance activities, is carried out before the onset of the flood. Then, floodwater from the west bank of the Gash River is diverted through seven diversion structures for an effective period of 60 to 70 days and conveyed by a system of canal networks to the prepared Misgas, in two turns. That is, the total Misgas to be irrigated in any one year are divided into two parts. These are known as the first and second turns, and irrigation does not begin on the second turn until that on the first turn has been completed [29]. In planning the Misga turns for the first and second turns, it has been assumed that the first flood (mid-July to mid-August) is stronger than the second flood (mid-August to mid-September), thus, 60% of the area to be irrigated that year is planned for the first turn. The portion of the prepared Misga that received water within the irrigation period (July–September) is referred to as the irrigated area. Usually, about one week after irrigation, when tractors can enter the field, seed sowing is done on the portion of the irrigated area that has been sufficiently watered and adequately prepared. The area that is actually sown with seed is referred to as the cultivated area [30]. The dominant crop in GDSIS is a sorghum variety called Aklamoy and is cultivated between September and February.

2.2. Data

The primary data used in this study were Landsat-7 ETM+ images. A total of 14 (seven for each season (2009–2010 and 2004–2005)) Landsat-7 ETM+ images were used for the algorithm development. In each season, one of the images was within the irrigation period, while the remaining six were within the crop growing period. These gap-present (scan line corrector (SLC)-off) images were downloaded from the United States Geological Survey (USGS) Global Visualization Viewer (http://glovis.usgs.gov). The selection of this data product was justified by the following reasons. First, Landsat-7 ETM+ has a relatively high spatial and temporal resolution. Second, the location of the case study was within the center of the image where the data quality was very similar to the SLC-on Landsat-7 ETM+ data. In addition, Google Earth was utilized to delineate block boundaries and to supplement the identification of sample points. Finally, field reported data on cultivated areas compiled by [30] were used to compare the results of the developed method.

3. Methods

The focus of this study was to develop a simple satellite remote-sensing based algorithm for mapping and estimating the cultivated area of GDSIS. The methodology consisted of four main steps: (1) the preparation and enhancement of satellite images; (2) the visual interpretation of the original and enhanced images to identify and select sample points; (3) building the threshold-based classification algorithm; and (4) testing the applicability of the methodology and the algorithm. Figure 2 depicts the flow chart of the process and the evolution of the satellite image through each process in this study.

3.1. Image Preparation and Enhancement

In this paper, image preparation dealt with (1) acquiring and organizing the satellite images, (2) selecting the bands and pre-processing, and (3) the visual and spectral enhancement of images. Owing to the availability of field data for verification of the results of the proposed methodology, the search for cloud free Landsat 7-ETM+ images was limited from 2004 to 2011. The 2009–2010 season was found to have many and almost evenly distributed images across the cropping season. Eventually, this season was selected to build the algorithm. A satellite image acquired on 28 July 2009 fell within the irrigation period, while all the other images were within the crop growing period. Bands 1 to 6 from the single-date image (28 July 2009) and 3, 4, and 6 from the other images were selected and processed for the production of visually and spectrally enhanced images. Landsat-7 ETM+ sensors capture reflected solar energy, convert these data to radiance, then rescale these data into an 8-bit digital number (DN) with a range between 0–255. The optical and thermal DN values were converted to the reflectance and land surface temperature (Ts), respectively, by applying the procedure explained in [31].

Spectral enhancement was performed by computing the relevant vegetation and water spectral indices. In the GDSIS, flood water from the Gash River is diverted to each prepared Misga for about three weeks within the irrigation period (mid-July to mid-September). This practice leaves the irrigated area covered with water or wet enough to be clearly detected on a satellite image acquired within this period. Hence, the Modified Normalized Difference Water Index (MNDWI) proposed by [32] was adopted to identify the typical irrigated area. It is the normalized ratio of green (G) and middle infrared (MIR) bands and is expressed as:

$$\mathrm{MNDWI}=\frac{\mathrm{G}-\left(\mathrm{MIR}\right)}{\mathrm{G}+\left(\mathrm{MIR}\right)}$$

(1)

This index is designed to (1) maximize the reflectance of water by using green wavelengths; (2) minimize the low reflectance of MIR by water features; and (3) take advantage of the high reflectance of MIR by vegetation and soil features. As a result, water features have positive values and are enhanced, while vegetation and soil features usually have zero or negative values and are, therefore, suppressed [32].

The vegetation spectral index used in this study was the Normalized Difference Vegetation Index (NDVI). It is the normalized ratio of near infrared (NIR) and red (R) bands, which are the most sensitive to the vegetation, and is expressed as:

$$\mathrm{NDVI}=\frac{\left(\mathrm{NIR}\right)-\mathrm{R}}{\left(\mathrm{NIR}\right)+\mathrm{R}}$$

(2)

3.2. Identification and Selection of Sample Points

Most studies that involve the classification of satellite images start with the collection of sample points for training the classification algorithm and validating the result. Generally, these sample points have to be collected by field survey at the same time as the satellite image acquisition [17]. This often constitutes the largest cost of the mapping project. However, when the mapped features are evidently identified from the background environment, photo-interpretation can substitute the costly field survey [33]. Hence, in this study, visual interpretation of the image was employed to identify irrigated, dry-bare, and forest areas. Due to the specific nature of the GDSIS, images during the irrigation period (or before sowing) were suitable for identifying the above features. Accordingly, the Landsat-7 ETM+ image on 28 July 2009 was selected for this purpose. The image was then processed to produce the standard false color, pseudo natural color, NDVI, and MNDWI maps. Finally, the above maps were visually interpreted to identify and select sample points of irrigated, bare-dry, and forest areas.

3.3. Building the Classification Algorithm

The goal of this study was to develop a simple land cover classification algorithm to extract the seasonal cultivated area of the GDSIS. First, NDVI and surface temperature (Ts) plots of all the sample points identified and selected by visual interpretation were produced to study the temporal dynamics of crop condition across a season. These plots, plus knowledge of the irrigation and cropping calendar, were then used to manually calibrate the classification algorithm. That is, the required variables and their corresponding thresholds of the classification were set based on the NDVI and Ts temporal evolutions of the selected sample points in conjunction with temporal information on crop planting, maturity, and harvest. Nemani and Running [23] established a conceptual diagram showing the seasonal trajectories of different land cover types in the Ts-NDVI space. In this study, we developed a similar Ts-NDVI space of the selected sample points.

3.4. Verification of Results

Due to the lack of ground truth data on the cultivated area, we relied on the field report data of the cultivated area for verification of the methodology. First, the developed algorithm was validated by comparing the cultivated area from the developed method with the field reported data. Second, the whole process was repeated for another season. Although the field reported data on the cultivated area were available for all six irrigation blocks, the most northern irrigation block (Hadaliya) was excluded from the analysis as it was very difficult to delineate it from the available irrigation network maps. Unlike the other five irrigation blocks, which are bordered by two main canals (one from the south and the other from the north), in the Hadaliya block, there is only one main canal to its south (Figure 1). The lack of a main canal to its north and the disappearing of the Gash River near the Hadaliya block made it difficult to clearly delineate it from Google Earth.

4. Results and Discussion

4.1. Identification and Selection of Sample Points

Visual interpretation of a conveniently selected image was done to confidently identify and select the pixels of interest. Figure 3 depicts the identification and selection of sample points of interest from the pseudo natural color, MNDWI, and NDVI maps. In the pseudo natural color map, the flooded area (flooded due to irrigation, hereafter referred as irrigated area) is shown as dark-blue, forest as green, and dry-bare area as brown. As expected, the MNDWI and NDVI enhanced the presence of water (flooded surface) and vegetation, respectively. In both images, white is high value and black is low value. The very sharp MNDWI image indicated a very high moisture difference in the GDSIS during the irrigation period. This proved the potential of employing MNDWI to identify typical irrigated areas. Although the pseudo natural color, MNDWI, and NDVI maps showed a distinction between the irrigated, bare-dry, and forest areas, information from standard false color and Google Earth were also used to increase the confidence of identifying and selecting the sample points of interest.

Considering the analysis of the visual interpretation, one could attempt to estimate or map the seasonal irrigated area. However, as the irrigation period of the GDSIS extends from July to September, irrigated areas irrigated after 28 July 2009 could not be identified from Figure 3. Therefore, extraction and mapping the seasonal irrigated area from only a single image acquired within the irrigation period would be incomplete. Moreover, the aim of this study was to map the seasonal cultivated area, which is usually less than the irrigated area [30], because no sowing of crops takes place on the portion of the irrigated area that: (i) either accumulated too little moisture to supply the total seasonal crop water requirement; or (ii) was not properly prepared for crop cultivation. In other words, in the GDSIS, irrigation is necessary, but not every irrigated area has the sufficient condition for cultivation. Considering the above situation, the selection of typical irrigated areas was done from areas with high MNDWI (high water) and low NDVI (no vegetation or no mesquite infestation). Hence, the selected typical irrigated areas are safely assumed as samples for cultivated areas. Accordingly, 150 typical points (50 each) from the irrigated, dry-bare, and forest area were selected and analyzed for the development of an algorithm that extracted the cultivated area. The selection of dry-bare points, as shown in Figure 3, was made outside, but near the GDSIS to avoid any possibility of being irrigated after 28 July 2009.

4.2. Building the Algorithm

Nemani and Running [23] established a conceptual diagram showing the seasonal trajectories of different land cover types in the Ts-NDVI space. In this study, we developed a similar Ts-NDVI space of the selected sample points as shown in Figure 4. As the purpose of this study was to extract the cultivated area in the GDSIS, Figure 4 is derived from the crop growing season in the GDSIS (September to February). Hence, the mean surface temperature and mean NDVI refer to the seasonal mean computed from satellite images within the crop growing period. However, the selection of sample points was done from an image from the irrigation period (in this case, 28 July 2009).

According to [23], forests maintain relatively lower temperature in comparison to crops. Forests—with their deep root systems—tend to dissipate more energy by transpiration through much of the growing season and, therefore, maintain canopy temperatures close to air temperature. Moreover, forests are aerodynamically rough and can dissipate energy efficiently as sensible heat, therefore, maintaining low surface temperatures. Furthermore, forests generally have higher canopy covers and low fraction of exposed soils, again contributing to low radiometric temperature [23]. On the other hand, crops show a higher surface temperature for two reasons: (1) under low-to-moderate wind speed, crops are aerodynamically smooth when compared to forests, therefore, the high aerodynamic resistances suppress sensible heat transfer, resulting in higher surface temperature; and (2) incomplete canopies through much of the growing season allow more radiation to penetrate and heat the underlying soil, thus contributing to high surface temperature. These phenomena can be clearly seen in Figure 4.

The selected typical dry-bare sample points exhibited a very low NDVI and high Ts. This result was consistent with the study by [23], where they set the threshold value of less than 0.15 mean NDVI and greater than 35 °C Ts for the barren surface. Therefore, these sample points must be representative of land surface that is neither cultivated nor covered with forest during that particular season.

To investigate the seasonal NDVI and Ts dynamics, the mean values of all the selected pixels in their respective features at each image date were computed and are depicted in Figure 5. As expected, Figure 5a shows that the forest pixels had relatively high NDVI value throughout, even before sowing (e.g., 28 July 2009) and after harvest (e.g., 05 February 2010). As clearly shown in Figure 5b, the forest pixels maintained a relatively low temperature (less than 32 °C) throughout the season. Moreover, the variability of NDVI for the forest was lower when compared to that of the crops (cultivated). Hence, the selected forest pixels must be evergreen forests, most likely mesquite trees, which are very dominant in the Gash Delta.

On the basis of the above analysis, the threshold-based rule to extract the cultivated area (shown in Figure 6) was developed. Different studies, such as [20,23], have employed many variables (mean, maximum, minimum, range, etc.) derived from many indices. In this study, however, primarily owing to the non-complex land cover of the GDSIS and in search of a simple algorithm, only two variables (mean Ts and NDVI) were found sufficient, as clearly seen in Figure 4. The threshold values of the selected variables were fixed manually from Figure 4. Caution was taken to maximize the inclusion of the cultivated points while minimizing the inclusion from other sample points (Dry-Bare and Forest). Accordingly, an NDVI of 0.2 to 0.4 and Ts of 30–35 °C were set to best define the region containing the cultivated sample points for the 2009–2010 season. This region was expressed in the form of a decision tree algorithm as shown in Figure 6. The classification decision tree was made to have two simple branches: the first branch, with variables and their corresponding threshold values, produces the cultivated area pixels; and the second branch lumps all the pixels excluded from the first branch as uncultivated.

In this study, our main purpose was to estimate (extract) the cultivated area, not to classify the land use/cover of the study area. At the same time, the selected land features in Figure 3 do not represent the whole land cover of the study area, but are sufficient for analyzing the seasonal NDVI and Ts dynamics to extract the cultivated area. Therefore, although only the NDVI threshold seemed to suffice in Figure 4, to ensure exclusion of any land cover other than those represented by the dry-bare and forest, adding the Ts threshold was necessary. The upper limit of the Ts, for example, is important in excluding short lived grasses grown in the insufficiently irrigated, hence uncultivated area. On the other hand, the lower limit of Ts was set to exclude, for example, well irrigated areas that were uncultivated due to poor land preparation (infested by trees).

4.3. Verification

The algorithm was set in ERDAS IMAGINE 8.4 to map and estimate the total cultivated area of the respective seasons in the GDSIS. The orientation of the cultivated fields (Figure 7) was found to match the orientation of Misgas (local name for irrigation unit) identified from Google Earth and the design layout of the irrigation system network.

To quantitatively assess the accuracy of the algorithm, we compared the extracted cultivated areas with field reported cultivated areas (

Table 1. Verification results of cultivated area for the 2009–2010 season.

Method of Estimation	Cultivated Area (ha)
	Kassala	Mekali	Degain	Tendelai	Metateib
Field report	5064	4876	5621	5371	2924
Remote sensing	4954	4648	5845	5330	3286
Error (ha)	−110	−228	224	−41	362
Error (%)	−2	−5	4	−1	12

). Although we acknowledged the inaccuracies associated with field data, we relied on the cultivated areas reported by [29], which seemed to be the most complete report (detailing the target irrigation area, actual irrigated area, and actual cultivated area for each block from 2004–2011), and hence was considered to be the best available data for verification.

As shown in Table 1, the cultivated areas estimated by the algorithm were very comparable to the field reported cultivated area. The algorithm slightly overestimated the cultivated area for the Degain irrigation block by 224 ha (4%) and for the Metateib irrigation block by 362 ha (12%). On the other hand, it marginally underestimated the cultivated areas of the other three irrigation blocks, Kassala by 110 ha (2%), Mekali by 228 ha (5%), and Tendelai by 41 ha (1%). These were in very good agreement, which demonstrates the potential of using remotely sensed data to estimate the cultivated area in the GDSIS.

To test its applicability, the developed approach was repeated for the 2004–2005 season. This season was selected merely on the basis of the availability of many cloud free Landsat-7 ETM+ images. Out of all these images, only the image dated 31 August 2004 fell in the irrigation period and was used for selecting the sample points. This date was near the end of the irrigation period, and as a result, sharp MNDWI, pseudo natural, and standard false color images could not be produced. This ultimately affected the quality of the sample data. This can be noted from Figure 8, where the boundary between the cultivated and forest points are not as clear as those in Figure 4 (2009–2010 season). Although the NDVI and Ts plots showed the same pattern, it also showed the necessity of modifying the threshold values slightly. Accordingly, the threshold values for the 2004–2005 season were set at a mean NDVI of 0.2 to 0.35 and mean Ts of 30–34 °C, as shown in Figure 9. The two cultivated sample points (blue, triangle mark) with a mean NDVI greater than 0.35 were neglected to avoid the inclusion of many sample points from the forest (green, rectangle mark). In other words, a mean NDVI upper limit of 0.35 was more representative of the cultivated sample points than a mean NDVI upper limit of 0.4.

The developed algorithm (Figure 9) was used to extract the cultivated area for the 2004–2005 season and was compared with the field reported cultivated areas as shown in

Table 2. Verification results of cultivated area from the 2004–2005 season.

Method of Estimation	Cultivated Area (ha)
	Kassala	Mekali	Degain	Tendelai	Metateib
Field report	4066	3459	5397	1806	2054
Remote sensing	4485	4156	4746	1471	1532
Error (ha)	419	697	−651	−335	−522
Error (%)	10	20	−12	−19	−25

. The algorithm overestimated the cultivated area of the Kassala irrigation block by 10% and Mekali irrigation block by 20%, whereas it underestimated those of the Degain irrigation block by 12%, Tendelai irrigation block by 19%, and Metateib irrigation block by 25%. Although the errors across all irrigation blocks of the 2004–2005 season were larger than those of 2009–2010, the authors still believe that these errors are reasonably acceptable. The overall good agreement of the estimated cultivated area with the field reported data demonstrates the potential of applying the proposed methodology in the GDSIS.

Assuming that other sources of error were the same for both seasons, the relative poor performance by the algorithm for 2004–2005 was primarily due to the quality of the sample points. As pointed out when setting the threshold values of the algorithm, the quality of the sample points was dependent on the date of the satellite image used for the identification and selection of the sample points. This implies the acquisition of a cloud free satellite image during the irrigation period, from which sufficiently watered areas could be clearly and confidently identified, was necessary for the accuracy of the developed algorithm.

Although the Ts-NDVI space for both seasons showed the same pattern, the resulting threshold values were slightly different. The influence of this small difference was examined by changing the NDVI threshold upper limit for the 2004–2005 season from 0.35 to 0.4. This change decreased the performance of the algorithm with an error in some of the irrigation blocks as high as 50%. Similar influence was also observed when the Ts threshold upper limit was set at 35 °C instead of 34 °C. This implies that while the model framework was simple and intuitive, the accuracy of the cultivated area mapping results was highly sensitive to the NDVI and Ts threshold values. Hence, fixing the threshold values using sample points as demonstrated by Figure 4 and Figure 8 was necessary to estimate the cultivated area of each season.

5. Conclusions

In this study, satellite remote sensing data were employed to develop a methodology for mapping and estimating the seasonal cultivated area in the GDSIS. Visual interpretations followed by NDVI and surface temperature (Ts) analysis of selected images within the crop period were performed to develop a threshold-based classification algorithm. Due to the specific nature of the GDSIS, the visual interpretations were found to be very effective in selecting sample points. The multi-date image analysis using NDVI and Ts plots of the selected sample points demonstrated the importance and capability of combining information from multiple bands of Landsat-7 ETM+ to uniquely characterize the cultivated area.

Verification of results by comparing the estimated cultivated area with the field reported area demonstrated the capability of the developed algorithm and, hence, the potential of our approach. Moreover, repeating the methodology in a different season proved that the seasonal cultivated area of the GDSIS could be mapped and estimated solely from Landsat-7 ETM+ images. The application of this simple methodology, which requires basic knowledge of remote sensing techniques and uses freely available satellite data, can help the GDSIS management office to regularly estimate the cultivated area at low cost. The main limitation of this study is the use of field reported cultivated area for the verification of results. Moreover, this study focused only on the cultivated area, however, estimation of the seasonal areas that are irrigated but uncultivated due to little water application and estimation of uncultivated areas due to poor land preparation are also very important to the management of the system. Therefore, future works should extend the current study to include the above areas and improve the verification process using ground truth data on cultivated areas.