*2.2. Equipment*

Remote sensing data was collected using a DJI Matrice 100 (M100) quadcopter UAV. The M100 was deployed at our case study locations with three types of camera payloads—visual, multispectral and infrared. These cameras include the DJI Zenmuse X3 visual (12 MP), Zenmuse X3 multispectral (Blue-Green-NIR 680–800 nm at 12 MP) and DJI Zenmuse XTR radiometric thermal (13 mm, 30 Hz and spectral bandwidth of 7–13 µm). Additionally, ground temperatures were validated using a Nubee NUB8380 Digital Infrared Thermometer.

#### *2.3. Data Collection Methods*

Two datasets were collected during the 2018 calendar year: (1) surface temperature measured at 12:00 PM across the entire year and (2) surface temperature measured on a diurnal cycle. To evaluate surface temperature across the entire year, fourteen flights in Milwaukee and one in El Paso were recorded between 26 February and 13 September 2018 (Table 2). To evaluate the diurnal cycle of temperature, three flights in Milwaukee and one in El Paso measured temperature throughout the day at 9:00 AM, 12:00 PM, 3:00 PM and 5:00 PM (Table 3). Weather data was collected at Marquette from a station on top of Engineering Hall and weather data at UTEP was collected from a weather station 10.5 km away at El Paso International Airport. Each station recorded air temperature, relative humidity, wind speed, wind direction, relative humidity, solar radiation and atmospheric pressure. Drone imagery was captured autonomously using a third-party photogrammetry software called Pix4Dcapture. Using this software, autonomous flight paths were programmed to the drone prior to each mission. Programmed flight path information included drone speed, altitude and image overlap. Drone speed was set at 54 km/h for visual and multispectral flights but set at a lower threshold of 30.6 km/h for thermal flights due to the difference in image capture speed between the two camera technologies. The flight altitude for each mission was set to the United States Federal Aviation Administration (FAA) maximum allowable limit of 120 m, which resulted in thermal imagery at a 13 cm pixel size. Finally, the image overlap was set to 85%, which provided reliable overlap for stitching an orthomosaic during data processing.

**Table 2.** Flight log and summary of meteorological variables recorded for Marquette and UTEP during fifteen noon flights.


**Table 3.** Flight log and summary of meteorological variables recorded for Marquette and UTEP during four diurnal flights.


#### *2.4. Thermal Data Processing*

After data collection in the field, a series of post-processing steps were performed using Pix4D and ESRI's ArcMap to stitch the drone thermal imagery into orthomosaics, correct temperature values for emissivity and extract surface temperature data for analysis. First, Pix4D was used to stitch the captured thermal images into orthomosaics, export the orthomosaics as a 32-bit TIFF and georeference them for application within ArcMap.

Once in ArcMap, an emissivity correction was applied to each thermal orthomosaic. Emissivity is a measure of how well a material can emit energy as thermal radiation and different materials have different values of emissivity depending on their surface properties [24]. Land use classifications that were previously delineated for each case study area were used to apply emissivity values to the target surfaces. The emissivity values for each land use classification used in this study are listed in Table 4 and are based upon a review of emissivity studies. These emissivity values were then applied in the following emissivity correction equation derived from Stefan-Boltzmann Law:

$$T\_{target} = \sqrt[4]{\frac{T\_{sensor}^4 - (1 - \varepsilon) \ast T\_{background}^4}{\varepsilon}} \tag{1}$$

where *Ttarget* is the actual temperature of the target surface [K], *Tsensor* is the temperature measured by the infrared camera [K], *Tbackground* is the recorded air temperature [K] and ε is the emissivity value of the target surface [25]. This equation was used to correct each surface type for their respective emissivity before performing spatial data analysis.


**Table 4.** Emissivity values for each surface type.

Once the thermal data were corrected for emissivity, spatial data analysis was performed in ArcMap. First, a land use feature map was created that categorized the surface types in each case study location. Then inconsistencies within these areas, such as a parked car within a parking lot, human traffic on a sidewalk or construction materials on the street, were clipped and removed for each flight. Once these inconsistencies were removed, zonal statistics was applied to compute summary statistics of each surface type such as mean and standard deviation of the land surface temperature. A complete flow-chart of the process from flight programming to developing summary statistics is shown in Figure 2. In total this process took about 3 h to complete for each flight.

