Impact of Land Use Land Cover Change on Land Surface Temperature

Impact of Land Use Land Cover Change on Land Surface Temperature– A Case Study of Surat City, India

1. Study Area

Surat is a city in the Indian state of Gujarat, administrative capital of the Surat district. It used to be a large seaport and is now a center for diamond cutting and polishing, also known as the Diamond City of India.It is the eighth largest city and ninth largest urban agglomeration in India. The city is located at 21°10′N 72°50′E and has an average elevation of 13 meters. It has a tropical monsoon climate with temperature in summer ranging from 37˚C to 44˚C and the average annual rainfall ranges from 1,000 to 1,200 millimetres. Surat is a major industrial and commercial center and port, and in recent years it has experienced one of the most rapid economic growth rates among cities in India. (Bhatt, 2013)

As per a study conducted by City Mayors Foundation, Surat was the world’s 4th-fastest growing city in 2016 (City Mayors Statistics, n.d.). It has also been selected as one of twenty Indian cities to be developed as a smart city under Smart Cities Mission. The present study focuses on Surat inner city (Figure 1) covering approximately 309.12 km² area. This area was selected as a case study owing to significant LULC changes as the reason of industrialization and urbanization. (Tran, 2017).

                                                 Figure 1: Study Area

2. Data and Methods
   1.  Data

The study employed multi-temporal satellite data procured from Landsat 5 TM (1998 and 2008) and Landsat 8 OLI/TIRS (2016) for Surat inner City to investigate land use land cover change and their impact on land surface temperature. The details regarding the data used and elaborate methodology is represented in Table 1 and Figure 2 respectively. The details regarding the adopted methods are discussed in the following sub sections.

Table 1: Details of landsat data used in the study

                                             Figure 2: Methodology

2.2. Method for land use/land cover classification

An unsupervised classification was performed on the landsat imageries using the ISODATA clustering method to classify the images into the desired classes of which four different classes were effectively identified. A thematic raster layer was generated using the ISODATA algorithm while running ArcGIS. The pixels were identified for each of the categories and they were grouped into land cover classes such as water, vegetation, built-up and bare soil (oyekola, 2018; Dafalla, 2013).

2.2.1.      Method for accuracy assessment of land use classification

Accuracy assessment is an important step in the processing of remote sensing data. It determines how closely the result associates with the true values and renders the qualitative information from the procured satellite data. To evaluate the accuracy assessment of unsupervised image classification, an error matrix of referenced data was generated to attain information such as overall accuracy, producer’s accuracy, user’s accuracy, omission error, commission error and Kappa coefficient. For calculation ArcGIS software was used and a total of 885 sample sites were selected from google earth to match them with the LULC map for verification (Choudhary, 2018; Rwanga, 2017).

The overall accuracy of the classified image determines the percentage of matched number of sites to the total number of sites and can be calculated using the following equation. (Choudhary, 2018; Rwanga, 2017).

                                         $\mathit{Overall\; Accuracy}= \frac{\sum \mathit{Diagonal\; Value}}{N}$

                                         (1)

1111

Accuracy for individual land use land cover classes can be determined using two approaches i.e. producer’s accuracy and user’s accuracy. Dividing the matched number of sites of an individual class to that of total sites of the same class, multiplied by 100, yield user’s accuracy (Story and Congalton, 1986). The commission error of user’s accuracy measures likelihood of a classified pixel matching the land use land cover type of its corresponding real-world location (Congalton, 1991; Jensen, 1996; Campbell, 2007). By using the following equation’s user’s accuracy and commission error can be evaluated.

                         $\mathit{Use}{r}^{‘}\mathit{s\; Accuracy}\left(\%\right)= \frac{\sum \mathit{Diagonal\; value\; of\; Row}}{\mathit{Row\; Total}}\mathit{X }100$

                              (2)

                                    $\mathit{Commission\; Error}= \frac{\sum \mathit{offDiagonal\; element\; of\; Row}}{\mathit{Row\; Total}}\mathit{X }100$

                        (3)   

