LAND SURFACE TEMPERATURE AS AN INDICATOR OF URBAN HEAT ISLAND EFFECT: A GOOGLE EARTH ENGINE BASED WEB-APP

At least 2 billion urban occupants will be concentrated in Asia and Africa, amounting to 70% of the global population by 2050. This rapid urbanization has caused an innate effect on the ecology and environment, which further results in intense temperature variations in urban and rural areas, especially in India. According to a recent IPCC report, 8 out of the 15 hottest cities in the world are situated in India. The rising industrial work, construction activities, type of material used for construction, and other factors have reduced thermal cooling and created temperature imbalance, thereby creating a vicious effect called ”urban heat island” (UHI) or ”surface urban heat island” (SUHI). Several researchers have also related it with climate change due to their contribution to the greenhouse effect and global warming. In this study, we have particularly emphasized northern India, including Punjab, Rajasthan, Haryana, and Delhi. We created a Google Earth Engine (GEE) based Web-App to assess the UHI intensity over the past 15 years (2003 – 2018). We are using Moderate Resolution Imaging Spectroradiometer (MODIS) images, Landsat 5, 7, and 8 data for studying UHI. The land surface temperature (LST) based UHI intensity (day and night time) will be available for major metropolitan cities with their respective clusters. With feasibility in SUHI monitoring, we can address an increasing need for resilient, sustainable, and safe urban planning of our cities as portrayed under the Sustainable Development Goals (SDG 11 highlighted by United Nations).