**Figure 2.** Flow chart of data collection and processing.

#### *2.5. Surface Parameters*

In addition to surface temperature, three other material properties were derived from visual and multispectral imagery, converted into spatial distribution rasters and averaged for each surface type. These include albedo (*S*), normalized difference vegetation index (*NDVI*) and apparent thermal inertia (*ATI*). Albedo, a measure of solar reflectance of a material, was derived from blue, green, red and near-IR image bands as shown in the following equation:

$$S = c\_b b\_k + c\_g g\_k + c\_r r\_k + c\_i i\_k \tag{2}$$

where *c<sup>b</sup>* = 0.17, *c<sup>g</sup>* = −0.13, *c<sup>r</sup>* = 0.33 and *c<sup>i</sup>* = 0.54 are derived constants and *b<sup>k</sup>* , *g<sup>k</sup>* , *r<sup>k</sup>* and *i<sup>k</sup>* are the band reflectance's for—blue, *b<sup>k</sup>* (420–492 nm); green, *g<sup>k</sup>* (533–587 nm); red, *r<sup>k</sup>* (604–664 nm); and near-IR, *i<sup>k</sup>* (833–920 nm) [34].

Visual and multispectral imagery were also used to derive *NDVI*, which is a measure of the degree of live vegetation and is commonly used to evaluate soil moisture dynamics, erosion potential and plant and crop health. As shown in Equation (3), *NDVI* is a function of near-IR and red band reflectance and is estimated on a scale of −1 to +1, with higher values indicating higher vegetative cover and greater plant health [35].

$$NDVI = \frac{(NIR - Red)}{(NIR + Red)}\tag{3}$$

Finally, *ATI* was derived for each surface type from albedo (*S*), solar correction (*SCR*) and the diurnal temperature amplitude (*DTA*) (Equation (4)). *ATI* is an estimation of thermal inertia from remotely sensed observations and can be estimated from diurnal changes in temperature. Specifically, *ATI* is derived from solar correction (*SCR*), albedo (*S*) and the diurnal temperature amplitude (*DTA*), where *DTA* is the difference between the maximum and minimum surface temperature recorded at the time the remote images were captured and *SCR* is the solar correction factor (Equation (5)), which is dependent on geographic location, the local latitude (θ) and the solar declination (ϕ) [36].

$$ATI = \frac{\text{SCR}(1-\text{S})}{\text{DTA}} \tag{4}$$

$$\text{SCR} = \sin\theta\sin\phi(1 - \left(\tan\theta\tan\phi\right)^2) + \cos\theta\cos\phi\,\arccos(-\tan\theta\tan\phi) \tag{5}$$

#### *2.6. Model Development*

Drone observations were applied to develop empirical models of land surface temperature. These include (1) a regression model to predict spatially averaged surface temperatures at 12:00 PM based upon meteorological variables and (2) a model to assess diurnal variability and predict surface temperatures throughout a given day.

Multi-variable regression models were developed to predict spatially averaged surface temperature of the fourteen Milwaukee and single El Paso 12:00 PM flights using MATLAB and the statistical software package JMP 13 [37]. Response screening was performed for each of the respective datasets to identify the strength of relationship between surface temperature (response) and meteorological parameters (predictors). Between the two case study locations, six surface types that were common to both locations were used as response variables: grass, canopy cover, concrete parking lot, concrete sidewalk, composite rooftop and road surface. Meteorological predictor variables included air temperature, relative humidity, preceding 24 h rainfall, wind speed, wind direction, barometric pressure and solar radiation. After response screening, stepwise linear regression was then performed to predict land surface temperature based upon meteorological parameters as represented in following equation:

$$y = \beta\_0 + \beta\_1 \mathbf{x}\_1 + \beta\_2 \mathbf{x}\_2 + \dots + \beta\_k \mathbf{x}\_k \tag{6}$$

where *y* is the response variable, β0, β1, . . . , β*<sup>k</sup>* are the regression coefficients and *x*1, *x*2, . . . , *x<sup>k</sup>* are the predictor variables for *k* predictors [38]. These models were developed using data from both the Milwaukee and El Paso flights; therefore, to evaluate the influence and leverage of the El Paso dataset we computed Cook's D influence and hat matrix leverage statistics [38].