Whereas, the producer’s accuracy is computed by dividing the number of matched sites to the total number of sites developed from referenced data multiplied by hundred. It measures how effectively an area has been classified. The omission error of the producer’s accuracy determines percentage of observed sites on ground, not categorized in referenced map. By using the following equation’s producer’s accuracy and omission error can be evaluated (Choudhary, 2018; Rwanga, 2017).

                                $\mathit{Produce}{r}^{‘}\mathit{sAccuracy}\left(\%\right)= \frac{\sum \mathit{Diagonal\; value\; of\; Column}}{\mathit{Column\; Total}}\mathit{X }100$

                   (4)

                       $\mathit{Omission\; Error}= \frac{\sum \mathit{off\; diagonal\; element\; of\; Column}}{\mathit{Column\; Total}}\mathit{X }100$

                      (5)

Another method of measuring the accuracy is the Kappa coefficient (K) (Foody, 1992; Ma and Redmond, 1995) In this study, Kappa coefficient were also computed for land use land cover maps of all three years. It is a discrete multivariate technique measuring agreement or accuracy. The value of Kappa varies from 0 to 1. Where, 0 (zero) represents the worst and 1 (one) represents the best. (Jensen, 1996; Congalton, 1991).

The value of Kappa is expressed as percentage (%). Kappa coefficient (K) is derived by the following formula:

                                                     $K= \frac{N\sum _{i=1}^r{x}_{\mathit{ii }}–\sum _{i=1}^r\left(x_i+{\mathit{Xx}}_{+1}\right)}{N^{2 }–\sum _{i=1}^r\left(x_{\mathit{ii}}X{x}_{+1}\right)}$

                                         (6)

where; r is number of rows and columns in error matrix, N is total number of observations (pixels), Xii is observation in row i and column i, Xi+ is marginal total of row i, X+i is marginal total of column i. 

A Kappa coefficient equal to 1 means perfect agreement where as a value close to zero means that the agreement is no better than would be expected by chance (Rwanga, 2017).

2.3. Method for extraction of NDVI from Landsat data

NDVI is a most widely used index for monitoring the growth status and spatial density distribution of vegetation (Sun et al. 1998). It is utilized as an indicator of greenness and biomass of the earth’s surface. (Chen and Brutsaert 1998). Moreover, it is also believed to be a noble indicator of surface radiant temperature (Lo et al. 1997). Vegetation is extremely sensitive to absorption and reflection of the red and infrared bands and hence NDVI values for

satellite images of all three years (1998, 2008, 2016) were calculated as the ratio between measured reflectance in the red (R) and near-infrared (NIR) bands based using the following equation (Xiao and Weng 2007).

                                                    $\mathit{NDVI}= \frac{\left(\mathit{NIR}–\mathit{Red}\right)}{\left(\mathit{NIR}+\mathit{Red}\right)}$

                                                             (7)

Red and NIR represents spectral reflectance measurements acquired in the red and near-infrared regions, respectively. The values of NDVI ranges from -1 to 1. The higher the NDVI value, healthier and denser the vegetation in that specific area. The index for normal healthy vegetation falls between the range of 0.1– 0.75, while NDVI values for rock and soil are close to zero, and negative for water bodies. (Sahana, 2016; Tan, 2010; Babalola, 2016).

2.4. Method for extracting LST from thermal band of Landsat data

Increase of temperature is a major concern of the present world. To investigate the relationship between increasing temperature and escalating urbanization, Land Surface Temperature (LST) was computed (Choudhary, 2018). In this study, Spectral radiance method was used to retrieve LST from the Landsat 5 TM (band 6) for year 1998, 2008. For Landsat 8 TIRS, year 2016 out of two thermal bands i.e. band 10 and 11, single band (band 10) was used since band 11 assumes large calibration uncertainty (USGS, 2014; Yu, Guo, & Wu, 2014). Besides that, a combination of method and formulas was exercised to retrieve LST (Sahana, 2016).

LST from Landsat TM 5 Image was estimated by following two steps process. Spectral

radiance was extracted by using the following equation: (Sahana, 2016)

                                                        L = M_L * QCAL + AL                                                  (8)

Where L represents Spectral Radiance, M_L is Band specific multiplicative recalling factor, AL is Band specific additive rescaling factor and Q_cal is thermal band.

