Urban climate and registration of thermal anomalies of territories using satellite imagery

This paper describes a simplified method of mapping of the ur-ban environment surface to obtain a map of the thermal anomalies distribu-tion and study the structure of the urban heat island. Those maps allow evaluating or planning the urban microclimate optimization methods and studying the effect of land cover type on the site temperature. The article discusses the processing of five satellite images for the summer and winter from 2002 to 2019. We propose a simpler and more automated processing of thermal images for Landsat 7 and Landsat 8. The stages of automatic atmospheric correction according to the DOS1 method and calculation of the emissivity with surface classification are considered. Image processing was carried out in the QGIS software package using the Semiautomatic Classification Plugin extension. As a result, thermal anomalies in Chelya-binsk were localized and a comparison of the thermal map for the specific region before and after urbanization was made.


Introduction
In recent decades, together with the rapid urbanization of territories, the problem of environmental comfort and safety of populated areas has become more acute. The so-called "urban heat island" (UHI) effect is becoming more pronounced in almost all major cities or industrial agglomerations, characterized by the appearance of zones of thermal anomalies with high urban environment surface temperatures. Thermal anomalies are most commonly concentrated in the most urbanized sectors of the area, and their occurrence mostly indicates cases of anthropogenic nature. Anthropogenic influence on the UHI effect can manifest in the presence of dense development and the land surface covering materials that actively absorb heat radiation, and also in reduction of green planting areas that lead to changes in surface radiative properties of the area (albedo and emissivity).
The influence of these factors leads to a change in the territory's energy balance, which leads to an increase in temperature in the city by more than 4 degrees Centigrade [1] and deterioration of the urban microclimate. Thus, the urban heat island effect is one of the important issues that arise when planning and developing urban areas, as it negatively affects the ecological and sanitary parameters of the environment. In addition, heat radiation from the surface can be a manifestation of phenomena and processes unseen in direct observation.
At the moment, the study of the phenomenon of urban heat island can be associated with more general tasks of studying and optimizing the urban microclimate. The microclimate may differ for different parts of the city, so it can be influenced locally, both positive-ly and negatively. For example, if you build a specialized planting or park area in a development zone, the average temperature in this area can significantly decrease.
Also, a significant influence on the temperature and environmental conditions of the area can produce changes in the aeration regime of this territory [2].
Thus, the study of the urban microclimate comes down to the compilation and processing of data on the spatial and quantitative distribution of thermal anomalies of the city for various sites and time intervals.
The city must conduct a complete mapping of the urban heat island over a wide time interval taking into account seasonal differences in order to comprise thermal anomalies in urban development planning and forecasting the dynamics of changes in the characteristics of the microclimate. We need such methods that would allow mapping UHI with adjacent territories in a sufficiently high spatial and temporal resolution to solve this problem.
Currently, there are satellite systems that allow carrying out a survey of the land surface in a very wide spectral range [3]. Special processing of these satellite images allows us to obtain sufficiently high-quality temperature maps of the earth's surface, so-called LSTmaps (land surface temperature maps); the images obtained from the newest satellites have such a high spatial and temporal resolution that they are quite suitable for research on a regional scale.
This paper describes simplified methods for determining the characteristics of heat radiation of the surface for Chelyabinsk from 2002 to 2019. We used multi-season thermal satellite images from ETM+ and TIRS sensors of the Landsat series using radiometric and atmospheric correction and surface type classification.
It seems clear that modern remote sensing data in the form of thermal images is a very valuable and rich source of information for the study of thermal anomalies in cities. Determining the nature and boundaries of thermal anomalies, as well as identifying patterns of spatial and temporal changes in the heat radiation of territories, will help to understand the causes of adverse environmental conditions in cities with a large number of industrial emissions, and to develop methods and recommendations for their minimization.

