ASSESSING THE MICRO-SCALE TEMPERATURE-HUMIDITY INDEX (THI) ESTIMATED FROM UNMANNED AERIAL SYSTEMS AND SATELLITE DATA

Direct heat and moisture conditions can lead to discomfort for humans and animals and can decrease health performance. The discomfort index or temperature-humidity index (THI) represents an important indicator that measures the heat sensed by humans for different climate conditions. In extreme situations, heatstroke may occur, which in unfortunate cases will lead to death. Many research studies have been conducted on the urban heat island (UHI) phenomenon, although a majority of such work focuses on regional-scale analyses and emphasizes the thermal trend through larger administrative units. Fewer micro-scale analyses have been performed at the local scale to detect the potential area for increased THI within a city. This work seeks to estimate the THI at the micro-scale level by utilizing the thermal camera on-board of unmanned aerial systems (UASs). Thermal information of the surface and visual images are collected by the UAS, while a thermohygrometer is used to collect the air temperature and the relative humidity at the ground surface for ground truth information. Solar radiation and wind exposure modeled from digital surface model (DSM) and normalized difference vegetation index (NDVI) data are used as explanatory variables, and a random forest machine learning method is implemented to model the spatial distribution of the THI. The results and discussion will provide future possibilities for micro-scale analyses of the UHI.


INTRODUCTION
In recent years, as an urban environmental problem, the urban heat island (UHI) phenomenon that accompanies urbanization has become prominent not only in Japan but also in other countries around the world. UHI is a cause of increased tropical nights and associated health hazards, such as heatstroke in urban areas, changes in ecosystems due to the overwintering of mosquitoes that carry infectious diseases, and extreme torrential rains (Deilami et al., 2018). It is a serious problem that greatly affects urban space environments and leads to associated health issues (Royé, 2017). Regarding global warming and climate change, irregular climatic conditions appear in various urban area locations. Japan is also facing irregular cases of climatic conditions (Imada et al., 2019). Extreme heat conditions can threaten health, and the temperature as well as humidity level can affect heat-related diseases (Fujibe et al., 2018;Ito et al., 2018). To understand the human sensations associated with different climatic conditions, various indices are prepared that consider both temperature and humidity. The discomfort index or temperature-humidity index (THI) represents how humans feel discomfort depending on the combination of low-high temperature and relative humidity. With the increasing development of urban areas leading to the expansion of hotter areas in combination with a changing climate, more threats to health will occur. Therefore, determining the heat hotspots or the distribution of THI in urban areas (i.e., where is the area exposed to higher discomfort) has become an issue. Factors in the formation of UHI include low albedo material, complex urban morphology, waste heat, low density vegetation in urban spaces, etc. (Aflaki et al., 2017). Many studies have been conducted on the heat in urban spaces (Kolokotroni et al., 2012). However, the results are tabulated to determine the trend for the whole city or district.
Such works estimate heat trends in a relatively large administrative unit (Imhoff et al., 2010), and fewer work has focused on urban heat in microscale units that form the city, such as individual buildings and roads. If there are solutions for observing and mapping potentially high THI areas, then strategic decisions may be made for mitigating various issues caused from UHIs. In this work, we focus on developing a new method for observing and estimating THI mapping at an extreme micro-scale unit by utilizing the unmanned aerial systems (UASs) and satellite remote sensing data.

UAS and Camera Specification
An unmanned aerial system (UAS) is utilized to collect remote sensing data of the site (Figure 2). This work utilizes the DJI Matrice 210 multicopter UAS (DJI, Shenzhen, China) for the aerial system, and the DJI Zenmuse-XT2 (DJI, Shenzhen, China) is the sensor that is equipped on board the UAS. Matrice 210 is an industrial application UAS improved for various flight performance and safety features. Zenmuse-XT2 has two sensors: visible and thermal. The camera can simultaneously collect both visible and thermal information; therefore, it is suitable when both data are needed.

GNSS and Thermohygrometer
The Reach (Emlid, Hong Kong, China) global navigation satellite system (GNSS) was utilized for recording the geographical coordinates ( Figure 3). The small and light-weight GNSS makes it simple to record positioning at a higher precision using multiple L1 GNSSs, e.g., GPS (Global Positioning System), QZSS (Quasi-Zenith Satellite System), GLONASS (Global Navigation Satellite System), Galileo, BeiDou, etc. The equipment is used together with the thermohygrometer to record the temperature and relative humidity at each geographical location. A MJ-UDL-20 thermohygrometer logger (SATOTECH, Kawasaki, Japan) can log the temperature and humidity with a minimum interval of 2 seconds. The temperature and humidity range can be logged from -20 to 70 ℃ and 5 to 95 %, respectively. The thermohygrometer is used to collect the ground truth data of the study site.