INTRODUCTION
Climate change, increase in urban-rural population, rise in population density, and poverty are few out of a huge lot leading to challenges gripping the world. The rise in temperature both air flux temperature and land surface temperature have been contributing factors in climate changes [Mustafa et al., 2020]. Several studies have shown that better social, economic and livelihood prospects lead to unprecedented urban migration, which may subsequently lead to temperature rise and cities behave as urban heat islands (UHI) [Li et al., 2013, Sheng et al., 2017.
India currently stands second just shy of China to become the most populous nation by 2027. Massive urban migration has become visible for a better life and work prospects in both formal and informal sectors, which may have caused changes in land-use patterns and fluctuations in regional temperature. A direct relationship can be established between the two [Chen and Zhang, 2017]. Research shows the positive effects of the presence of greencover, agriculture, and water-bodies helping to regulate the UHI effect, and simultaneously the negative consequences of built-up area and uncovered landcover accelerate UHI at an alarming rate [Amiri et al., 2009, Chen and Zhang, 2017, Song et al., 2014.
Remote sensing is very beneficial for UHI and LST studies. Multispectral images obtained from Landsat series sensors i.e. TM (thematic mapper), ETM+ (enhanced thematic mapper), and TIRS (thermal infraRed sensor) are * Corresponding author beneficial for estimating LST and determining the appropriateness of 30m spatial resolution for UHI studies [Guha et al., 2018]. The Landsat series is available at a global scale at 30m spatial resolution but with poorer temporal resolution at 16 days. For carrying out UHI and LST studies, this research uses thermal infrared bands of Landsat-5 and Landsat-8, available at 30 m spatial resolutions [Sobrino et al., 2004, Jiménez-Muñoz et al., 2014. Landsat series with help of new emerging methods such as monowindow algorithm [Qin et al., 2001], single-channel algorithm Sobrino, 2003, Jiménez-Muñoz andSobrino, 2009] have made UHI and LST estimations easy.
Global surface UHI explorer app 1 , developed by Yale University helps to monitor UHI intensities at a global scale using Google Earth Engine (GEE). The UHI dataset was based on a simplified urban-extent algorithm (SUE) using MODIS Terra and Aqua dataset at 1 km spatial resolution [Chakraborty and Lee, 2019]. The app shows the annual daytime and nighttime temperatures of 10,000 big cities around the globe and contains a year-wise slider from 2003 until 2018. It also generates two charts for the selected cluster, where the first one displays the change in the UHI intensity from 2003 to 2018, the second chart shows the monthly variation of the UHI intensity derived from 15 years of MODIS data.
In this paper we have developed a GEE-based Webapp for determining LST and analyzing the UHI effect at 30 m spatial resolution using Landsat 5 and Landsat 8 images.
2. STUDY AREA AND DATA UHI and LST estimation was carried out for Punjab state in the northern part of India. Punjab is bordered by the mountainous region and union territory of Jammu and Kashmir (J & K) to the north, Chandigarh city to the east, Himachal Pradesh to the northeast, Haryana to the southeast, and Rajasthan to the southwest. The state adds an area of 50,362 sq. kilometers, accounting for 1.53% of India's total geographical area. The state lies between 29.54°N to 32.57°N latitudes and 73.88°E to 76.94°E longitudes.
In this study, we have used thermal infrared (TIR) bands available at 120m and 100m from Landsat 5 and Landsat-8 satellites respectively. Emerging technologies and algorithms were used or LST estimations. We have also used MODIS (or Moderate Resolution Imaging Spectroradiometer) on-board Terra and Aqua satellites with global coverage of Earth at a temporal resolution of 1 day. For the development and validation of a global, interactive Earth system, MODIS is playing a vital role thus able to predict global change accurately enough to assist policymakers in making sound decisions concerning the protection of our environment. Even at coarser resolutions such as 250m, 500m, and 1km due to global and day-night time coverage, MODIS datasets are extremely fruitful in UHI and LST estimations (https://modis.gsfc.nasa.gov/ about/).

METHODOLOGY
Google Earth Engine (GEE) has petabytes of freely available datasets from the United States Geological Survey (USGS). This research uses various datasets available in GEE such as Landsat-5, Landsat-8, MODIS etc. For computation of LST from Landsat-8, either of the two bands 10 or 11 can be used in split-window method; however band 11 has some quality issues and large calibration constraint. Hence, band 10 is preferable for thermal studies in determining brightness temperature. The digital number (DN) image can be converted to top of atmosphere (TOA) reflectance, which can be further used to compute surface reflectance (SR). The SR is used to determine normalized difference vegetation index (NDVI) using NIR and RED band of the images. SR reflectance data is readily available from GEE datasets derived from Land Surface Reflectance Code (LaSRC) where the coastal band is used for aerosol inversion, and MODIS data for auxiliary atmospheric and atmospheric correction performed with radiative transfer model [Vermote and Kotchenova, 2008]. Quality assessment band (BQA), can be used to retrieve cloud coverage masking, which includes cloud shadowing for which a function was created in GEE [Ermida et al., 2020]. The NDVI equation is mentioned as:

Surface Emission
Using the NDVI thresholds, min, max values, surface emissivity (ε) was estimated using the equations given as the NDVI thresholds method [Sobrino et al., 2004, Sobrino et al., 2001. Also, the fractional vegetation (Fv), of each pixel was determined from the NDVI threshold information using the following equation [Carlson and Ripley, 1997]: Using the F v parameter, (ε) can be determined using equation [Guha et al., 2018]: The (ε) is derived from Landsat bands as follows [Ermida et al., 2020]: • F v is derived from ASTER using the NDVI equation.
• Using original ASTER emissivity, bare ground emissivity and corresponding F v with a prescribed value of ε b,veg =0.99 • ASTER emissivity for bare ground can be used for each TIR band 10,11 of Landsat-8 with spectral fixing.
• Corresponding NDVI values for F v determination.
• Fractional vegetative cover is fundamentally used for calculating emissivity in each TIR band combination.

CHIRPS Daily: Atmospheric Data
Climate Hazards Group InfraRed Precipitation Station (CHIRPS) gridded data is a 30+ years dataset representing quasiglobal rainfall/precipitation. It incorporates 0.05°arc resolution satellite images along with in-situ/field-based station data to create tile/grid-based rainfall time series for trend analysis and seasonal drought forecasting. The CHIRPS dataset is available at a global scale from 1981 -until now. Precipitation band is the only band available with its raster values as a function of mm/day.

Landuse Landcover
Landuse classification for the given study area can be carried out mainly into 4 classes i.e. Agriculture, Forest, Urban-city, and Waterbody. There are two methods to carry out classification i.e. unsupervised (without the presence of class labels) and supervised (with the presence of class labels). General steps to carry out a supervised based classification are (developers.google.com/earth-engine/guides/ classification): • Store training data, store class labels, identify specific properties of each class label as a numeric value to define part of predictors.
• Combine all these labels with the specific property as a feature collection.
• Identify or create a classifier, hyper-tune its parameters.
• Train the classifier with set parameters with defined training data.
• Classify image or feature collection with defined labels.
• Validate with test dataset, create error-matrix.
• Compute the overall accuracy and classified map.
• Iterate the process after tuning the hyper-parameters until desired overall accuracy isn't obtained.
With property storing class labels and properties storing predictor variables, training data is considered a feature collection. Class labels need to be continuous, starting from 0 and stored as integers. When it is necessary to convert class values to continuous integers, use remap() function. The prediction labels have to be numeric values.
Validation, testing, and/or training dataset can be made from multiple sources. Preparation of training dataset can be intensely carried out in GEE, with the help of geometry drawing tool. Alternatively, a predefined training set can be imported from a GEE as an asset. A classifier can be defined from one of the constructor tabs as an ee.Classifier and can be trained using the classifier.train() option. The current research uses a statistical machine intelligence and learning engine (SMILE) Classification and Regression Trees (CART) classifier [Breiman et al., 1984] to easily predict more than two classes.

WORKFLOW
The entire workflow can be segregated into 3 phases Figure 2: Workflow starting with input where we have Landsat-8 time series data, layer stacked. Cloud-cover mast at less than 10% cover. In the Data processing section, here I am showing single snippets using Landsat 8 imagery. I am determining the vegetation cover with normalized difference vegetation index. Using the thermal band of Landsat 8 image collection and brightness temperature determine the land surface temperature and emissivity factor. For land cover classification I am using the smile cart supervised classification algorithm. Using the climate hazards group infrared precipitation station data abbreviated to CHIRPS which gives a global daily precipitation data at 0.05 arc resolution I determine the precipitation. Also, MODIS TERRA and Aqua dataset is used for surface heat Island estimation. In the output section, I am showing the Indian outline map and my highlighted study area Punjab state. in the bottom corner, shown is the LST time series chart variation for the year 2016.
MODIS Terra and Aqua datasets are used to witness and analyze the annual summer, winter day-time, and nighttime variations across the years from 2003 until 2018. Overlaying city-wide clusters for the study area to identify The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIV-M-3-2021 ASPRS 2021 Annual Conference, 29 March-2 April 2021, virtual the UHI and LST variations. Annual max and min for seasonal variation across the day-time and nighttime are also added as additional drop-down options. The advantages are the availability of Terra as a daytime sensor and Aqua as a nighttime sensor for assessing thermal temperature variations. The temporal resolution of 1 day is also an additional advantage. The drawbacks, however, include coarser resolution of 1km and not taking into account the other factors for example precipitation variations, landcover changes, vegetation cover which have a tremendous impact in controlling the UHI effects [Chakraborty and Lee, 2019]. These factors were covered in the proposed research using Landsat-8 OLI/TIRS (LS), CHIRPS (precipitation) datasets. LS series have been beneficial to cover the UHI and LST estimates from 2013 -until now, at 30m fine resolution with one drawback being the poor temporal resolution of 16 days. CHIRPS daily global precipitation since 1981 until now provides data at fine resolution for annual precipitation in (mm/hr) was incorporated for the study. Landuse landcover was done for the entire study area using the SMILE CART algorithm. GEE user interface platform (https://code.earthengine.google.com/) is used for the implementation of all these factors. Finally, a GEE-based web app is also developed to monitor the LST estimation which soon will be available to use for the general public.

RESULTS AND DISCUSSION
Now for the results section, in the figure-4 top left cor- Figure 3: LST, Land Cover Classified Image, Precipitation (mm/hr) ner, we have the precipitation map using CHIRPS data set measured in mm/day. The red portion signifies Higher precipitation cover followed by mild precipitation in the green region and lower precipitation in the blue region. In the top right corner, the land cover classification was carried out using a smile cart supervised algorithm, there are four classes water body agriculture City, and forest. the classification results overlap and match with the ground truth labels used. the water body is signified with blue color, agriculture is signified with the red region as Punjab state is a hub for agriculture yield. The city is one with yellow color and forest with a dark green cover. In the bottom corner, I have the LST map estimated from the brightness temperature equation, the blue color predominates in the Northern and North-Eastern regions of Punjab signifying the lower temperatures. The urban settlement has red and yellowish red color along with the South Western region signifying the higher temperature.
In figure-5 of the result top left corner signifies images for Landsat 8 image collection visualized as true color composite (TCC). The top left corner signifies the vegetation cover determined using normalized difference with index dark green color showcases forest cover, light green color is predominantly agriculture and Orange color is for non-vegetation cover. The bottom left corner indicates the thermal band of Landsat-8, the light blue portion has a lower temperature compared to White and light green portions which indicate high temperature. Emissivity shown in the bottom right corner determined from the LST brightness temperature equation indicates with higher temperatures in the red and yellow region and lower temperature is indicated by light green and darker green color which is dominated by forest cover and Agriculture. In figure-6 here, annual precipitation as a function of mm/day for the year 2020 is shown as a time series for this study area month of March till April we have average precipitation during months of May and June due to Lu condition and absence of tropical winds, the temperature is much higher thus precipitation is lower. the month of mid-July till early October there is a sudden spike in precipitation due to the onset of southwest monsoon winds. December to January slight rainfall is occurring due to the presence of North-Western disturbance.
In figure-7 here, the global urban heat island explorer app is created by the center of earth observation, Yale University (https://yceo. yale.edu/research/global-surface-uhi-explorer). In, yellow is a cluster of 10,000 cities available globally. To find the surface UHI intensity it has been implemented using MODIS Terra Aqua data set from the year 2003 to the year 2018. The options include annual day and night time, summer day and night time, winter day and night time.
Finally, figure-10 shows the interface on GEE Explorer App that showcases the algorithm to create land surface temperature estimation with Landsat-8 OLI/TIRS sensors. This tool maps land surface temperature in the Punjab region from 2013-onwards using a brightness temperature (BT) parameter from Landsat-8, Vegetation cover (NDVI), Landcover classification using Smile Cart Classification derived from Landsat imagery. The landcover classification had an overall accuracy of 90.17%. One can use the tools below to explore changes in LST estimation, vegetation cover, land cover changes, and precipitation changes. The Dropdown menu has been made available to the user to filter the satellite images for different dates. Selection of images as part of image collection can be car-The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIV-M-3-2021 ASPRS 2021 Annual Conference, 29 March-2 April 2021, virtual In figure-8 here, LST variation for Chandigarh city is shown for the state of Punjab, the annual temperature variation can be seen for the year 2018 where light glow signifies the annual day temperature varying less than -1.5 degree Celsius while light brown region indicates annual day temperature varying from 0 to 1.5 degree Celsius.

CONCLUSION
In this paper, we have shown the UHI and LST estimation for the state of Punjab, India using both Landsat-8 OLI/TIRS sensor datasets and MODIS Terra and Aqua datasets. We can conclude that MODIS Terra and Aqua datasets are crucial for LST and UHI determination due to mainly two reasons. Firstly, availability of 1km global daily data and secondly, availability of both daytime and nighttime mean, max and min temperature variations for 10,000 cluster cities mainly can be considered the relevant reasons to use MODIS data. Using Landsat-8 OLI/TIRS sensor-based imagery is beneficial mainly for its fine spatial resolution at 30m and optimizing with other factors (NDVI, F v , (ε) and Thermal information (Band 10,11)). Punjab is an agricultural hub, thus the NDVI values are more than 0 throughout the year ranging between (0.1-1.0). The annual precipitation derived using the CHIRPS dataset from 1981-on wards, shows that Punjab is predominately affected by SW and NW trade winds. The LST, thermal information also shed the same concept that the rural regions in SE Punjab (Abhor) along Rajasthan state have shown a large heatwave condition due to the presence of loo during summers, lack of vegetation cover, and water scarcity. Emissivity variations for SE Punjab are  higher than in other areas. Precipitation (mm/hr) annually has receded throughout the years. These, all factors are reasoning's behind the high-temperature surge. Incorporating drought parameters can help establish the rising UHI and Surface Heat Island (SHI) effects.