Assessing the Link between Human Modiﬁcation and Changes in Land Surface Temperature in Hainan, China Using Image Archives from Google Earth Engine

: In many areas of the world, population growth and land development have increased demand for land and other natural resources. Coastal areas are particularly susceptible since they are conducive for marine transportation, energy production, aquaculture, marine tourism and other activities. Anthropogenic activities in the coastal areas have triggered unprecedented land use change, depletion of coastal wetlands, loss of biodiversity, and degradation of other vital ecosystem services. The changes can be particularly drastic for small coastal islands with rich biodiversity. In this study, the inﬂuence of human modiﬁcation on land surface temperature (LST) for the coastal island Hainan in Southern China was investigated. We hypothesize that for this island, footprints of human activities are linked to the variation of land surface temperature, which could indicate environmental degradation. To test this hypothesis, we estimated LST changes between 2000 and 2016 and computed the spatio-temporal correlation between LST and human modiﬁcation. Speciﬁcally, we classiﬁed temperature data for the four years 2000, 2006, 2012 and 2016 into 5 temperature zones based on their respective mean and standard deviation values. We then assessed the correlation between each temperature zone and a human modiﬁcation index computed for the year 2016. Apart from this, we estimated mean, maximum and the standard deviation of annual temperature for each pixel in the 17 years to assess the links with human modiﬁcation. The results showed that: (1) The mean LST temperature in Hainan Island increased with ﬂuctuations from 2000 to 2016. (2) The moderate temperature zones were dominant in the island during the four years included in this study. (3) A strong positive correlation of 0.72 between human modiﬁcation index and mean and maximum LST temperature indicated a potential link between human modiﬁcation and mean and maximum LST temperatures over the 17 years of analysis. (4) The mean value of human modiﬁcation index in the temperature zones in 2016 showed a progressive rise with 0.24 in the low temperature zone, 0.33 in the secondary moderate, 0.45 in the moderate, 0.54 in the secondary high and 0.61 in the high temperature zones. This work highlighted the potential value of using large and multi-temporal earth observation datasets from cloud platforms to assess the inﬂuence of human activities in sensitive ecosystems. The results could contribute to the development of sustainable management and coastal ecosystems conservation plans.


