Study of the Urban Heat Island (UHI) Using Remote Sensing Data/Techniques: A Systematic Review

: Urban Heat Islands (UHI) consist of the occurrence of higher temperatures in urbanized areas when compared to rural areas. During the warmer seasons, this effect can lead to thermal discomfort, higher energy consumption, and aggravated pollution effects. The application of Re ‐ mote Sensing (RS) data/techniques using thermal sensors onboard satellites, drones, or aircraft, al ‐ low for the estimation of Land Surface Temperature (LST). This article presents a systematic review of publications in Scopus and Web of Science (WOS) on UHI analysis using RS data/techniques and LST, from 2000 to 2020. The selection of articles considered keywords, title, abstract, and when deemed necessary, the full text. The process was conducted by two independent researchers and 579 articles, published in English, were selected. Qualitative and quantitative analyses were per ‐ formed. Cfa climate areas are the most represented, as the Northern Hemisphere concentrates the most studied areas, especially in Asia (69.94%); Landsat products were the most applied to estimates LST (68.39%) and LULC (55.96%); ArcGIS (30.74%) was most used software for data treatment, and correlation (38.69%) was the most applied statistic technique. There is an increasing number of pub ‐ lications, especially from 2016, and the transversality of UHI studies corroborates the relevance of this topic. processes, whose soils can present higher thermal signatures when revegetated [132]. in


Introduction
Urbanization is an anthropic alteration that generates modifications in surface materials due to vegetation suppression, albedo variation, and soil sealing, influencing the local energy balance, contributing to the formation of the Urban Heat Island (UHI) [1]. This effect is a consequence of the greater absorption of electromagnetic energy and the slow cooling of urbanized surfaces compared to surrounding areas with the presence of vegetation [2,3].
This local dynamic contributes to an increase in surface temperature, with a subsequent reduction in relative humidity and latent heat, and an intensification of sensible heat. The main causes of UHI formation are: (i) the ability of building materials to store heat; (ii) anthropogenic heat production; (iii) alteration and minimization of wind speed as a function of surface roughness; and (iv) increased absorption of solar radiation from lower albedo surfaces, among others [4][5][6][7].
Due to the relevance of the topic, studies focusing on understanding the causes and effects of UHI are increasing, including the analysis of the variability in air temperature in urbanized areas and its underlying mechanisms [8][9][10][11]. Voogt and Oke (2003) [12] suggested a sub-classification analogous to UHI, called Surface Urban Heat Island (SUHI) which, in addition to studying surface temperature differences between urban and rural  [3]).
In addition, there is an influence regarding the thermal properties of the constituent materials and their arrangement, which can generate shading effects and corroborate to the local heat exchange characteristics: the more heterogeneous their composition, the more complex the local heating and cooling relationship (as in the case of civil constructions). In urban areas, SUHI gains intensity around noon, and during the night, urban areas retain more heat due to the thermal behavior of buildings, surface materials, and areas that have reduced view of the skydome/hemisphere, as a function of lower Sky View Factor (SVF) (Figure 1) [19][20][21][22].
There are several techniques to study the thermal behavior of a site, such as Remote Sensing (RS) [23], data from fixed meteorological stations, data from qualified entities, in situ campaigns with portable thermal cameras, etc. Some studies adopt more than one of the aforementioned data sources to complement the information, both at the atmospheric and surface-level (UHI and SUHI, respectively). This variety of sources allows the access to diurnal and/or nocturnal data and in different seasons of the year, although the effect is more intense in summer and winter [19][20][21][22].
From the RS data, it is possible to compute the Land Surface Temperature (LST), a relevant variable that can be used to determine the radiative load of the earth's surface. It conducts longwave radiation and turbulent heat fluxes between the earth and the atmosphere, i.e., it is a variable in the physics of local and global surface processes, linked to radiative, latent, and sensible heat fluxes at the interface with the surface, allowing analyses of warming trends at the earth's surface [24][25][26][27].
The LST, also known as radiometric or "skin" temperature, refers to the direct measurement of the earth's surface temperature. Unlike the measurements taken by meteorological stations that record the temperature near the surface, the LST allows for a more detailed scale of analysis: in areas of dense vegetation it will represent the temperature of the leaves of the canopy; in areas of sparse vegetation the temperature will correspond to the whole canopy, subsurface (limbs, branches, etc.), and the ground surface; and on the bare ground it will correspond to the temperature of the top (few micrometers) from the ground surface [25].
In RS, the aggregated radiometric temperature is represented per pixel, according to the sensor display field. The LST is recovered by estimating the emitted surface radiance, obtained by the atmospheric correction in the radiance sensor, and inventing the Planck function, considering the effects of emissivity variation [25].

Satellites More Used to Estimate Land Surface Temperature
Considering the use of RS data for studying UHI and SUHI, the most commonly used sensors in RS concerning the works evaluated in this review were thermal sensors, which detect emitted and/or reflected terrestrial radiation [2,24]. Thermal sensors present challenges in recording the collection, especially when clouds are present, which can affect the reliability of the results and the temporal analysis of a site [36,37].
Another technique is the adoption of passive microwave sensors [36,38] which, although they do not present as detailed spatial scale as thermal sensors, overcome the challenge of atmospheric influence and allow for a thermal analysis of the atmospheric boundary layer, as in the case of Moscow, where profiles up to 600 m high were created for the study of UHI [38].
It is emphasized that there are studies that use data from both collection methodologies, to complement the information in the study area [36]. Details regarding the operation of each methodology will be presented in sequence (Sections 1.1.1 and 1.1.2).

Thermal Sensors
Thermal sensors operate in the 8-14 μm range, known as one atmospheric window. This range includes the emission part of the earth's spectrum, making heat detection by RS sensors possible. The emitted radiation is recorded as Digital Number (DN) and subsequently converted into temperature data [2,26].
Each sensor has its spectral resolution that defines its ability and sensitivity to generate data within the electromagnetic spectrum, which includes both visible and non-visible zones. In practice, the sensors are cooled to temperatures close to absolute zero so that their eventual emissions do not influence the measurements and the temperature data records of the targets. A comparison is made between the measured radiation temperature (target) and the internal reference temperatures of the sensors (related to the absolute radiation temperature) [2]. The main advantage of RS sensors is their continuous spatial coverage and temporal repeatability of the study area [39].
The most applied satellite/sensors to study LST are Landsat, Moderate Resolution Imaging Spectroradiometer (MODIS), and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER). Some studies used data obtained by polar and/or geostationary satellites, which prioritize the temporal mapping of the study area [40]. Table 1 presents a summary of the thermal sensors most applied in UHI studies. It is emphasized that this list is not intended to end or encourage that only these satellites be used to obtain information about the LST (the OSCAR-Observing Systems Capability Analysis and Review Tool, offers a query base for users, who can filter the choices of satellites, sensors, or products) [41]. In the next subsections, a brief description will be given considering these sensors. 