Flight Design
On September 9 th , a flight campaign was conducted at study site 1 from 10:58 AM to 11:03 AM and at site 2 from 14:13 PM to 14:18 PM. The UAS flew at an altitude of 75 m in a single grid path ( Figure 4) and collected multiple aerial photos during the flight. The ground sampling distance (GSD) at the 75 m altitude was approximately 1.75 cm and 6.7 cm for the visible and thermal sensor, respectively. Forward (side) overlap was set to approximately 80 % (75 %), which are the parameters of focus for the thermal camera. This setting will correspond to approximately 87 % (85 %) of the forward (side) overlaps for the visible sensor.

Collecting Ground Truth Data for the THI
The Reach GNSS rover and thermohygrometer was attached on top of a safety helmet, and another GNSS was set on the ground nearby using a tripod or equipment fixed to the ground as the base station ( Figure 5). During the aerial survey of the UAS, a crew equipped the helmet and walked around the study site to record both the geographical position, temperature and humidity from each devices. The Reach GNSS was set to observe GNSS signals from GPS, QZSS, Galileo and Beidou at a logging frequency of 1 Hz, and the thermohygrometer was set to a 2 second interval for the recording. Furthermore, the GNSS signal data from the base station was used together with the rover data for the Post Processing Kinematic (PPK) to improve the precision of the positioning. Then, the timestamp of the GNSS log and the timestamp from the thermohygrometer were matched to generate point vector data of the temperature and humidity at each geographical location. When the point data were generated, an Ordinary Kriging interpolation was performed on the temperature and humidity data to develop a spatial map. Using equation 1 (Yoo and Chung, 2018), the Temperature-Humidity Index (THI) was computed for the study site. The THI represents how much discomfort a person feels from the combination of temperature and humidity. According to Yoo and Chung (2018), a THI below 68 indicates comfort, a THI from 68-75 indicates when discomfort begins to occur, a THI from 75-80 will make 50 % of the people feel discomfort and a THI >80 will lead all people to feel discomfort (critical condition affecting health). The index is used to facilitate computing the THI, and the function is similar to what is often utilized in the study country.  ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume V-3-2020, 2020 XXIV ISPRS Congress (2020 edition)

Study Area
The study area is located at a closed elementary school and junior high school in Maebashi City, Gunma Prefecture, Japan ( Figure  1) located near the central area of Maebashi City. The climatic conditions of the area vary within a year. Our field work was carried out on September 9, 2019. The daily average temperature on September 2019 was 24.6 ℃, where the maximum reached 37.0 ℃. It was an unusual month with hotter days. The total precipitation was 75.5 mm, and the average relative humidity was 70 %, which indicate very hot and moist conditions and higher discomfort. Figure 1. Overview of the study area (Gunma Prefecture, Maebashi City, Japan).

Aerial Image Mosaics and Digital Surface Model:
Using the collected visible and thermal aerial images, the Structure from Motion (SfM) technique is implemented using Metashape Pro ver. 1.6.0. (Agisoft, St. Petersburg, Russia) to generate a mosaicked image of the whole study area. Various settings are available in the software for reconstructing the 3D model of the scenery. For the visible image, the alignment setting is set to "high" with default tie and key points, "high" for dense point generation and further processing to generate a digital surface model (DSM) with a dense point cloud, and using the DSM for orthorectifying the image. The DSM is generated by the inversed distance weighting (IDW) interpolation algorithm. For the thermal image, the "highest" setting is for the alignment process with default tie and key points and the "ultra-high" setting for dense cloud generation. Instead of creating a DSM, mesh data are built from the dense points and a mosaicked thermal image is generated with the mesh data.

Solar Radiation and Wind Exposure:
Using the generated DSM, the total solar radiation, direct solar radiation and diffuse solar radiation are modeled using the ArcGIS Pro ver. 2.4. (ESRI, Redlands, USA) for the same date and time as the observation. For site 1, the solar radiation is modeled from before sunrise to 11:00 AM. For site 2, it is modeled from before sunrise up to 2:00 PM. System for Automated Geoscientific Analyses (SAGA) GIS (Conrad et al., 2015) ver. 6.0.0. is used to generate a wind exposition index (Boehner and Antonic, 2009) to determine the degree of exposed areas. The local area had some degree of wind, which will affect the local heat exchange related to evaporation (Lei et al., 2018). It is difficult to map the wind velocity for each area; however, by modeling the magnitude of wind effects, it is expected to work as an alternative.

