Dense Temperature Mapping and Heat Wave Risk Analysis Based on Multisource Remote Sensing Data

As high temperature and heat wave have become great threats to human survival, social stability, and ecological safety, it is of great significance to master the spatial and temporal dynamic changes of temperature to prevent high temperature and heat wave risks. The meteorological station can provide accurate near-ground temperature, but only within a specific space and time. In order to meet the needs of large-scale research, spatial interpolation methods were widely used to obtain spatially continuous temperature maps. However, these methods often ignore the influence of external factors on temperature, such as land cover, height, etc., and neglect to supplement temporal-wise information. To deal with these issues, a joint spatio-temporal method is proposed to obtain dense temperature mapping from multisource remote sensing data, which combining a geographically weighted regression model and a polynomial fitting model. Besides, a heat wave risk model is also built based on the dense temperature maps and population data, in order to evaluate the heat wave risk of different areas. Accuracy evaluations and experiments have verified the effectiveness of the proposed methods. Case study on the four cities of Zhejiang Province, China have demonstrated that areas with higher degree of urbanization are often accompanied by higher heat wave risks, such as the northern part of the study area. The study also found that the heat wave risks have presented a centralized distribution and spatial autocorrelation characteristics in the study area.


I. INTRODUCTION
E XTREME heat events are growing more and more frequent as a result of long-term climate change [1]. Heat waves, or extended periods of extreme heat, can be harmful to human wellbeing [2], agricultural production [3], and ecological stability [4]. According to recent studies, over 200 urban regions around the world experienced excessive heat during the past 40 years, from 1973 to 2012 [5]. This indicates that the occurrence of heat waves is rising, which poses a serious threat to global ecosystems The authors are with the Guangdong Provincial Key Laboratory for Urbanization and Geo-simulation, School of Geography and Planning, Sun Yat-sen University, Guangzhou 510275, China (e-mail: liumx23@mail2.sysu. edu.cn; lixzh23@mail2.sysu.edu.cn; chaizhq@mail2.sysu.edu.cn; chenanq29@ mail2.sysu.edu.cn; zhangyy388@mail2.sysu.edu.cn; zhanggiszsu@163.com).
Digital Object Identifier 10.1109/JSTARS.2023.3260467 [6]. Therefore, it is crucial to monitor the dynamics and analyze the impact of extremely hot weather and heat waves [7], [8].
As a disastrous weather, the research on climate change in recent years has paid more and more attention to the impact of high temperature climate and heat waves [9], [10], [11]. Anderson and Bell [12] found that, compared with nonheat wave periods, the mortality rate could be increased by 3.74% during heat waves. Thereinto, the death risk will increase by 2.49% for every 1°F increase of heat wave intensity. Zhu et al. [13] verified that hot weather will reduce the global average planting frequency, since the drought caused by heat waves will interfere normal water supply and harvest window of many crops. Keenan and Riley [14] believed that continuous warming will lead to a significant reduction of temperature-sensitive vegetations. It can be seen that high temperature and heat waves have very detrimental effects on human health, agriculture, and ecosystems [15], [16]. Therefore, intensive and dynamic monitoring of nearground temperatures is necessary to minimize the occurrence of potential heat hazards and heat wave threats.
The near-ground temperature is a popular and practical temperature indicator, which can reflect the temperature that the human body can perceive on the ground [17], [18], [19]. Usually, meteorological station data is the first-hand data to obtain the near-ground temperature distribution [20], [21]. Although meteorological station data can provide the temperature of the point where the station is located, accurate and precise heat wave analysis requires continuous temperature distribution [22]. Therefore, previous studies have concentrated on spatial interpolation to obtain polygon-wise surface temperature from point-wise meteorological station data [23], [24]. For example, based on the meteorological station data, Wu and Li [25] refined the temperature maps of the study area through the spatial Kriging interpolation method. Li and Heap [26] explored the performance of three commonly used spatial interpolation methods in temperature interpolation, including nongeostatistical interpolation methods, geostatistical interpolation methods, and combination methods. These studies have proved that spatial interpolation techniques enable the quick acquisition of spatially continuous temperature maps and provide fundamental data for the study of extremely heat weather.
However, there are still a few deficiencies with current temperature mapping methods. First, the simple spatial interpolation only considers the spatial distance differences, and ignores the influences of external factors to near-ground temperature, such as land cover and terrain height [27], [28]. Second, current This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ approaches are usually limited to retrieving annual, monthly, or daily time-series temperature data, which is difficult to satisfy in-depth analysis relying on dense temperature distribution (like hourly) [29]. Therefore, it is a pressing issue to figure out how to use multisource remote sensing data more comprehensively to obtain dense temperature maps.
In order to grasp the distribution of temperature more precisely and forecast the risk of heat waves, this study brings forward a dense temperature mapping method and a heat wave risk analysis method based on multisource remote sensing data. The main contributions of the study can be summarized as follows.
1) A joint spatio-temporal method is proposed, which contains a geographically weighted regression (GWR) model [30] and a polynomial fitting model [31], for dense temperature mapping from multisource remote sensing data. 2) A heat wave risk model is built based on the dense temperature maps and population data, in order to evaluate the heat wave risk of different areas. 3) A case study on four cities of Zhejiang Province, China using the proposed methods has been conducted. Based on the hourly temperature maps from January 1, 2011 to December 31, 2014 in the study area, experiments and analysis have fully demonstrated the validity of the joint spatio-temporal method and the heat wave risk model. The rest of this article is organized as follows. Section II describes the study area and involved materials, while Section III demonstrates the methodology. The experimental results will be evaluated and analyzed in Section IV. Finally, Section V concludes this article.