Landsat mission
The Landsat mission began in 1972, but only in the following decade, with Landsat 4, thermal data began to be collected, due to the recording capabilities of the following sensors: in Landsat 4 and 5, the Thematic Mapper (TM) [42,43]; in Landsat 7, the Enhanced Thematic Mapper Plus (ETM+) [44]; and in Landsat 8, the Thermal Infrared Sensor (TIRS) and Operational Land Imager (OLI) [45]. The Landsat family operates in heliosynchronous and near-polar orbit with a nominal altitude of 705 km, an orbital inclination of 98.2°, and a revision time of 16 days [42][43][44][45].
The Thematic Mapper (TM), present on Landsat 4 (operation from July 1982 to June 2001) and Landsat 5 (operation from March 1984 to January 2013), is a multispectral sensor, with seven spectral bands that detect information, simultaneously. Thermal infrared radiation is detected by band 6, whose records occur only at night. The Instantaneous Field of View (IFOV) of this sensor is 30 × 30 m in bands 1-5 and 7, and 120 × 120 m in band 6 [42,43]. ETM+ of Landsat 7 began operation on 15 April 1999, and since June 2003 has acquired and provided data with gaps caused by the Scan Line Corrector (SLC) failure [44]. It consists of a fixed scanning, multispectral radiometer with eight spectral bands, including a panchromatic and a thermal band. The spectral resolution of the bands is in the range of 0.45 to 12.5 μm and it has a spatial resolution of 30 m (except for band 6, which is 60 m, resampled to 30 m, and the panchromatic band, which is 15 m). The radiometric resolution, which is the sensor's ability to discriminate small variations in energy, which in images refers to the number of distinguishable shades of gray between black and white, is related to the sensor's bits, and is 8 bits and records Digital Number (DN) values between 0 and 255 gray levels [44].
Concerning Landsat 8, the OLI consists of one panchromatic band and eight multispectral bands, with a resolution of 15 and 30 m, respectively. In the TIRS sensor, there are two thermal bands with a resolution of 100 m, whose data are resampled to 30 m [45]. Landsat 9 is planned to be launched at the end of 2021. There is a similarity regarding the operation and characteristics of the Landsat 8 thermal sensors, renamed to TIRS 2 and OLI 2 [46].
Landsat has a long and uninterrupted observation program [47] and its heliosynchronous orbit favors temporal studies of the same place, since its revisit occurs at the same point and UTM time, with an interval of 16 days. Its data are released free of charge by the United States Geological Survey (USGS), in levels 1 and 2 (with absence and presence of pretreatment, respectively).
To compute the LST, it is necessary to process the two levels of data: NIR and RED bands from level 1 and thermal band from level 2 [48]. The standard processing will be presented below [27].
The conversion from Digital Number (DN) to Top Of Atmosphere (TOA) spectral radiance (Lλ)-physical quantity defined as the radiant flux in each direction, considering a surface normalized concerning surface area and unit solid angle [26]. It is calculated according to Equation (1) [49,50], using the thermal band. (1) where Lλ is the TOA spectral radiance (Watts/(m2 * srad * μm)); band-specific multiplicative rescaling factor from the metadata (radiance_mult_band_x, where x is the band number) is the ML; the band-specific additive rescaling factor from the metadata (radi-ance_add_band_x, where x is the band number) is represented by AL; and the quantized and calibrated standard product pixel values (DN) is Qcal, i.e., the thermal band.
where TK is the Top Of Atmosphere brightness temperature (K); Lλ is TOA spectral radiance (Watts/(m 2 * srad * μm)); K1 and K2 are band-specific thermal conversion constants from the metadata (K1 or K2 _constant_band_x, where x is the thermal band number).