Finally, we explored the variation in surface temperatures as they change throughout the day (9 AM, 12 PM, 3 PM and 5 PM) and evaluated if this variation could be explained by any meteorological parameters. In addition to exploring diurnal changes in variability, we applied the data to develop a model to predict land surface temperatures throughout the day for the six land use types common to each location. To do so, we applied the drone data collected on the four diurnal flight missions to estimate land surface temperatures based upon the solar radiation and the difference between the air and land surface temperatures, which have been found to be statistically significant predictors for diurnal estimates of pavement temperatures [39].

First, we computed a parameter (*g*) based upon the drone-derived mean land surface temperature and measured air temperature and solar radiation:

$$\mathbf{g} = \left(\overline{T\_s} - T\_a\right) \mathbf{\*} \mathbf{S} \tag{7}$$

where *T<sup>s</sup>* is the mean surface temperature of the land use, *T<sup>a</sup>* is the measured air temperature and *S* is the measured solar radiation (kW). Next, *g* at a given hour *i* was estimated using a Gaussian peak model given by the following:

$$g\_i = a \ast e^{-0.5 \ast \left(\frac{i-b}{c}\right)^2} \tag{8}$$

where *g<sup>i</sup>* , is the parameter *g* at hour *i*, *a* is the peak value, *b* is the critical point and *c* is the growth rate [40]. Using this model, the mean land surface temperature can be predicted based upon air temperature and solar radiation for any time of day using the following:

$$T\_{s,i} = T\_{a,i} + (g\_i/S\_i) \tag{9}$$

where *Ts,i* is the estimated surface temperature at hour *i* and *Ta,i* and *S<sup>i</sup>* are the air temperature and solar radiation at hour *i*. Taken as a whole, these models test both the suitability of predicting drone-derived mean land surface temperatures based upon meteorological variables, as well as the generalizability of our findings by including data from sites in two different geomorphologic and climatic regions.

#### **3. Results**

#### *3.1. Surface Temperature Variability*

We evaluated the land surface temperature variability of each flight across common land use types and generally found that green surfaces had a greater degree of variability than gray surfaces, with the exception being the rubber rooftop. As an example, the distribution of surface temperature data (1,986,543 total data points) is shown in Figure 3 for a flight recorded on July 11, 2018. The six gray surfaces recorded a smaller distribution of temperature on average but had more extreme values than green surfaces (Figure 3a). Gray surfaces retain more heat from the sun because of their high emissivity and ATI and therefore typically have higher surface temperatures. Additionally, non-normal behavior was identified for both canopy cover and rubber rooftop (Figure 3b). Canopy cover exhibits a left skew while the rubber rooftop exhibits a right skew. The canopy cover had a variation of tree types and therefore a variation of leaf area indices (LAI), which may be a reason for the skew in the temperature data. The rubber rooftop also exhibited a strong right skew, which may be due to small materials on the roof surface, such as ventilation pipes and drainage grates, that were difficult to detect and may not have been removed from the dataset. Therefore, this caused a distribution of lower temperatures to be recorded.

௦, = , + (

௦, ,

/ )

**Figure 3.** Boxplot distribution of surface temperature (**a**); and histogram of surface temperature (**b**). Data from flight recorded on July 11, 2018. Note GRS = grass; SM = shrub/mulch; CPY = canopy; PL = parking lot; SW = sidewalk; RTC = composite rooftop; RTR = rubber rooftop; RD = road; SLR = solar.

We then summarized the average temperature, standard deviation and coefficient of variation of each surface type for the fourteen recorded flights in Milwaukee, WI and the single flight in El Paso, TX (Table 5). Generally, gray surfaces exhibited higher temperatures throughout the year than green surfaces. In El Paso, the asphalt parking lot exhibited the highest average temperature (51.7 ◦C) and grass exhibited the lowest (41.6 ◦C), while in Milwaukee the black rubber rooftop exhibited the highest average temperature (57.4 ◦C) and canopy cover exhibited the lowest (30.4 ◦C). In terms of variation, the lowest degrees of variation typically occurred in the parking lots and grass. However, there is a noted difference in the variation between the two locations; the road in Milwaukee had the highest coefficient of variation of 0.32, while the road in El Paso had the lowest at 0.04. This may be due to a difference in traffic on the days that flights were conducted. The location in Milwaukee is located near the city center and is subject to heavy and constant vehicular traffic, while the location in El Paso is in a restricted traffic area and experienced very low vehicle activity on the weekend that the flight was conducted.