A. Study Area
Four cities, including Huzhou, Jiaxing, Hangzhou, and Ningbo, in Zhejiang Province, China are selected as the study area. Located in the south wing of the Yangtze River Delta on the southeast coast of China, Zhejiang Province (27°02'-31°11'N, 118°01'-123°10'E) has a monsoonal humid climate and usually has hot weather in summer. The average annual temperature of Zhejiang Province is between 15°C and 18°C, while the annual sunshine hours of that are between 1100 and 2200 h, and the average annual precipitation of that is between 1100 and 2000 mm. Therefore, as shown in Fig. 1, the four cities in northern Zhejiang are chosen as the study area in this article. The brief introduction of each city is as follows.
1) Huzhou: Huzhou City locates in the north of Zhejiang Province, with an area of 5820 square kilometers. Its terrain roughly inclines from southwest to northeast, with mountains in the west and plain water system in the east, belonging to the north subtropical monsoon climate zone. 2) Jiaxing: Jiaxing City is located in the northeastern part of Zhejiang Province, with an area of 3915 square kilometers. It is one of the core cities in the Yangtze River Delta.

3) Hangzhou: Hangzhou City is the capital of Zhejiang
Province, located in the south of the Yangtze River Delta and the Qiantang River basin, is an important central city in the Yangtze River Delta and the transportation hub of southeastern China. It has a subtropical monsoon climate. 4) Ningbo: It is an important port city on the southeast coast of China, located at the southern wing of the Yangtze River Delta. The topography of Ningbo is high in the southwest and low in the northeast. It also has a subtropical monsoon climate.

B. Datasets and Preprocessing
1) Near-Ground Temperature Data: The basic data for this research area is the near-ground temperature data from 600 meteorological stations in the study area. All these meteorological stations are provided and maintained by the Zhejiang Meteorological Bureau. According to Fig. 1(b), the meteorological stations in the study area are densely and uniformly distributed. Each station has recorded daily temperature of 0-3, 9-15, and 21-23 h, including hourly temperature observations, date, time, and latitude and longitude of each station. All of these data have undergone rigorous quality control to eliminate the influence of random errors, with a maximum error of 0.1°C for 91% meteorological stations. Therefore, temperature information from these stations between January 1, 2011 and December 31, 2014 are employed for subsequent study.
2) MOD13A2 NDVI Product: The MOD13A2 NDVI product [32], which is built based on MODIS satellite imagery, with a spatial resolution of 1 km and a temporal resolution of 16 days, is introduced in this study to provide information of vegetation coverage. Considering vegetation has a high reflectance in the NIR band and a significant absorption in the RED band, the NDVI is designed to be calculated as follows: where ρ NIR represents the reflectance in the near-infrared band and ρ red represents the reflectance in the visible red band in MODIS imagery. It is a frequent and practical method to use NDVI to denote vegetation distribution. Therefore, the MOD13A2 NDVI product is utilized in this study to obtain vegetation coverage in the study area, as shown in Fig. 1(c).
3) SRTM3 DEM Product: In order to obtain DEM information of the study area, the Shuttle Radar Topography Mission (SRTM) topographic product [33] is utilized, which is one of the most accurate and high-resolution global topographic data set. The SRTM product has two resolutions of about 30 m and about 90 m, with an absolute horizontal and elevation accuracy of 20 m and 16 m, respectively. Limited by data availability, the 90 m SRTM3 DEM data is used, as shown in Fig. 1(d).