Materials and methods for studying thermal anomalies
Obtaining long-term and uniform data series on the structure and characteristics of the urban microclimate and thermal anomalies can be very useful or even necessary for the work of architects and designers of populated areas. A number of modern urban climatologists point out the importance of the issue of transferring knowledge about the urban climate into the daily practice of urban planning and design . It is also important to crystallize the obtained observations and methods into a set of modern comprehensive recommendations for urban planning and planning of urbanized territories.
An important aspect of urban climate optimization is an integrated approach to assessing several factors that affect the climate parameters of a territory with a wide spatial and temporal data coverage. To make the right design decisions for planning urban development on the early stages, it is necessary to assess the effect of the proposed measures to account for or optimize the existing urban microclimate as accurately as possible. For these purposes, researchers are increasingly resorting to various numerical methods, as well as mathematical modeling packages and finite element (or finite volume) analysis. However, for a more complete representation of the source data model, it is necessary to draw up a temperature map of the studied territory as a starting point for modeling and optimization.
There are several ways to obtain data on the earth's surface temperature; among those, there are two main directions: direct temperature measurements by sensors and remote measurements at different scales (aerial photography or satellite observations). Urban climate research currently uses observations from stationary weather stations more often [4], and much less frequently at high-altitude observation points, or using microwave profilers [5]. It is worth noting that the high requirements for quality, uniformity, and stability of surface temperature data favorably distinguish methods of remote sensing of the earth among other methods. They allow to get large amounts of data without resorting to manual measurements or expensive experiments. This is why, in recent decades, researchers have increasingly used remote sensing data to solve applied problems and to describe the surface heat island, and diagnose thermal anomalies [6,7]. This work is aimed at generalizing modern methods of processing satellite images and searching for the most optimal and simple method for studying the distribution of temperature anomalies in the central part of the city of Chelyabinsk. The obtained methods aim at researchers who are not professionals in the field of remote sensing data processing, such as specialists of urban planning profile.
The city of Chelyabinsk considered in this work is a large industrial megapolis with high rates of industrial development and urbanization, whereas the environmental situation in the city according to data for 2016-2019 is quite critical [8]. High-density development, the proximity of industrial enterprises to the center, flat terrain, and relatively low wind speeds lead to the accumulation of harmful impurities over the city center and the frequent appearance of smog. The situation may further worsen due to the stagnation of air masses in the vicinity of new neighborhoods under construction, due to their suboptimal orientation and location [9]. As the initial data for analyzing the thermal structure of the city of Chelyabinsk it was decided to use a number of images for the summer and winter period in several years, with an emphasis on control zones with sharp changes in the type of the earth's surface during the urbanization of the territory.

