Sub-pixel Area Calculation Methods for Estimating Irrigated Areas.

The goal of this paper was to develop and demonstrate practical methods forcomputing sub-pixel areas (SPAs) from coarse-resolution satellite sensor data. Themethods were tested and verified using: (a) global irrigated area map (GIAM) at 10-kmresolution based, primarily, on AVHRR data, and (b) irrigated area map for India at 500-mbased, primarily, on MODIS data. The sub-pixel irrigated areas (SPIAs) from coarse-resolution satellite sensor data were estimated by multiplying the full pixel irrigated areas(FPIAs) with irrigated area fractions (IAFs). Three methods were presented for IAFcomputation: (a) Google Earth Estimate (IAF-GEE); (b) High resolution imagery (IAF-HRI); and (c) Sub-pixel de-composition technique (IAF-SPDT). The IAF-GEE involvedthe use of "zoom-in-views" of sub-meter to 4-meter very high resolution imagery (VHRI)from Google Earth and helped determine total area available for irrigation (TAAI) or netirrigated areas that does not consider intensity or seasonality of irrigation. The IAF-HRI isa well known method that uses finer-resolution data to determine SPAs of the coarser-resolution imagery. The IAF-SPDT is a unique and innovative method wherein SPAs aredetermined based on the precise location of every pixel of a class in 2-dimensionalbrightness-greenness-wetness (BGW) feature-space plot of red band versus near-infraredband spectral reflectivity. The SPIAs computed using IAF-SPDT for the GIAM was within2 % of the SPIA computed using well known IAF-HRI. Further the fractions from the 2 methods were significantly correlated. The IAF-HRI and IAF-SPDT help to determine annualized or gross irrigated areas (AIA) that does consider intensity or seasonality (e.g., sum of areas from season 1, season 2, and continuous year-round crops). The national census based irrigated areas for the top 40 irrigated nations (which covers about 90% of global irrigation) was significantly better related (and had lesser uncertainties and errors) when compared to SPIAs than FPIAs derived using 10-km and 500-m data. The SPIAs were closer to actual areas whereas FPIAs grossly over-estimate areas. The research clearly demonstrated the value and the importance of sub-pixel areas as opposed to full pixel areas and presented 3 innovative methods for computing the same.


