International Journal of Applied Earth Observation and Geoinformation 23 (2013) 95–108
Contents lists available at SciVerse ScienceDirect
International Journal of Applied Earth Observation and Geoinformation journal homepage: www.elsevier.com/locate/jag
Downscaling of thermal images over urban areas using the land surface temperature–impervious percentage relationship W. Essa a , J. van der Kwast b,c,∗ , B. Verbeiren a , O. Batelaan a,d a
Department of Hydrology and Hydraulic Engineering, Vrije Universiteit Brussel, Pleinlaan 2, BE-1050 Brussels, Belgium Flemish Institute for Technological Research (VITO), Boeretang 200, BE-2400 Mol, Belgium UNESCO-IHE Institute for Water Education, Westvest 7, 2611 AX Delft, The Netherlands d Department of Earth and Environmental Sciences, K.U.Leuven, Celestijnenlaan 200e - bus 2410, 3001 Heverlee, Belgium b c
a r t i c l e
i n f o
Article history: Received 6 August 2012 Accepted 14 December 2012 Keywords: Land surface temperature Urban Thermal sharpening Impervious percentage DisTrad Landsat 7 ETM+
a b s t r a c t Intensive expansion and densification of urban areas decreases environmental quality and quality of urban life as exemplified by the urban heat island effect. For this reason, thermal information is becoming an increasingly important data source for integration in urban studies. It is expected that future spaceborne thermal sensors will provide data at appropriate spatial and temporal resolutions for urban studies. Until they become operational, research has to rely on downscaling algorithms increasing the spatial resolution of relatively coarse resolution thermal images albeit having a high temporal resolution. Existing downscaling algorithms, however, have been developed for sharpening images over rural and natural areas, resulting in large errors when applied to urban areas. The objective of this study is to adapt the DisTrad method for downscaling land surface temperature (LST) over urban areas using the relationship between LST and impervious percentage. The proposed approach is evaluated by sharpening aggregated LST derived from Landsat 7 ETM+ imagery collected over the city of Dublin on May 24th 2001. The new approach shows improved downscaling results over urban areas for all evaluated resolutions, especially in an environment with mixed land cover. The adapted DisTrad approach was most successful at a resolution of 480 m, resulting in a correlation of R2 = 0.84 with an observed image at the same resolution. Furthermore, sharpening using the adapted DisTrad approach was able to preserve the spatial autocorrelation present in urban environments. The unmixing performance of the adapted DisTrad approach improves with decreasing resolution due to the fact that the functional relationship between LST and impervious percentage was defined at coarse resolutions. © 2012 Elsevier B.V. All rights reserved.
1. Introduction In 1950 30% of the world’s population was living in urban areas, while this value will increase to 60% in 2030 (Golden, 2004). The growing population in urbanized areas causes an increase of anthropogenic heat. As a result the core of cities in temperate climate zones becomes warmer than their periphery, thus forming an urban heat island (UHI) (Voogt and Oke, 2003), which has an impact on the microscale and mesoscale climate. Heat-related morbidity and mortality are among the primary health concerns and are expected to increase as a result of climate change that increases the UHI effect (Johnson and Wilson, 2009).
∗ Corresponding author at: UNESCO-IHE Institute for Water Education, Westvest 7, 2611 AX Delft, The Netherlands. Tel.: +31 152151877. E-mail addresses:
[email protected],
[email protected] (J. van der Kwast). 0303-2434/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jag.2012.12.007
The UHI intensity is influenced by local and synoptic weather, season, time of day, size of the city and its geographical location, urban morphology, and anthropogenic heat (Kim and Baik, 2005). Therefore, a range of spatial and temporal scales should be considered in UHI studies. One of the great challenges for current UHI studies is how to support this research with techniques to retrieve land surface temperature data with high temporal and spatial resolution. This data has to support the scientific and policy makers’ community in understanding impacts, both positive and negative, to set policies for mitigating those impacts or to develop indicators to address the UHI effect within a sustainable development context. Therefore land surface temperature is a physical parameter of prime importance for a wide variety of scientific urban studies. In urban climatology, land surface temperature is the key input parameter used in algorithms estimating net radiation, sensible and latent heat flux, thermal inertia, surface urban heat island intensity, perceptible water, evapotranspiration, urban-induced surface runoff and surface moisture (Stathopoulou and Cartalis, 2009). Various spaceborne thermal sensors provide
96
W. Essa et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 95–108
useful information for urban studies. The choice of an appropriate sensor is currently a trade-off between spatial and temporal resolution: sensors with a high temporal resolution have a low spatial resolution and vice versa. Thermal sharpening is a promising technique for enhancing the spatial resolution of mixed pixels acquired by coarse resolution sensors. During the last ten years various thermal sharpening techniques for unmixing land surface temperature have been developed based on land cover. Thermal sharpening methods in remote sensing can be grouped into two categories. The first method, known as the physical downscaling method, is based on establishing a functional relationship between land surface temperature and ancillary data, which satisfies two conditions: (1) the functional relationship is physically meaningful and holds across different resolutions; and (2) the ancillary data can be easily scaled (Zhukov et al., 1999; Kustas et al., 2001; Liu et al., 2006; Agam et al., 2007a,b; Liu and Pu, 2008). The second method, known as the statistical downscaling method is based on the direct estimation of subpixel land surface temperature from higher spatial resolution ancillary data without referring to physical models at different resolutions (Ottlé et al., 2008; Stathopoulou and Cartalis, 2009). Testing of the latter method is still in a preliminary stage for applications in urban areas. In a previous study Essa et al. (2012) evaluated the relationship between land surface temperature of urban pixels and 15 indices and concluded that the relationship with impervious percentage gives the best results. This conclusion is also in agreement with findings from Yuan and Bauer (2007), who suggest that impervious percentage accounts for most of the land surface temperature variation in urban areas. Furthermore, Essa et al. (2012) and Yuan and Bauer (2007) conclude that there is a strong linear relationship between land surface temperature and impervious percentage, independent of the season of image acquisition. In contrast, the relationship between land surface temperature and NDVI varies with the seasons (Kaufmann, 2003). Therefore, Essa et al. (2012) suggested to use impervious percentage instead of NDVI in the DisTrad sharpening algorithm (Kustas et al., 2001) when applied to urban areas. They demonstrated the improvement for sharpening of a Landsat 7 ETM+ image of Dublin. Due to the low temporal resolution of Landsat images, the current study evaluates the adapted DisTrad algorithm developed by Essa et al. (2012) for sharpening images with a spatial resolution comparable with daily MODIS land surface temperature products. While in Essa et al. (2012) only 25% of the pixels with the lowest coefficient of variance (CV) within a coarse resolution pixel were used in the regression, this study evaluates different threshold values. Furthermore, a detailed analysis of the spatial autocorrelation of the unsharpened and sharpened products is presented. The images sharpened with the original and the adapted DisTrad algorithm are compared with each other and with non-sharpened images using visual, statistical, and geostatistical analysis methods.
2. Study area and data This study uses the 60 m resolution land surface temperature image that was derived by Essa et al. (2012) as a basis for deriving coarse resolution images and reference images for evaluation of the sharpening. The image was derived from a Landsat 7 ETM+ scene from May 24th 2001 covering 2700 km2 , bounded by the longitudes 6◦ 00 and 6◦ 48 W and latitudes 53◦ 07 and 53◦ 35 N. It covers the urbanized Mid-East of the Greater Dublin region (Fig. 1). This area is characterized by heterogeneous urban landscape patterns, which makes the study area ideal for evaluating the adapted DisTrad procedure.
Fig. 1. Ireland with inset figure indicating the study area: the Greater Dublin region.
A 30 m resolution impervious percentage map derived by Van de Voorde et al. (2010) will be used in the regression analysis. More details on the derivation of the impervious percentage map can be found in Van de Voorde et al. (2010). 3. Methodology 3.1. Sharpening experiment Ideally, the adapted DisTrad approach applied to a coarse resolution land surface temperature product needs to be evaluated with high resolution land surface temperature data acquired at the same moment. Image products from different sensors, however, differ in acquisition time and date, data level, processing algorithms, temporal integration, geometric accuracy, etc. Therefore, in this study a Landsat 7 ETM+ image was aggregated to 960 m resolution and evaluated for downscaling using the approach proposed in Essa et al. (2012). In this way the performance of the algorithm can be optimized and verified in a best case scenario, independent of errors introduced by differences between sensors. Landsat 7 ETM+ has both VNIR and TIR data at a relatively high resolution acquired at the same date and time. This makes this sensor ideal for evaluating the downscaling, without having to take care of differences between TIR and VNIR channels from different platforms. The performance of the adapted DisTrad approach using impervious percentage to derive high resolution land surface temperature data over urban areas was evaluated by comparing the sharpened results with the original DisTrad approach using NDVI (Kustas et al., 2001). In addition, a baseline case (UniTrad), representing no sharpening, i.e. no incorporation of high-resolution information, was also determined by assigning each sub-pixel the land surface temperature of the corresponding coarse resolution pixel. The Landsat 7 ETM+ land surface temperature image, impervious percentage image, and NDVI image were aggregated by a spatial averaging function to images with a resolution of 60, 90, 120, 240, 480, and 960 m, which will be referred to as “observed images”. Aggregation using spatial averaging is often used to simulate remote sensing images. Dominguez et al. (2011), for example, aggregated higher resolution ATLAS data (albedo, NDVI and land surface temperature observed at 10 m resolution) for evaluation of the High-resolution Urban Thermal Sharpener (HUTS) method. Agam et al. (2007b) aggregated aircraft data (NDVI and land surface temperature observed at 10 m resolution) for evaluation of the TsHARP method. Both studies used data aggregated to 90 m resolution to simulate ASTER satellite data. In this study the aggregated
W. Essa et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 95–108
97
Semivariogram analysis was used to compare the spatial autocorrelation within the land surface temperature images produced by DisTrad and within the observed Landsat 7 ETM+ land surface temperature images. Similarity in spatial variability of land surface temperature for the observed and sharpened products can give an indication for downscaling accuracy. The experimental semivariogram represents the spatial structures and the spatial variation in a graphical way (Isaaks and Srivastava, 1989; Goovaerts, 1997). Sharpened and observed land surface temperature images at all targeted resolutions were analysed by calculating the semivariogram at each resolution. A detailed discussion on the use of semivariograms in remote sensing can be found in van der Meer (2012). The spatial autocorrelation of all sharpened and observed land surface temperature images was evaluated by fitting a semivariogram model to the full data set of each land surface temperature image.
4. Results and discussion Fig. 2. Flowchart showing the procedure of comparing ‘observed’ images, spatially averaged from Landsat 7 ETM+ , with sharpened images.
960 m resolution land surface temperature image was used to be downscaled to the other five higher resolutions. The sharpened and observed images were compared (Fig. 2). 3.2. DisTrad sharpening procedure DisTrad is a sharpening technique for downscaling land surface temperature developed by Kustas et al. (2001) and adapted by Agam et al. (2007a,b). DisTrad has been developed to sharpen land surface temperature over agricultural areas in order to improve evapotranspiration estimation for crop monitoring purposes. Therefore, the original DisTrad procedure relies on vegetation information estimated using indices such as NDVI (Kustas et al., 2001) or fractional vegetation cover (Agam et al., 2007a,b) that are derived from shortwave bands available at higher resolutions than the thermal bands. The linear regression of those vegetation indices with land surface temperature is used as an estimator of land surface temperature at multiple resolutions. It is assumed that the regression equation is resolution independent. Vegetation cover information is typically derived at spatial resolutions 4–16 times finer than the resolution of the thermal band (Agam et al., 2007b). For the application of DisTrad, it is important that a wide range of land surface temperatures and impervious percentages are present within the image scene, in order to develop a representative regression function (Kustas et al., 2001; Agam et al., 2007b). Therefore, DisTrad is most effective when applied to heterogeneous areas with images acquired between mid-morning and mid-afternoon. The goal of this study is to improve the sharpening of coarse resolution daily land surface temperature products over urban areas using DisTrad based on the assumption that a unique, scale independent linear relationship exists between land surface temperature (dependent variable) and impervious percentage (independent variable). For a detailed description of the DisTrad procedure the reader is referred to Kustas et al. (2001) and Essa et al. (2012). 3.3. Geostatistical analysis The level of agreement between the observed and sharpened land surface temperatures was evaluated by matching land surface temperature histograms of both images, using the normal distribution model. Mean-absolute-error (MAE) and root mean square error (RMSE) were used as goodness-of-fit measures.
4.1. Regression analysis The selection of a subset of the most uniform high resolution pixels within a coarse pixel using the CV, as suggested by Kustas et al. (2001), was evaluated for two objectives. The first objective was to evaluate the applicability of using 25% instead of all pixels at the five observed image resolutions (60, 90, 120, 240 and 480 m). Only pixels within the urban mask were considered. The second objective was to determine the optimal index pixel fraction to be used in the downscaling application by evaluating the use of 100%, 75%, 50% or 25% of the pixels with the lowest CV values within the simulated 960 m image. The analysis was done for both impervious percentage (Fig. 3) and NDVI (Fig. 4). Fig. 3 indicates that the correlation between land surface temperature and impervious percentage is higher when 25% of the pixels are selected based on the lowest CV values instead of using all pixels. The selection procedure removed the redundant land surface temperature pixels and increased the inter-variability of impervious percentage and land surface temperature values. The regression parameters and statistical variables e.g. mean, standard deviation and median were fairly insensitive to a change of resolution (Table 1) over all classes, supporting the assumption that the land surface temperature–impervious percentage relationship holds for one sensor scene at multiple spatial resolutions. The original DisTrad method, using NDVI, was compared with the adapted method using impervious percentage for deriving sub-pixel land surface temperature over an urban environment. A higher correlation between NDVI and land surface temperature was expected for the acquisition date, because of the high vegetation cover in May, which strongly contributes to the land surface temperature variability. The results, however, show that land surface temperature has a higher correlation with impervious percentage (Fig. 3) than NDVI for all studied resolutions, except 480 m (Figs. 3 and 4), with all pixels used in the regression. The selection of 25% of the NDVI pixels with the lowest CV does not improve the correlation for NDVI (Fig. 4), compared to impervious percentage (Fig. 3). It even results in a worse correlation with land surface temperature for the 60–120 m resolutions and does not improve much the correlation for the 240–480 m resolutions, compared to using 100% of the NDVI pixels. For this reason, it is expected that it is preferable to use impervious percentage instead of NDVI even for seasons with less vegetation cover in the urban environment. Selection of uniform pixels of land cover based on a fraction of pixels with the lowest CV is more effective when the land cover is more heterogeneous within its landscape. Consequently, impervious percentage, with its larger value distribution
98
W. Essa et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 95–108
Fig. 3. Relationship between Landsat 7 ETM+ land surface temperature (LST) (◦ C) and impervious percentage (I%) for all pixels within the urban mask (right) and the same relationship for 25% of the urban pixels with the lowest CV (left). The red arrows show the mixing temperature problem in urban areas, where the same I% has a wide range of land surface temperatures.
in the urban mask, shows better results than NDVI. NDVI images show a wider value range in the natural and agricultural environment, while impervious percentage shows more heterogeneity in the urban environment. The wide range of impervious percentage values in the urban mask, which reflects the wide heterogeneity range of the manmade urban fabrics, correlated to a wide range of land surface temperature values and was an essential precondition for using DisTrad for downscaling land surface temperatures. NDVI shows a short value range in the urban environment, resulting in
worse correlations with land surface temperature, which confirms the conclusion of Agam et al. (2007a). The relationship between the indices and land surface temperatures are plotted in Figs. 5 and 6 respectively. The correlation between impervious percentage and land surface temperature at 960 m resolution (Fig. 5) was highest (R2 = 0.78 and a RMSE = 1.28 ◦ C) when 50% of the pixels with the lowest CV were used in the regression analysis. However, the other scenarios of sharpening with DisTrad, with thresholds of 100%, 75% and 25%,
W. Essa et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 95–108
99
Fig. 4. Relationship between Landsat 7 ETM+ land surface temperature (LST) and NDVI for all pixels within the urban mask (right) and the same relationship for 25% of the urban pixels with the lowest CV (left).
Table 1 Correlation coefficients and goodness of fit for the correlation between observed land surface temperature (LST) and impervious surface percentage (I%) calculated for all classes (urban, vegetation and bare soil). Resolution (m)
60 90 120 240 480
Statistics
Correlation LST and I%
Mean
Std. Dev.
Median
a0
a1
R2
RMSE
30.63 30.63 30.63 30.63 30.63
3.18 3.06 3.00 2.75 2.45
30.00 30.07 30.07 30.17 30.23
29.50 29.17 29.10 28.99 29.05
7.479 8.357 8.736 9.197 9.075
0.5876 0.5285 0.5470 0.5079 0.5408
6.505 4.080 3.491 2.603 1.858
100
W. Essa et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 95–108
Fig. 5. Relationship between observed Landsat 7 ETM+ land surface temperatures versus different selections of impervious percentage pixels at the coarse resolution (960 m).
also show an acceptable correlation, with respectively R2 = 0.67, 0.69, and 0.76. The correlation between NDVI and land surface temperature was higher (R2 = 0.65 and 0.63 respectively and the RMSE = 1.54 ◦ C and 1.49 ◦ C respectively) when 100% and 75% of the pixels with the
lowest CV were used in the regression analysis (Fig. 6) compared to other scenarios where 50% and 25% (R2 = 0.40 and 0.34 respectively and the RMSE = 1.59 ◦ C and 1.84 ◦ C respectively) were used as a threshold. These correlations, however, are still lower than all results for impervious surface percentage.
Fig. 6. Relationship between observed Landsat 7 ETM+ land surface temperatures versus different selections of NDVI pixels at the coarse resolution (960 m).
W. Essa et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 95–108
To conclude, the impervious percentage proved to be a useful parameter for the DisTrad downscaling algorithm since it has a higher correlation with land surface temperature than NDVI. Based on this analysis, the impervious percentage index shows the highest correlation with land surface temperature at 960 m resolution, when a subset of 50% of the pixels with the lowest CV is used. The following regression equation was obtained for this scenario and will be used in this study: Tˆ960 = 29.76 + 8.716I960
(1)
An appropriate decision about the degree of the statistical fit between land surface temperature and the index value is important in order to get a proper accuracy of the sharpening. Essa et al. (2012) showed that the correlation between impervious percentage and land surface temperature improved when using 25% of the pixels with the lowest CV in the regression. In this study a threshold of 50% of the pixels with the lowest CV within a coarse resolution resulted in the highest correlation. The threshold value can be different for other cities depending on their structure. Further studies should investigate better selection techniques for urban pixels, which result in higher correlations with land surface temperature. The correlation should also be improved for pixels with different materials having different thermo-physical and thermochemical properties but the same impervious percentage for which the general scene-specific land surface temperature–impervious percentage relationship is not valid. This was found to be a shortcoming in the proposed DisTrad procedure. Future studies have to develop classification methods with the aim to classify impervious percentage into classes based on their thermo-physical properties, and then downscale those classified pixels separately based on a different linear regression with land surface temperature. Using a land-cover map that is related to physical thermal properties using lookup tables could be an important improvement to physical downscaling algorithms and techniques. Another possibility is to develop new indices that are able to estimate the thermal weight of each urban pixel. Such indices can replace DisTrad or be integrated within DisTrad to be used for downscaling land surface temperature. Multiple regression analysis using different remote sensing indices as independent variables is another option for evaluating downscaling of urban land surface temperatures. 4.2. Scene characteristics It is assumed that the correlation between NDVI or impervious percentage and land surface temperature is higher when the regression equations are specifically calculated for particular land-cover classes. Tables 1 and 2 show two scenarios using masks of different land cover. In the first scenario (Table 1) the regression equation is based on impervious percentage within all classes, while in the second scenario (Table 2) the regression equation is calculated for impervious percentage only within the urban mask. Results show that the first scenario has a higher error (RMSE) than the second scenario for all resolutions. The correlation (R2 ), however, is higher for the first scenario for 60, 90 and 120 m resolutions, while the second scenario shows a higher correlation at 240 and 480 m resolution. Essa et al. (2012) show that using regression coefficients of the relationship based on pixels in all classes results in an underestimation of land surface temperatures in urban areas compared to the observed land surface temperature images. Furthermore, the first scenario shows a systematically lower land surface temperature (∼−3.17 ◦ C) in all classes, averaged over all resolutions, compared to the second scenario, which reflects that the estimation of land surface temperature in urban areas is biased because of the cooling effect by evapotranspiration in the vegetation class in the first scenario. Therefore, the second scenario that applies the DisTrad downscaling methodology to the urban mask of Dublin is preferred.
101
Results show that it is instructive to examine the DisTrad procedure in urban areas based on the relationship of land surface temperature and impervious percentage within the context of an urban mask, as each land-cover type has different thermal properties, which lead to a different correlation of the land surface temperature. Applying DisTrad to pixels with a mixed land cover will lead to a bias of the land surface temperature estimation. Meaningful regression parameters for the vegetation class could not be derived, therefore the sharpening algorithm did not retrieve additional useful subpixel information over areas of the Landsat 7 ETM+ scene composed of natural vegetation. In contrast, in this study, significant high-resolution land surface temperature data could be successfully retrieved from its relationship with impervious percentage over the urban areas within the scene. 4.3. Evaluation of sharpening The sharpening products are evaluated by comparing the results of sharpened land surface temperature images with the observed land surface temperature images at the corresponding resolutions (Table 3). In addition, the land surface temperature–impervious percentage correlations before and after sharpening are compared, which is summarized in Tables 2 and 3. The sharpening was evaluated for the five investigated resolutions (60, 90, 120, 240, and 480 m) by comparing the regression parameters and coefficients of determination. Statistical measures within the urban mask were also calculated. The mean and median land surface temperatures of the impervious surface within the urban mask before and after sharpening were very similar (∼34 ◦ C) with a small overall underestimation of land surface temperature of about −0.30 ◦ C averaged over all resolutions. Furthermore, the standard deviations of the sharpened images were systematically decreased with an average of only −0.28 ◦ C averaged over all resolutions, compared to UniTrad images. This indicates a good performance of DisTrad to downscale land surface temperature over urban areas. The correlation (R2 ) of the land surface temperature with impervious percentage before and after sharpening is systematically increased by 0.33 for all resolutions due to the regression procedure. DisTrad preserves the resolution dependant variation, which is indicated by the regression coefficients and goodness-of-fit measures with their similar values and trends especially at low resolutions. Moreover, the regression error determined by the root mean square error after sharpening shows a significant systematic decrease of 1.18 ◦ C averaged over all resolutions. The goodness-of-fit measures (R2 and RMSE) obtained from the correlation analysis between sharpened and observed land surface temperature, are shown in Table 3. Based on these measures it can be concluded that DisTrad was able to retrieve reliable high resolution land surface temperature information for the urban areas within the scene at all resolutions. At a resolution of 60 m, which has the highest subpixel heterogeneity, the sharpened land surface temperature variability could still be resolved (R2 = 0.40), while correlation is significantly increased proportionally to coarser resolutions, which are more homogeneous (e.g. R2 = 0.84 at 480 m). The regression parameters were fairly insensitive to resolution, supporting the assumption that the relationship between land surface temperature and impervious percentage is scale independent in urban areas. RMSE is found to be higher at higher resolutions (3.8 ◦ C) and lower at lower resolutions (1.0 ◦ C), which results from a higher sub-pixel heterogeneity at higher resolutions of impervious percentage. Fig. 7 shows the sharpening result when Eq. (1) is used to sharpen the 960 m land surface temperature image to 60 m, the highest resolution considered in this study. The sharpened image is smoothed compared to the observed 60 m land surface
102
W. Essa et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 95–108
Table 2 Correlation coefficients and goodness of fit for the correlation between observed land surface temperature (LST) and impervious surface percentage (I%) calculated within the urban class. Resolution (m)
Statistics
60 90 120 240 480
Correlation LST and I%
Mean
Std. Dev.
Median
a0
a1
R2
RMSE
34.02 33.81 33.94 33.79 33.45
3.40 3.99 3.14 2.83 2.63
34.10 34.08 34.08 33.92 33.54
30.31 29.81 29.57 29.47 29.38
6.571 7.786 8.431 8.957 9.286
0.1708 0.3569 0.4449 0.5479 0.6326
4.733 2.955 2.558 1.967 1.612
Table 3 Correlation coefficients and goodness-of-fit for the correlation between observed and sharpened land surface temperature (Corr. observed–sharpened) using the regression equation of the correlation between land surface temperature and impervious percentage (Corr. LST–I%) within the urban class. Resolution (m)
60 90 120 240 480
Statistics
Corr. LST–I%
Corr. observed–sharpened
Mean
Std. Dev.
Median
a0
a1
R
RMSE
a0
a1
R2
RMSE
34.94 33.78 33.65 33.34 32.99
3.02 2.95 2.93 2.87 2.82
34.2 34.0 33.9 33.5 33.1
28.88 28.89 28.80 28.76 28.62
9.307 9.632 9.815 10.000 10.380
0.6211 0.7948 0.7979 0.8001 0.8023
2.677 1.374 1.335 1.287 1.257
8.923 7.580 6.547 5.773 4.768
0.7401 0.7818 0.8141 0.8392 0.8673
0.3967 0.5286 0.6013 0.7293 0.8405
3.809 2.518 2.162 1.521 1.060
temperature image. Fig. 7 also shows the scatterplot of all urban pixels in the observed and the sharpened image. 4.4. Descriptive analysis The sharpened images are compared with the observed land surface temperature images at the corresponding resolutions (60, 90, 120, 240 and 480 m) by calculating the error difference (observed image minus sharpened image). Absolute error differences of 2 ◦ C represents 63%, 68%, 73%, 81%, and 88% of the area of the downscaled images of respectively 60, 90, 120, 240 and 480 m resolution. Large error differences are caused by mixed pixels located at the class boundary between urban areas and vegetation and/or between urban areas and water bodies. Fig. 8 shows two examples (a and b) with extreme differences between sharpened and observed land surface temperatures. The first example (Fig. 8a ) represents misclassified pixels located at the boundary between vegetation plots and water bodies. Such pixels (blue pixels in Fig. 8a ) were cooler in the observed land surface temperature image than in the downscaled images, due to the fact that for those pixels the impervious percentage was overestimated by the unmixing procedure. Consequently, the estimated
2
land surface temperatures were overestimated. Positive land surface temperature errors (red pixels in Fig. 8a ) are an indication of an underestimation of downscaled land surface temperatures due to the fact that downscaling was here applied to impervious percentage pixels with a high thermal emission caused by their thermo-physical properties (Section 4.6) and/or by additional sources of anthropogenic heat in various forms, which can be explicitly incorporated into the energy balance (e.g. heat from transportation, industry, households, etc.). Those pixels were found to be randomly distributed in this example and in the whole central part of the city of Dublin. The second example (Fig. 8b ) shows Dublin port. The red pixels are underestimated land surface temperature values derived with DisTrad. DisTrad was not able to detect the anthropogenic heat released from the intensive oil industry, transportation and port businesses characterizing Dublin port, because these pixels are outliers in the regression between land surface temperature and impervious percentage. Furthermore, DisTrad underestimated pixels concentrated in bigger clusters along the coastal edge, due to mixed pixels containing water and land. This is clearly present at all resolutions. Due to this problem urban coastal edge pixels need to be masked.
Fig. 7. Observed versus sharpened land surface temperature at 60 m resolution.
W. Essa et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 95–108
103
Fig. 8. (a) and (b) are the LANDSAT 7 ETM+ land surface temperature at 90 m resolution for two sub-urban areas of the study area, representing respectively examples of sub-urban areas with a boundary with vegetated areas and port of Dublin; (a ) and (b ) are the standard deviations of the sharpened products from the observed land surface temperatures respectively.
4.5. Geostatistical analysis In this study, all semivariograms for the observed and sharpened land surface temperature products (Fig. 9) can be described by combinations of two permissible models (Goovaerts, 1997), i.e. (1) an exponential model (Eq. (2)) and (2) a Gaussian model (Eq. (3)):
3h
(h) = c 1 − exp −
(h) = c 1 − exp
−
R 3h2 R2
(2)
(3)
where c is the sill and R the (practical) range. An exponential semivariogram corresponds to a spatial structure that quickly diminishes with distance as a decreasing exponential, while a Gaussian model indicates a much more persistent spatial structure which only diminishes gradually after a certain distance, similar to a Gaussian bell shaped curve. When a spatial structure is described using a combination of both basic models the total sill of the spatial variable is equal to the sum of the sills of the basic models. The relative contribution
of each basic structure can be quantified by the so called Q-value, which is the sill of the basic model divided by the total sill. The results indicate that the original impervious percentage map with a resolution of 30 m has a complex spatial structure that can be described by three basic models: two exponential and one Gaussian. The first exponential model has a very short range of only 60 m, but a high sill of 0.041 so that it accounts for 53% of the total spatial variation. This indicates a spatial structure of short distance, more precisely only two cells wide, which likely expresses impervious areas that are connected to neighbouring impervious areas as neighbouring buildings along streets. The second exponential model has a larger range of 600 m and accounts for 24% of the overall spatial variation. This probably reflects whole streets or housing blocks, with a typical size of 600 m. The Gaussian structure has a very large range of 5000 m, a sill of 0.018, and accounts for 23% of the overall spatial variation. This structure is persistent over large distances and probably reflects suburbs, which in Dublin are distributed along the coast and mayor traffic roads and are about 5 km in size. When the impervious area image is repeatedly smoothed by averaging, the short range structures are affected while the long range structure persists. The exponential structures are gradually smeared out, which implies that the sill
104
W. Essa et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 95–108
Fig. 9. Estimated semivariograms for the sharpened land surface temperature images and for the observed Landsat 7 ETM+ land surface temperature images at resolutions 60, 90, 120, 240, and 480 m.
decreases step by step while the range increases with each new image resolution. Hence, when the resolution becomes 960 m the first exponential structure has completely disappeared, while the second exponential structure has a much wider range of 2460 m
(i.e. 600 + 60 + 120 + 240 + 480 + 960 m), a sill of only 0.011, and a Qvalue of 38%. On the other hand, the Gaussian structure still has a sill of 0.018 and a range of 5000 m, and now accounts for 62% of the spatial structure.
W. Essa et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 95–108
The observed land surface temperature images have similar characteristics as the impervious percentage images, although the spatial structure is somewhat more blurred. For instance, the Landsat 7 ETM+ land surface temperature image with a resolution of 60 m has a semivariogram that can be described by one exponential model with a range of 600 m and one Gaussian model with a range of 6000 m. Hence, the first exponential model seems to combine the 60 m and 600 m exponential structures of the 30 m resolution impervious area image, and accounts for 69% of the overall spatial structure, while the Gaussian structure has now a range of 6000 m and accounts for the remaining 31%. So the land surface temperature image is also strongly conditioned by urbanization, but seemingly the spatial variation of the land surface temperature has somewhat dissipated. This becomes more evident when the Landsat 7 ETM+ land surface temperature images are smoothed by averaging. The exponential structure is gradually smoothed so that, at a resolution of 960 m, the exponential structure has a range of 2400 m and accounts for only 43% of the spatial structure, while the Gaussian structure still has a range of 6000 m and accounts for the remaining 57% of the spatial structure. This spatial structure very much resembles the spatial structure of the impervious area image with the same 960 m resolution. The sharpened land surface temperature images are obtained from the combination of impervious percentage images and the 960 m resolution residual image and, hence, will show characteristics of both images. The residual image is a combination of the 960 m resolution observed LST images and impervious percentage images, which have very similar structures, and consequently these spatial structures will also be present in the residual image. The spatial structure of the residual image can be described by an exponential model with a range of 2430 m and a sill of 1.9, which accounts for 40% of the spatial structure, and a Gaussian model with a range of 5500 m and a sill of 2.9, which represents 60% of the spatial structure. For the sharpened images it becomes more complicated, because these images are a combination of images of different resolutions. As a result three exponential structures can be identified in the sharpened images and one Gaussian structure. Two exponential models result from the impervious percentage image and have ranges that increase as a result of the impervious percentage image averaging. One exponential model results from the residual image and, hence, has a fixed range and sill in all images. Gradually the first two exponential structures decrease in importance as they are removed by the averaging of the impervious percentage images. The Gaussian structures of the impervious percentage images and residual image blend together, so that in all sharpened images one Gaussian structure can be identified with a fixed range of 6000 m and a constant sill of 3.4. This Gaussian structure is the most dominant feature of all LST images. In this case, where the sharpened images are compared to the observed Landsat 7 ETM+ images with semivariograms calculated for the same resolution (Fig. 9), the following conclusions can be drawn: the observed Landsat 7 ETM+ images are much more complex in their structure than the sharpened land surface temperature images, which is clearly reflected in more land surface temperature variation than in the sharpened images. The sharpened images have smaller total sills especially at 60, 90 and 120 m resolutions, while the shape and the range of all variograms are similar to the observed Landsat 7 ETM+ image at all resolutions. The Gaussian structure dominates the spatial variation in both the observed Landsat 7 ETM+ and sharpened images, especially for the lower resolutions. The next important spatial structure in both sharpened and observed Landsat 7 ETM+ images is an exponential structure with a range that increases from 600 m to about 2400 m, depending upon the resolution. Also this structure becomes more and more important for the low resolutions as other exponential structures with smaller ranges are smoothened out.
105
Table 4 Comparison of mean-absolute-error (MAE, ◦ C) for each sharpened land surface temperature image versus each observed land surface temperature image of the same resolution using UniTrad and the adapted DisTrad approach. Resolution (m)
DisTrad
UniTrad
60 90 120 240 480
1.98 1.73 1.58 1.28 0.97
2.28 2.25 2.00 1.60 1.10
The semivariogram analysis shows the potential of geostatistical techniques to be used for sharpening land surface temperature images. For example cokriging can be applied to low resolution images with a high resolution impervious percentage image as secondary variable. Alternatively, kriging of a low resolution image can be applied with an external drift, i.e. the mean of the spatial land surface temperature variation is considered to be a linear function of the high resolution impervious percentage image. Such techniques would make full use of the spatial autocorrelations of the images and the correlations between images, and also enable to obtain error bands for the estimates. It should however be noted that isotropic behaviour was assumed in the semivariogram analysis. Further research is needed to investigate if anisotropic behaviour exists and how significant this is for the different resolutions. The distribution of land surface temperature values in the scene was examined by plotting the observed Landsat 7 ETM+ land surface temperature distribution, the uniform land surface temperature distribution (UniTrad), and the sharpened land surface temperature distribution resulting from the adapted DisTrad approach (Fig. 10a). The 2.5% tails from the range of the input land surface temperature maps were excluded in order to analyse the 95% range. The distribution plots show that the distributions from results of the adapted DisTrad approach follow the observed land surface temperature and the mean value much closer than UniTrad for 60–240 m resolution images, while both the adapted DisTrad approach and UniTrad show similar results at 480 m resolution. Fig. 10b shows the normal distribution of the errors estimated by subtracting results from the adapted DisTrad approach and UniTrad images from the observed land surface temperature images. The errors for the results of the adapted DisTrad approach are symmetrically distributed around zero. The adapted DisTrad results show a higher fraction of pixels with a close to zero error, compared to UniTrad errors which are asymmetrically distributed at resolutions from 60 to 240 m. Results from the adapted DisTrad approach have a lower error than UniTrad results at all resolutions, which further demonstrates the feasibility of using the adapted DisTrad approach to sharpen land surface temperature images over urban areas, especially at the higher resolutions. The mean-absolute-error (MAE) was used to assess the magnitude of error divergence of results of the adapted DisTrad approach and UniTrad from the observed land surface temperature images. Results from the adapted DisTrad approach show a significant systematic decrease in MAE by −0.4 ◦ C averaged over the 60–240 m resolutions (Table 4). However, at the lowest resolution of sharpening (480 m) the adapted DisTrad approach results only a slight decrease in MAE compared to UniTrad (−0.13 ◦ C), because at coarser resolutions UniTrad is more similar to the simulated 960 m resolution image. The adapted DisTrad approach improved the results compared to UniTrad, which is visible in the regression coefficients and goodness-of-fit measures of the adapted DisTrad and UniTrad results plotted against the observed land surface temperature pixels (Table 5). The slopes of the regression lines of land surface temperature sharpened with the adapted DisTrad approach versus
106
W. Essa et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 95–108
Fig. 10. The left side (a) shows the distribution of the observed, the sharpened (DisTrad), and the uniformly disaggregated (UniTrad) land surface temperatures. The right side (b) shows the error distributions resulting from subtracting the observed land surface temperature from DisTrad and UniTrad. (a) and (b) are calculated at different resolutions (60, 90, 120, 240, and 480 m).
observed land surface temperature at resolutions from 90 to 240 m are closer to the 1:1 perfect agreement line than those of UniTrad (0.78–0.84 versus 0.76–0.81 respectively), while at the 60 and 480 m resolutions the slope of the DisTrad regression line shows a slightly worse fit with the 1:1 line than UniTrad (0.74 and 0.87 versus 0.75 and 0.88 respectively). The coefficients of determination (R2 ) are also higher for the results of the adapted DisTrad approach compared to the UniTrad results with a systematic increase of +0.13 ◦ C at resolutions from 90 to 240 m, and
increases to +0.08 ◦ C at resolutions 60 and 480 m. Moreover, the root mean square errors (RMSE) are also significantly lower for the adapted DisTrad results compared to UniTrad: −0.41 ◦ C at 60 m resolution, −0.33 ◦ C at resolutions from 90 to 240 m, and −0.22 ◦ C at 480 m resolution. From this it can be concluded that the adapted DisTrad sharpening algorithm considerably improves the land surface temperature downscaling over urban areas, compared to UniTrad, especially for the higher to middle resolutions (90–240 m).
W. Essa et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 95–108 Table 5 Regression parameters and goodness-of-fit measures for the correlation between observed land surface temperatures and DisTrad or UniTrad sharpened images at different resolutions. Resolution (m)
a0
a1
R2
RMSE
Slope
DisTrad 60 90 120 240 480
8.92 7.58 6.55 5.77 4.77
0.74 0.78 0.81 0.84 0.87
0.40 0.53 0.60 0.73 0.84
3.81 2.52 2.16 1.52 1.06
0.74 0.78 0.81 0.84 0.87
UniTrad 60 90 120 240 480
9.45 8.92 8.52 7.19 4.47
0.75 0.76 0.77 0.81 0.88
0.32 0.41 0.46 0.60 0.77
4.22 2.84 2.52 1.84 1.28
0.75 0.76 0.77 0.81 0.88
4.6. Error sources in the DisTrad algorithm Scatterplots of observed land surface temperature versus impervious percentage (Fig. 3) show a wider scatter in land surface temperature at high impervious percentages, and narrower scatter in land surface temperature at lower impervious percentages. Although this is also visible in the scatterplots where all impervious percentage pixels were used, the scatterplots using 25% of the pixels with the lowest CV emphasize this effect. This is caused by the presence of thermally mixed pixels at all resolutions, indicated by a red arrow in Fig. 3. Thermally mixed pixels have a large range of temperatures for the same impervious percentage. Although these pixels have the same impervious percentage, their materials differ in terms of thermo-physical and thermo-chemical properties such as heat capacity, thermal conductivity and thermal inertia. These pixels cause outliers in the land surface temperature–impervious percentage relationship. The intensity of this phenomenon was distinct at high impervious percentage (I% ∼ 100%). Pixels with a low impervious percentage show a smaller range in temperatures, because the variability in impervious surface materials is less in these pixels. 5. Conclusions The adapted DisTrad algorithm proposed by Essa et al. (2012) was evaluated for its feasibility to sharpen images with a spatial resolution comparable with daily MODIS land surface temperature products. In addition, the appropriate threshold for the selection of pixels with the lowest CV that are used in the regression analysis was determined. Although sharpening results were not evaluated using real images, the adapted DisTrad approach does appear to have potential for generating useful land surface temperature information over urban areas from thermal images at a high temporal and improved spatial resolution (e.g. when applied to MODIS land surface temperature products). DisTrad using the impervious percentage has a higher correlation with land surface temperature compared to NDVI at all investigated resolutions, even though the image was acquired in a season with high vegetation cover. Moreover, the correlation between land surface temperature and the index improved more for impervious percentage than for NDVI when only pixels with a low CV are used in the regression analysis. The best results were obtained when a subset of 50% of the pixels with the lowest CV are used in the regression analysis. High resolution land surface temperature information can be retrieved from low resolution sensors using DisTrad based on impervious percentage relationship with the urban land surface temperature. Downscaled land surface temperature information,
107
at higher resolutions i.e. 90–480 m, showed a higher correlation (R2 > 0.50) and smaller error (MAE <1.75 ◦ C) when compared with observed land surface temperature. The adapted DisTrad algorithm was able to reconstruct patterns of land surface temperature in urban areas which have a clearly defined land surface temperature–impervious percentage relationship, although R2 were generally low. Where impervious surfaces border the coast or vegetated areas, mixed pixels exist with large sub-pixel differences in thermal characteristics. The mixed pixel problem increases with pixel size. In areas where this is a significant problem, DisTrad cannot be a replacement for high-resolution thermal imagery. Nevertheless, in the absence of high resolution thermal imagery, DisTrad seems to be able to enhance the resolution of images over urban areas, and provides higher resolution land surface temperature distributions. Acknowledgements This research was part of a PhD study, which was supported by an ERASMUS MUNDUS Scholarship based on the bilateral agreement between Al-AQSA University (Palestine) and the EU. The research presented in this paper was supported by the Belgian Science Policy Office in the frame of the STEREO II programme – project SR/00/105. We are indebted to the MAMUD project for providing this research with the basic images, maps, and information. The authors of this research acknowledge the assistance of Prof. Dr. ir. F. De Smedt for calculating and discussing the semivariogram models. References Agam, N., Kustas, W., Anderson, M., Li, F., Colaizzi, P., 2007a. Utility of thermal sharpening over Texas high plains irrigated agricultural fields. Journal of Geophysical Research 112, 19110–19120. Agam, N., Kustas, W., Anderson, M., Li, F., Neale, C., 2007b. A vegetation index based technique for spatial sharpening of thermal imagery. Remote Sensing of Environment 107, 545–558. Dominguez, A., Kleissl, J., Luvall, J., Rickman, R., 2011. High-resolution urban thermal sharpener (HUTS). Remote Sensing of Environment 115, 1772–1780. Essa, W., Verbeiren, B., Van der Kwast, J., Batelaan, O., 2012. Evaluation of DisTrad thermal sharpening for urban areas. International Journal of Applied Earth Observation and Geoinformation 19, 163–172. Golden, J., 2004. The built environment induced urban heat island effect in rapidly urbanizing arid regions – a sustainable urban engineering complexity. Environmental Sciences 1, 321–349. Goovaerts, G., 1997. Geostatistics for Natural Resources Evaluation. Oxford University Press, New York, 483 pp. Isaaks, E.H., Srivastava, R.M., 1989. An Introduction to Applied Geostatistics. Oxford University Press, Oxford. Johnson, D.P., Wilson, J.S., 2009. The socio-spatial dynamics of extreme urban heat events: the case of heat-related deaths in Philadelphia. Applied Geography 29, 419–434. Kaufmann, R.K., 2003. The effect if vegetation on surface temperature: a statistical analysis of NDVI and climate data. Geophysical Research Letters 30, 2147. Kim, Y.H., Baik, J.J., 2005. Spatial and temporal structure of the urban heat island in Seoul. Journal of Applied Meteorology 44, 591–605. Kustas, W., Norman, J., Anderson, M., French, A., 2001. Estimating sub-pixel surface temperatures and energy fluxes from the vegetation index–radiometric temperature relationship. Remote Sensing of Environment 85, 429–440. Liu, Y., Hiyama, T., Yamaguchi, Y., 2006. Scaling of land surface temperature using satellite data: a case examination on ASTER and MODIS products over a heterogeneous terrain area. Remote Sensing of Environment 105, 115–128. Liu, D., Pu, R., 2008. Downscaling thermal infrared radiance for subpixel land surface temperature retrieval. Sensors 8, 2695–2706. Ottlé, C., Kallel, A., Monteil, G., LeHgarat, S., Coudert, B., 2008. Subpixel temperature estimation from low resolution thermal infrared remote sensing. In: IEEE International Geoscience and Remote Sensing Symposium, pp. 403–406. Stathopoulou, M., Cartalis, C., 2009. Downscaling AVHRR land surface temperatures for improved surface urban heat island intensity estimation. Remote Sensing of Environment 113, 2592–2605. van der Meer, F., 2012. Remote-sensing image analysis and geostatistics. International Journal of Remote Sensing 33, 5644–5676.
108
W. Essa et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 95–108
Van de Voorde, T., van der Kwast, J., Uljee, I., Engelen, G., Canters, F., 2010. Improving the calibration of the MOLAND urban growth model with land-use information derived from a time-series of medium resolution remote sensing data. In: Taniar, D., Gervasi, O., Murgante, B., Pardede, E., Apduhan, B.O. (Eds.), Computational Science and Its Applications – ICCSA 2010. Springer, Berlin, pp. 89–104. Voogt, J.A., Oke, T.R., 2003. Thermal remote sensing of urban climates. Remote Sensing of Environment 86, 370–384.
Yuan, F., Bauer, M.E., 2007. Comparison of impervious surface area and normalized difference vegetation index as indicators of surface urban heat island effects in Landsat imagery. Remote Sensing of Environment 106, 375–386. Zhukov, B., Oertel, D., Lanzl, F., Reinhackel, G., 1999. Unmixing-based multisensor multiresolution image fusion. IEEE Transactions on Geoscience and Remote Sensing 37, 1212–1226.