NDVI:
The Sentinel-2B Multispectral Instrument (MSI) Level 2A (L2A) product from the observation date of October 10, 2019 was collected via the webpage. The L2A product is processed with atmospheric correction and derived from the associated level-1C products and then converted into the surface reflectance (bottom of atmosphere (BOA) reflectance). Using the red (band 4) and near infrared (band 8) bands, the normalized difference vegetation index (NDVI) was computed for the study site and clipped for the study area. The NDVI data were resampled to the same resolution as the thermal image using the bilinear method for further processing purposes.

Modeling the THI with the Random Forest Algorithm
Using the thermal data (total, direct and diffuse solar radiation; wind exposition index; and NDVI data) the random forest machine learning method is implemented to model the THI of the study site. R Studio (Integrated Development for R. RStudio, Inc., Boston, MA) and the random forest library (library "randomForest") were installed to process the modeling. A grid search method (Garcia et al., 2011) is implemented to train various numbers of models with slightly different parameters: node size = 10, 20…50; mtry = 1, 2…5; and ntree = 500, 600…1500. Bootstrap resampling is implemented with a sample size of 90 % and iteration = 5. From the multiple generated models, the final tuned model is chosen from the lowest errors, and it uses the parameter of node size = 10, mtry = 4, ntree = 700. The modeled THI is validated via comparison with the reference set. To evaluate the model accuracy, the correlation of determination (R 2 ), Mean Absolute Error (MAE) and the Root Mean Square Error (RMSE) are computed. We would like to emphasize that the R 2 here is considered as how well the predicted model fits to the 1:1 line. To avoid confusion, further it will be denoted as fitting R 2 (equation 2).
, recommending a fitting R 2 > 0.6 indicates high correlation, thus succeed in model development. Figure 6 shows the successfully developed ortho mosaic imagery of the visible and thermal images for site 1 and site 2. For the thermal image, at site 1, the trend shows that bare soil areas are warmer than vegetation areas (e.g., trees). However, the building roofs on the east side show rapid heating, which is possibly from the material and the exposure to direct solar radiation. The northern buildings are composed of white material that has higher albedo and lead to less heating. Buildings with darker rooftops are also cooler compared to buildings with hotter roofs, which seems to be due to the moss grown on the roof. At site 2, the thermal trend is similar to that of site 1 for bare soil and vegetation, and exposed asphalt has been heated and white roof buildings on the south are also heated but show some heat gradient because the roofs have a dome-like shape, which leads to differences in exposure to direct solar radiation. The thermal conditions are clearly visible through UAS sensing.  Figure 7 shows the mapping result from the interpolation of the temperature and humidity vector points collected through the ground survey. The approach for using a light-weight GNSS with the thermohygrometer to log the spatial location has succeeded, and the kriging interpolation lets us interpret the trend of the THI on the ground. However, the data were collected by walking the ground surface; therefore, the area of the building rooftops is not represented correctly. The modeling process excludes those areas for the training samples. At site 1, bare soil areas showed higher THI and the northern part between one tree and the building on the north show small patches of higher THI. Near the tree, a lower THI is observed, and the area between the buildings and behind the buildings where direct solar radiation is blocked shows a lower THI. During the survey, it was confirmed that under trees and areas with shadows are much cooler than directly exposed areas. The ground THI trend seems to delineate such areas successfully. At site 2, it was clear that the bare soil area and sparse vegetation areas show lower THI, while the surrounding areas clearly show higher THI, which is due to the condition of the land use characterized by materials such as asphalt. The northeast area shows a gradient change from the center position, which leads to a slightly higher THI. The northeast area has a smaller fraction of grassy vegetation cover, while the lower THI area has a higher fraction. The temperature seems to cool down in those vegetative areas compared to exposed bare soil areas.  Figure 8 shows all the variables generated from the DSM and the NDVI collected from the satellite. The ultra-high resolution UAS imageries provide a good solution for developing a detailed solar radiation model and wind exposure of the sites. Few studies have collected such detail modeling results at building to street units, and the use of the UAS made it possible to enhance the scale to what we usually face during daily life. Just by adjusting the accumulation of the solar radiation time, the area can be interpreted for higher radiation or shadow areas. For the NDVI data, although the Sentinel-2 satellite is at the original 10 m resolution, the NDVI trend for each study site can be distinguished well. It is clear that at site 1, building areas have lower NDVI values, trees show higher NDVI, and bare soil areas show lower NDVI. At site 2, the trees along the southeast to the north direction show a higher NDVI, the bare soil area at northeast side shows a lower NDVI, and the southwest part shows a slight increase of the NDVI. Therefore, although the resolution is 10 m, the slight green changes on the playground can be sensed. These data are used for the explanatory variables in the THI modeling process. Solar radiation shown here is only the total radiation for the visual purpose. Note that site 1 is modeled up to 11:00 AM, while site 2 is modeled until 2:00 PM for solar radiation. Figure 9 shows the final result of mapping the THI at a microscale unit. On site 1, since the time of observation was 11:00 AM, most of the area shows a lower THI. A higher trend is observed for rooftops; however, those areas were not considered during the model training process. Thus, the actual THI condition is uncertain. However, some possible trends can be observed since those areas are exposed to more heat; thus, the higher THI could be expected. The vegetation area seems to show a lower THI than the bare soil area, which again seems to reflect the actual environment. For site 2, at 2:00 PM, most of the asphalt area shows high THI values. The bare soil area shows a lower THI compared to the asphalt area, and the bare soil area with a greater fraction of vegetation also shows a lower THI. Some areas behind buildings show a low THI even if the land use is asphalt. Such areas are cast in shadow; therefore, even if the land use includes highly heat absorbing material, the area would show lower THI compared to directly exposed areas. The approach for utilizing the UAS has succeeded in the modeling. From the validation samples, the modeled THI is compared with the reference THI (Figure 7). The fitting R 2 showed a good fit of 0.9596, and the RMSE was 0.2371 and the MAE was 0.1526 ( Figure 10). The larger residual above a THI of 83.3 was associated with sensor lag during ground surveying. Entering a shadowed area after a highly exposed THI area for only few seconds was not enough for the sensor to stabilize to the actual temperature and humidity; thus, the reference sample shows a higher THI than the modeled THI, and vice versa. These models can enhance our understanding of potential areas with higher THIs, which can further contribute to decision making for sustainable city planning. Hotspots of high THI areas can be identified to so that announcements can be made to prevent heatstroke. These issues can be further analyzed in-depth if the thermal conditions at different temporal times are collected throughout a day for understanding the THI transition through space and time. Figure 11 shows the contribution of the variables to the developed model. The figure indicates the mean decrease in node impurity or the Gini index. Diffuse solar radiation has the greatest effect, followed by thermal information based on the UAS, modeled total solar radiation, NDVI, wind exposure and modeled direct solar radiation. This finding indicates that directly observed thermal information corresponds to the THI to some degree but will not always reflect such data correctly; moreover, using the modeled solar radiation is effective for retrieving the THI condition at each site. Figure 9. THI modeling result for site 1 and site 2 implementing the random forest regression (Note that site 1 is at 11:00 AM and site 2 is at 2:00 PM). The scaling and the color bar is synchronized for the interpretation purpose. Figure 10. Evaluation between the reference set and the modeled result. The R 2 here indicates how well the model fits to the 1:1 line. Figure 11. Importance of each variable explained as the mean increase in node purity (i.e., decrease in node impurity: Gini Index). The upper area shows more influential variables used in the random forest model for estimating the THI.

