Temperature of coastal waters and of watercourses from ASTER images

This paper presents an algorithm to improve the spatial resolution, from 90 m to 30 m, of the thermal mapping of small bodies of water or near coasts obtained from the ASTER satellite sensor. The entire procedure is based only on ASTER images. The first part of the work deals with the physical and mathematical basis on which the algorithm was constructed, the schema of the main steps and the methods of validation of the algorithm. In the second part two applications of the algorithm are shown, the first on the area of the delta of the Po River (Italy), the second on the lagoon of Venice (Italy).


Introduction
The temperature of coastal and river waters is a physical quantity very important in various areas of environmental concern. Applications that make use of this quantity are many and among them can be mentioned the study of river, lake and sea ecosystems, the fish farming and shellfish, the study of eutrophication, the study of mass and energy budgets and the monitoring and characterization of the heat discharged from industrial plants in water bodies [Richter et al., 1986;Giardino et al., 2001;Ritchie et al., 2003;Cherkauer et al., 2005;Tcherepanov et al., 2005;Handcock et al., 2006]. This last case also directly involves regulatory issues, for the Italian law, for example, the references are contained in Legislative Decree no. 152/06. Remote sensing by satellites and aircrafts, given its well known ability to obtain thematic maps, could be a very valuable tool in this area. However, it is important to remember that remote sensing provides skin temperatures, while bulk water temperature is generally required. One of the main challenges from an environmental monitoring perspective in using satellite derived skin temperatures is understanding bulk water temperatures from the water skin temperatures [Imberger and Patterson, 1990;Hook et al., 2003;Setlinger, 2005]. Airborne sensors, such as MIVIS (Multispectral Infrared and Visible Imaging Spectrometer) and DAIS (Digital Airborne Imaging Spectrometer), have the most appropriate technical characteristics for the purpose; nevertheless their employment is limited by cost and organization of airborne campaigns. As for the satellite remote sensing, the main limiting factor is undoubtedly the low spatial resolution (hereafter also indicated by SR) of the sensors currently in use. Indeed, most of the sensors that operate in the thermal infrared (TIR) installed on board of satellites have pixel size greater than or equal to 1 km, such as MODIS (Multispectral Infrared and Visible Imaging Spectrometer, AVHRR (Advanced Very High Resolution Radiometer) and AATSR (Advanced Along-Track Scanning Radiometer). Exceptions in this sense are represented by the ETM+ (Enhanced Thematic Mapper Plus) sensor and the ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) sensor. Unfortunately, ETM+ experienced a failure of its Scan Line Corrector mechanism in May 2003, as a consequence, at the moment its use is problematic. These considerations about sensors comparison are focused on the aim of this work, which is spatial resolution of thermal mapping. Concerning other aspects, as temporal resolution, sensors like AVHRRR or MODIS offer better performances than sensor like ASTER and ETM+. Thus their well known usefulness for environmental applications is not diminished. ASTER [NASA, 2004] acquires images in five bands of TIR at a spatial resolution of 90 m and in ten bands distributed between 500 nm and 2500 nm at spatial resolutions varying between 15 m and 30 m. Despite the better spatial resolution, in many cases a spatial detail of 90 m can be still inadequate. This is evident, for example, considering the typical transverse dimensions of the rivers or of the discharge canals of industrial plants. Given the above considerations, the main purpose of the work presented in this paper is the construction of an algorithm to improve the spatial resolution, from 90 to 30 m, of the thermal mapping obtained from ASTER on coastal and river waters. Techniques for improving the spatial resolution of remote sensing images are generally based on algorithms for image fusion [Pohl and Van Gender, 1988;Chavez et al., 1991;Wald, 1999]. In most cases these techniques are aimed at improving the spatial resolution of multi-spectral images by "injecting" into them the spatial details extracted from an image acquired with a better spatial resolution, as those acquired by panchromatic sensors. In this case, the two basic requirements are: 1) the acquired area must be the same; 2) the images at high and low spatial resolution must have enough spectral overlapping to ensure a good correlation between them. The most famous family of these techniques is Pan Sharpening. In other cases the information at high spatial resolution is extracted from other sources, such as images acquired by other sensors on different dates or thematic maps in raster or vector format. An example of this approach, very significant also for this paper, is the work of Setlinger et al. [2008]. This work deals with sub-pixel water temperature estimation using a priori knowledge of water boundaries derived from vectorized water features. The algorithm presented in this paper falls into this last category. It has been developed for the improvement of ASTER images acquired in the TIR based on information on the type of coverage of the area extracted from ASTER images acquired in the VNIR (Visible and Near Infrared) and with spatial resolution varying between 15 m and 30 m. This algorithm is still in testing phase and currently the greatest efforts are aimed at verifying the results. Below a summary description of the structure of the algorithm followed by two applications on the coastal areas of the lagoon of Venice and the Po River delta are given.