Methods for obtaining satellite images
The method for determining the temperature of the earth's surface from satellites is to obtain multispectral images, in different bands of which information about a certain radiation spectrum is stored. Information about the temperature of the earth's surface is obtained by capturing invisible to the eye radiation in the thermal infrared range by the satellite sensor. The resulting data is stored as a monochrome bitmap with high bit depth. Each pixel of this raster stores data about radiation received directly by the satellite's sensor, and therefore may contain incorrect values affected by cloud cover and other distortions. The theoretical instrumental accuracy of surface temperature estimation is about 0.5 °C, but the parameters of the atmosphere and aerosol clouds in the air can underestimate the values by several degrees.
In this regard, it is important not only to get the necessary image but also to process or decode it, to extract information, translate it into a visible image and simultaneously correct and compensate all the distortions and disturbances persisting in the image bands.
Currently, several publications describe in detail a large number of different methods for processing satellite images for research at the global and regional level. The oldest and most developed methods related to research at the global level. These include research on oceans, global climate change, and natural disasters or major forest fires. Studies of this scale usually use images of low or ultra-low spatial resolution (more than 1 km per pixel), but with very high temporal resolution (less than an hour between images of the same area). Examples of such systems are MODIS and AVHRR [10]. Despite all the advantages of these systems, the low spatial resolution of images makes them unsuitable for regional studies on a city or agglomeration scale.
Modern studies of regional scale surface temperature maps for various megacities, including Moscow [11], often refer to the ETM+/Landsat-7 system as an example of a system with high spatial and temporal resolution. Information about this satellite, as well as other current remote sensing systems, is shown in Table 1. From this table, it is obvious that images from the ETM+ sensor are the most suitable for studying the internal thermal structure of cities due to their high resolution (60m per pixel), but unfortunately, since the end of May 2003, due to the failure of the satellite longitudinal motion compensation unit, the ETM+ radiometer provides data with a defect in the form of missing data bands in the images. Although there are methods for circumventing these defects by summing up data with other images due to their superabundance, using this satellite for research after 2003 has not been practical in the context of this material.
Of the two remaining satellites, only two systems have a resolution suitable for our task: ASTER-Terra and Landsat 8. A number of studies [12] pointed out the advantage of the ASTER system in terms of the number of spectral capture bands and radiometric resolution, which can be useful for a more detailed classification of types of the earth's surface, as well as for conducting additional studies related to tracking vegetation parameters (  However, by the middle of 2020, it can be stated that data from this satellite for a fixed point is received with a low frequency (periodicity is over a month); therefore, no uniformity in time resolution is can be observed. In addition, mid-infrared data ceased to be transmitted in early 2008 due to overheating of sensors in this range [13]. Since this is the only thermal infrared range imaging system, it is now advisable to use data from the Landsat 8 satellite. However, it is important to note that in studies covering time periods up to 2013, it is acceptable to use data from the ASTER system, if available, as well as from the ETM+ sensor for the data up to 2003.
To observe seasonal changes in the distribution of an urban heat island, as well as to analyze changes in the structure of thermal anomalies over time, it is necessary to obtain images for the maximum possible number of time intervals for similar seasonal periods.
For this work, images from August 2002 and 2019 were used as a summer pair, as well as a number of images from the winter of 2001, 2013, and 2019. During the search for various sources of satellite images, numerous services were analyzed, which differ in the ability to preview and sort images by date and coordinates. The most well-known and convenient ones are shown in Table 2. It is important to note that it is necessary to look for images with the least cloud cover in the area of interest. It is also useful to familiarize yourself with the global notation system for Landsat data, i.e., the Worldwide Reference System (WRS). This system allows the user to request satellite images anywhere in the world by specifying the conditional center of the scene of interest with two parameters, namely, PATH and ROW. The first parameter defines the number of the vertical column of images or the orbit, and the second parameter indicates the image row or the number of the image from the North Pole to the South. This way, you can uniquely set the number of images included in the area of interest, taking into account the overlap of neighboring sections along with the patch and row, which will significantly speed up and simplify the search for the desired image. If the shooting periodicity is once every 16 days, the satellite can cross the zone of interest several times in one day through neighboring paths. This creates redundancy in the data, which can greatly help the researcher if there are defects in the images or if clouds cover the desired area.

The need to process and correct images
As a rule, images received from Landsat 7 and 8 systems are provided in the form of socalled digital sensor values (Digital Number, DN), a pre-adjusted form for the visible range with radiometric and geometric correction. This correction allowed us to take into account the curvature of the Earth's surface and the terrain features of the territory being photographed.
Additional corrections were usually required to obtain information from the thermal bands required to produce a thermal map of the surface of the so-called LST (land surface temperature) map [14].
Currently, images from on the Earth Explorer portal are distributed as sets of preprocessed images (collections) designated as Level 1 data and marked L1TP in the name of each snapshot. Not only these products have an atmospheric correction but can be independently processed to the L2 level based on metadata from images by using any suitable atmospheric correction algorithm.
Atmospheric image correction is the process of reducing the effect on the image of the atmosphere and converting the values of radiation reaching the satellite sensors (Top of Atmosphere (TOA) radiance) to the values of the actually reflected from the earth's surface spectral solar radiation.
There are several ways to correct the influence of the atmosphere on satellite images that do not require additional information about the state of the atmosphere at the time of the shooting. For example, the empirical line methods, as well as the Dark Object Subtraction (DOS) method [15]. The choice of a specific method depends on the type of application and is selected based on the content of the snapshot and the zoom level.
Collections of images of the second level of L2 processing, representing the albedo of the earth's surface with correction for the angle of the sun and the influence of the atmosphere, are also available on demand. However, these collections contain no thermal band and not suitable for our tasks.
To get a map of the earth's surface temperature, you need to upload an image from the L1TP collection for the date of interest with the lowest possible cloud cover; then, it is advisable to perform atmospheric correction of the image to more accurately plot the temperature. However, it is important to note that researchers are still discussing the need for the atmospheric correction of images, since there is still no reliable way to determine the amount of moisture in the atmosphere. In addition, the process of obtaining atmospheric indicators for the site and the time of shooting is quite time-consuming. Therefore, as long as there is a clear confirmation of no need for atmospheric correction for thermal bands, it is advisable to further consider the simplest automatic correction methods, provided that the most cloudless images are used.
An example of uploading source data for WRS Patch: 163 and Row: 21 from the Earth Explorer portal is shown in Figure 2. This image was obtained for the date of 02.08.2002 at 11:44 a.m. local time.