The next step is conversion of spectral radiance to surface temperature (T_B) using thermal constants provided in the metadata file. 

                                                            $T_B=\frac{K_2}{\ln \left(\frac{K_1}{L}+1\right)}–273.15$

                                               (9)                                             

Where K₁ and K₂ are band specific thermal conversion constant from the metadata, T_B is the surface temperature. For obtaining the result in Celsius, the temperature was revised by adding absolute zero (approx. -273.15) (Xu, 2004).

Furthermore, Land Surface Temperature (LST) from Landsat 8 TIRS (band 10) was derived using combination of methods and formulas in the algorithm. The values from above two equations (8 and 9) along with fraction of vegetation and land surface emissivity were taken into account for LST estimation. The surface emissivity, ε, was estimated using the NDVI thresholds method (Sobrino, Munoz, & Paolini, 2004; Sobrino, Raissouni, & Li, 2001). The fraction of vegetation, F_V, of each pixel was reckoned from the NDVI using the following equation (Carlson & Ripley, 1997):

                                                        $F_v={\frac{\mathit{NDVI}– {\mathit{NDVI}}_{\mathit{min}}}{{\mathit{NDVI}}_{\mathit{max}}–{\mathit{NDVI}}_{\mathit{min}}}}^2$

                                                   (10)

Where, NDVI_min stands for minimum NDVI value where pixels are considered as bare soil and NDVI_max stands for maximum NDVI value where pixels are considered as healthy vegetation (Guha, 2018).         

dε is the effect of the geometrical distribution of the natural surfaces and internal reflections. For heterogeneous and undulating surfaces, the value of dε may be 2%.

                                                            dε = (1- ε_s) (1 – F_v) Fε_v                                              (11)

where ε_v is vegetation emissivity, ε_s is soil emissivity, F_V is fractional vegetation and F is a shape factor whose mean is 0.55 (Sobrino et al., 2004).

                                                            ε= ε_v F_{v +} ε_s (1 – F_v) + dε                                                   (12)

where ε is emissivity. From Equations (11) and (12), ε may be determined by the following equation:

                                                              ε = 0.004 * F_v + 0.986                                             (13)

Finally, the Land Surface Temperature (LST) was derived using the following equation (Weng, Lu, & Schubring, 2004).                      

                                                              $\mathit{LST}= \frac{T_B}{1+\left(\frac{\mathit{\lambda \sigma }{T}_B}{\mathit{hc}}\right)\ln \epsilon }$

                                                  (14)                         

where h is Plank’s constant (6.626 × 10−34 Js), c is the velocity of light in a vacuum (2.998 × 10−8 m/sec), λ is the effective wavelength (10.9 mm for band 10 in Landsat 8 data), σ is Boltzmann constant (1.38 × 10−23 J/K), and ε is emissivity.

2.4.1.      LST Validation

One of the major LST validation model is through near-surface air temperature (Li, 2013; Srivastava, 2009) which was adapted in this study. The final retrieved LST results were validated using, the mean near-surface air temperature (Liu, 2011) and considering not just the mean temperature but also the actual temperature in the given pixel at the instant of the satellite passing over the area of meteorological station. The comparison was done with air temperature, which is disparate and can sometimes result in extensive differences since the resolution of LANDSAT 8 for the thermal band used for LST retrieval is 100m and 30m for the red and NIR bands. The LST was calculated and extracted for the pixel in which the meteorological station fell. Occasionally, the differences can be huge depending on the weather condition and other factors (Gallo, 2019).

Daily summary surface data was collected from the National Oceanic and Atmospheric Administration (NOAA) website for the meteorological station existing in the study area. The near surface air temperature data coinciding with the date of acquisition of satellite images for all the three years i.e. 1998, 2008 and 2016 were taken into account. The differences between the retrieved LSTs and the air temperatures and details on the stations are presented in Table 2. From Tables 2 it can be concluded that, the difference between the LST retrieved from the the satellite images and the near-air temperature varied from 4.06 (˚C) to 6.87 (˚C). 

                                  Table 2: Details and difference between LST and air temperature

2.5. Statistical Analysis

