1. Introduction
The wide coverage, cost efficiency, and time series of observations over the entire globe are the main advantage of remote sensing data, which allows the usage of land surface temperature (LST) as a variable to study rapid urbanization, destruction of vegetated areas, and climate change at either local or global scale [
1]. LST increased in urban areas of many countries of the world, due to land cover changes from natural environments into urban impervious surfaces (UIS), resulting in the creation of new surface urban heat islands (SUHI) [
2]. SUHIs induce thermal stress on human bodies especially in summer seasons, and increase the intensity and duration of heat waves [
3], which could result in rising death rates during heat waves in coastal areas [
4]. Hence, monitoring surface temperature in urban areas might help planners abate and mitigate potential thermal stress threats to people living and working in those areas.
Monitoring temporal LST and SUHI changes is challenging due to the seasonal variance in temperature, and even daily changes within the same season. Some studies revealed LST changes based on only a single image per period [
5,
6], raising doubts on their results reliability. While some other studies sacrificed the time series analysis by compose images of many years—sometimes 30 years divided into two periods—to do a reliable average of LST values [
7]. Other efforts have been made to fuse coarse spatial resolution and high temporal resolution data such as MODIS (Moderate Resolution Imaging Spectroradiometer), with moderate spatial resolution and low temporal resolution data such as Landsat and ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) data, to overcome cloud cover and long time span of moderate spatial resolution data problems. However, this solution in many cases is inapplicable as it is time intensive demanding large computing resources. While, taking the advantage of using MODIS high temporal resolution data alone might be useful to select appropriate images for LST processing in large study areas, but it will not be helpful in studying small areas of interest such as coastal reclaimed areas because of the data coarse spatial resolution.
Land reclamation is a well-known solution for land augmentation in coastal areas under pressure from population increases, rapid urbanization, and economic development [
8,
9]. Land reclamation, however, can result in unanticipated outcomes as it changes the thermodynamic condition of the coastal area’s surface, leading to the increase in LST, and results in the creation of SUHIs [
10]. As the reclaimed land is a dynamic environment, the study of spatial and temporal changes of SUHI on these areas faces additional difficulties due to continuous land cover change that alters land surface temperature rapidly over time. Hence, reclaimed lands in regular cloudy weather need an LST temporal comparison method, independent of seasonal surface temperature variations for SUHIs detection. Therefore, we have developed a new normalized method to minimize the temporal variance of LST in coastal areas—especially on coastal reclaimed areas—for LST monitoring and SUHI extraction.
Several researchers used normalization methods—especially the min–max feature scaling—to represent LST changes over time to support their findings [
5,
11]. Nevertheless, to the best of our knowledge, previous studies seldom used any normalization method for cross-season comparisons or SUHI extraction. We call our method the land surface temperature normalization method (LSTn), which is based on the average of water surface temperature for time-varying comparisons and SUHI extraction. We selected the average temperature of the water surface for LST normalization because water was the original land cover in the area of interest before the reclamation process. Therefore, we estimate the land surface temperature changes based on the average of the original land cover LST.
We selected reclaimed areas of Lingding Bay, Southern China, to apply the proposed method, because of the rapid reclamation processes and high urbanization rates occurring there since the late 1970s. Besides, the area of Lingding Bay has a dense cloud cover during the year with rare available moderate spatial resolution images suitable for LST processing. Previous studies estimated LST in the Lingding Bay coastal area as a part of the Pearl River Estuary (PRE) region, depending on single image per period due to the lack of data in this study area [
4,
12,
13,
14]. These studies revealed a positive correlation between LST and urban impervious surfaces (UIS), and negative relationship between LST and vegetation. Since these findings have been confirmed in the literature, we correlated the land cover type and UIS results from our proposed method to validate the applicability of LSTn for temporal monitoring and SUHI extraction in the study area.
The rest of this paper is as follows,
Section 2 describes the methodology and background of the proposed model.
Section 3 presents the results of the applied method in Lingding Bay and shows advantages and disadvantages of the LSTn method.
Section 4 discusses the results of the method in Lingding Bay and its correlations to the literature. Conclusions are drawn in
Section 5.
2. Materials and Methods
Figure 1 illustrates the workflow of the LSTn method. In the first step, we define the area of interest (AOI), by applying the normalized difference water index (NDWI) on the first available Landsat image of the target study area. Reclaimed areas are extracted by subtracting the NDWI of a recent image from the NDWI of the first Landsat image. A single free cloud image is selected for each studied year to classify land cover and extract the urban impervious surfaces (UIS) fraction using the multiple endmember spectral mixture analysis (MESMA) technique, based on the visible near-infrared and short wave infrared (VNIR/SWIR) wavelength bands.
In the second step, we estimate the land surface temperature (LST) for all available images after isolating cloud cover pixels from the scenes. Results for each season in a studied year are averaged to map seasonal LST. Another averaging estimates the annual LST from seasonal results. The third step consists of a normalization of seasonal and annual LST averages (LSTn) based on the average of open water surface temperature. The fourth and last step correlates LSTn to land cover and the UIS fraction validates the LSTn results for temporal monitoring of surface temperature and SUHIs extraction.
In this method, we did not use LST outputs from passive microwave sensors, because these products are still at a coarse resolution unlike those products extracted from thermal infrared remote sensing data [
15], which are inapplicable in the relatively small coastal reclaimed areas. We therefore, used thermal remote sensing data available at various spatial scales for assessing SUHI in reclaimed areas and other relatively small coastal environments. Details of the LSTn method application in the Lingding Bay are presented in the following subsections.
2.1. Study Area and Data Collection
The Lingding Bay is the main outlet of the Pearl River to the South China Sea. It is located in the south of Guangdong Province, an area of more than 1000 km
2, with a sub-tropical monsoon climate. The coastal area is under an intense pressure from land reclamation; many parts are considered artificial environments given the extent of recent economic development. Land reclamation in this area was a solution for the local governments to maintain the high population density of the Pearl River Estuary (PRE) region, exceeding 1200/km
2, related to the rapid evolution of industrialization in the Pearl River Delta since the late 1970s [
16].
Figure 2 presents the study area consisting of the coastal zone of Shenzhen, Dongguan, Guangzhou, Zhongshan, and the north of Zhuhai (hereafter, it will be referred to as Zhuhai) districts in Guangdong Province.
Figure 2 illustrates the area of interest (AOI) in light blue, mainland Chinese cities have red borders, Hong Kong have electron gold borders, and Macau within green borders. To define the reclaimed lands in Lingding Bay as the area of interest (AOI), we used the first available Landsat Multispectral Scanner System MSS image in the PRE, acquired on December 23, 1973, to extract the water body as AOI in 1973.
Table 1 illustrates all Landsat images used in this study with their respective seasons and data types. Among all images, we employed four cloud free Landsat images separated by about a decade, defined the areas reclaimed since 1973 (starred in
Table 1).
As shown in
Table 1, each period has data within a year buffer before and after each year period. Because, cloud free Landsat data are rare in this study area during a single year period. That decision was an effort to extract a very reliable LST average of all seasons during each period, rather than depending on a few sets of images during a single year. We combine the data of spring and autumn seasons together in one set for the same reason mentioned above. Winter season images are those acquired in December, January, and February. Images sensed in June, July, and August are corresponding to the summer season. While, spring/autumn images are those acquired in March, April, and May representing the Spring season, and September, October, and November images of Autumn season. Images from less than 60% cloud cover were collected for this case study. However, for the summer of 1987 period and winter of 1997 period there was only one image per season for the three corresponding years of the study period, as shown in
Table 1.
2.2. Reclaimed Areas Extraction
We used the normalized difference water index (NDWI) on the Landsat MSS image to extract the Lingding Bay water body in 1973. This method uses the fact that water shows high reflectance in the green light wavelength and its low reflectance in the near-infrared wavelength using the expression [
18]
where G is a band that reflects green light and NIR is the near-infrared band. In the next step, we converted the extracted water surface of Lingding Bay to a vector layer in a base map, and used the Landsat 8 image collected on 7 February 2016 to extract the reclaimed area extent since 1973. The extracted reclaimed lands were defined as the study area of interest for monitoring LST and SUHI in Lingding Bay.
The Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) module was applied to optical bands of the four Landsat images, starred in
Table 1, to correct atmospheric effects. Then, we used VIPERTools add-on in the ENVI 5.1 software environment to classify land cover in the study area based on multiple endmember spectral mixture analysis (MESMA) technique following the guideline of the add-on tool [
19]. MESMA has the advantage of decomposing each pixel using a combination of representational endmembers [
20], and it applies by running several models for each pixel to select the best fit one with minimum RMSE and has successfully been tested to map the majority of image pixels with two-endmember models only [
19].
VIPERTools is a series of steps started from collecting spectra from the image to end by mapping the land cover based on single or multi-spectral mixture analysis technique. In short, this process started by collecting different pixel spectra for five land cover types (water, shallow water, vegetation, bare soil, and impervious surfaces) to choose the purest pixels in each group as endmembers. VIPERTools has the advantage of using different approaches to select the optimal endmembers for each chosen feature (the Count-based Endmember Selection (CoB), the Endmember Average RMSE (EAR), and the Minimum Average Spectral Angle (MASA)). The optimum selection of endmembers had the highest in-CoB value (the total number of spectra modeled within the class) and the lowest EAR and MASA [
19]. Selected endmembers for each land cover classification were applied iteratively in a linear spectral unmixing analysis process using VIPERTools to define the best-fit spectra to map land cover types based on the multi-spectral analysis technique of MESMA. We processed the MESMA analysis two times; first, to map the five land cover types mentioned above. Second, we processed it with two-endmember models—one of UIS endmember spectra, and other of vegetation endmember spectra—to extract UIS composition inside each pixel representing the reclaimed areas. Spectra presented in
Figure 3 for UIS and vegetation were used to produce the two multi-endmembers classification image for each studied period.
Figure 3a, shows vegetation endmember spectra used in the second MESMA processing with high reflectance in the NIR wavelength, in comparison to other Landsat wavelengths; especially the red band. Variant spectra reflected by UIS endmembers represent different UIS materials inside urban areas, seen in
Figure 3b. Results of the second MESMA procedure include values of UIS and vegetation combined with shade effects inside each pixel. Hence, Equation (2) was used to calculate urban impervious surfaces composition inside each pixel by isolating the shade fraction results
where, U is the urban impervious surfaces fraction, V is the vegetation fraction, and C is the composition of the impervious surface inside each pixel. This study used land cover classification result in the first MESMA procedure to analyze the performance of LSTn in different types of land cover. Additionally, we used land cover classification and UIS composition to evaluate our results by measuring the correlation between the two results and LSTn, based on previous studies that noted significant correlations between UHIs, land cover type, and UIS composition [
4,
12,
13,
14].
2.3. Land Surface Temperature Calculation
To estimate the LST based on Landsat thermal infrared data, the digital numbers (DNs) in the thermal bands of Landsat 5 (Band 6) and Landsat 8 (bands 10 and 11) were converted to radiance data using the formula in (Equation (3)) [
21,
22]
where the Lλ is the radiance value, Qcalmin is the minimum quantized calibrated pixel value, and Qcalmax is the maximum quantized calibrated pixel value; Lmin is spectral radiance scales to Qcalmin, Lmax is spectral radiance scales to Qcalmax.
The spectral radiance of each image was converted to brightness temperature values as [
21,
22]
where Tk is the brightness temperature in Kelvin, Lλ is the spectral radiance value; K1 and K2 are the constants of Landsat calibrations. To convert the brightness temperature from Kelvin to Celsius, expressed as [
23]
where Tc is the brightness temperature in Celsius and Tk is the brightness temperature in Kelvin. Thirdly, LST (in Celsius) was calculated by using the formula in (Equation (6)) [
21,
22]
where Ts is the LST; Tc is the brightness temperature in Celsius; λ is the emitted radiance wavelength; P is the result of (h*c/b) in which h is the Planck’s constant, c is the velocity of light, and b is the Boltzmann constant; and Ɛ is the land surface emissivity (LSE).
The only missing value in this formula is the LSE. The estimation of LSE starts by calculating the NDVI for reclaimed areas using (Equation (7)) [
14,
24]
where NIR is the near infrared band, and RED is the red band in each image. Then, the NDVI results were used to calculate the proportion of vegetation by applying the expression in (Equation (8)) [
22,
24]
where NDVI
min is the NDVI minimum value, and NDVI
max is the NDVI maximum value. The LSE was estimated as [
22,
25]
After estimating the LST of all images, seasonal averages were calculated, followed by the overall average of the period LST to be ready for normalization procedure.
2.4. LST Normalization (LSTn)
After estimating LST average for seasonal and over all period data, we applied the LST normalization based on the average of water surface temperature. This normalization process is expressed in the formula
where
is the normalized LST (LSTn), X is the LST value, W
avr is the average value of bay’s water surface temperature, LST
max is the maximum value of LST, and LST
min is the LST minimum value within the study area, and the result ratio is varying from −1 to 1. We applied this equation to the estimated average seasonal data for each period, winter, spring/autumn, and summer. The LSTn of all seasons were average to extract LSTn results for each period.
4. Discussion
LST monitoring supports governmental human and environmental health initiatives and decision-making. SUHI affects the human body, soil moisture, and energy consumption in urban areas and surrounding peri-urban zones [
26]. A method for LST comparisons over time, specifically for rapidly changing coastal areas under regular cloudy conditions is urgently needed, despite the limited availability of thermal data due to cloud cover. LST variation over time is always problematic in monitoring applications over time [
27]. In this study, we propose a new normalized method to minimize seasonal changes in the LST for a flexible image selection to monitor LST and SUHI over time. LSTn overcomes seasonal and other temporal challenges by normalizing and scaling LST values estimated from Landsat images acquired at different dates, enabling temporal comparisons of surface temperature and extraction of SUHIs over coastal reclaimed areas.
LST varies widely over space and time; hence, it is difficult to evaluate the accuracy of the LST estimation obtained from Landsat images due to a lack of high-resolution thermal infrared data with identical times of image acquisition. Errors stemming from atmospheric corrections, sensor noise, false land surface emissivity estimation, aerosols and other gaseous absorbers, algorithm processing, sensor view angle effects, and wavelength uncertainty [
28,
29,
30,
31], could cascade errors and uncertainty in the LST retrievals. This study used two sets of Landsat sensor data; Landsat5 and Landsat8. On-board calibration of Landsat5 sensor had a consistent offset error of about 0.7 K over the lifetime of the instrument [
32]. The uncertainty of Landsat8 TIR laboratory calibration before launch was estimated by about <0.3% at all calibration temperatures for the two TIR bands, with maximum allowed radiance error W/(m
2 sr µm) of 0.059 (0.4 K at 300 K) for Band 10 and 0.049 (0.4 K at 300 K) for Band 11 [
33]. However, the stray light from far out-of-field has influenced the absolute calibration of Landsat8 thermal data since launch. Stray light refers to unwanted radiance from outside the field-of-view recorded by the sensor in the optical system of the scene and could reach rise to about 9 K [
34]. Since standard calibration techniques failed to correct for this error, efforts have been made to develop a method for stray light correction in Landsat8 TIR data. The first attempt to fix this error succeeded in reducing the stray light error by 0.29 W/m
2/sr/um (~2.1 K) in Band 10 and 0.51 W/m
2/sr/um (~4.4 K) in Band 11. The root-mean-square (RMS) variability was 0.12 W/m
2/sr/um (~0.8 K) for Band 10 and 0.2 W/m
2/sr/um (~1.75 K) for Band 11. [
34,
35]. In this correction effort, stray light reduced extreme errors from 9 K to only about 2 K [
34]. In early 2017, Gerace and Montanaro [
36], reduced the stray light error from 2.1 K at 300 K to 0.3 K for Band 10 and from 4.4 K to 0.19 K for Band 11, with less variability of 0.52 K at 300 K for Band 10 and 0.91 K at 300 K for Band 11. This final calibration was used in our Landsat8 images, as it was implemented in all Landsat8 products provided by USGS since early 2017 [
36].
Jimenez-Muoz and Sobrino [
28], concluded the uncertainty of atmospheric effects on thermal remote sensing LST retrieval between 0.2–0.7 K. While estimated the uncertainty of land surface emissivity between 0.2–0.4 K, and concluded a total error of thermal LST ranges between 0.3–0.8 K. Studies with in situ measurements have lower uncertainty values; since we do not have in situ available data, we considered the higher values of uncertainty in our results. Yu et al. [
29] compared three different approaches for LST inversion from 41 Landsat8 TIRS images: the radiative transfer equation-based method, the split-window algorithm and the single channel method. The results from this study revealed a small RMSE difference of 0.122 between all processed methods, with an advantage to the radiative transfer equation-based method by RMSE 0.903 K when applied on Landsat8 Band 10. The highest RMSE was 1.67 K for the single channel method when applied on Band 11 of Landsat8 data. Considering all error sources including noise error, bandpass effects, and wavelength indetermination, with no available in situ measurements, we took 0.9 °C as the LST retrieval accuracy of our analysis based on literature [
28]. In future work, we plan to measure LST uncertainty and influence on the proposed LSTn method results; and estimate the correlation between uncertainty values and normalized LST values in LSTn. In the current study however, we assume a correlation between LSTn results and land cover, as well as a correlation to the UIS composition as reported in the literature.
Previous studies revealed a positive correlation between land cover types, mainly UIS and LST [
6,
7,
8,
20,
21,
37]. Therefore, we considered the validation of our proposed method by correlating LSTn to land cover types, and UIS fraction. Land cover types were sorted based on ascending rank as we concluded from the literature, ranked on ordinal values for water (1), vegetation (2), soil (3), and urban impervious surface (4). The Pearson’s correlation coefficient results are revealed in
Table 4. All correlations were significant at the 0.01 level for annual and seasonal LSTn results in all studied years. Highest correlation coefficient value for the annual average was presented in 1997 for both land cover classification and UIS fraction. The highest correlation coefficient for seasonal LSTn was varied between winter seasons in 1997 and 2007 for both correlation tests and 1987 UIS fraction test, transitional seasons for both tests in 1997, and land cover classification correlation of 1987.
To the best of our knowledge, no previous method was conducted to eliminate or reduce the seasonal influence of estimated LST for temporal comparison determination. Our LSTn method yielded a decrease in seasonal variations of surface temperature allowing for improved temporal comparisons independently of seasonal effects. In comparison to the regular LST variations method to detect SUHI, and its correlation to buildings and impervious surfaces [
38], the proposed LSTn method revealed sharper edges of surface temperatures between UIS and other land cover types that make SUHIs be more distinguishable in the area of study.
As mentioned, normalization methods were used previously for surface temperature comparison. Sirous et al. [
11] calculated a min-max feature scaling normalized land surface temperature (NLST) to compare their SUHI intensity results between different images that used a different normalization method, but did not consider a specific land cover type for normalization estimation. Amiri et al. [
5] used the same normalization method to compare the surface temperature overtime in transitional land cover areas, mainly from vegetation and bare soil to urban impervious surfaces.
In contrast, our proposed method might consider as a new LST intensity method since it depends on a specific land cover type out of the urban area, as previously estimated in literature based on forest, agriculture areas, and water surfaces [
27]. However, this method overcome the uncertainties of urban and surround reference land cover definition as it depends on the average of open water surface temperature as a reference. Nevertheless, this method has some limitations, as it decreases the seasonal influence on the temporal comparison and time series analysis, not removing it permanently, and results could be influenced significantly in case of extreme temperature events such as heat waves in the selected cloud free images.
5. Conclusions
For a long time, the study of temporal LST changes using thermal infrared remote sensing was a matter of identifying similarities during the same season over multiple years. Cloud cover in optical remote sensing added uncertainty and limited the availability of data usage in many areas worldwide. In this paper, we propose a new method for surface temperature monitoring on newly reclaimed areas over time, which is independent from seasonal effects for more flexible optical remote sensing dataset selection over time by normalizing LST data based on the average open water temperature. We used Landsat data at 10-year intervals from 1987 to 2017 to compare LST normalized values (LSTn) over reclaimed areas of Lingding Bay, Southern China.
The LSTn method revealed a relative decrease in seasonal surface temperature variability compared to the original LST results. Hence, the proposed method provides a more robust means to monitor surface temperature over time. LSTn trends had a relatively lower variability average (0.01) to the overall images annual mean in comparison to the LST variability average (3.9 °C). Therefore, detecting of SUHI is easier when using LSTn because it increases urban impervious surface temperature values relative to other land cover types.
Experiments are required to evaluate the original LST uncertainty influence on LSTn results. Caution is needed when selecting the cloud free images for surface temperature monitoring in cases of extreme temperature events such as strong heat waves in the summer. Additionally, more tests are needed to examine this method for winter surface temperature on soil and water covered zones, as the winter season reveals relatively higher difference to the overall mean images in some studied years in comparison to other studied seasons.
This method can be applied for coastal urban areas in general because of its dependency on the average of open water temperature for calculating the normalization average. Future work will focus on the influence of estimated LST uncertainty on the normalized results in LSTn, and the applicability of this method in surface temperature comparisons of coastal cities in different climatic zones, based on the nearby open water temperature average.