4) Global Annual Urban Dynamics (GAUD) Data:
The GAUD data [34] is used to acquire the impervious surface coverage of the study area. The GAUD data is the first continuous time-series (1985-2015) urban expansion dataset with 30-m resolution, which is obtained from Landsat remote sensing image data based on the Google Earth Engine (GEE) platform. Among them, the GAUD data of 2011-2014 were retrieved to support the research in this article. As can be seen in Fig. 1(e), the impervious surface in the study area mostly concentrates in northeastern region.

5) MCD12Q1 Land Cover Product:
To verify the influence of different land cover types to temperature, the MCD12Q1 Land Cover product [35] is adopted in this study, which has a spatial resolution of 500 m. According to the definition of the International Geosphere-Biosphere Program, the MCD12Q1 Land Cover product includes 17 major land cover types, including water, evergreen coniferous forest, evergreen broad-leaved forest, deciduous coniferous forest, deciduous broad-leaved forest, mixed forests, dense shrubs, sparse thickets, woody savanna, savanna, the grass, permanent wetlands, agricultural land, cities and buildings, agricultural lands/natural vegetation, snow and ice, sparse vegetation and uninhabited areas classification. Based on the MCD12Q1 data, the influence of land type on the temperature fitting function can be verified more comprehensively.

6) Population Data:
To study the risk of heat waves, we introduce population data to construct a heat wave risk model. The population data used come from China's 2010 census data, which includes the total population (total, male, female), population of different ages (0-14 years old, 15-64 years old, 65 years old and above), number of households, local population, and other information. Based on the needs of this research, the total population and population by age are used. Fig. 2. Flowchart of the proposed method for temperature mapping and heat wave risk analysis. The joint spatio-temporal employs a GWR model and a polynomial fitting model to obtain dense temperature maps from multisource remote sensing data, on which a heat wave risk model will be applied to evaluate heat wave risks.

A. Overview
The overview diagram of the proposed method is demonstrated in Fig. 2, which consists of the following steps.
1) First, a GWR method is used to obtain spatially continuous temperature data, in which the meteorological station data, MOD13A2 NDVI Product, SRTM topographic product, and GAUD data are utilized, in order to provide multidimensional information of internal and external factors. 2) Second, a polynomial fitting function is applied to obtain dense temperature maps both spatial-and temporalwisely, based on the spatially interpolated temperature maps at first step. 3) Given the dense temperature maps obtained by the proposed joint spatio-temporal method, a heat wave risk model is applied to analyze the risk of high temperature and heat waves to different populations.

B. GWR Model
Relevant studies have shown that there is a close relationship between temperature change and factors such as NDVI [36], DEM [37], and land cover types [27], [38]. However, traditional temperature mapping method based on simple spatial interpolation failed to consider the influence of various external factors. Therefore, in this section, both NDVI, DEM, and impervious surface will be visualized as independent variables, and a geographically weighted model will be adopted to obtain temperature maps of the study area. Among them, NDVI, DEM, and impervious surface ratio are extracted from MOD13A2 product, SRTM data, and GAUD data, respectively.
Specifically, the original 600 meteorological station data will be randomly divided according to the ratio of 9:1 into training set and validation set, as shown in Fig. 3(a). Data of the training set will be used as the input of the GWR model to optimize model parameters and obtain spatially constant temperature The process of using GWR to get temperature maps can be expressed as follows: where (u i , v i ) denotes the location of the ith regression analysis point, , respectively, represent the three regression coefficients.