Statistical analysis such as zonal statistics, correlation and regression analysis were performed in this study. Zonal statistics using ArcGIS was exercised to summarize the values of NDVI and LST within each LULC class while correlation and regression was used to determine the relationship between LST and NDVI. 

 3. Results and Discussion

 3.1. LULC Change

In this study, the land use land cover classification was designated into four different classes viz. Water body, vegetation, built up area and bare soil for the year 1998, 2008 and 2016 (Figure 3). The accuracy assessment of classified images was analyzed and represented in the form of confusion matrix (Table 3). The overall accuracy (%) of the classified images for the year 1998 was 89.49, for the year 2008 was 91.96 and for 2016 was 94.20. The reliability of the result was derived with the help of Kappa coefficient and the values for year 1998, 2008 and 2016 was 0.85, 0.88 and 0.92 respectively (Table 3).

Besides that, total area of each land use land cover category and percentage of each class between 1998 and 2016 were calculated and are presented in Table 4. It is revealed from the table that over the last 18 years there was a dramatic increase in built up area and noticeable decrease in water body, vegetation and bare soil. One of the most conspicuous changes was noticed in built up area which has gone up from 19.10% to 26.41% at the rate of 38.28% in the first phase of study (1998-2008) and from 26.41 to 35.45% at the rate of 34.23% in the later phase of the study (2008– 2016). This change can be attributed to urban expansion particularly in the central part of the city. Similarly, water body showed decreased in area (0.36% during 1998-2008 and 0.21 during 2008-2016). While, vegetation and bare soil showed increase in one phase and decrease in other phase of study. In case of vegetation, first phase of study showed increase by (10.19% during 1998–2008) and decrease in second phase (12.92% during 2008–2016). Likewise, bare soil showed decrease in first phase (17.14% during 1998-2008) and increase in the second phase of study (4.08% during 2008-2016). 

   Figure 3: Spatial distribution of land use land cover classes for year 1998, 2008 and 2016

Table 3: Confusion matrix

                           Table 4: Land use land cover distribution over the study area

3.2. NDVI Change

Spatial variation of NDVI depends on various factors such as slope, topography, radiation availability and so on. (Liu, 2004). NDVI is a commonly used index to measure land surface greenness based on the inference that it is positively proportional to the volume of green vegetation in an image pixel area (Liu, 2004). The NDVI values of the pixel vary between -1 and 1. Extreme negative values represents water, values close to 0 indicates bare soil while values close to 1 indicates green vegetation. Comparing the NDVI maps of 1998, 2008 and 2016, changes were identified. The spatial distribution of NDVI over Surat inner City is shown in Figure 4. In 1998 the NDVI values varied between -0.37 and 0.66, which gradually reduced to ranging between -0.25 and 0.59 in 2008 and further reduced to the range of -0.14 to 0.56. Moreover, the region exhibiting high mean NDVI values in 1998 i.e. vegetation showed substantial reduction through 1998 to 2016 while the region with low mean NDVI values after water body i.e. built up area demonstrated gradual increase through 1998 to 2016 (Table 4) (Table 5).   

Table 5: NDVI values for different land use land cover categories

                      Figure 4: Spatial distribution of NDVI for year 1998, 2008 and 2016

3.3. LST Change

Satellite images for all the three years i.e. 1998, 2008 and 2016 were investigated for the characteristics of land surface temperature for each land use land cover class. This helped in determining the impact of LULC change on LST. The mean surface temperature of the study area showed increase at the rate of 2.42 C per decade. Table 7 shows that all the land use land cover classes identified recorded increase in the surface temperature over the timeline of the study. Spatio-temporal dispersion of surface temperature reveals that built up area has recorded the highest mean temperature 28.46 C in 1998, 32.46 C in 2008 and 32.62 in 2016. This implies that increase in urbanization does elevate the surface temperature by replacing vegetation with non-evaporating surface. Moreover, on analyzing the pattern of LST for all the three years it was observed that highest LST was recorded in the central part of the study area which is mostly built up. Overall, the study shows a positive correlation between built up area and LST.