Introduction
Humans are consistently modifying land surfaces through expansion of settlements, intensification of agriculture, and construction of infrastructure for transport, energy development, and exploration of natural resources. Anthropogenic modifications contribute to ecosystem changes, including vegetation loss, degradation of biodiversity, change of land surface and atmospheric temperatures, and pollution among other losses of ecosystem services. Coastal areas, particularly small offshore islands are among the vulnerable ecosystems [1][2][3][4] and can thus be easily decimated by uncontrolled human activities. It is therefore necessary to analyze, quantify and to mitigate the influence of human modification on the integrity of critical coastal ecosystems. Studies on land surface temperature have mainly focused on understanding urban heat island (UHI), a phenomenon first illustrated in 1818, reinforcing the knowledge that urban areas are characterized by higher air and surface temperature when compared with the temperatures in rural areas [1,[5][6][7][8][9][10]. The UHI phenomenon affects not only the health of urban dwellers, the livability of urban areas, but also the climatic and ecological processes [11][12][13]. Anthropogenic heat release, surface cover, climatic conditions, air pollution and other factors are considered to be numerous causes of UHI [12,14]. Surface UHI is characterized as the relative difference between the land surface temperature (LST) in urban areas and the temperature of surrounding areas and is influenced by urban development [12,[15][16][17][18]. UHI is generally linked to urban development and land cover transformations associated with urban land uses [19]. Land cover transformation in urban areas include the conversion of natural vegetation and agriculture lands to impervious surfaces [19]. While studying the link between land surface temperature and vegetation, the Normalized Difference Index (NDVI), has commonly been used as a proxy for vegetation abundance [20]. Understanding the contribution of human modification on changes of land surface temperature and on UHI remains an open research question [21].
Land surface temperatures (LST) vary strongly across the globe. Unusual changes in LST-for instance within one or two decades-can usually be credited to heightened degree human activities. Studies on the influence of human activities on LST have focused on the contribution of land use/land cover, human settlement, urbanization, and social economic factors. Sahana et al. investigated the temperature response to land use/land cover using Landsat 5 TM and Landsat 8 OLI data and observed that increase in non-evaporating surfaces and decrease in vegetation had increased the surface temperature in Sundarban Biosphere Reserve, India [22]. Moreover, artificial urban surfaces have been found to have a positive exponential relationship with LST [23]. Quantitative evidence of the effects of land use change on the LST has also been documented [24]. Xiao and Weng studied the impact of land use and land cover changes on land surface temperature for the Karst area in China and concluded that LST changes were mainly associated with changes in construction materials [25]. There are also other studies that focus on the impacts of human-induced land cover changes on the dynamics of surface parameters, such as impacts on surface roughness and fractional vegetation coverage in East Asia [26]. Fu and Weng investigated the influence of urbanization induced land use/land cover change on LST in the Atlanta metropolitan area in the US [27]. Additional studies have analyzed the mutual relationships and linkages between land surface temperature, UHI, land cover characteristics, social-ecological and social-economic variables in different cities [28,29]. Regional relationships between surface temperature, vegetation and human settlement was also assessed in the Phoenix, USA, which is a very rapidly urbanizing area [30].
From the ensuing research, it is apparent that it is not possible to conclusively understand the complexity of the influence of human activities on landscapes from limited variables, such as from the extent of the natural and converted land. It is important to account for the cumulative impacts from multiple human activities both at the local level and at the global scale. However, information on cumulative effects of human modifications are usually sparse. This is particularly true when studying the influence of human activities in coastal islands. Technological advances in the fields of earth observation and computation, coupled with the availability of other rich archives of data, provide an enviable resource for assessing and evaluating the extent of human modification and Remote Sens. 2020, 12, 888 3 of 17 the potential influence of this modification on the environmental and climatic characteristics of our planet. In this study, archived data on human modification and land surface temperatures were adopted to evaluate the spatial and temporal link between human activities and changes in land surface temperatures. The Global Human Modification (GHM) dataset provides the cumulative human modification of terrestrial lands and their estimated impacts using global datasets for the median year of 2016 [31]. For the medium-term changes in land surface temperatures, we adopted MODIS 8-day Land Surface Temperature data. Both data sets are archived in Google Earth Engine (GEE). As GEE is a cloud-based platform for geospatial analysis which consists of multi-source data catalog with massive computational capabilities to bear on a variety of societal and environmental issues, it was adopted as the tool of data exporting and analyzing [32]. To pursue the objectives of this study, Hainan Island in China was selected as the study site. Being the most popular tourism and in-migration destination from other province and the Special Economic Zone in China [33,34], the development of Hainan Island has attracted a variety of human activities in the past decades. With the large influx of residents and tourists and the high-speed of urbanization, the anthropogenic activities have been linked to the thermal environment causing elevated temperature, loss of the valuable, temperature-sensitive onshore habitats (coral reefs, mangroves), and ecosystem vulnerability [33,[35][36][37]. Making Hainan Island, a candidate site for investigating the spatio-temporal links between human activities and LST in Hainan.

Study Area
Hainan is the second largest Chinese island in the South China Sea and is approximately 24 km away from the Chinese mainland, separated by the Qiongzhou Strait ( Figure 1). The island is located the southernmost province of China within latitudes 18 • 10 -20 • 10 and longitudes 108 • 37 -111 • 03 E. The spatial extent of the island covers an area of approximately 35,354 km 2 and is home to a population of approximately 9.34 million people in 2018 from the Hainan Statistics Yearbook [38]. The island is characterized by tropical climate throughout the whole year. Hainan, being an important tourist destination and having been designated as China's largest Special Economic Zone since 1988 [33], underwent rapid economic growth that was mainly triggered by tourism. The developments in the island have come with high-speed urbanization, rapid land-use changes and a large influx of residents from the mainland China [33]. The number of tourists from China and other countries visiting one single city, Sanya, was increasing from 2.4 million in 2001 to 5.38 million in 2007 [34] and the total tourists on the whole island reached 76.27 million in 2018 from Hainan Statistics Yearbook [38]. Consequently, natural resources in the island have been explored and exploited extensively to support the vibrant tourism and other sectors. This has particularly been encouraged by the fact that Hainan and its surrounding coastal ecosystem traditionally encompass a high diversity of natural habitat, including coral reef, mangroves, and various types of shores for tourism [35]. Hainan hosts one third of the total mangrove forest areas of China [36], hence, is an important inter-tidal estuarine ecosystem that is currently very vulnerable to human influence. The rapid increase of economic activities in the area has already led to the farmland loss, conversion of forests to orchard, conversion of farmlands to artificial surfaces, pollution, and degradation onshore habitats [37]. The reduction of plants and psammolittoral organisms (coral reefs, mangroves, and seaweed) will further cause ecosystem vulnerability and exposure to further degradation, hindering sustainable local natural resource management and development. Assessing and documenting the impact of human activities on environmental and climate sustainability is therefore vital for sustainable development policy for the local ecosystem. Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 17

Data source
The main data sets in this work included MODIS Land Surface Temperature (LST) products and the Global Human Modification layer, both of which are available in Google Earth Engine (GEE). In addition, we used administration data for referencing. The main datasets are outlined in Table 1.

Data Source
The main data sets in this work included MODIS Land Surface Temperature (LST) products and the Global Human Modification layer, both of which are available in Google Earth Engine (GEE). In addition, we used administration data for referencing. The main datasets are outlined in Table 1. The MODIS 8-day temperature products in GEE were used to assess the temperature change between 2000 and 2016. This time scale was selected to coincide with the GHM data, which was produced from the year 2016. For the human modification, we used the Human Modification index (HM) as a proxy for human footprints from the GHM dataset provided by GEE. The GHM dataset is developed by the scientists at the Nature Conservancy and Conservation Science Partners and is available on GEE platform. HM provides a cumulative measure of global human modification of terrestrial lands at 1 km spatial resolution. In this dataset, the global intensity of modification is based on modeling the physical extents of 13 anthropogenic stressors and their estimated impacts using spatially explicit global datasets for the year 2016. The HM values range from 0 to 1 according to the degree of land modification, with 0 representing the no evidence of human modification and 1 representing the highest degree of modification. HM is a composite index that is computed from five major categories of thirteen environmental stressors. The five categories of stressors are: (a) human settlement (population density, built-up areas), (b) agriculture (cropland, livestock), (c) transportation (major roads, minor roads, two tracks, railroads), (d) mining and energy production (mining, oil wells, wind turbines), and (e) electrical infrastructure (powerlines, nighttime lights) [33].

Retrieval of Images from GEE
Data for GHM and land surface temperature (LST) were retrieved from GEE. The main methodological steps in this study are outlined in Figure 2. The data pre-processing step involved using the outlines of the area of study to filter out the MODIS LST data for the area of study from 2000 to 2016, calculation of annual LST statistics from the series of 8-Day data, estimating changes in LST and extracting HM data for the area of study. Specifically, a region of interest (ROI) was built from the shapefiles of the boundaries of the Hainan Island. The region of interest was used to filter and to clip the datasets from GEE. As GHM data was only available for only one year (2016), the data was clipped to the ROI and exported for further analysis and visualization. The LST data that used in this study was derived from the MODIS Terra land surface temperature and emissivity 8-day global dataset in GEE. This data is available from 2000 and contains an average of 8-day land surface temperature in a 1200*1200 km grid at 1 km resolution. We used the LST_Day_1km band with 1 km resolution for the daytime LST analysis. Each pixel value is a simple average of all the corresponding pixels collected within that 8 days period. We aggregated the images for a single year to create a single raster representing the annual mean temperature for that particular year. We also calculated mean, maximum and standard deviation statistics from the 17 annual LST datasets. The 17 annual datasets represented individual raster for each year from 2000 to 2016. The annual images for the entire time scale were used for the comparison with the global human modification.
As the temperature units for the MODIS LST data is provided in Kelvin, we developed a function using the mathematical operations provided in GEE to retrieve annual mean LST images in degrees centigrade. The output from GEE was then exported to a GIS platform for further analysis. In order to estimate the temperature changes, we subtracted the mean temperature for the baseline year (2000) from the temperature of the year 2016. The annual temperature images were used to calculate and plot the mean annual LST changes and for the classification of the temperature zones. Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 17

Land surface temperature pattern classifications
Values of HM from the global human modification data ranged from 0 to 1 according to the proportion of landscape modification which is a fundamental indicator of ecological function [31,39,40]. In generating the HM data, the human modification model [41] was used to derive the potential magnitude of impact on terrestrial lands based on the 13 human stressors with 5 categories from the most recent datasets. The 13 stressors were categorized into 5 categories including: (a) human settlement (population density, built-up areas), (b) agriculture (cropland, livestock), (c) transportation (major roads, minor roads, two tracks, railroads), (d) mining and energy production (mining, oil wells, wind turbines), and (e) electrical infrastructure (powerlines and nighttime lights). Specifically, the 2015 UN-adjusted, Gridded Population of the World dataset [42] and the Global Human Settlements Layer [43] were used for the human settlement category acquisition; the Unified Cropland Layer [44] was used to estimate cropland, while livestock densities were identified from the Gridded Livestock of the World v2 database [45]. OpenStreetMap (OSM) was the source of the transportation and mining categories [46]. For electrical infrastructure, OSM was used to derive the above-ground powerlines and nighttime lights was from the recent version of Defense Meteorological Satellite nighttime lights [47]. Based on the above most recent datasets, the physical extent of 13 human stressors was mapped. Then the spatial extents were weighted by their respective intensity levels to maintain different land use activities. The result of HM is considered to be a comprehensive and continuously scaled spatial assessment ranging from 0 to 1 and represents the proportion of landscape modified by human activities. On the basis of the multiple advantages of HM, providing an update of extent where human activities have modified, using the most recent global datasets and more anthropogenic drivers for data sources, having a high impact on biodiversity [48], the dataset was used in this research as a proxy for human modification. For human modification classification, HM data was clipped by the ROI layer and categorized into four classes. Based on the global, non-normal, distribution of HM [31] and the manner consistent with literatures [31,[49][50], the HM was binned into four classes: low (0 ≤ HM ≤ 0.1), moderate (0.1 < HM ≤ 0.4), high (0.4 < HM ≤ 0.7), and very high (0.7 < HM ≤ 0.91). In this area, the no presence of stressor (value equal to 0) was merged into the low class.
In order to analyze the spatial patterns of the thermal environment from 2000 to 2016, we quantified the temperature variations during the research period by reclassifying the LST in the four years of 2000, 2006, 2012, and 2016. For this reclassification, we adopted the mean-standard deviation method which was used for thermal environmental research for the urban expansion [51], or land cover response to the LST in the Nanjing metropolitan region [52]. The mean-standard deviation method and the variables for reclassification are outlined in Table 2 to show how we classify the LST

Land Surface Temperature Pattern Classifications
Values of HM from the global human modification data ranged from 0 to 1 according to the proportion of landscape modification which is a fundamental indicator of ecological function [31,39,40]. In generating the HM data, the human modification model [41] was used to derive the potential magnitude of impact on terrestrial lands based on the 13 human stressors with 5 categories from the most recent datasets. The 13 stressors were categorized into 5 categories including: (a) human settlement (population density, built-up areas), (b) agriculture (cropland, livestock), (c) transportation (major roads, minor roads, two tracks, railroads), (d) mining and energy production (mining, oil wells, wind turbines), and (e) electrical infrastructure (powerlines and nighttime lights). Specifically, the 2015 UN-adjusted, Gridded Population of the World dataset [42] and the Global Human Settlements Layer [43] were used for the human settlement category acquisition; the Unified Cropland Layer [44] was used to estimate cropland, while livestock densities were identified from the Gridded Livestock of the World v2 database [45]. OpenStreetMap (OSM) was the source of the transportation and mining categories [46]. For electrical infrastructure, OSM was used to derive the above-ground powerlines and nighttime lights was from the recent version of Defense Meteorological Satellite nighttime lights [47]. Based on the above most recent datasets, the physical extent of 13 human stressors was mapped. Then the spatial extents were weighted by their respective intensity levels to maintain different land use activities. The result of HM is considered to be a comprehensive and continuously scaled spatial assessment ranging from 0 to 1 and represents the proportion of landscape modified by human activities. On the basis of the multiple advantages of HM, providing an update of extent where human activities have modified, using the most recent global datasets and more anthropogenic drivers for data sources, having a high impact on biodiversity [48], the dataset was used in this research as a proxy for human modification. For human modification classification, HM data was clipped by the ROI layer and categorized into four classes. Based on the global, non-normal, distribution of HM [31] and the manner consistent with literatures [31,49,50], the HM was binned into four classes: low (0 ≤ HM ≤ 0.1), moderate (0.1 < HM ≤ 0.4), high (0.4 < HM ≤ 0.7), and very high (0.7 < HM ≤ 0.91). In this area, the no presence of stressor (value equal to 0) was merged into the low class.
In order to analyze the spatial patterns of the thermal environment from 2000 to 2016, we quantified the temperature variations during the research period by reclassifying the LST in the four years of 2000, 2006, 2012, and 2016. For this reclassification, we adopted the mean-standard deviation method which was used for thermal environmental research for the urban expansion [51], or land cover response to the LST in the Nanjing metropolitan region [52]. The mean-standard deviation method and the variables for reclassification are outlined in Table 2 to show how we classify the LST images into five LST zones. There are two variables needed for the calculation in the mean-standard deviation method: the mean LST of the study area (µ) and the standard deviation of the LST (std). The mean-standard deviation method can reflect the spatial and temporal variation for the land surface temperature in detail through reclassifying the LST into five temperature zones based on the LST range of each zone: extreme low-temperature zone, low-temperature zone, moderate-temperature zone, high-temperature zone, extreme high-temperature zone ( Table 2). We applied the mean-standard deviation method into the mean LST images for the years 2000, 2006, 2012, and 2016 (Table 3), respectively, and reclassified the LST images into the five temperature zones to quantify the temporal and spatial variations for the thermal environment. For the year 2000, the mean-standard deviation method yielded mean LST of 28.35 • C and standard deviation of 2.01 • C. While for the year 2016, the mean-standard deviation method yielded mean LST of 28.73 • C and standard deviation of 1.85 • C. The classifications of five surface temperature zones for the LST in the four years were shown in Table 3. The LST in 2000 was subtracted from the LST in 2016 to calculate the LST changes and was classified afterwards. We used zero, the maximum temperature change, the minimum temperature change, and approximately 0.5 times of the maximum, respectively, to classify the temperature changes and to show the temperature change ranges over the past 17 years.

Spatio-Temporal Variation of Land Surface Temperature
Apart from assessing the pattern changes in temperature zones at the end points, that is in 2000, 2006, 2012 and in 2016 as has been described in Section 3.2, we estimated the mean, maximum and standard deviation of annual temperatures in each pixel over the 17 years. Specifically, we computed cell statistics for the 17 annual temperature raster files using ArcGIS. The results of mean, maximum and standard deviation of temperature for each pixel was used to compute the correlation between temperature statistics over the 17 years with the human modification in the study area.

Correlation of the LST Variations and the Human Modification
In order to analyze the relationship of human modification and LST variations, we exported two outputs. The first output was the distribution of human modification in each temperature zones in 2016 while the second output was the correlations with human modification and the temperature statistics (illustrated in Section 3.3) over the 17 years separately.
Random samples of data of human modification values in the five temperature zones as reclassified from the 2016 temperature data were extracted. From the sampled, we plotted box plots to highlight the distribution of human modification indices in each temperature zones. To create a representative sample for the correlation analysis, we randomly generated 1500 points, with approximately 300 points for each temperature zones and used the points to extract values from both the temperature and human modification data. We also used zonal statistics method in ArcGIS to calculate the statistics of human modification index in each temperature zone as at the year 2016. This provided the basis for the plots to visually demonstrate the variation of human modification in the various temperature zones. To furtherly interpret the human modification activities from the perspective of land use, we also adopted 30m resolution land use/land cover data [53] to assess the mean temperature and human modification distribution in each type of land use/land cover (Section 4.4) on Hainan Island.
Except for assessing the pattern changes in temperature zones at the four years has been described in Section 3.2, we estimated the mean, maximum and standard deviation of annual temperatures in each pixel over the 17 years in the analysis. Specifically, we developed a Python code in ArcGIS to compute cell statistics for the 17 annual temperature raster files. The results of mean, maximum and standard deviation of temperature for each pixel was used to compute the correlation between temperature statistics over the 17 years against the human modification in the study area.

The LST Pattern Changes of the Thermal Environment
From the analysis of temperature data, we analyzed the spatial (Figure 3a) and temporal annual mean temperature changes in Figure 3b. Figure Table 4 shows the statistics of the five temperature zones changing in the four years. Visually from Figure 4, the spatial pattern of distribution of the temperature zones in 2006, 2012 and 2012 are similar to the scenario in 2000 although the areas of temperature zones are either expanding or shrinking and the changing ranges are different. The moderate temperature areas were distributed widely in the island and the low temperature zones were mainly located in the south-west and central parts of Hainan. The high temperature zones were in the western and southern parts of the island, while the secondary high temperature areas were mainly located in the northern island for both four years. The moderate temperature zone dominated the whole Hainan in the four years, with the percentage of pixel in this zone being approximately 43.6%. 40.6%, 42.73% and 41.3% (Table 4) in the four years separately. The secondary high temperature is expanding in the four years and the Low temperature zone is increasing slightly. The other two zones are with fluctuations. Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 17  Table 4 shows the statistics of the five temperature zones changing in the four years. Visually from Figure 4, the spatial pattern of distribution of the temperature zones in 2006, 2012 and 2012 are similar to the scenario in 2000 although the areas of temperature zones are either expanding or shrinking and the changing ranges are different. The moderate temperature areas were distributed widely in the island and the low temperature zones were mainly located in the south-west and central parts of Hainan. The high temperature zones were in the western and southern parts of the island, while the secondary high temperature areas were mainly located in the northern island for both four years. The moderate temperature zone dominated the whole Hainan in the four years, with the percentage of pixel in this zone being approximately 43.6%. 40.6%, 42.73% and 41.3% (Table 4) in the four years separately. The secondary high temperature is expanding in the  four years and the Low temperature zone is increasing slightly. The other two zones are with fluctuations.

The tempo-spatial variations of the thermal environment and human modification classification
To visualize the relationship between human modification and land surface temperature, we mapped human modification against mean, maximum and standard deviation of temperature as computed from the 17 years of temperature data. Figure 5 presents maps of the statistics of the 17 yearly land surface temperature data and the geographic distribution of human modification index data in the study. Specifically, Figure 5(a) is a map of the mean temperature from the 17 years, Figure  5 (b) is the map of standard deviation while Figure 5 (c) is the map of the maximum temperature per pixel from the 17-yearly data. When visually compared against the human modification index (Figure  5d), we observed that for all the temperature statistics data, the northern edges of the island appeared

The Tempo-Spatial Variations of the Thermal Environment and Human Modification Classification
To visualize the relationship between human modification and land surface temperature, we mapped human modification against mean, maximum and standard deviation of temperature as computed from the 17 years of temperature data. Figure 5 presents maps of the statistics of the 17 yearly land surface temperature data and the geographic distribution of human modification index data in the study. Specifically, Figure 5a is a map of the mean temperature from the 17 years, Figure 5b is the map of standard deviation while Figure 5c is the map of the maximum temperature per pixel from the 17-yearly data. When visually compared against the human modification index (Figure 5d), we observed that for all the temperature statistics data, the northern edges of the island appeared to have relatively higher mean, maximum and standard deviation in the temperature data. A similar pattern was evident in the human modification data.

The Correlation of Human Modification with the LST Variations
To analyze the correlation between human modification and land surface temperature we used two outputs. Firstly, we extracted random samples of data of human modification values against the five temperature zones as reclassified from the 2016 temperature data (Figure 3b). From the sampled, we plotted box plots to highlight the distribution of human modification index in each temperature zones. In Figure 6, while the mean value of the human modification in the low temperature zone was 0.24, the mean of human modification in the high temperature zone was 0.61. With a progressive rise of the mean human modification value with 0.24 in the low, 0.33 in the secondary moderate, 0.45 in the moderate, 0.54 in the secondary high and 0.61 in the high temperature zones, thus Figure 6 shows a potential link between human modification and temperature zones as at the year 2016.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 17 to have relatively higher mean, maximum and standard deviation in the temperature data. A similar pattern was evident in the human modification data.

The correlation of human modification with the LST variations
To analyze the correlation between human modification and land surface temperature we used two outputs. Firstly, we extracted random samples of data of human modification values against the five temperature zones as reclassified from the 2016 temperature data (Figure 3d). From the sampled, we plotted box plots to highlight the distribution of human modification index in each temperature zones. In Figure 6, while the mean value of the human modification in the low temperature zone was 0.24, the mean of human modification in the high temperature zone was 0.61. With a progressive rise of the mean human modification value with 0.24 in the low, 0.33 in the secondary moderate, 0.45 in the moderate, 0.54 in the secondary high and 0.61 in the high temperature zones, thus Figure 6 shows a potential link between human modification and temperature zones as at the year 2016. Secondly, from the sampled points, the values from the mean, maximum and standard deviation raster files were extracted and compared with the values of human modification at the same sampled locations. Specifically, the values at 1800 random sampling points in the study area were used. We then plotted the mean, maximum and standard deviations of temperature against the human modification and computed Pearson's correlation in each case. Figure 7 presents the scatter plots of mean, maximum and standard deviation of temperature respectively against the human modification index. There was a weak positive linear correlation (0.33) between the standard deviation of temperature and the human modification index. For the mean and maximum temperature, there was a strong positive correlation of 0.72 with the human modification index, indicating a potential link Secondly, from the sampled points, the values from the mean, maximum and standard deviation raster files were extracted and compared with the values of human modification at the same sampled locations. Specifically, the values at 1800 random sampling points in the study area were used. We then plotted the mean, maximum and standard deviations of temperature against the human modification and computed Pearson's correlation in each case. Figure 7 presents the scatter plots of mean, maximum and standard deviation of temperature respectively against the human modification index. There was a weak positive linear correlation (0.33) between the standard deviation of temperature and the human modification index. For the mean and maximum temperature, there was a strong positive correlation of 0.72 with the human modification index, indicating a potential link between human modification and mean and maximum temperatures over the 17 years of analysis.
Secondly, from the sampled points, the values from the mean, maximum and standard deviation raster files were extracted and compared with the values of human modification at the same sampled locations. Specifically, the values at 1800 random sampling points in the study area were used. We then plotted the mean, maximum and standard deviations of temperature against the human modification and computed Pearson's correlation in each case. Figure 7 presents the scatter plots of mean, maximum and standard deviation of temperature respectively against the human modification index. There was a weak positive linear correlation (0.33) between the standard deviation of temperature and the human modification index. For the mean and maximum temperature, there was a strong positive correlation of 0.72 with the human modification index, indicating a potential link between human modification and mean and maximum temperatures over the 17 years of analysis.

Relationship between Land Cover Classes, Human Modification and Mean Temperature
Using the land cover characteristics in Hainan Island ( Figure 8) as zones, we extracted values from human modification index and from the mean temperature rasters for each of the land cover zones. We then plotted the variation of human modification and mean temperatures for each land cover zone.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 17 Using the land cover characteristics in Hainan Island ( Figure 8) as zones, we extracted values from human modification index and from the mean temperature rasters for each of the land cover zones. We then plotted the variation of human modification and mean temperatures for each land cover zone. When comparing human modification against the land use categories (Figure 9), impervious surfaces had the highest average human modification at about 0.65 while the forest land cover classes had the lowest average human modification index at 0.37. Shrubland, water bodies, grassland, wetlands, bareland and croplands had average human modification values of 0.47, 0.48, 0.5, 0.53, 0.55, and 0.56 respectively, indicating a potential link between land use characteristics and human When comparing human modification against the land use categories (Figure 9), impervious surfaces had the highest average human modification at about 0.65 while the forest land cover classes had the lowest average human modification index at 0.37. Shrubland, water bodies, grassland, wetlands, bareland and croplands had average human modification values of 0.47, 0.48, 0.5, 0.53, 0.55, and 0.56 respectively, indicating a potential link between land use characteristics and human modification index with forest cover classes having low human modification indices while impervious built-up areas having the highest human modification indices. In terms of mean temperature, bare surfaces had the highest mean temperature values at about 30.3 • C followed by impervious surfaces with mean temperature values of about 29.8 • C. The land cover class with the lowest relative mean temperature value was forest cover class at 27.6 • C. Showing the forested areas which were the least modified by humans were at least 2 • C cooler than the highly modified impervious surfaces and bareland.

Discussion and Conclusions
The focus of this work was on using archives of satellite-derived products to assess the potential link between human modification and land surface temperature variation in Hainan Island in China. For the land surface temperature, MODIS Terra Land Surface Temperature from 2000 to 2016 was used, while for human modification index from global human modification data as at 2016 was used. In analyzing the land surface temperatures, the differences in the per-pixel temperature between the baseline data from 2000 and in the follow up data from 2016 was considered. Also, the medium-term descriptive statistics, including mean, maximum and standard deviation for the 17 years of analysis was considered. Further, the correlation coefficients between human modification indices and the various temperature characteristics were calculated. Mean temperature and human modification were plotted against land cover classes to visualize the relation between temperature, human modification and land cover characteristics. This analysis showed that there was a strong positive correlation of about 0.72 between mean temperature (from the 17-yearly temperature) and human modification index. Similarly, there was a strong positive correlation (0.72) between maximum temperature and the human modification. The results from this work are consistent with other studies which have found positive correlation between human-induced land cover changes and land surface temperature [54]. Similarly, it has been recorded that there is a very high likelihood that human influences can explain some of the incidences

Discussion and Conclusions
The focus of this work was on using archives of satellite-derived products to assess the potential link between human modification and land surface temperature variation in Hainan Island in China. For the land surface temperature, MODIS Terra Land Surface Temperature from 2000 to 2016 was used, while for human modification index from global human modification data as at 2016 was used. In analyzing the land surface temperatures, the differences in the per-pixel temperature between the baseline data from 2000 and in the follow up data from 2016 was considered. Also, the medium-term descriptive statistics, including mean, maximum and standard deviation for the 17 years of analysis was considered. Further, the correlation coefficients between human modification indices and the various temperature characteristics were calculated. Mean temperature and human modification were plotted against land cover classes to visualize the relation between temperature, human modification and land cover characteristics. This analysis showed that there was a strong positive correlation of about 0.72 between mean temperature (from the 17-yearly temperature) and human modification index. Similarly, there was a strong positive correlation (0.72) between maximum temperature and the human modification. The results from this work are consistent with other studies which have found positive correlation between human-induced land cover changes and land surface temperature [54]. Similarly, it has been recorded that there is a very high likelihood that human influences can explain some of the incidences of heatwaves [55]. For the comparison of mean temperature, human modification and land cover characteristics, highest values of human modification were observed in the areas with impervious surfaces. Similarly, the two land cover classes with the highest mean temperature values were barelands and areas of impervious surfaces (Figure 9), reaffirming studies that have linked urbanization to increased land surface temperature and urban heat island.
One limitation of this work was the inadequate temporal coverage of the human modification data which was only available for the year 2016. This limited the evaluation of annual link of human modification against temperature changes. Having access to human modification data with longer coverage and higher temporal resolution may yield better results and may help in pinpointing particular human activities that contribute significantly to environmental change. Since the human modification index is an accumulative assessment, we used the mean, maximum, standard deviation of annual temperature in 17 years to assess the links between LST variations and human modification index, mitigating the influence of the lack of annual human modification data. In future, annual human modification should be considered for assess the annual link.
In recent years, earth observation datasets as well as derived information covering large spatial scales has become more easily accessible and usable through the creation of cloud platforms like GEE. This study demonstrates the potential contribution of multi-source earth observation datasets from cloud platforms to the development of sustainable management plans and policy for the management and conservation of coastal ecosystems. Therefore, future studies could focus on including other variables in analyzing the influence of human modification on coastal ecosystems. For instance, additional data may include annual changes in vegetation characteristics, land cover patterns, and the various policy and management characteristics and regulations for the sites under investigation. Secondly, the methods could be up-scaled and applied in detecting the correlation between human modification and temperature changes, thus highlighting hotspots of anthropogenic influence. Furthermore, the methods could be extended to provide evidence-based support for conservation and protection of sensitive ecosystems and habitats from human disturbance.
Author Contributions: L.C. and F.O. conceived the study, analyzed the data, visualized the results, and developed the manuscript. H.B. contributed to the analysis of the data, wrote, and edited the manuscript. T.B. reviewed, wrote, edited the manuscript, and supervised the work. All authors have read and agreed to the published version of the manuscript.

Funding:
The work is partially funded by the Austrian Science Fund (FWF) through the Doctoral College GIScience at the University of Salzburg (DK W1237-N23).