A generalized algorithm for determining the earth's surface temperature
After receiving the necessary images, their initial processing, and conducting atmospheric correction, you can proceed to direct decoding of the images to obtain the final LST map or map of the temperature of the earth's surface layer.
The initial data for determining the temperature are the values of the radiation intensity arriving at the satellite sensor and registered by the corresponding heat band. These values are denoted as DN-digital numbers. There are several approaches to further converting the original numerical values of image brightness to surface brightness temperatures in degrees Kelvin and Celsius. If we generalize the main stages of transformations [16], we can get the following algorithm for the transition from the DN sensor values to the LST values: 1. Getting the original image in the thermal range (Band 6 and 10 for ETM+ and TIRS) 2. Transition from DN to TOA Spectral Radiance according to the equation (1) 3. Transition from TOA radiance to TOA brightness temperature according to equation (2) 4. Defining the emissivity parameter for equation (3) 5. Calculation of the final LST map based on the received data, equation (3) The following equations are listed with explanations: Lλ=RADIANCE_MULT(TBand) * DN(TBand) + RADIANCE_ADD_(TBand). (1) where: Lλ = TOA Spectral Radiance (Watts/(m2*srad*µm)); TBand = thermal band of the image from the satellite of interest (6 for Landsat 7 and 10 for Landsat 8); RADIANCE_MULT, and RADIANCE_ADD = multiplicative and additive conversion coefficients.
Switching to TOA brightness temperature: where: TB=TOA brightness temperature(K); Lλ = TOA spectral radiance; K1 = TBandspecific thermal conversion constant from the metadata for thermal band; K2 = TBandspecific thermal conversion constant from the metadata for our thermal band The immediate transition to calculating the LST value using equation 3 is complicated by an additional factor: the need to calculate the emissivity parameter, i.e., the reflecting capacity of the earth's surface [17]. There are several different ways to account for this parameter since it is dynamic for each pixel of the image and significantly affects the final result. For example, this parameter can be used as a value based on NDVI [18] or classification of the earth's surface by macro classes, or be a constant in the case of a homogeneous area. This factor largely determines the existing set of approaches to the assessment of the LST in terms of using additional sources of information about the state of the earth's surface.
Existing studies [19] confirm a fairly high convergence of the obtained temperature data for the two main methods of calculating emissivity: the ones based on the NDVI parameter and surface classification. Using supervised semi-automatic surface classification methods can be a more flexible option for temperature studies on a larger scale since it allows to manually edit surface classes for known landmarks and then save this classification map for a whole series of images. The choice of method depends on the task type and size. In this study, the method was chosen based on the supervised semi-automatic classification due to its simplicity, clarity, and more precise control.
After reviewing the general aspects of image acquisition and processing, we can conclude that in order to obtain LST maps of the earth's surface from Landsat data, in general, it is necessary to take into account only two main factors in all the variety of calculation methods: the factor of atmospheric image correction and the factor of taking into account the emissivity of the earth's surface.
The choice of a particular method for each factor depends on the scale of the task, the survey area, as well as on the uniformity of the surface and it cannot be determined unambiguously. One of the tasks of this work is to find the easiest way to obtain an LST map with an emphasis on solving the problems of studying the urban heat island in conditions of dense development and high variability of the earth's surface parameters.
Anyway, L1TP images with geometric and radiometric correction are already available for us. All we need is to perform atmospheric correction and calculate the emissivity parameter using the simplest automated method.