C. Polynomial Fitting Model
Since the original meteorological station data only provides temperatures at 0-3, 9-15, and 21-23 h every day, it is impossible to obtain dense time-series temperature data for more precise research simply by the GWR model. In order to fill in the temporal gaps, the polynomial fitting model is applied to the spatially constant temperature maps obtained by the GWR model, so as to get the hourly dense temperature map of the study area of total four years from January 1, 2011 to December 31, 2014.
The polynomial fitting model aims to fit the missing data from known variables with a polynomial function, where the least square method will be used to optimize the function. In order to find the appropriate terms of the polynomial function, experiments on the hourly temperature data of the China Meteorological Information Center (CMIC) were conducted first. Comparative experiments were conducted on the CMIC data using the least squares loss function. Finally, a six-degree polynomial fitting function is finally adopted to fit missing values of the hourly temperature of each day. The loss function and the polynomial fitting function can be expressed as follows:

D. Heat Wave Risk Model
In order to explore the risk of heat wave weather of different areas, a heat wave risk model is built based on the spatio-temporal dense temperature map obtained by the proposed method. Previous studies have proved that high temperature has different effects on the human health of different ages. In particular, teenagers and the elderly usually have lower tolerance to heat waves. Therefore, when measuring regional heat wave risk, it is very important to take special groups into considerations. Therefore, taking the street as the basic unit, a heat wave risk index (HRI) is constructed by combining heat wave days and populations of different ages. The HRI to calculate yearly heat wave risks can be denoted as follows: where n 1 is the total number of people of a street, n 2 is the number of young people aged 0-14 of a street, and n 3 is the elderly population aged above 65 of a street. days refers to the number of days with daily maximum temperature higher than 35°C each year. In addition, in order to test whether the heat wave risk has autocorrelation spatially, the global Moran's I [39], [40] was introduced to conduct spatial autocorrelation analysis on the yearly HRI results, which can be expressed as follows: In the formula, i = j, n is the number of observations of the variable x, i, j represent two observation positions, x is the observed value,x is the average value of all observed values, and w ij represents the spatial weight matrix.
The value of Moran's I ranges from −1 to 1. When I>0, it means that the attribute values of all regions have positive correlation spatially; when I = 0, it means that the regions are randomly distributed without spatial correlation. When I<0, it indicates that the attribute values of all locales are negatively correlated spatially.
In addition, the results of global Moran's I were tested using the z score and p-value. The z score is calculated by the expection and the variance, and their formulas can be represented as follows:   6. Accuracy evaluation on polynomial fitting results. "ABE" means "absolute error," and "STD" means "standard deviation." V ar (I) = n 2 (n−1) 2 When performing a significance test on Moran's I value, at the 5% significant level, when the z score is greater than 1.96, it means that a certain phenomenon has a significant correlation distribution within the research range; if the z score is between [−1.96, 1.96], it means that a phenomenon has no significant correlation distribution within the research range, and the spatial correlation is also weak; if the z score is less than −1.96, it means that a certain phenomenon has a negative spatial autocorrelation distribution within the research range.

A. Accuracy Evaluation on Temperature Maps
The temperature mapping results by the GWR model are tested by the meteorological data of the validation set in Fig. 3(a).
The accuracy evaluation results from 2011 to 2014 are all depicted in the box plots in Fig. 4, where the average value of the absolute value error in each year is calculated. It can be seen that the errors mainly range between 0.5 K and 0.9 K. The year of 2013 has the highest error but also within 1 K. In general, the accuracy verification results show that the temperature map predicted by the GWR model can meet the accuracy requirements of subsequent research.
In order to verify the effectiveness of the GWR model, comparisons of the mean absolute error (MAE) of the results by the GWR model and traditional Kriging model in 2014 are provided. As shown in Fig. 5, between 150 and 237 days in 2014, the MAE values by Kriging model were mostly higher than those of GWR model. Specifically, the average MAE of GWR model is 0.7898, while that of Kriging model is 0.8415, which is 0.0517 higher. The comparative results demonstrate the utility of the GWR model to get spatially constant temperature maps.
To determine the optimal terms of the polynomial fitting function, comparative experiments based on the temperature data from the CMIC were conducted, which includes hourly temperature data from September 9, 2020 to September 20, 2020. Then, the MCD12Q1 land cover product is introduced to test the accuracy of the polynomial fitting model, and the impact of different land types to the fitting results. The final accuracy of polynomial fitting results is sorted out in Fig. 6. It can be seen that the fitting accuracy of different land types has little difference. The absolute value error basically falls between 0.5 K and 0.8 K, and the variance is within 0.4 K. It indicates that the underlying surface has little influence on the polynomial fitting accuracy, and the fitting results can meet the judgment of heat wave weather.
For better perception of the extracted dense temperature maps, the number of hours exceeding 35°C in each day of 2013 (horizontal axis) is further counted according to the basic unit of different streets in the study area (vertical axis). Due to space limitation, we only show the results from the 150th to 270th day of the year, that is, the days when the heat wave mainly occurs in the study area, as shown in Fig. 7. The results show that  there are multiple peaks of heat waves within a year, and there are large spatial differences among regions. The visualization results show the temporal and spatial differences of heat waves as a whole, which further reflects the effectiveness of dense temperature maps.