The lowest average temperature was recorded in water body is 23.42 C in 1998, 27.23 C in 2008 and 28.23 C in 2016. (Figure 5 and Table 6). Vegetation divulged a low radiant temperature in all the three years since vegetation can reduce the amount of heat stored in the soil surface through transpiration (Omran, 2012). Surface temperature of water body was low compared to other classes, but the rate of increase was high because the data acquisition date corresponds to end of October and November and the radiant reflected from the water body is considered to be lower than other object during winter season. Rate of increase of surface temperature throughout 1998 to 2016 was highest over water body (4.81 C) followed by bare soil (4.33C), vegetation (4.18C) and built up area (4.16C) (Table 7)

                           Table 6: LST values over different land use land cover categories

                      Table 7: Change in LST values over different land use land cover categories

                        Figure 5: Spatial distribution of LST for year 1998, 2008 and 2016

3.4. LST and NDVI relationship

To assess the relationship between LST and NDVI, regression and correlation analysis were performed. Analysis based on linear regression showed that the coefficient of determination (r2) ranges from 0.965, 0.969 and 0.88 in the year 1998, 2008 and 2016 respectively (Figure 6). Likewise, the correlation coefficient (r) was observed to be -0.982, -0.984 and -0.938 for the year 1998, 2008 and 2016 respectively (Table 8).   

This observed relationship depicts that there is a strong negative relationship between LST and NDVI. From the analysis, it was observed that the vegetation covers had shown noticeably low mean temperature during the timeline of the study i.e. 1998 to 2016. This can be attributed to the fact that dense vegetation tends to minimize the portion of heat stored in the soil and surface structures by means of transpiration.

        Figure 6: Regression analysis between LST and NDVI for year 1998, 2008 and 2016

         Table 8: Correlation analysis between LST and NDVI for year 1998, 2008 and 2016

References

• Babalola, OS., Akinsanola, AA. (2016). Change Detection in Land Surface Temperature and Land Use Land Cover over Lagos Metropolis, Nigeria. J Remote Sensing & GIS, 5, 171.
• Bhat, G.K., Karanth, A., Dashora, L., Rajasekar, U. (2013). Addressing flooding in the city of Surat beyond its boundaries. Environment & Urbanization, 25(2), 429–441.
• Campbell, J.B. (2007) Introduction to Remote Sensing. 4th Edition, The Guilford Press, New York.
• Carlson, T.N., & Ripley, D.A. (1997). On the relation between NDVI, fractional vegetation cover, and leaf area index. Remote Sensing of Environment, 62, 241–252.
• Chen D., Brutsaert, W. (1998). Satellite-sensed distribution and spatial patterns of vegetation parameters over a tallgrass prairie. J Atmos Sci, 55, 1225–1238.
• City Mayors Statistics: World’s fastest growing urban areas. http://www.citymayors.com/statistics/urban_growth1.html. Retrieved 20 April 2019.  
• Choudhury, D., Das, K., Das, A. (2018). The Egyptian Journal of Remote Sensing and Space Sciences, https://doi.org/10.1016/j.ejrs.2018.05.004.
• Congalton, R.G. (1991). A Review of Assessing the Accuracy of Classifications of Remotely Sensed Data. Remote Sensing of Environment, 37, 35-46.