Introduction
Pixel size plays an important role in area computations, especially when coarser-resolution data are used. For example, a single Advanced Very High Resolution Radiometer (AVHRR) 10-kilometer pixel constitutes a full pixel area (FPA) of 10,000 hectares and Moderate Resolution Imaging Spectroradiometer (MODIS) 500-meter constitutes an FPA of 25 hectares. So, in many cases, only a fraction of a coarse resolution pixel falls under a particular land cover. The sub-pixel areas (SPAs) represent actual areas. The SPAs are computed by multiplying FPAs with irrigated area fraction (IAF). A comparative study for china [1] in estimating the areas derived from AVHRR Version 2.0 International Geosphere-Biospere Programme (IGBP) DIScover [2] dataset showed that about half of the DIScover cropland pixels had less than 60 % fractional cropland cover within a pixel size of 1-km. The pixel was named "irrigated" because it has certain percentage of area within the pixel which is irrigated-which can, typically, vary from a nominal 10 % to 100 %. It is, thereby, obvious that counting whole pixels can lead to over estimation of actual areas [3]. In the AVHRR 10-km, this issue becomes even more critical since every pixel encompasses 10,000 hectares. The implication of using FPAs in place where SPAs need to be reported is of significant importance in many applications such as water use calculations, food production estimates, and global scenario modeling.
Nevertheless, the coarser resolution imagery still remains the only practical data for global or regional studies. However, they invariably result in significant errors as a result of mixed land cover composition [4]. A number of methods have been used for un-mixing the sub-pixel in the coarser resolution imagery. These methods include [5]: (a) artificial neural network (ANN), (b) mixture modeling, and (c) fuzzy c-means classification. Even though ANN is the most accurate method; accurate co-registration and the availability of a training data set are real problem [5]. The other methods of sub-pixel area estimations include [6][7][8][9]: (a) regression based approaches, (b) high resolution imagery (HRI) use, and (c) groundtruth fractions. Fang et. al. (1998) estimated rice areas from AVHRR data by using higher resolution Landsat data by linking them through a linear statistical model [10]. Comparing the cropland area [11] on a smaller landscape-scale study using SPOT VGT and Landsat TM showed an increase in the Landsat TM computed area due to the type of the landscape. Gallego et. al. (1993) and Gonzalez et. al. (1997) used high-resolution satellite images to evaluate crop areas through regression estimator in an area frame survey [6,7]. DeFries et. al. (1996) derived percent forest cover in each of the AVHRR 8-kilometer grid cell by two different methods over central Africa and found that the percent cover estimated by classified MSS scenes were more accurate than the multiple linear regression and regression method [12]. Quarmby et. al. (1992) used linear mixture modeling for estimating crop area for a region of 2500 km 2 using multi-temporal AVHRR and recommends testing this technique to estimate crop area at national or continental scale [13]. All these studies prove that the sub-pixel fraction is essential in estimating the percentage area of a particular land cover and land use.
There are two, contrasting, inferences on the relationship between area and resolution. This is, mainly, as a result of whether the areas are calculated using SPAs or FPAs. This is discussed below taking irrigated cropland areas as an example. First, Ozdogan and Woodcock (2006) imply that finer the resolution of the imagery lesser is the cultivated area [14]. This is because, at finer resolution one can separate non-agricultural areas such as roads, settlements, barren areas, and fallow areas from cultivated agricultural areas. These areas can be significant and often as high as 30 to 40 %. In coarser resolution, the pixel will be considered irrigated when they are dominated by and\or significantly occupied by irrigated areas; but not necessarily entirely occupied by irrigated areas. So, in reality coarse resolution irrigated pixels, many times, consists of irrigated cropland areas as well as noncropland areas (e.g., roads, permanent fallows, settlements), leading to over-estimation of irrigated cropland areas if we consider the entire FPA as the actual area. This will lead to over-estimation of irrigated areas.
Second, finer the spatial resolution of the imagery greater is the area [8] because at finer resolution one can capture all or most of the fragmented and scattered cultivated areas where as at coarserresolution there is a good likelihood that one may miss fragmented irrigated areas completely unless the fragments within the pixels are highly significant and the irrigated cropland area fractions are determined accurately through sub-pixel de-composition. In this scenario, the coarse resolution FPA will not map fragmented cropland areas as croplands, thus under-estimating the cropland areas.
Based on the above two contrasting scenarios, the need for further investigation on the relationship between the areas and the resolutions are quite obvious. This will require building more reliable, robust, and practical methods of determining area fractions (AFs) that will lead to determining accurate areas or sub-pixel areas (SPAs).
Given the above discussions, the overarching goal of this paper was to develop methods for establishing IAFs that will lead to determining SPAs. In order to practically illustrate the methodology development, we took the global irrigated area map (GIAM; http://www.iwmigiam.org). Three distinct and unique methods of irrigated area fractions (IAFs) are discussed and illustrated. The sub-pixel irrigated areas (SPIAs) are then calculated by multiplying the full pixel irrigated areas (FPIAs) with IAFs. The study will compare the FPIA and SPIA computed from the GIAM map with actual irrigated areas obtained from the National systems.

Methods
The sub-pixel irrigated area (SPIA) computation methodology was illustrated using global irrigated area map at 10-km resolution (GIAM 10-km; Figure 1) and an irrigated area map for India derived from MODIS 500-m resolution. First, the GIAM which was produced using nominal 10-km resolution remote sensing data in conjunction with a number of other secondary data ( Figure 1, http://www.iwmigiam.org; 8) was used to compute SPIAs and compare it with areas determined using the national statistics: Where, IAF is irrigated area fraction and FPIA is full pixel irrigated area. The FPIA's are computed directly from the digital images using any commonly used remote sensing and\or geographic information systems (GIS) software packages.
Second, MODIS 500-m data was used to determine SPIAs and compare them with national statistics.
The methods for determining the IAFs remain the same for various resolutions; illustrated, mostly, for only 10-km resolution below.

IAF-GEE
Google Earth (http://earth.google.com/) provides a large volume of sub-meter to 4 meter, very high resolution (VHR), imagery. In addition, the ability to instantly "zoom into" any spot on earth is an attractive feature for "viewing" landscape and its details in fine resolution. When zoom into a particular class (Figure 2a; already labeled as irrigated in GIAM map as shown in Figure 1) in Google Earth VHR imagery, the fraction of area actually irrigated can be determined (Figure 2b). IAF-GEE is defined as the irrigated area fraction (IAF) of areas that are irrigated at any given point of time plus areas that are left follow but are "equipped" for irrigation at the same point of time.
The irrigated area fractions using Google Earth Estimates (IAF-GEE) were determined for each of the 28 classes ( Figure 1). The procedure involved distributing 30 to 50 well distributed random points (4 points illustrated in Figure 2a) for every class and "zooming into" these points (Figure 2a) in Google Earth to determine the IAF-GEE (Figure 2b) for each point. It is important to take note that the area has already been labeled irrigated in GIAM [8] and hence the entire area was considered irrigated in FPIA. However, the actual area irrigated was represented by SPIAs. Since the resolution of GIAM was 10-kilometers, the interpretation window in Google Earth was also fixed to 10-km by 10-km ( Figure 2a). Within the interpretation window of 10-km by 10-km, the percentage area irrigated was estimated by visual interpretation (Figure 2a). Once we have the IAF-GEE for every point, they are averaged to get one IAF-GEE for the class (illustrated taking 4 points in Figure 2b). Adequate numbers were selected to represent the class area and distribution in different parts of the world. Larger the size of the class, greater was the number of points. Also, when class was distributed in different parts of the world, as far as possible, points were selected to represent various locations in the world. The same approach is then repeated for all 28 classes (Figure 1b) leading to IAF-GEE values for these classes (Table 1). When we randomly zoom into points for a class where sub-meter to 4-meter data is absent in Google Earth, we select alternative points.
This method was earlier illustrated for establishing sub-pixel areas of forest canopies mapped using AVHRR 8-km by DeFries et. al., (1996) using Landsat MSS [12].  taking images that coincide with peak vegetation growing dates of the season (Figure 3c).

IAF-SPDT
The irrigated area fraction using sub-pixel de-composition technique (IAF-SPDT) was based on the RED-band versus near-infrared (NIR-band) reflectivity of the pixels in a class and their location in the brightness-greenness-wetness (BGW) plots ( Figure 4). The figure 4 shows the spectral reflectivity of every pixel of class 4 of the GIAM 28 class map for the season 1 (June-September) by plotting AVHRR band 1 min (absorption maxima) versus AVHRR band 2 max (reflection maxima). Through this effort, every pixel in the class was assigned a particular IAF percentage ( Figure 4) based on where it occurs in the BGW plot. The general rule is that the IAF percentages are highest in greenness zone, lower in brightness and wetness zones, and lowest near soil lines. But, there are significant exceptions to it that needs to be noted. One such example was the flooded rice pady, the significantly lower reflectivity in AVHRR band 1 and band 2 as a result of background water even though 100% pixel area was irrigated. The pixels will then cluster in an area between wetness and greenness zones. When assigning percentage area irrigated, we will retain 100% for this intermediate zone as well as the peak greenness zone in SPDT; since both have 100% area irrigated even when reflectivity differ significantly. The exact irrigated area percentages of pixels were determined based on observing the composition of pixels falling on different portions of SPDT in: A. groundtruth data and digital photos, B. high-resolution images, C. extensive literature review showing relationships between spectral indices and percent cover [17][18][19][20][21][22][23], and D. relative positioning of the RED and NIR reflectivity of pixels of a class in 2-dimensional feature space (2-d FS) SPDT [8,15].
Greater the understanding one has between the percent irrigated area versus band reflectivity, greater the reliability in assigning IAF percentages and the associated sub-pixel areas calculations. Depending on where the pixel falls in SPDT plot and IAF is determined.

Total area available for irrigation (TAAI)
The TAAI does not consider the intensity of irrigation. The TAAI is defined as the area irrigated at any given point of time and the area left fallow at the same point of time.

Annualized irrigated area (AIA)
The AIA considers the intensity of irrigation. The AIA is defined as the area irrigated during different seasons: season 1 + season 2 + continuous year-round irrigated crops.
For each of the GIAM 28 classes, the seasonality or intensity of (e.g., single crop, double crop, and year-round crop) irrigation was established by plotting the NDVI time-series ( Figure 5) which was a proxy for crop calendar and crop growth dynamics.

Results and Discussion
The irrigated area fractions (IAFs; Table 1) from the 3 methods were determined for each of the 28 GIAM classes. The IAF-HRI and IAF-SPDT were determined for season 1, season 2, and continuous year round crops (Table 1 and 2).  Table 1 and 2.
The sub-pixel irrigated areas (SPIAs) for the GIAM 28-classes ( Figure 1) were computed (Table 3) by multiplying full pixel irrigated areas (FPIAs) with irrigated area fractions (IAFs) ( Table 1). The total area available for irrigation (TAAI), or the net irrigated area, was determined by taking IAF-GEE of each class and multiplying it with FPIAs of the classes. The annualized irrigated areas (AIAs), or the gross area, was determined by taking the average of the IAF-HRI and IAF-SPDT of each class and multiplying them with FPIAs of the classes for season 1, season 2, and year-round continuous crops (see crop calendar of classes in Table 2) and summing their areas ( Table 3).
The SPIA and FPIA of the TAAI are represented in the Table 3. The FPIAs of the world was 589 million hectares or Mha (rounded off to nearest million) when compared with the SPIAs of 399 Mha, a difference of 190 Mha ( Table 3). The SPIA of the AIAs was 480 Mha ( Table 3). The AIAs considered cropping intensity (Table 2) and was calculated by multiplying FPIA with average IAF-HRI and IAF-SPDT. The fractions estimated by two independent methods ( Figure 6) were essential in developing confidence in IAFs and in refining them when one sees significant difference between two methods for any class.
The validity of these remote sensing based numbers were compared with the reported irrigated areas from the national statistics synthesized by the United Nations Food and Agricultural Organization and the University of Frankfurt (FAO/UF) [24,25]. The FAO/UF reported "area equipped for irrigation" (equivalent of GIAM TAAI) as 278 Mha. The differences between the FAO/UF and IWMI GIAM SPIAs (Figure 7b and 8b) were significantly smaller than the differences between the FAO/UF and IWMI GIAM FPIAs (Figure 7a and 8a). This was due to inadequate accounting and/or complete absence of informal irrigation (e.g., ground water, small reservoir, tanks) statistics in the national census, uncertainties in IAFs of GIAM, and methodological and definition differences in mapping irrigated areas. Nevertheless, it was obvious that SPIAs (Figure 7b and 8b) provide significantly improved estimates of areas when compared with FPIAs (Figure 7a and 8a).

FPIAs and SPIAs at AVHRR 10-km resolution versus National statistics (FAO Aquastat)
For the purpose 40 leading irrigated area countries that together encompass about 90% of all global irrigation were considered. The FPIA and SPIAs from the GIAM 10-km were plotted against the FAO/UF as shown in the figure 7a and figure 7b respectively. The slope of the 1:1 line improved from a poor 0.36 (Figure 7a) to a decent 0.54 (Figure 7b). There was still a considerable difference between the national statistics compiled by FAO/UF and the remote sensing based results of IWMI GIAM. This was mainly as a result of inadequate accounting of informal irrigation (e.g., ground water, small reservoirs, and tanks) in the national irrigated area statistics [8]. The FPIAs provide gross over-estimation of areas [8,26,27] requiring SPIAs to determine actual irrigated areas.

FPAs and SPAs at MODIS 500-m resolution versus National statistics (FAO Aquastat)
The census-based national irrigated area statistics from the Ministry of Water Resources (MoWR), Govt. of India [27] were compared with the MODIS 500-m derived FPIAs (Figure 8a) and SPIAs (Figure 7b). The 1:1 line shows that the MoWR statistics was about 75% of the SPIA (Figure 8b) and only 52% of FPIA (Figure 8a), this clearly indicates the significantly better relationship with SPIA than FPIA. However, still there was a significant gap in areas between the remote sensing based estimates and the national statistics (Figure 8b and 7b). The causes for which are described in section 3.0. The results also showed an improved relationship between the national statistics with 500-m data (Figure 8b and 8a) when compared to national statistics and 10-km data (Figure 7b and 7a). A MODIS 500-m provides an area of 25 hectares per pixel compared to 10,000 hectare per pixel from an 10-kilometer data, an improvement of 400 times. So, it was expected that the MODIS areas were more precise and nearer to the actual areas compared to AVHRR based estimates.

Uncertainties and errors in SPA estimates
Uncertainties and errors can creep in the SPIA estimates if sufficient care was not taken in IAF-GEE, IAF-HRI, and IAF-SPDT estimates. The global SPIA by summing season 1, season 2, and continuous year round cropping when computed by multiplying the FPIA with average of IAF-HRI and IAF-SPDT was 480 Mha (Table 3). The SPIA was: (a) 484 Mha when only IAF-HRI was used, and (b) 476 Mha only when IAF-SPDT was used. The difference in area calculated from the 2 methods was about 8 Mha; which accounts to less than 2% of the total GIAM area. This was a very low margin of error. We did re-visit the fractions for every class that had significant difference between two methods and found that these differences were inevitable as a result of the differences in methods. However, the use of two or more methods help us to compare IAFs obtained from different methods for a class and when there were significant differences in IAFs of different methods, which assist in re-evaluate the fractions for the those class leading to improved IAFs. The resultant IAFs reduced the uncertainties and errors in SPIA estimates. Multiple methods improved the robustness of the IAFs. The IAF from the two methods ( Figure 6) provide significant correlation. For IAF-SPDT, further improvement may be possible if we have a better understanding of the percentages for every class in SPDT plot. This will require better groundtruth knowledge of the class in consideration and building relationships between spectral reflectivity and/or NDVI with IAF. Even a small error in IAF estimation can accelerate errors in areas. However, use of 3 methods should provide sufficient robust estimates. Uncertainities can be further reduced if we have greater field knowledge of IAFs in a given area.
The IAF-GEE found to be powerful and easy to use and has the ability to "zoom into" any spot in the world instantaneously. The synoptic view provided by Google makes the IAF-GEE estimates very reliable. The IAF-GEE was ideal to calculate the total area available for irrigation (TAAI) as it assisted in determining the area that was irrigated at the time of estimation as well as area left fallow but equipped for irrigation at the same point of time. There were three significant limitations in IAF-GEE method. First, the absence of very high resolution imagery (VHRI) for every spot in the world. Indeed, for large portions of the world VHRI were absent. Second, the absence of multidate images. Third, the presence of images of varying dates. Third, there is an inherent assumption that the areas "equipped for irrigation" is always not irrigated at particular point of time but will be irrigated at other times. The likelihood that these areas are under permanent fallow has not been fully investigated. IAF-HRI too has the same limitations as second and third points mentioned under IAF-GEE. The IAF-SPDT found to be potentially most powerful method considering the limitation in other toiw methods. Its strength were lies in the: (a) ability to determine areas for every pixel, (b) scale the areas logically in a 2-dimensional brightness-greenness-wetness (BGW) plot, (c) make corrections to areas as we learn more about the class reflectivity and/or NDVI versus percent cover that would aid to assign IAF percentages better in the SPDT plot ( Figure 4).

Conclusions
The study established the importance and the need for sub-pixel area (SPA) computation for determining actual areas from coarse resolution remote sensing data and espoused practical methods for computing the same. The full pixel areas (FPAs) obtained from coarse-resolution remote sensing data significantly over-estimate areas. Three unique methods of SPA estimation were developed and illustrated by taking a global irrigated area map (GIAM). The sub-pixel irrigated area (SPIA) was computed by multiplying full pixel irrigated area (FPIA) with irrigated area fractions (IAFs). The uncertainties and errors in SPIA computation were directly proportional to errors in irrigated area fractions (IAFs), given that the FPIA remains constant. The three distinct methods of IAFs computations were: (a) Google Earth estimate (IAF-GEE); (b) High resolution imagery (IAF-HRI); and (c) Sub-pixel decomposition technique (IAF-SPDT).
The IAF-HRI and IAF-SPDT fractions were found to be useful in computing SPIAs and account for seasonality or intensity (e.g., first crop, second crop, continuous crop). Thus, they provide net as well gross irrigated areas. The significant correlation was found between IAF-HRI and IAF-SPDT. The overall, the areas from all classes, determined by these 2 methods differed by less than 2%. However, there were individual classes that showed much higher significant differences in areas between 2 methods. Given this fact, using more than one IAF was recommended in order to: (a) reduce uncertainties and errors in areas, and (b) provide more robust estimate of areas. The IAF-GEE stands on its own and was useful for computing net irrigated areas without considering intensity or seasonality.
The SPIAs provided significantly better relationships with the national statistics than FPIAs. The FPIAs were also shown gross over-estimate of areas. The paper: (a) highlighted the importance of computing sub-pixel areas for determining accurate areas, and (b) developed and demonstrated 3 unique and practical methods of computing sub-pixel areas.