B. Risk Analysis on Heat Waves
To further analyze the heat wave risk at the street level, the number of heat wave days per street needs to be extracted first. According to the regulations of the China Meteorological Administration, when the local maximum temperature is equal to or higher than 35°C for 3 days, it can be called as a mild heat wave; when the maximum temperature is equal to or higher than 35°C for 5 days, it can be defined as a moderate heat wave; when the maximum temperature is equal to or higher than 38°C for 3 days, it is a severe heat wave. Therefore, we take streets as the basic unit and extract their number of heat wave days per year for the subsequent heat wave risk analysis. Taking the results of 2013 as an example, the distribution and degree of heat waves in different areas are shown in Fig. 8. It can be seen from the results that, moderate heat waves of the study area were occurred from the 183rd to 198th day of the year (about June), while severe heat waves were happened from the 204th to 230th day of the year, that is, about July. Based on the geographical location of the study area, we can see that the occurrence time of the moderate heat waves and severe heat waves correspond to the early summer and the midsummer of the study area, respectively, which further verifies the accuracy and reliability of the results.
Then, based on the number of heat wave days per year obtained by the above steps and combined with the population data, the HRI is calculated. Fig. 9 shows the HRI results for each street in the study area in 2011. It can be seen that the area along the river with a higher degree of urbanization usually has a higher HRI. For example, there is a higher risk of heat waves in the north of Yuyao City. In addition, the heat wave risk in the study area shows a centralized distribution, with the middle as the high point and gradually decreasing outwards.
Finally, spatial autocorrelation analysis on the HRI results from 2011 to 2014 in each street is conducted based on global Moran's I. The results are summarized in Table I. It can be seen that for each year from 2011 to 2014, Moran's I index is higher than 0.4, the z-value is higher than 1.96, and p-value is close to 0. This shows that the heat wave risk presents a strong positive correlation and aggregation distribution spatially. Combined with the above visualization results, we can find that the heat wave risks are clustered in the north of Yuyao City, the northeast of Hangzhou City, and the west of Jiaxing City.

V. CONCLUSION
In this article, in order to obtain more intensive temperature data to meet the needs of large-scale and time-series research, a joint spatio-temporal method is proposed for dense temperature mapping. Specifically, a GWR model is adopted to invert spatial continuous temperature maps from multisource remote sensing data, while a polynomial fitting model is applied to fill in the missing data at temporal-wise. With the aim to evaluate the heat wave risk of different areas, a heat wave risk model is built based on the dense temperature maps and population data. The accuracy evaluations have verified the effectiveness of the proposed joint spatio-temporal method for dense temperature mapping. And the case study on the study area, four cities of Zhejiang Province, China, have also demonstrated the utility of the heat wave risk model. The experimental results and analysis also found that the areas with higher degree of urbanization often have higher heat wave risks, which also show strong spatial autocorrelation characteristics.
Mengxi Liu (Graduate Student Member, IEEE) received the B.S. degree in geographic information science in 2019 from Sun Yat-sen University, Guangzhou, China, where she is currently working toward the Ph.D. degree in cartography and geographic information system with the School of Geography and Planning.
Her research interests include intelligent understanding of remote sensing images, change detection, and domain adaptation. Xuezhang Li received the B.S. degree in geographic information science and the M.S. degree in cartography and geographic information system from Sun Yat-sen University, Guangzhou, China, in 2020 and 2022, respectively.
His research interests include machine learning and deep learning in phenology and urbanization.
Zhuoqun Chai is currently working toward the B.S. degree in geographic information science with the School of Geography and Planning, Sun Yat-sen University, Guangzhou, China.
His research interests include remote sensing image classification and change detection with deep learning.