The procedure of semi-automatic image processing to obtain a thermal map of the surface
The QGIS software was chosen as the main tool for working with raster maps for this task since it is a free and open-source tool that supports writing third-party extensions and scripts in python.
Since this article describes the method of simplified semi-automatic calculation, the QGIS Semi-Automatic Classification Plugin (SCP) [20] was chosen as a tool for simplified image processing and automation of a significant part of the image decoding process. This extension allows to perform all necessary corrections, as well as calculate the emissivity coefficient using supervised signature classification. Some sources [21] also mention the possibility of fully automatic calculation of the LST with atmospheric correction based on the calculated atmospheric profile with the use of the QGIS Land Surface Temperature plugin, but this plugin does not imply the use of signature classification and its support is currently over.
The SCP plugin is an extension for the QGIS software with its own panel-based graphical interface that has the main tools for processing images, as well as subsequent classification of the earth's surface. Working in the plugin should be described as a sequence of actions performed with a set of different spectral ranges (bands) of a single satellite shot. However, the images are supported not only from the Landsat system but also from AS-TER, which provides for convenient comparison of the images from different systems, as well as their validation. This plugin allows to preprocess snapshots just after downloading or after importing the bands into the program and creating a band set.
After receiving images in DN format, bands 2-7 are automatically converted to surface reflectance using DOS1 atmospheric correction method. For thermal bands, it converts the values to TOA brightness temperature through equation 1 and 2. In general, the existing set of tools and scripts built-in into this extension allows to automate the following processing steps: 1. Getting a snapshot at the specified coordinates using the Earth Explorer details 2. Atmospheric correction by DOS1 method 3. Transition to TOA brightness temperature for the thermal bands 4. Cropping the entire set of bands for the region of interest 5. Semi-automatic signature-based classification 6. Calculation of the resulting LST using equation (3) in a raster calculator. In order to test the applicability of this extension for the faster production of LST maps, the problem of mapping the distribution of thermal anomalies in the city of Chelyabinsk was considered. This example also involves the study of the distribution of the heat island effect for this city for the winter and summer periods in different years. Below are the results of processing different-season images of the city using the SCP extension for different satellites and seasons.
A couple of shots of the city for summer 2002 and 2019 was selected as the initial data for the experiment. These shots were used to compare the distribution of temperature anomalies and also to study the impact of urbanization on the urban heat island effect. The initial data for the summer pair of images are shown in Table 3.  All acquired images were processed in the QGIS software package using the SCP plugin for processing and calculating the emission coefficient using the supervised surface classification method. The first step in determining the emission factor is to determine the type of surface. The entire area is divided into four classes: water, open soil, vegetation, and buildings. Afterimage classification: each class is assigned an emission coefficient value corresponding to the surface type: Water -0.98; Built-up area -0.94, Vegetation -0.98, Bare soil -0.93.
The SCP plugin allows us to perform semi-automatic supervised signature classification of the surface or so-called classification with training. By getting the values of spectral characteristics in the selected image, we can classify image pixels with similar values to a particular class. Obtaining a raster map of surface classification can be a useful side effect since this map can be used in other applications, for example, for territorial zoning or estimating the distribution of different types of land use. As an aid to sample selection, the extension provides visualization of the RGB image composite in false colors using satellite bands in a modified order. This method allows highlighting areas of vegetation in red tones to facilitate classification.
A graphic representation of a color-based signature classification map using this extension is shown in Figure 3. The surface classification raster map obtained after this operation was used to determine the emissivity coefficient for the transition to the LST map.
The last stage of image processing is the calculation of the final LST value using equation (3) in a raster calculator with the use of the resulting emissivity map as an operand. The map was then converted to degrees Celsius by subtracting the 273.15 value from the resulting map.
These operations were performed for all images, after which they were assigned a color gradient; the histogram equalization procedure [22] was also performed to increase the image contrast. The final map rasters were combined with the image of the city map using the multiply overlay mode.