Structure of the algorithm
The algorithm is applicable to "coastal" pixels (CP) of the images acquired by ASTER in the TIR bands. A TIR pixel is classified as "coastal" if the sub-pixels, with a size of 30 m that compose it, are partially on water and partially on land (Fig. 1). The algorithm is based on two assumptions: Hp1: There is a good degree of correlation between the radiance emitted in the TIR by the pixels of the image with a set of variables obtainable from ASTER VIS-NIR images at spatial resolutions of 15 m or 30 m; Hp2: The correlation assumed in Hp1 is invariant to small changes of scale, such as for the ratio 3:1 (e.g. 90 m: 30 m). The hypothesis Hp1 leads to the basic relationship of the algorithm, according to which the radiance (L) emitted by a given pixel in a TIR band can be expressed using the following linear combination: where f w and f s are the water and non-vegetated soil cover fractions respectively, z v is a variable that describes the presence of vegetation. a 0 , a w , a s e a v are the coefficients of the linear combination and are computed for each CP considered. The variables that can be used to describe the vegetation and retrievable by VIS-NIR images are many. The algorithm developed is able to use the fraction of vegetated cover soil (f v ) and the vegetation indices NDVI (Normalized Difference Vegetation Index), PVI (Perpendicular Vegetation Index) and SAVI (Soil Adjusted Vegetation Index). All the mentioned indexes are well known in literature and are defined as function of surface reflectances at about 0.8 µm (NIR, ASTER band 3) and at about 0.65 µm (Red, ASTER band 2) respectively. For each CP, it is considered the set of pixels that surround it (neighbourhood U ). Two kind of data are extracted from U: the TIR radiances at the spatial resolution of 90 m (L 90 ) and the values of the variables f w , f s and z v , calculated at the spatial resolution of the VIS-NIR images (SR = 15 m) and subsequently upscaled at 90 m (average value). The local values (i.e. only valid for that CP) of the coefficients a 0 , a w , a s and a v are estimated from this set of variables using a multiple linear regression method: The neighbour of the considered CP (U ) is a square box centred on that pixel. The width of the box is one of the parameters of the algorithm. In this work a width of 5 pixels were used. , the superscript "reg" indicates the estimated values) is less than 0.15 W m -2 μm -1 sr -1 , which roughly corresponds to a temperature variation of 0.9 K, at the reference wavelength of 10 μm and at the reference temperature of 300 K. The value of 0.9 K has been chosen considering the 1 K temperature accuracy in the 270 K -240 K range of ASTER [Abrams et al., 2002].
The previous steps of the algorithm are aimed at assessing the radiance (L 30 ) of the "coastal" pixels for which the hypothesis Hp1 is verified. The next step is to construct the thermal map (radiance and brightness temperature) for all the water surfaces (pure water pixels). This phase is divided into the three main steps summarized below and sketched in Figure  2: 1) For all the CPs for which the hypothesis Hp1 is verified the radiances L 30 of the sub-pixels (SR = 30 m) entirely covered by water are reported on the map. As an example, in Figure 2 three CPs of this kind, located at the raw-columns R1C3, R2C2 and R3C2, and the corresponding L 30 values are shown. For further statistics, these sub-pixels are labelled as of K 1 class. 2) The original radiance values at SR = 90 m (L 90 ) are assigned to the sub-pixels fully covered by water (not CP, as the pixels R2C3 and R3C3 in Fig. 2).
3) The values of each sub-pixel obtained in the previous step and completely surrounded by pure water sub-pixels are replaced with a new value computed by bilinear interpolation of the values in a 3x3 sub-pixels box cantered on the sub-pixel (L 30,int in Fig. 2). This operation is done to make more realistic the spatial patterns of the thermal map and is based on the reasonable consideration that on water surfaces there are not strong spatial discontinuity [Sentlinger et al., 2008]. It is important to realize that, in this way, the values of some of the sub-pixels considered in this step (marked with a star in Fig. 2) are affected by the L 30 values obtained in step 1. These last sub-pixels are labelled as of K 2 class.