273.15
(3) where TC is Top Of Atmosphere brightness temperature (°C). The thermal band includes the Emissivity (E) of the soil and vegetation. Although it is possible to deduce radiance and temperature values, it is necessary to include emissivity in the computation for the generation of the LST, using the NIR and RED bands. For that purpose, first, the Normalized Difference Vegetation Index (NDVI) is calculated (according to Equation (4)). NDVI was first proposed by Rouse et al. (1974) [52].
where NIR is the near-infrared satellite band and RED is the red satellite band. This result is used to obtain the Vegetation Proportion (PV) (Equation (5) [49,53].
where NDVImin and NDVImax are the minimum and maximum values obtained in the NDVI calculation, respectively. In Equation (6), the emissivity (E) is estimated [54].
where LST is the temperature, with correction by emissivity (°C); TC is the temperature of the brightness at the sensor (°C); λ is the wavelength of the emitted radiance; E is the emissivity; P is which is deduced from Equation (8) [23,49]. The standard process is presented in Figure 2.  [55]. With a revisit time of every one or two days, its spatial resolution is low: bands 1 and 2 are recorded with 250 m, bands 3 to 7 with 500 m, and bands 8 to 36, with 1 km. MODIS sensors can measure the physical and biological properties of the oceans and land and the physical properties of the atmosphere [55].
Thermal data has a resolution of 1 km, which makes its products applicable to larger regions. There are several pre-processed products for users, such as ocean surface temperature, ice, snow, evaporation, precipitation [56], LST with daytime and nighttime data, and emissivity, from MOD11C3, MOD11A1, and MOD11A2 products, respectively.
In MODIS products, the LST value is retrieved by applying the generalized splitwindow algorithm [57]. In MODIS LST Collection-6 products, the MODIS Land Surface Temperature and Emissivity (LST&E) product (MOD21), the LST is calculated by the AS-TER Temperature Emissivity Separation (TES) algorithm and are more sensitive to land cover changes when compared to other emissivity products [55,58].
TES applies physics concepts to dynamically and simultaneously retrieve both LST and spectral emissivity for the MODIS thermal infrared bands (29, 31, and 32) through simulations of radiative transfers for atmospheric correction, called Water Vapor Scaling (WVS), and a model of the variation of surface radiation data, resulting in accuracies at the 1 K level for various surfaces (deserts, water, and vegetation) [59].
Studies using products from sensors onboard satellites may be useful to show differentiation in the daytime and nighttime LST levels when there is not much cloud concentration. For this purpose, a study used a sensor onboard the Aqua satellite to analyze heat wave episodes in Cyprus [60]. The combined use of Terra and Aqua MODIS data can also provide valuable results, especially when complemented with daytime data, as it increases the spatial, temporal, and angular coverage of clear-sky data [61].  ASTER Advanced Spaceborne Thermal Emission Reflection Radiometer (ASTER) is a sensor onboard NASA's Terra satellite, launched in December 1999. It operates with three independent imaging subsystems and collects data in different regions of the electromagnetic spectrum: in the Visible Spectrum (VIS) and NIR region there are three spectral bands with 15 m spatial resolution (bands 1 to 3); in the shortwave infrared there are six bands with 30 m resolution (bands 4 to 9), and in the thermal infrared region there are five bands with 90 m resolution (bands 10 to 14) [62][63][64][65].
Additionally, there is the VIS-NIR system, consisting of two telescopes. One of them operates with rear-view, with a difference of seconds relative to the nadir view, making it possible to generate stereoscopic pairs in band 3. The revisit time is 16 days or less since the subsystems have a side view of ±24 degrees beyond nadir. Each image covers a 60 × 60 km area [62][63][64][65].
The data are available in different treatment levels. Level 1A, unsuitable for most users, refers to raw data, without calibration or noise correction pretreatment, and with non-co-recorded bands. Level 1B, a product suitable for many applications, has geometric and radiometric calibration and co-registered bands. Level 2 includes more specific products, identified by acronyms, the most used of which are AST04-surface temperature; AST05-surface emissivity; AST06-decorrelation-enhanced color compositions; AST07-surface reflectance; AST08-kinetic temperature (in which it is possible to obtain the LST, with an accuracy of 1.5 K and there may be nighttime records of the study area, which is interesting for UHI studies); and AST09-surface radiance. Level 3 products are more sophisticated and include, for instance, digital elevation models (AST14). Level 2 and 3 products are freely available on the Internet, although they were charged from 1999 to April 2016 (resulting in fewer search results in this period) [58,[62][63][64][65][66]. 

Polar and geostationary satellites
Polar and geostationary satellites can also be used for LST estimation. Their fixed positioning at a point on earth provides wide temporal coverage of local data, even allowing for a more adequate characterization of the diurnal cycle compared to other orbits [40]. However, the spatial resolution and coverage area do not present much detail, which can compromise the identification of land surfaces in a more refined way [37,67], and are not suitable for analyses of UHI patterns and dynamics about heat-associated health risk at urban scales [67].
NPP consists of a mission to monitor biological productivity and long-term climate trends. It acts between the measurements gauged by EOS Terra and Aqua, complementing the data between NASA's EOS missions and the Joint Polar Satellite System (JPSS) [70].
MGS operates in geostationary orbit, 36,000 km above the equator and can be directed for data collection at different locations, and is currently operated over Europe, Africa, and the Indian Ocean. Its sensor, SEVIRI, has 12 spectral channels and provides data on the atmosphere, assisting in the generation of predictive weather models. It has eight channels in the thermal infrared region, including High-Resolution Visible (HRV), which performs sampling with a nadir distance of 1 km, with daytime cycle coverage every five minutes, due to its fast scanning [68,71].
The GOES mission is controlled by the National Aeronautics and Space Administration (NASA) and operated by the National Oceanic and Atmospheric Administration (NOAA). Its orbit is geosynchronous equatorial, 35,800 km above the earth, and scans the continental United States and its surroundings (Atlantic and Pacific oceans), Central and South America, and southern Canada [67]. The GOES LST product is based on the splitwindow technique, which corrects atmospheric absorption and surface emissivity. It has a spatial resolution between 4 and 8 km [72].
Considering that most of the materials analyzed in this review do not use data from polar and geostationary sensors, the treatment processes for estimating the LST will not be detailed.

Passive Microwave Sensors
LST can also be computed from processing data from passive microwave radiometers. Passive sensors are microwave instruments that receive and measure natural emissions, which are produced by the earth's atmospheric and surface constituents. The power associated with the measurement is related to the composition and roughness of the surface being measured, and the physical temperature. The frequency bands are determined by fixed physical properties of what is being measured and cannot be duplicated in other bands [73]. They generate products capable of measuring biogeophysical variables, such as air temperature [36,74], even at night [36]. Compared to thermal sensors, passive microwave sensors have a greater ability to penetrate clouds, snow, rain, arid surfaces, and vegetation, complementing the data for UHI studies [75].
The limitations of these data relate mainly to the coarseness of spatial resolution (low) and the inability to recover temperatures on frozen surfaces [36].
The first attempts to obtain data from passive sensors were performed in the 1990s, i.e., compared to thermal sensors, there are limitations of historical and coverage areas [76][77][78].
The Advanced Microwave Scanning Radiometer-Earth observing system (AMSR-E) on the Aqua satellite was responsible for collecting data at high temporal resolution and contributed to the mapping of global near-surface air temperatures on a near-daily basis and for over a decade, until there was a failure of its antenna in October 2011 [36]. A combination of data from AMSR-E and AMSR2 (onboard the Japan Aerospace eXploration Agency's Global Change Observation Mission-Water1 (JAXA's GCOM-W1)) was released in May 2012 [36,79]. The data gap was filled with data obtained by the Microwave Radiation Imager (MWRI), onboard on the Chinese satellite FengYun-3B (FY3B) [36]. An example application is the use of data from these sensors to obtain information about the threedimensional structure of urban temperature. The Atmospheric Boundary Layer (ABL) acts as a buffer zone and accumulates heat, pollutants, and moisture. The ABL influences the intensity of vertical exchange, and the more unstable it is, the more favorable the removal of pollutants in the lower layer [38].
Moreover, the information from passive sensors can be combined with thermal ones for the complementarity information in locations with a higher presence of clouds or obstacles that influence the collected results. It is noteworthy that the methodology employing the use of passive sensors was not recurrent in the analyzed papers.

SUHI and LST Studies Applying Different Techniques
Information on the surface temperature of a location can be useful in different fields of investigation, either to assess climate change of natural and/or anthropogenic origin or to predict climatic behavior based on statistical models and time series. LST prediction can be an important tool for mitigating possible social impacts and promoting improvements in the quality of life of the populations.
Although LST is widely applied for studies that aim to analyze SUHI, this parameter is cross-cutting and can corroborate for the analysis of other topics, for example: environmental analyses with research on aquifers/water bodies [80], biomes [2], forest dynamics [81], and natural hazards such as volcanic activity studies [82]; air pollution [83] and determination of particulate concentration in urbanized areas [84]; architectural assessment to identify the building area and building layout with the least impact to the local microclimate [85] and identification of Local Climate Zones (LCZ) [86], from the identification of regions with uniform characteristics, from hundreds of meters to several kilometers, with respect to land cover, surface structure, building materials and human activities [10], use of hybrid data, including physical and mobile equipment for climate characterization of a site [87], socioeconomic and cultural analysis, including studies on biophysical and socioeconomic impacts [88][89][90], and social events [91]; and public health and mortality associated with HWs [92].

Study Objectives
Given the relevance of the studies in this area, this work aimed to identify the publications from 2000 to 2020 that address the UHI effect using RS data, namely the LST, and to perform a qualitative and quantitative analysis considering the following parameters: (i) number of publications in the considered period; (ii) geographic location and climate zone of the study areas; (iii) keywords; (iv) general objectives and application areas; (v) equipment/sensor used to obtain information about the thermal data; (vi) information sources used to obtain Land Use Land Cover (LULC) and additional data; (vii) software employed to process the data; (viii) Vegetation Index (VI) and additional data used for the comparative analyses; (ix) statistical methods employed to evaluate the relationship between indices and LST; (x) authors' main conclusions (quantitative and qualitative analysis); (xi) bibliographical references; (xii) authors; and (xiii) citations.

Bibliographic Database
The databases selected for this study were Scopus and Web of Science (WOS). Scopus has been available since 2004. Its publications are peer-reviewed, and it is one of the largest databases of literature abstracts and citations worldwide, with more than 25,100 active titles from over 5000 international publishers, divided into the fields of science, technology, medicine, social sciences, arts, and humanities. In 2020, it included 23,452 peer-reviewed journals. Its tools for information retrieval and aggregation are efficient, allowing data extraction in several formats, which facilitates the analytical process [93]. WOS is a global database of various publishers. It has consistent indexing in its main base, which is approximately 50 years old. In the WOS platform version, there are more than 34,586 publications, including journals, books, proceedings, patents, and datasets. Like Scopus, peer review occurs independently, and it is possible to filter by area or publication type and to issue customized reports in different extensions, which allows the analysis of the publications available in the database [94].

Search Strategy and Validity
Considering the relevance of reviews to provide a broad picture on a given subject from the contextualization of a topic, its methodological approaches, the direction of research, and the challenges and points for improvement discussed [95], the selection process of the evaluated papers is important.
The following steps were performed: firstly, the papers were screened to select those applying RS in UHI studies; secondly, the data from the selected papers were systematized, identified, and filtered, considering the titles, abstracts, methodologies, and other parts that could correspond to the 13 criteria analyzed (described in Section 1.3); subsequently, the resulting data were identified and filtered, descriptive statistics were applied, and analytical graphs were plotted; finally, the results were discussed, by class.
Not all steps contained in a meta-analysis were adopted in this review, since it was not intended to assign weight or relevance to any of the criteria evaluated in the papers, such as the number of citations or year of publication [96]. The intention was to select and evaluate all papers quantitatively and qualitatively, from a manual review, using synthesis study techniques similar to those adopted in other reviews in the area of RS [97].
The details of each step are presented below. The search in both databases considered publications that occurred from 2000 to 2020. The keywords ((UHI*) OR (satellite image*)) existing in the article title, abstract, or keywords fields were searched. Search results were saved in a list so that the base was static and Excel software was used for data tabulation and manipulation.
In the first stage, only the publications classified as "papers" were considered, to maintain in the database as many unpublished documents as possible, since some research may be first presented at conferences and later submitted as journal papers. Thus, 222 (17.9%) publications classified as "Conference/Proceedings Paper" and 39 (3.2%) as "Reviews" were discarded. The second step to refine the database consisted of the identification of duplicated documents based on the criteria of similarity of the Digital Object Identifier (DOI), title, and abstract, resulting in the elimination of 136 papers (13.91%). For each repeated title, one was kept, and the other was deleted.
Next, the title, abstract, and keywords were evaluated to check if the UHI studies used RS data (satellite, drone, or aerial imagery), to estimate LST. In cases where the fields did not spell out which data sources were used, the researchers performed a full reading of the paper methodology. At this stage, 43 papers (5.10%) were excluded, resulting in 799 publications (64.49% of the initial selection).
To refine the database and allow a more effective qualitative analysis of the remaining papers, a third analysis was performed, considering the papers published in English and without the "in press" status. Thus, of the 799 approved publications, 756 (94.62%) were in English; 36 in Chinese (4.5%); 3 in Spanish (0.38%), and the others correspond to 1 publication in other languages: German, French, Croatian, and Portuguese (totaling 0.50%). Thus, 43 papers (5.38%) were discarded. Additionally, 14 papers with "in press" status were also discarded for the next selection step.
Thus, 742 papers were selected for the next stage of analysis, which consisted of a thorough reading regarding the adherence of their objectives and methods regarding this review. In this step, 163 papers were discarded because they did not use RS data as the main sources of information or because they did not present in a clean form the main objective (UHI effect).
As each researcher performed this process separately, the titles approved by both were included in the final database. Finally, 579 papers (46.73% of the initial total), were selected. All these papers applied remote sensing techniques for UHI and SUHI analysis. Figure 3 presents a summary of the steps and the respective pass/fail criteria at different stages. For the presentation of the results, in addition to graphs generated in Excel, the VOSviewer tool [98] and ArcGIS version 10.7 software [99] were used to present the data in grouping and map format, respectively.

Results
In this item, the most relevant results identified during the development of this work will be presented.

Publication Period
The evolution in the number of publications concerning the study of UHI with RS data/techniques in the last 20 years can be observed in Figure 4. From 2000 to 2008, this methodological combination resulted in fewer publications when compared to the following years. In 2001 and 2002, the database search did not show any publications. In this case, a second search was performed, using only the term "UHI" as a criterion, to assess whether there had been any inconsistency. As a result, 26 papers published in this period were identified (adding the two databases), but none met the selection criteria of this paper (either using data only from meteorological stations or applying the acronym UHI to other phenomena, such as in health areas). It is possible to verify a significant increase of publications in the last 20 years, especially in the last five years, the period which concentrated 399 (68.91%) of the analyzed papers. This increase reinforces the relevance of the topic in the scientific community, as well as the applicability of RS to analyze the UHI effect.

Geographic Location and Climate Zone of the Study Areas
Regarding the study areas, the most studied continents or subcontinents were Asia (405-69.94%), Europe (65-11.23%), and North America (59-10.19%), totaling 529 (91.36%) papers. Besides the papers with continental delimitation, 14 papers were produced with data from at least two continents, representing 2.42% of the analyzed papers. Table 2 shows the number of publications and their respective percentages. Considering the countries, China was the most studied, with 233 (40.24%) papers, followed by India with 69 (11.92%), the United States (50-8.64%), and Iran (19-3.28%). These four countries are accountable for 64.08% (371) of the material analyzed. Figure 5 shows the spatial distribution of the papers by country. Studies that were conducted on more than one continent (3-0.52%) were not considered in this map in order not to generate duplication of information. Regarding climatic zones, in this work, we adopted the Köppen-Geiger global classification, proposed in 1900 by climatologist Wladimir Köppen and improved in 1918Köppen and improved in , 1927Köppen and improved in , and 1936 with the help of Rudolf Geiger to separate the papers analyzed. According to the concepts of phytosociology and ecology, the regional natural vegetation correlates with the predominant climate of the area. In this classification, these climate types are designated considering seasonality, and annual and monthly mean values of precipitation and air temperature [100,101].
In the nominal representation of the climate type, there are uppercase and lowercase letters that symbolize the types and subtypes considered, respectively. Table 3 shows the division of the study sites by climate zones. Regions with a Cfa climate were the most studied, followed by Dwa and Cwa. It is noteworthy that, of the 579 papers evaluated, 41 (7.08%) addressed more than one climate zone in the same study. In addition, 31 papers analyzed different climates in the same study. As expected, there is a correspondence between the continents with most case studies and the climate of such studies.

Keywords
From the data generated by the VOSviewer software, a total of 3050 keywords were identified, corresponding to an average of five keywords per selected article (n = 579). Figure 6 shows the 338 most used keywords in the papers. The diameter of the circle in the image is proportional to the number of repetitions of each term.   All keywords are associated with topics and methodologies that can be employed for the study of UHI, except for the word "China", which stands out among the 15 most cited because of the high number of publications in the area.

General Objectives and Application Areas
The qualitative analysis of the objectives of the papers resulted in the identification of seven main areas addressed, presented in Table 5. The most recurrent objectives/areas in the papers focused on studying the UHI effect, using a comparison between the temperature amplitude of anthropically altered areas versus vegetated ones (416 papers-71.85%). Secondly, studies focused on analyzing the environmental characteristics, conditions, and impacts of UHI, considering biophysical, biochemical parameters, and data on water bodies and wetlands, totaling 61 (10.53%) papers.
It is important to emphasize the subjectivity of this classification, especially in the case of papers that worked with interdisciplinary topics, which were framed in the subclass that presented more affinity.

Equipment/Sensor Used to Obtain Information about the Thermal Data
Regarding the satellite/sensor selection to obtain the thermal data, 446 papers (77%) used RS data from a single source, 114 (19.69%) combined at least two sources of information, 16 (2.76%) used three data sources, 2 (0.35%) used data from four sources, and 1 (0.17%) combined data from five sources.
The most used satellites/sensors were Landsat, with 396 papers (68.39%), considering those that used complementary sources of data and 317 (54.75%) that used only Landsat products; MODIS, with application in 185 (31.95%) papers, being the exclusive source of data in 103 (17.79%) of the papers; and ASTER, which figures as the third most used, with 36 papers (6.22%), being the exclusive source in 15 (2.59). Considering the sum of the papers that cited the three satellites as exclusive sources, there are 435 (75.13%) of all the papers studied. It is emphasized that since most of the papers aimed a more detailed spatial monitoring of the study area, the authors adopted to use data from heliosynchronous orbits, rather than increasing the temporal resolution. Table 6 presents the list of five satellites most used and the number of papers that cited them. Considering data from complementary sources, the use of meteorological networks was cited in 41 (7.08%) of the papers, followed by in situ measurements in 14 (2.42%) publications. The data collection by drones, helicopters, and airborne were cited in 13 publications (2.25%).
The combination RS data with complementary information, such different resolution [102], in situ measurement data, especially in daylight [69] and model information, can be a good strategy to obtain data with better temporal resolution (due to the possibility of collection in the interval between satellite passes), as well as to improve spatial and spectral scales.
Considering RS sources, 324 papers (55.96%) used Landsat data, with 211 papers (36.44%) using this data exclusively. MODIS was the second most used, with 95 (16.41%) citations, where 51 (8.81%) were exclusive sources. Next, ASTER was cited in 32 (5.53%) papers, of which 11 (1.90%) were exclusive. Table 7 shows the five main satellites/sensors used in LULC computation. In addition to the derivation of LULC through satellite data, eight papers (1.38%) used information recorded by drones and sensors onboard helicopters and airborne, and two (0.35%) used exclusively these sources. One of the advantages of using data from these sources is the improvement in spatial and temporal resolution when compared to satellites, making it possible to record information at more detailed levels and, in areas of considerable heterogeneity, the distinction between land-use types can be useful for the researchers. However, the cost associated with this technology can be a limiting factor, given the small number of publications found.

Software Employed
Regarding RS data processing, 208 (35.92%) papers did not specify the name of the software used for processing the data. Considering those that did, 193 (33.33%) processed the data in a single software program, 111 (19.17%) used two software packages, and the other 67 papers (11.58%) used three to six different packages. Depending on the complexity of the study, a single software program may not answer all the processing needs, which justifies the combination and complementary use of different software. Figure 7 presents the number of papers concerning the number of software programs used.   ArcGIS and ENVI were the most used software programs in the papers analyzed. ArcGIS was used as the exclusive tool in 60 articles (10.36%) and 178 (30.74%) combined with other software. For ENVI, 28 articles (4.84%) used it exclusively and 93 (16.06%) combined with other software.
In addition to spatial image processing and analysis software, other software/systems were used to supplement information. For instance, FRAGSTATS, used for calculating landscape metrics, was used as complementary software in 39 papers (6.74%) and SPSS, applied for statistical processing, was used in 33 papers (5.70%).

Vegetation Index and Extra Data
Several complementary/auxiliary indices and indicators have been used in addition to the LST analysis. The ten indices/extra data most used are shown in Table 9.  Most indicators are associated with different topics, such as urban morphology and construction data, and indices obtained through RS data processing, such as NDVI, Normalized Difference Built-Up Index (NDBI)-an index that allows the identification mapping of built-up areas [103], census (demographic), environmental, social, and financial data. NDVI figures as the most cited indicator (190), followed by indices that use vegetation data such as EVI, FVC, and photosynthesis (105). It should be noted that the classification of the indicators considered the adoption of identical names by the authors.

Applied Statistical Methods
As for the statistical methods, 36 papers (6.22%) did not explicitly mention how the results were evaluated, 232 (40.07%) relied on one type of analysis, and 195 (33.68%) used two methods for data analysis. Table 10 presents the proportion of papers concerning the methods employed. The combination of methods can be a strategy when it comes to evaluating the bestfitting model for the data set, as well as complementing analytical information. Table 11 shows the list of the most applied statistics. It is noteworthy that, since there are papers that applied different statistical techniques, the sum of the percentages exceeds 100%. The correlation technique was the most applied in the studies selected, corresponding to 224 papers (38.69%), followed by the regression technique (210 papers, 36.27%), and descriptive statistics (128 papers, 22.11%). The data most applied in the correlations were the NDVI with the LST values.

Authors' Main Conclusions
This section highlights the main conclusions of the analyzed in this work, considering quantitative data, and highlighting some works for qualitative analysis, that were chosen to exemplify the study of UHI and SUHI in different contexts/areas, aiming for a systemic understanding of the causes identified for the interpretation of this effect. For that, included was the discussion regarding the methodologies employed and statistical analysis, suggestions, and mitigation strategies.

Quantitative Analysis
Considering the main authors' findings, in all studies selected, the UHI effect was evidenced. In addition to this finding, some papers discussed additional measures, highlighted in Table 12. Most of the papers identified the presence of the UHI in the studied areas and recognized the importance of local studies to identify the most appropriate mitigation measures. Two hundred and twenty-nine papers (39.55%) recommended mitigation measures, including such actions as interventions on the built-up structures, and maintenance/increase of water bodies and vegetation areas.

Qualitative Analysis
In this item, some UHI and SUHI studies will be presented to explore possible causes, consequences, methodologies, and mitigations. The selection of the papers was made considering the abstract and the diversity prioritized in the approach, without intending to emphasize any research or recommend standardized methodologies.
Each site presents its specificities, both in relation to LULC and Urban Functional Zones (UFZs), which refers to divisions of the urban surface, considering socioeconomic and landscape data [104].
To propose appropriate mitigation measures, it is necessary to understand these structures, as well as the climatic characteristics and local soil type/LULC [105], which vegetation would best fit the study area, the socioeconomic data, technical resources, and other additional indicators, relevant to public management.
The UHI and SUHI were associated with the formation factors shown in Table 13. Table 13. Highlight of some reasons for UHI and SUSHI´s formation.

Civil Constructions
Topographical features, building standards, building cooling, landscape configuration, urban morphology, spatial patterns of land use and occupancy [106][107][108], texture and spatial distribution [109], rural/urban areas with low vegetation density, in bare soil [110][111][112][113]. Three-dimensional urban design (such as height, volume, and surface area of the buildings, and the footprint and shadow volume) showed greater influence relative to two-dimensional ones (such as building rooftop area, street, and vegetation) [114]. Heterogeneity in the composition of built surfaces and areas [115].
Urban Size Influence on surface heat flux, soil/air temperature, humidity, and wind circulation, resulting in differences in urban-rural energy balance [117].

Environment
Wind [123], wind speed, and land surface conditions, especially in the morning, corroborate the advance of warm air [124], aerosols can influence/limit diurnal surface heating in dry seasons [125]. Table 14 presents some of the elements that can help in understanding the dynamics of UHI and SUHI.

Indicator
IV NDVI, NDBI [126], PV (that can show better results in the regressions applied when compared to NDVI) [127], landscape effects (that may have more impact than socioenvironmental data) [89].

Construction and soil
Characteristics and types of soils, LULC [128,129], imperviousness [130], albedo-exposed land areas [131], presence of surface contaminants, as hydrocarbons used in bioremediation processes, whose soils can present higher thermal signatures when revegetated [132].
As for data analysis, several comparative and correlative statistical tools were applied. Among the papers analyzed, some compared two of these tools and presented critical conclusions regarding the applicability and reliability of the results, namely: the Ordinary Least Squares (OLS) and the Geographically Weighted Regression (GWR).
The OLS is based on the assumption of data independence, ignoring the spatial effect between spatial items [146,147], whereas the GWR assumes that relationships between variables are influenced and change according to geographic distribution, allowing estimation of regression coefficients for each site [148][149][150][151].
GWR was found to be effective [151], more than OLS [152][153][154][155], and provides an analysis of the multiple spatial data layers with autocorrelated structures [156], predicts the behavior of variables with better accuracy, and can describe spatial non-stationarity [157]. It is recommended to apply GWR at finer scales than 480 m and in flat cities, since at larger spatial scales (720-1200 m), the model does not show significant differences from OLS [158].
As for the satellites used in the studies, most use Landsat [159,160], and there are considerations pointed out in some of the analyzed papers. MODIS is pointed out as interesting for research, but some studies concluded that temperature data are overestimated, especially during the daytime, and underestimated at the nighttime when compared to in situ data [161,162], or showed higher LST results in winter when compared to in situ data [163] and variability in measured temperature [164], especially in very dry or very wet weather conditions [165]. Landsat was found to be more efficient than ASTER in differentiating targets and temperature records. In areas with clear seasonality and thermal anisotropy, LST acquisition is challenging [166].
Regarding the consequences, Table 15 shows some examples. Table 15. Main consequences of the UHI and SUHI effect.

Economic
Regulation of the values practiced in the real estate market, whose properties most affected by the effect are cheaper [138], and preference in housing green areas and open spaces by high-income and educated families [88], which suggests understanding and studies regarding social vulnerability associated with residential segregation [172]. Increased electricity consumption in summer (due to the need for cooling/use of air conditioning), and minimized in winter (where the need for heating was not so necessary) [173], especially in areas inhabited by the high-income population, which have the means to regulate temperature, promoting thermal comfort [140]. It is noteworthy that the use of air conditioning can corroborate the increase of primary pollutant emissions and ozone generation [174], and increases the financial costs [175].
As mitigation measures, the analyzed papers suggest a series of actions that range from behavioral changes to structural changes in the analyzed areas, as presented in Table  16. Table 16. Mitigation of the UHI and SUHI effect.

Construction and LULC
Conversion of vacant land into green areas to minimize the effects of thermal discomfort [139], knowledge of local parameters, urban plantation, decentralization of urban areas [176][177][178][179], minimizing impervious areas and bare soils (such as abandoned land) [176], and diversification and balance of LULC elements, combining urban and vegetated areas, optimizing spatial configuration [178]; industrial relocation, buildings demolition, and brownfield redevelopment [180]; optimization of green land cover configuration [181] (in relatively small cities) or controlling sprawl (in larger cities) reducing the ratio of compaction in urban sprawl [182]; keeping the balance between urban and non-urban uses [183]; light-colored skyscrapers with glass curtain walls systems (which showed relatively low LST) [184].

Environmental
Preservation of water bodies [80,174,185,186] with simple geometry. The area of the waterbody is also a factor that influences temperature variation and, with high surface moisture, are less efficient than a water body (such as a river) [187] but can attenuate surface temperature better than vegetation cover, especially in dry seasons [188]. Increased proportion of vegetation at the expense of impermeable areas [189,190], distribution of trees and taking advantage of their shade effects [191], green surfaces/roofs [192][193][194][195][196][197][198][199][200]. The relationship between the presence of vegetation and cooling is not linear [201], and additional studies are needed to evaluate local specificities and needs. There are cases, for example, that associate the use of large trees with negative effects on nighttime cooling [202], due to higher heat retention as a consequence of a lower SVF [203]. Therefore, their characteristics must be evaluated to identify the most suitable and favorable vegetal species, arrangements, and strategies for heat exchange in a region. In addition to increasing the volume of vegetation areas at a site [204], evaluating the spatial pattern, shape, canopy cover [205], leaf density, and area of influence of heat regulation/shade effects [206][207][208]. Large parks did not present advantages over small ones [209].

Social, economic, and health
Encouraging the use of public transportation to reduce the circulation of cars and the emission of atmospheric pollutants [137]. Development of social measures, such as subsidies for electricity [173], and creation of heat warning and health surveillance systems, which alert users to the occurrence of heat waves [141].
As for punctual events, a study conducted in Dehradun, India, identified the increase in the number of hot spots and the minimization of the thermal comfort level associated with the post-lockdown period (COVID-19) considered in the paper on 14 April 2020, compared to data from 28 April 2019, 25 April 2018, and 8 May 2017 [145]. About the adaptations made for the 2008 Beijing Olympic Games, mitigations were carried out, such as increasing water-permeable areas, with the inclusion of green parks distributed among the central cities, minimizing the effects of the UHI [91].

Bibliographical References
The VOSviewer software was used to verify the similarities regarding the materials used for the theoretical basis of the papers. Table 17 shows the eight most cited references. It is important to note that the software performs combinations considering the similarity of the work cited; that is, the identical texts.

Authors
With the support of the VOSviewer software, 1475 authors were identified, with the number of publications ranging from 1 to 18. Table 18 shows the relationship between the number of publications, the number of authors, and the respective percentage. It was found that 1136 (77.02%) of the identified authors are involved in the publication of one article. The sum between the authors who published from one to three articles corresponds to 1.395 (94.58%), which concentrates most of the authors.
Of the 1475 authors, 719 (48.75%) showed interrelationship in publications. Using the software's association strength tool, the authors' interrelationships were identified and separated into 37 classes, distinguished by colors in Figure 8.  Table 19 presents the ten most cited papers in the analyzed papers.  Table 20 shows a scale of citations and the respective papers, in absolute numbers and percentages. A relatively small number of the analyzed articles were not cited (8.29%). Most articles had between 1 and 25 citations, with 345 papers (59.59%) in total, corresponding to more than half of the articles.

Discussion and Conclusions
The growing number of papers published on the UHI effect, especially since 2016, from 48 to 118 publications at the time, emphasizes the interest by the scientific community in the dissemination of this issue, considering its causes and consequences in several dimensions, such as environmental, social, and economic. This scientific relevance should add to the need for the promotion of an integrated urban climate planning [217], using valuable data as a primary input [218], as a way to tackle the challenges posed by this complex effect.
In the evaluated papers, it was observed that the methodological approaches and the impacts analyzed concerning the UHI were diverse and transversal, which adds to the complexity of this issue and corroborates with the need for researchers to know, as much as possible, the specificities inherent to the areas of study.
Most of the studies were concentrated in Asia (405-69.94%) which can be explained by interest from the scientific community in understanding the impacts of rapid urban growth in some parts of this continent that is often associated with the generation of intense UHI [219]. China was the most studied country.
As for climate regions, Cfa was the most recurrent. In this climatic region, the minimum temperatures range from 0 to −3 °C, there are at least four months with an average temperature of 10 °C, and at least one month with an average temperature above 22 °C. Precipitation does not vary significantly between seasons, i.e., there is no extremely dry month [220]. Considering that there are several areas in Asia and North America with this classification, it is natural that the number of publications regarding this particular climate region is higher (Figure 9, with Cfa in focus). The most recurrent keywords were "land surface temperature", "urban heat island", and "heat island". Considering the words used for the initial selection of papers for analysis, it was expected that "urban heat island" and "heat island" are recurrent, because they are associated with the acronym (UHI*) and "land surface temperature" is a product that can be derived even from RS, associating the keyword (satellite image*).
There are several techniques to obtain LST data, from thermal sensors present in heliosynchronous, geostationary, polar satellites, passive microwave radiometers, and data processed with statistical modeling. The choice of the most suitable methodology should consider the available resources (whether financial, data processing capacity, hardware, software, human resources, etc.).
The most used RS data came from Landsat missions, both for LULC identification and LST calculation. This preference can be explained by the free availability of historical data, the easy and intuitive access to the platform (USGS) [48], and the greater detail of the thermal band, especially of the transition and heterogeneous areas, compared to MODIS.
LULC data were mainly obtained by satellites, followed by models available from competent agencies such as the European Environment Agency (EEA) [222]. Techniques combining in situ data collection and/or UAV use were also applied and present advantages in data validation, providing greater spatial detail, especially in small and heterogeneous areas, and a cost-benefit analysis is indicated to assess the feasibility of their implementation.
NDVI was the most used index, especially for its effectiveness in distinguishing areas with the presence of vegetation, water bodies, and anthropic zones. Considering that the computation of the indices is performed from sensor data collected at the same instant when the thermals are obtained, as in the case of Landsat, MODIS, and ASTER, there is temporal concomitance of collection in the same study area.
Regarding the software used for data processing, among the papers that explained, ArcGIS stood out as the most recurrent. Although it is not free software, its interface is intuitive, and its data processing is effective. Moreover, some universities have licenses for research purposes, which can corroborate the adoption of this software in the methodology.
Correlation and regression models were the most employed to evaluate the quantitative relationship between indicators, such as LST, LULC, and additional data (socioeconomic, climatic, and from weather stations). The choice of these models is due to the methodological need to assess whether an action causes some reaction in another indicator, such as whether the increase in urban density causes changes in the local microclimate. In this sense, researchers must evaluate and compare relevant subjects, in a sample number adequate for the research objective, to avoid possible biases in the results obtained.
Studies that related the UHI effect to surrounding areas, especially vegetated surfaces, were the most recurrent, which converges with the publication by Voogt and Oke (2003), which mentions that this comparison between anthropic structures with the anisotropic thermal behavior of vegetated surfaces is more modeled and documented [12]. It is worth mentioning that there are publications focused on improving UHI methodologies and analyses, such as applications of mathematical and computational modeling to understand the behavior of the effect and the consequences for the microclimate, in addition to identifying gaps in mitigating measures, such as evaluating suitable vegetation to mitigate UHI influences at the site, pointing out the basic requirements regarding canopy, leaf density, and physical distribution at a site, etc.
The most cited bibliographies in the papers reviewed were from Voogt and Oke (2003) [12] and Yuan and Bauer (2007) [210]. The abstract of both publications is adherent to the UHI study and both methodologies employ RS data. Voogt and Oke (2003) [12] cite that despite the progressive scenario, research must transcend the application of simple correlations in data processing, as well as a qualitative description of urban surface land use and include indicators such as emissivity and RS techniques, both to analyze surface radiative parameters and to better describe the urban surface and ensure the effectiveness of the atmospheric models used in the studies, respectively. The authors also recommend advances in the spectral and spatial resolution of current sensors and the next-generation satellites to optimize the details of urban surfaces, as well as the availability of affordable high-resolution portable thermal scanners.
As for the consequences, the analyzed papers highlighted several areas that can be influenced and affected by UHI, and that are intrinsically associated with local specificities and anthropic activities developed in the region. In the atmosphere, for example, UHI can influence the displacement of air masses and the transport of particles harmful to health. From an economic perspective, some studies have addressed the influence of UHI on real estate market values, since, in regions where the phenomenon is more evident, costs can be reduced. As for the population's health, thermal comfort was pointed out as a concern, especially when considering low-income families, who may have restrictions on the acquisition and maintenance of equipment such as air conditioners. The risks of death associated with cardiovascular diseases, development of depression, and sleep complaints [142] were also addressed.
Concerning UHI and SUHI mitigation measures, existent papers can be used as a source of information for local managers and representatives to identify the effect, understand its dynamics, and minimize its impacts. Understanding the local microclimate, urban growth, LULC, and complementary data (economic, socioenvironmental, and public health), can be key elements to structure a strategy.
Regarding the possible applications and perspectives of UHI studies, it was possible to identify multiple approaches, areas, and objects of study, which reinforces that the applicability transcends the identification of areas affected by the microclimate, in terms of temperature, but that encompasses social, economic, and health factors and well-being of the population.
Moreover, most of the papers analyzed presented solutions for minimizing and mitigating the UHI effect, which can compose advisory material for managers. At this point, the need to understand the local specificities is reinforced, so that more assertive decisions can be made for each case. Publications aimed at improving the methodology of UHI studies, data processing, which can occur in software that operates in clouds, as in the case of Google Earth Engine [102], and the inclusion of artificial intelligence. Machine learning algorithms, such as artificial neural networks [223], are of great value to the research and can assist in the creation of new techniques. In addition, the combination of data from different sources can contribute to minimizing the technical limitations of collection equipment or techniques. A more systemic look at the UHI effect will favor the advancement in studies in the area, already identified as growing from the results presented in this review.
The design of a local bioclimatic model can be an interesting tool, as it is based on the strategic use of urban elements to minimize heat storage [224,225]. Some strategies that can be incorporated include: (i) analyze and adjust the distribution of built-up areas, reducing continuous high rise and dense urban development; (ii) add humid or vegetated areas to urban areas, corroborating for heat exchange; (iii) preserve and revitalize parks, water bodies, and vegetation surfaces, considering the extension, density, and size of the canopy aiming to maximize heat exchange capacity concerning the surroundings; and (iv) evaluate civil construction materials, prioritizing those with lower heat retention.
It should be noted that UHI and SUHI are a consequence of specific factors on the study sites, such as climatic, anthropogenic, and environmental characteristics, corroborating so that the choice of methodology, scale, data sources, study objectives, diagnosis, and possible mitigation proposals consider these particularities, transcending a standard methodology.
Despite this, there is recurring information and methodological approaches in the studies reviewed. Figure 10 presents a suggestion of the steps necessary for the study of IHU and SUHI, which does not disregard the need to assess local characteristics and identify approaches specific to the study area. We believe that this work can contribute to a comprehensive overview of the current state of the art regarding the UHI effect using RS data, considering the last 20 years.
Author Contributions: C.R.d.A. participated in the manuscript design, study selection, analysis and selection of study data, data processing in software, and systematic review. A.C.T. participated in the manuscript design, study selection, analysis and selection of study data, and systematic review. A.G. participated in the manuscript design and systematic review. All authors have read and agreed to the published version of the manuscript.
Funding: This work was funded by National Funds through the FCT-Foundation for Science and Technology and FEDER, under the projects UIDB/04683/2020 and PT2020 Program for financial support to CIMO UIDB/00690/2020.