Analysis of summer images
To analyze the distribution of thermal anomalies, as well as to increase the visibility and simplicity of localization of anomalies, the resulting rasters were combined with the image of the city map in the multiply overlay mode. The resulting images were analyzed for the presence of striking differences in distribution in Figure 5 It is obvious that while there is a certain difference in absolute temperatures, the overall distribution is not too different. In the presence of common temperature maximums in the zones of industrial enterprises and urban landfills, one can note a larger size of heat spots on later photographs, as well as a larger number of warm zones. This fact can be explained by the appearance of buildings in areas previously covered with trees and shrubs.
For a more detailed analysis of the images, individual fragments with a larger scale were selected. The resulting rasters were passed through the histogram equalization procedure to enhance local contrast and saved in GeoTIFF format. Then these fragments were combined in the Google Maps software with satellite images of the city at the time of 2005 and 2019. After combining the images, an area of interest with a sharp change in land cover between 2002 and 2019 was identified in the zone of Shershnevsky Bor. By enlarging the fragment of interest and switching the image to the multiply overlay mode, we can see a clear change in the thermal structure of the site and the appearance of a bright thermal anomaly in the zone of new development (Fig. 6). Based on the obtained images, we can confidently state a strong change in the local temperature in the area of development relative to the period with a clear and green landscape. In this case, you can see that the color spots contour with the borders of the new neighborhood.
Without any doubt, these observations confirm the conclusion that the process of urbanization and compaction of buildings, accompanied by an increase in the number of cars and asphalt surfaces, inevitably leads not only to an increase in the overall average temperature in the city center, but also to the appearance of new temperature peaks and anomalies.

Analysis of winter images
Let's look at the results of observations for the winter series of images. These images were obtained by a similar method with the classification of the earth's surface with subsequent histogram equalization to increase contrast and make the comparison more visible. Results of calculations of LST and RGB color composite for December 12, 2001 are shown in Figure 7. An interesting observation was the presence of pronounced contrasting lines of the Miass River for the 2001 image. Given the relatively low average temperature in the image, we can conclude that the lack of contrast on reservoirs is caused by their freezing and being covered with a thick layer of snow, while the presence of a bright temperature anomaly after the dam indicates that water is being discharged, and flows of uncovered water have a positive temperature. In these images, you can clearly distinguish the change in the size of the city heat island's halo as the average temperature increases. It is obvious that with a decrease in air temperature, the thermal contrast of the city relative to suburban areas multiply increases. However, despite the approach of the main heat spot closer to the city center, the main thermal anomalies remain located in the areas of thermal power plants, industrial enterprises, and urban landfills.

Conclusions and discussion
After the images obtained from different satellites have been analyzed, we can draw conclusions about a sufficiently high spatial and temporal resolution of existing sensors of Landsat 7-8 satellites. It was confirmed that the data from the OLI sensor's thermal bands can be used on a city and district scale. Based on the obtained images, we can point out the expected prevalence of higher temperatures in the urban center zone, as well as higher temperatures of buildings and bare soil relative to areas with abundant landscaping.
The result of this article is a description of a method for estimating and mapping the surface temperature distribution of the urban environment based on satellite images of the Landsat 7-8 system in semi-automatic mode with the use of open-source software. The re-sulting method allows obtaining a large amount of data on seasonal and long-term changes in the temperature of the land surface in the outskirts of the town with the minimal amount of time and effort.
The obtained method allows us to acquire a large amount of data on seasonal and longterm changes in the temperature of the earth's surface in the vicinity of the city with minimal time and effort and to map the temperature anomalies of the city. The resulting method can also be further optimized by simplifying the calculation of the emissivity parameter using signature classification only for seasonal control points during one year. Application of this assumption will significantly reduce the amount of work and optimize the process of bulk processing of temperature maps.
According to the results of the study, it was noted that the temperature maxima on the thermal map coincide with the coordinates of industrial facilities in the city, as well as technological waste scrapheaps and dumps. Based on the obtained images of Chelyabinsk, a number of temperature anomalies were localized, and time intervals for the occurrence of new temperature maxima were determined, primarily caused by the process of urbanization and the displacement of plant cover.
The results and methods obtained in this work can be useful for solving applied problems of urban planning, as well as in dealing with problems of urban climatology and environment. Further development of this work may be the study of the applicability of the data obtained for localization of environmental violations and eco-monitoring of territories. After demonstrating the principle of obtaining data and their applicability for monitoring on a city scale, we can mention that these methods can be applied in the field of regional programs for a specific city. The use of Landsat observation data may be useful for several regional programs for the city of Chelyabinsk; among them, we can mark areas related to monitoring air and water pollution, as well as the program of Integrated environmental improvement of the area of the "Green City" settlements in the Chelyabinsk region.