CONCLUSIONS AND FUTURE WORK
This work has focused on modeling the temperature-humidity index (THI) at a micro-scale unit by utilizing various remote sensing data obtained from both unmanned aerial system (UAS) and satellite data. The UAS was utilized to collect ultra-high resolution aerial imageries for both visible and thermal data, and mosaic imagery was produced for the scene. The structure from motion (SfM) method was used to develop a 3D model of the scene, and a digital surface model (DSM) of the area was extracted. The DSM was used to model the solar radiation and wind exposure, and the normalized difference vegetation index (NDVI) was collected via Sentinel-2B satellite. All of these data were used as the explanatory variables for the THI modeling using the random forest machine learning method. The results indicate that the micro details of THI trends, which cannot be seen from conventional approaches by using low resolution data, were successfully identified, thus revealing potential areas that might be exposed to higher THI and represent potential threat areas that could lead to health-related issues. Diffuse solar radiation was the most influential variable for the developed random forest model. The challenge of this study was integrating remote sensing data from multiple platforms at different resolutions or acquisition times. Indices such as the NDVI were carefully selected at the nearest day to the UAS observation, and the authors acknowledge that it would be the best if all data are collected on the same day to improve the modeling. For future work, a broader spatial extent can be considered, such as for the whole city, to identify the temporal and spatial characteristics of the high THI environment and the potential locations where such risks lie. Such work was not possible in the current study due to complex regulations and concerns for UAS flights within the city center. In addition, similar challenges may occur for where the impact of THI is considered to be more serious than Maebashi City, i.e., larger cities than Maebashi City, such as Tokyo, and subtropical or tropical cities outside of Japan, where more severe weather conditions are observed. The progress will be shared with the city in the future for supporting the development and planning of smart cities.