Validation
In general, the quantitative assessment of the goodness of algorithms is a step required but very difficult. The most reliable methods of validation are certainly those based on comparison with measured values. Often, however, as in the cases presented in this paper, ground truth data are not available, so two alternative methods were used. The first method is based on the test of the hypothesis Hp1. The second is a simulation method and is heavily based on the hypothesis Hp2.
The test of Hp1 consists in: -Computation of the number of CP (N CP ) and of the number of CP that verify Hp1 (N CP * ); -Computation of the number of sub-pixels of K 1 class (N K1 ) and of K 2 class (N K2 ). These values are indicators of the algorithm efficiency. -Evaluation of the goodness of the regression for the CPs that verify Hp1. The goodness of the regression is estimated by considering the mean values, over the N CP * pixels, of: the coefficient of determination <R 2 >; the coefficient of multiple linear correlation <r M >; the standard error of the fit <SE>. In many cases, when the radiances L 90 are nearly constant in the regression neighbour U, the first two statistical parameters may not make much sense. For this reason the mean values are computed using only the CPs, between the N CP * obtained, for which the standard deviation of L 90 in U is greater than 0.1 W m -2 sr -1 µm -1 . The number of the CPs that verify this condition are indicated by N CP,stat . The second test follows an approach of upscaling-reconstruction-comparison already used in literature [Teggi et al., 2003;Wald et al., 1997], structured in the following steps:

Case studies
Below two applications of the algorithm are shown: the first on an area of the lagoon of Venice, the second on an area of the delta of the River Po. The ASTER image used for the Po River was acquired on August 22, 2008 at 10:10 GMT. The studied area, shown in Figure 3.a has an extension of 44010 m (x direction) per 17010 m (y direction). This area is characterized by the presence of long stretches of coast, with many lagoons and coves, often swampy, and long branches of various watercourses which width is very variable. The main branch of the Po has a width varying from 200 m to 600 m, while the main branches forming the delta have widths as low as a few tens of meters. The ASTER image used for the lagoon of Venice was acquired on September 5, 2007 at 10:07 GMT. The studied area, shown in Figure 4.a has an extension of 18540 m (x direction) per 7200 m (y direction). This area includes a large part of the Venice lagoon, where there are wetlands, islands (including the urbanized island of Venice) and the industrial area of Porto Marghera, characterized by wide canals. Both the ASTER images, before be furnished in input to the algorithm, passed through different elaboration steps, the most important of which were atmospheric corrections, classification and computation of the images of f w , f s , f v , NDVI, PVI and SAVI. The atmospheric corrections were done using a procedure based on the radiative transfer models MODTRAN and 6S [Gillespie at al., 1998;Vermote et al., 1997;Bogliolo et al., 2004] and local meteorological data as input for the models. The procedure requires in input the radiance measured by ASTER and gives in output the surface reflectance (VIS-NIR) or the surface brightness temperature (TIR). In the VIS-NIR range the procedure uses MODTRAN to model gas absorption and 6S to compute atmospheric scattering contributions (the reason of the use of both the radiative transfer models is purely linked to software aspects, more details can be found in the work of Bogliolo et al. [2004]). In the TIR range only MODTRAN is used. The three vegetation indexes were computed using the well known relations [Richardson and Wiegand, 1977;Huete, 1988 where ρ 2 and ρ 3 are the surface reflectance in the ASTER bands 2 and 3 (0.66 µm and 0.82 µ); α and β are the slope and the intercept of the "soil line" in the (ρ 2 , ρ 3 ) reflectance space; L is a constant that minimize the soil brightness influence. α and β were computed using an automated soil line identification procedure as that described by Fox et al. [2004].
Assuming intermediate vegetation densities, L was set to 0.5 according to Huete [1988]. The classification step consisted in assigning each pixel of the studied areas to one of the three classes non vegetated soil, vegetated soil and water. Pixels covered by water were individuated using the maximum likelihood method applied to ρ 2 , ρ 3 and NDVI images. The remaining pixels were divided in vegetated and non vegetated using an NDVI threshold value of 0.4. Having the classification images, the calculation of the cover fractions images at different spatial resolution was trivial. The algorithm presented above was applied to each of the TIR bands of ASTER. Since no significant differences were observed, only the results obtained using the 14th band (11.3 µm) of ASTER are shown in this work. Figure 3b and 4b show the temperature maps obtained for the Po River and for the lagoon of Venice respectively. The spatial enhancement of the temperature maps produced by the algorithm with respect to the analogous maps obtained using only the original TIR data acquired by ASTER is evident by comparing Figures 3b and 4b with Figures 3c and 4.c. These last two figures show the temperature maps obtained using only the data of the starting TIR radiance images, that is to say, applying only steps 2 and 3 of the algorithm.
In the first couple of maps many more river branches, coves and canals can be clearly distinguished and the temperature patterns follow better and more realistically the coastline with respect to the second couple of images. it is interesting to observe that the Canal Grande canal, which flows through the island of Venice, is not visible in either thermal mapping obtained. This canal has a width of 40 m -50 m (only in a few points it reaches 60 m). This implies that the neighbourhood U of the CP forming the canal are mostly composed by non water pixels and by a few fraction of mixed (water and non water) pixels. It has been verified, that in these cases, the hypothesis Hp1 is rarely verified and consequently, the results of the algorithm are almost always rejected. Table 1 reports the statistical parameters obtained with the first method of validation for each studied area and for each kind of variable used to quantify the vegetation. The first consideration that can be carried out from this table is that all the z v variables produced very similar results.    The variable f v performed slightly worse than the other variables but the differences can be considered unimportant. The hypothesis Hp1 is verified for most of the CPs: more than 80% in the Po image and more than 70% in the Venice image. The independent variables (regressors) and the predicted variable (radiance) used for the local multiple linear regressions shown a mean (computed over the N CP,stat CPs) multiple linear correlation coefficient close to 0.9. The local multiple linear regressions explained (average value) more than the 80%, for the Po case, and more than the 70%, in the Venice case, of the local variability. In all the cases, the mean value of the error of the fits was less than, or equal to 0.1 W m -2 µm -1 sr -1 , corresponding to 0.5 K. With respect to the temperature map obtained using only the original radiance data (that is to say, applying only steps 2 and 3 of the algorithm) the "higher spatial resolution" temperature map obtained using the algorithm contains more than 23000 (N K1 ) additional pixels (SR = 30m). Moreover, the temperature value of more than 10500 (N K2 ) pixels is affected by the new pixels added. Similar considerations can be made for the temperature map of Venice, in this case the pixels of class K 1 are more than 6700 and those of class K 2 are more than 3800. The results obtained using the second validation method are summarised in Table 2. Given the considerations made in the discussion of the previous validation method, only the case z v = NDVI is reported. For both the areas the comparisons between the original L 90 image and the reconstructed L 90,sim image produced similar results. Considering the pixels of K 1 class the image are practically unbiased, well correlated (r > 0.8) and with RMSD less than 0.1 W m -2 µm -1 sr -1 (0.5 K). The statistics does not change significantly when considering pixels of both class K 1 and class K 2 . These results show that the reconstruction of the original image can be considered satisfactory.

Conclusions
In this work an algorithm for improving the spatial detail of the maps of coastal waters and watercourses temperature derived from ASTER images has been shown. This algorithm produces a downscaling of the ASTER TIR images, from 90 m to 30 m, spatially more complete than that obtained by a classic interpolation technique. In the two case studies presented, the temperature maps obtained by the algorithm were compared with those obtained directly from ASTER TIR images using a classic interpolation method. This comparison showed that the former are much more detailed than the latter. In particular, there is a better definition of narrow watercourses, of small coves and of waters in rugged coastlines. The results achieved so far, which includes the two cases presented in this work, are satisfactory, although the experimental verification of the results is still weak due to the lack of ground-truth temperature measurements. Moreover, these performances indicate the possibility of using the algorithm for the downscaling from 90 m to 15 m of TIR ASTER images or for images acquired by other sensors.