• Duy, X. Tran, Filiberto Pla, Pedro Latorre-Carmona, Soe W. Myint, Mario Caetano, Hoan V. Kieu. (2017). Characterizing the relationship between land use land cover change and land surface temperature. ISPRS Journal of Photogrammetry and Remote Sensing, 124, 119-132.
• Foody, G.M. (1992). On the compensation for chance agreement in image classification accuracy assessment. Photogramm. Eng. Remote Sens. 58, 1459– 1460.
• Gallo, K., Hale, R., Tarpley, D., Yu, Y. (2011). Evaluation of the relationship between air and land surface temperature under clear and cloudy-sky conditions. Journal of Applied Meteorology and Climatology, 50(3), 767–775.
• Guha, S., Govil, H., Dey, A., Gill, N. (2018). Analytical study of land surface temperature with NDVI and NDBI using Landsat 8 OLI and TIRS data in Florence and Naples city, Italy, European Journal of Remote Sensing, 51(1), 667-678.
• Jensen, J.R. (1996). Introductory Digital Image Processing: A Remote Sensing Perspective. 2nd Edition, Prentice Hall, Inc., Upper Saddle River, NJ
• Liu, J.G., Mason, P.J., Clerici, N., Chen, S., Davis, A., et al. (2004). Landslide hazard assessment in the Three Gorges area of the Yangtze river using ASTER imagery: Zigui-Badong. Geomorphology, 61, 171-187.
• Lo CP, Quattrochi D, Luvall J. (1997). Application of high resolution thermal infrared remote sensing and GIS to assess the urban heat island effect. Int J Remote Sens 18:287–304
• Ma, Z., Redmond, R.L. (1995). Tau coefficients for accuracy assessment of classification of remote sensing data. Photogramm. Eng. Remote Sens. 61, 435–439
• Omran, ESE. (2012). Detection of land-use and surface temperature change at different resolutions. J Geograph Inf Syst, 4, 189–203.
• Oyekola, M.A., Adewuy, G.K. (2018). Unsupervised Classification in Land Cover Types Using Remote Sensing and GIS Techniques. International Journal of Science and Engineering Investigations, 7(72), 11-18.
• Li, Z. -L., Tang, B. -H., Wu, H., et al. (2013). Satellite-derived land surface temperature: current status and perspectives,” Remote Sensing of Environment, 131, 14–37.
• Liu, L., Zhang, Y.Z. (2011). Urban heat island analysis using the landsat TM data and ASTER data: a case study in Hong Kong,” Remote Sensing, 3(7), 1535–1552.
• Rwanga, S.S. and Ndambuki, J.M. (2017). Accuracy Assessment of Land Use/Land Cover Classification Using Remote Sensing and GIS. International Journal of Geosciences, 8, 611-622.
• Sahana, M., Ahmed, R., Sajjad, H. (2016). Analyzing land surface temperature distribution in response to land use/land cover change using split window algorithm and spectral radiance model in Sundarban Biosphere Reserve, India. Model. Earth Syst. Environ., 2(81), 1-11.
• Srivastava, P. K., Majumdar, T.J., Bhattacharya, A.K. (2009). Surface temperature estimation in Singhbhum Shear Zone of India using Landsat-7 ETM+ thermal infrared data. Advances in Space Research, 43(10), 1563–1574.
• Story, M., Congalton, R.G. (1986). Accuracy assessment: a user’s perspective. Photogramm. Eng. Remote Sens. 52, 397–399.
• Sun, H., Wang, C., Niu, Z. (1998). Analysis of the vegetation cover change and the relationship between NDVI and environment factors by using NOAA time series data. J Remote Sens 2(3), 205–210.
• Sobrino, J.A., Munoz, J.C., & Paolini, L. (2004). Land surface temperature retrieval from Landsat TM5. Remote Sensing of the Environment, 9, 434–440.
• Sobrino, J.A., Raissouni, N., & Li, Z.L. (2001). A comparative study of land surface emissivity retrieval from NOAA data. Remote Sensing of the Environment, 75(2), 256–266.
• Sun, Q., Wu, Z., Tan, J. (2012). The relationship between land surface temperature and land use/land cover in Guangzhou, China. Environ Earth Sci, 65, 1687–1694. 
• Tan, K.C., Lim, H.S., MatJafri, M.Z., Abdullah, K. (2010). Landsat data to evaluate urban expansion and determine land use/land cover changes in Penang Island, Malaysia. Environ Earth Sci, 60, 1509–1521.
• USGS. (2014). USGS-earthexplorer. Retrieved from http://earthexplorer.usgs.gov/
• Weng, Q., Lu, D., Schubring, J. (2004). Estimation of land surface temperature-vegetation abundance relationship for urban heat island studies. Remote Sensing of Environment, 89(4), 467–483.
• Xiao H, Weng Q. (2007). The impact of land use and land cover changes on land surface temperature in a karst area of China. J Environ Manag, 85, 245–257.
• Xu, H.Q., Chen, B.Q. (2004). “Remote sensing of the urban heat island and its changes in Xiamen City of SE China,” Journal of Environmental Sciences, 16(2), 276–281,
• Yu, X., Guo, X., & Wu, Z. (2014). Land surface temperature retrieval from Landsat 8 TIRS comparison between radiative transfer equation based method, split-window algorithm and single channel method. Remote Sensing, 6(10), 9829–9852.