**Table 5.** Average temperature, standard deviation and coefficient of variation of nine surface types From 14 recorded flights on Marquette University campus and from one flight recorded on UTEP's campus.

The distribution of the average temperature, standard deviation and coefficient of variation for the fourteen Milwaukee, WI flights is further illustrated in Figure 4. The shrub/mulch, composite rooftop and solar panels have the most consistent variability among the land use types as shown in the boxplot distribution of their coefficient of variation, while the greatest spread in variation occurred in the road and sidewalk. This may indicate that areas that are not subject to human traffic (e.g., shrub/mulch flower beds, rooftops and solar panels) have more consistent variability in their temperatures, while other areas that are subject to intermittent human traffic (e.g., roads and sidewalks) have inconsistent temperature variabilities. We also evaluated if the variability in land surface temperature correlated with any meteorological parameters but found no statistically significant predictors.

**Figure 4.** Boxplot distribution of average temperature, standard deviation and coefficient of variation from 14 recorded flights in Milwaukee, WI.

We evaluated the variation in surface temperatures throughout the day and found that the highest degree of variation occurred at noon. This is demonstrated in the Figure 5, which shows box plots of the standard deviation for six land use types: grass, canopy, parking lot, sidewalk, composite roof and road for data from both MU and UTEP. As illustrated, all land use types have the greatest standard deviation in temperatures during 12:00 PM, with lower levels of deviation in the morning and late afternoon. This trend suggests that as surfaces heat up, they do so at different rates, which contributes to more variability during mid-day.

#### *3.2. Impact of the Built Environment*

0

5 10 15 20 Standard Deviation (°C) Grass Canopy Parking Lot Sidewalk Rooftop Road We also evaluated the spatial distribution of surface temperature to locate and identify factors of the built environment that contribute to temperature variability. Figure 5 illustrates the spatial distribution of surface temperatures for a flight on July 8th, 2018. One factor of variability is the reflectance and shaded cover from nearby buildings. For example, sidewalks in close proximity to Engineering Hall exhibited higher temperatures, most likely due to the sun's reflectance off its glass paneling. Two similarly sized sidewalk areas were compared and results show the average temperature was 4.7 ◦C hotter for the location closer to the building than the one farther away. In comparison to the sidewalk, parking lot land uses had more consistent variability, perhaps because there were fewer nearby buildings or large trees to exacerbate (glass reflectance) or reduce (shaded cover) their temperature. This indicates proximity to nearby buildings or other structures can be a significant factor of uncertainty in predicting surface temperatures.

Time (hr)

20 40 60

0.2 0.4 0.6 0.8

Grass Shrub/mulch Canopy Parking Lot Sidewalk RT Composite RT Rubber Road Solar

**Standard Deviation**

Grass Shrub/mulch Canopy Parking Lot Sidewalk RT Composite RT Rubber Road Solar

**Coefficient of Variation**

Grass Shrub/mulch Canopy Parking Lot Sidewalk RT Composite RT Rubber Road Solar

**Average Temperature**

**Figure 5.** Standard deviation distributions for six land use types at hour 9, 12, 15 and 17 at both the MU and UTEP locations. °

Other sources of land surface temperature uncertainty are traffic and parked cars. Traffic flow along a roadway intermittently blocks the suns radiation, thereby impacting the surface temperatures of the roadway pavement below. This creates a concentrated pocket of cooler surface temperatures called a heat shadow, which results in variations in surface temperatures across the pavement. This is especially pronounced in pavement lots with parked cars as illustrated in Figure 6b, which shows the distribution of surface temperatures within a parking lot. In this figure a parked car rooftop, pavement surface and heat shadow recorded temperatures of 69.6 ◦C, 47.8 ◦C and 49.0 ◦C, respectively, all within a space of ~50 m<sup>2</sup> .

**Figure 6.** Spatial distribution of temperature from a flight recorded on 11 July 2018 (**a**) and zoomed in spatial distribution of temperature for the concrete parking lot from a flight recorded on 11 July 2018 (**b**). The hotter surfaces (red) in the right image are parked cars and the cooler surfaces (blue) are heat shadows visible after parked cars leave.
