ThermStats: An R package for quantifying surface thermal heterogeneity in assessments of microclimates

Variation in temperature at a fine spatial scale creates critically important microclimates for many organisms. Quantifying thermal heterogeneity at this scale is challenging and, until recently, has been largely restricted to the use of dataloggers to record air temperature. Thermography is becoming an increasingly viable alternative. A single photo from a thermal imaging camera contains thousands of spatially explicit surface temperature measurements, making thermal cameras ideal for rapidly assessing temperature variation at fine scale. Here, we present an R package—ThermStats—for processing thermal images and other gridded temperature data. The package addresses current constraints on applying thermography in ecology, by speeding up and simplifying the extraction of data from thermal images, and by facilitating the calculation of different metrics of thermal heterogeneity. The metrics capture both the frequency distribution and spatial patterns of temperature, and the package functions are designed to accommodate different sampling strategies and data in either matrix or raster format. We demonstrate how ThermStats can be used to capture temperature variation at fine spatial scales in structurally complex habitat, such as tropical rainforest. Using thermal images collected in the field (~0.5 cm2 resolution), we found that thermal hetereogeneity varied little between primary and logged forest, but did vary with time of day. Comparing temperature extremes in a microclimate layer estimated from LIght Detection And Ranging data (LIDAR; 2,500 m2 resolution) , we found that both hot and cold extremes were hotter inside oil palm plantations than in the neighbouring forest. Our package simplifies the processing of thermal data, and our metrics capture key spatiotemporal temperature trends that underpin physiological, behavioural and demographic responses to environmental change. As such, ThermStats can advance a wide range of studies requiring fine‐scale surface temperature data for microclimate investigations.


| INTRODUC TI ON
Temperature variation over spatial scales of mm to m has a central role in the ecology of many species. For small organisms, this is the scale at which nearly all thermal experiences will occur (Potter, Arthur Woods, & Pincebourde, 2013). Mobile species utilize finescale thermal heterogeneity on a daily basis to avoid or exploit extremes of heat. Over longer time periods, climate at this scale ('microclimate') can also maximize fitness and thus influence the fine-scale distribution of less mobile species (Maclean, Suggitt, Wilson, Duffy, & Bennie, 2017). Species can persist within 'microrefugia' where climatic conditions are otherwise unsuitable, potentially maintaining higher overall biodiversity in thermally heterogenous habitats (Hannah et al., 2014). The same mechanisms underlying fine-scale thermoregulation and thermal microrefugia could temper species' exposure to global climate change (Scheffers, Evans, Williams, & Edwards, 2014;Suggitt et al., 2018), especially in structurally complex habitats such as tropical rainforests (Scheffers et al., 2017). To understand and predict species' responses to temperature change we must therefore be able to efficiently capture thermal heterogeneity at fine spatial scales, but there are challenges associated with acquiring and analysing temperature data at this resolution.
While dataloggers remain instrumental in the field of thermal ecology (Bramer et al., 2018), technological advances in recent years have made thermal cameras an increasingly affordable and practical addition (Faye, Rebaudo, Yánez-Cajo, Cauvy-Fraunié, & Dangles, 2016;Scheffers et al., 2017). With one click a single thermal image provides thousands of spatially explicit surface temperature measurements at the mm-cm scale. To make the most of the wealth of data provided by thermal imagery, ecologists need some understanding of infrared technology and its limitations, as well as simple tools that enable a clear and reproducible workflow from potentially thousands of raw images through to simple summary statistics relevant to the biological questions of interest.
Extracting data from thermal images is rarely straightforward.
Raw data can be difficult to access, the processing steps unclear and parameters such as emissivity can strongly influence the accuracy of measurements, but may not be well understood by the novice user (Bramer et al., 2018). The extraction of data from FLIR thermal cameras specifically (one of the most commonly used brands) can be achieved using the freely available FLIR Tools software (https :// www.flir.com/produ cts/flir-tools ; cf. Scheffers et al., 2017), but cannot be done in batch and the automatic conversion of raw values to temperature lacks transparency. The process was simplified by the development of the Thermimage package (Tattersall, 2017) in R (R Core Team, 2018). However, Thermimage still does not directly facilitate batch processing, nor does it calculate (or suggest) what metrics best capture spatiotemporal temperature variation.
We present the package ThermStats to extend the functionality of Thermimage to meet the needs of ecologists: allowing hundreds of images to be processed with ease, and summarizing data in terms of variables that matter to individual organisms.
The most appropriate metrics to quantify thermal heterogeneity will depend on the taxonomic group and research questions of interest. Temperature varies across time and space in many ways that can easily be captured by thermal images; it is important to exploit the information provided without becoming overwhelmed. Both Shi, Wen, Paull, and Guo (2016) and Faye et al. (2016) provide a useful summary of some important metrics, and Faye et al. (2016) introduce a spatial component by borrowing metrics from landscape ecology, such as Shape Index and Cohesion Index (McGarigal, Cushman, & Ene, 2012).
Extending this approach reveals other relevant techniques, such as hot spot analysis (Getis & Ord, 1996) and connectivity (McGuire, Lawler, McRae, & Theobald, 2016). Together this suite of metrics provides a powerful way for ecologists to characterize temperature variation.
ThermStats integrates ideas, techniques and metrics from previous work into one simple framework for quantifying thermal hetereogeneity. The package alleviates two key constraints on the application of thermography in ecology: 1. Automates and speeds up the extraction and conversion of raw data from (FLIR) thermal images.
2. Provides a variety of tools to calculate metrics of thermal heterogeneity using any gridded temperature data.
We illustrate the utility of our package using temperature data for the island of Borneo, including thermal images collected in the field using a FLIR thermal camera, and a coarser resolution microclimate layer estimated using LIDAR data (Jucker et al., 2018).

| Step 1: thermal image collection
Surface temperature measurements at high spatial resolution can easily be collected in the field using a handheld thermal camera.
We used a FLIR model E40 camera, which costs ~US$4,000, weighs 825 g, and takes 19,200 measurements (160 × 120 pixels) in a single photo (FLIR, 2016;Scheffers et al., 2017). Various other models exist, including the smaller and more affordable FLIR ONE smartphone attachment at ~US$300, 34.5 g and a resolution of 80 × 60 pixels. While FLIR cameras appear to be the most widely used it should be noted that other brands are available, such as Optris and Testo. Currently ThermStats cannot directly extract data from other camera models, but readers are welcome to request enhancements via GitHub (https ://github.com/rasen ior/Therm Stats/ issues).
As with any field study, the sampling design should aim to achieve sufficient coverage over the study area and over time, such that the images are representative samples of the treatments of interest. For example a single image of the ground from 1 m away encompasses an area of 0.9 × 1.1 m (1 pixel ≅ 0.5 cm 2 ) using a FLIR E40 camera (FLIR, 2016), and so it may be necessary to take multiple photos in different cardinal directions and at different times of day to effectively represent the temperature of a study plot (Scheffers et al., 2017; Senior, Hill, Benedick, & Edwards, 2018) ( Figure 1).
Before any data are collected, we recommend users familiarize themselves with the technology (for a thorough review, see Tattersall, 2016). There are various sources of the infrared radiation detected by a thermal camera, but we want to focus only on the radiation emitted by the object or scene of interest, which is a function of its temperature. The amount of radiation emitted by a particular object, for a given temperature, depends on its emissivity. A perfect blackbody has an emissivity of 1, whereas surfaces that an ecologist is likely to photograph typically have an emissivity from 0.92 (for dry, bare soil; FLIR, 2016) to 0.99 (for green broadleaf forest; Snyder, Wan, Zhang, & Feng, 1998). Additionally, the temperature and relative humidity of the atmosphere and the distance between the object and the camera will affect (a) the amount of emitted radiation that is absorbed by the atmosphere and (b) the amount of background radiation, with some of this also being reflected by the object (reflected apparent temperature).
To accurately quantify surface temperature, environmental parameters (emissivity, reflected apparent temperature, atmospheric temperature, atmospheric relative humidity and object distance) can be set in the camera or defined during data processing (see 'Step 3: Conversion of raw data'). In the latter approach the user can measure atmospheric temperature and relative humidity concurrently with thermal image collection, and set these parameters during image processing. Object distance should be minimized and kept constant where possible. Emissivity can either be estimated from the literature (cf. Scheffers et al., 2017) or sampled in the field (FLIR, 2016).
Temperature measurements will be more accurate if there is less variation in surface emissivity within each image. That said, common components of the ground surface, such as leaves and soil, have very similar, high emissivity (> 0.9; Snyder et al., 1998) and thus the impact of emissivity variation on temperature measurements is relatively low.
Likewise, reflected apparent temperature can be sampled in the field (FLIR, 2016), but for high emissivities and short object distances, relatively little radiation is reflected and thus apparent temperature can be assumed to equal the atmospheric temperature (Tattersall, 2017).
It is recommended that thermal cameras are regularly calibrated (FLIR Systems suggest doing so once per year).

| Step 2: data extraction
A single thermal photo from a FLIR model E40 camera comprises 160 × 120 pixels, each of which is a unique measurement of received infrared radiation encoded as a raw 16-bit value. Images are extracted in batch by the function batch _ extract, which requires only the path to a directory of FLIR images, and the external software ExifTool.
Instructions for installing ExifTool can be found here: https ://sno.phy. queen su.ca/~phil/exift ool/insta ll.html. Windows users may find it easier to download ExifTool to a custom location, specified in the argument exiftoolpath of batch _ extract.
All the code and examples that follow were run on a 64-bit Windows 10 machine with an Intel Core i7-4600U @ 2.10 GHz processor and 8GB RAM. The function batch _ extract took 12 min to process 290 images on this machine, and 8 min on another F I G U R E 1 Overview of the main ThermStats functions and a typical workflow. Dashed lines indicate optional steps. Note that the precise functions used in Step 4 will depend on the study design machine (64-bit Windows 10 machine with an Intel Core i7-7820X CPU @ 3.60 GHz processor and 32GB RAM). Manual extraction using FLIR Tools, for comparison, took 62 min.

| Step 3: conversion of raw data
The raw values embedded in a FLIR thermal image can be converted to temperature in °C using equations from infrared thermography (FLIR, 2016;Tattersall, 2017) in the function batch _ convert. Temperature conversion will be more accurate when the user defines the environmental parameters, described in Step 1. Notably, default emissivity is 1, but should usually take a value between 0.95 and 0.97 (Tattersall, 2017), and relative humidity is highly dependent on habitat and weather conditions. Conversion of raw data also requires various calibration constants that are specific to each camera; these are retrieved automatically in batch _ extract to be passed directly to batch _ convert.

| Pixel-level statistics
The most relevant metrics to quantify thermal heterogeneity depend on the particular research questions. The function get _ stats takes a single thermal dataset, in the form of a matrix or raster , and calculates user-defined summary statistics across all pixels.
Standard summary statistics could include measures such as mean and standard deviation, but may also include metrics such as thermal richness (the number of unique temperature values) and thermal diversity indices (cf. Faye et al., 2017). Several helper functions are available to implement less standard summary statistics. On the basis of discussions in Faye et al. (2017) and Shi et al. (2016), we recommend some statistics for use by ecologists in Table 1.
For organisms seeking cooler microclimates under climate warming, the extent to which pixels in an image are connected to cooler pixels could be a proxy for the efficiency of thermoregulatory movements. From a given starting pixel, the greater maximum temperature difference that can be achieved by traversing from hotter to cooler pixels, the more effective local movement is likely to be as a mechanism to avoid excessively high temperatures. This concept of thermal connectivity is based on the landscape-level climate connectivity measure of McGuire et al. (2016), and can be calculated using the function connectivity. In absolute terms this metric is more suitable for coarser spatial scales that encompass the species' range, but could also be useful in relative terms to compare thermal connectivity between sites.

| Patch-level statistics
The function get _ stats identifies hot and cold spots in thermal images using a standalone function get _ patches. Hot and cold spots are based on the G* variant of the Getis-Ord local statistic (Getis & Ord, 1996), calculated using the spdep package (Bivand & Piras, 2015). The statistic is calculated for individual pixels by comparing the local weighted average of the pixel and its neighbours to the global average. High positive values exceeding the Z-value threshold (defined according to the sample size; Getis & Ord, 1996) are classified as hot spots, and low negative values as cold spots. Several spatial statistics are then calculated to characterize the hot and cold spots (Table 2; cf. Faye et al., 2016). There is an option to return patch outlines as a SpatialPolygonsDataFrame, which can be plotted on the temperature data using plot _ patches alongside an (optional) histogram of the temperature distribution (Figure 2).

Summary statistic Description
Average temperature Provides context for all other statistics, and could be used as a measure of the macroclimate for small, surface-dwelling organisms. The median is more robust than the mean to spurious extreme values that can sometimes arise in thermal images.
Temperature extremes While more rarely encountered, extreme values can be more significant to organisms, for example by exceeding upper thermal limits or by providing cool refugia from average conditions. The difference between extremes provides a measure of thermal diversity/stability (Shi et al., 2016), whereas the difference between extremes and average temperature provides a measure of the potential for thermal buffering. Again, we suggest the 5th and 95th percentiles are more robust to spurious extreme values than the minimum and maximum (respectively).
Temperature variability Over space and time, the standard deviation or coefficient of variation in temperature represents another measure of thermal stability (Shi et al., 2016), which may be particularly significant for organisms requiring constant temperatures, for example juveniles with a lower capacity for thermoregulatory behaviours. In contrast, for other mobile organisms-particularly ectotherms-high thermal diversity is likely to maximize opportunities for thermoregulation.
Thermal diversity indices Capture both the richness and evenness of different temperatures. Similar to temperature variability, the biological relevance of this measure is through its influence on the necessity and potential for thermoregulation. As discussed by Faye et al. (2016), Shannon's thermal diversity index quantifies how reliably one can predict the temperature of a pixel sampled at random from the temperature data. Simpson's thermal diversity index is similar; it captures the likelihood of two pixels being the same temperature (or temperature class) when sampled at random.
Thermal connectivity Calculated for each pixel as the potential for temperature change, which is the maximum temperature difference that can be achieved by traversing a gradient of hotter to cooler pixels. Adapted from McGuire et al. (2016).
TA B L E 1 Suggested summary statistics that can be applied by get _ stats

| Multiple images
We assume that for most users the spatial unit of replication will comprise multiple thermal images. In this case, the user can specify a grouping variable in stats _ by _ group. Matrices from each group will be bound together and get _ stats applied over the combined matrix. Processing times increase with more or larger images, although this can be managed by disabling the calculation of hot and cold spots More irregular microclimates could be less thermally stable but easier to locate, because of the greater proportion of edge.
Aggregation Index Degree of clustering in space of hot/cold spots, with zero representing no clustering (He, DeZonia, & Mladenoff, 2000) Dispersed, non-clustered microclimates are easier for animals to locate (Sears et al., 2016).

Patch Cohesion Index
Physical connectedness of hot/cold spots (Schumaker, 1996) More connected microclimates may facilitate more efficient travel.  (Getis & Ord, 1996), with high positive values assigned to hot spots (outlined and hatched in pink) and low negative values to cold spots (outlined and hatched in blue) or of thermal connectivity (patches = FALSE and connectivity = FALSE), since these are the most computationally intensive statistics. We assume that images are not adjacent in space, and therefore pad matrices with NA values before binding. Table S2 gives an example of the output from stats _ by _ group. Users can also average across repeated samples of the same area using average _ by _ group, the output of which can be passed directly to stats _ by _ group.
Raster stacks are generally quicker to process, and can be created from a list of matrices using stack _ imgs.

| Thermal heterogeneity over time in logged and unlogged forests
We demonstrate a typical workflow using fine-scale temperature data collected in the field with a FLIR thermal camera. To investigate how thermal heterogeneity varies over time and with selective logging, we sampled surface temperature in a large area of contiguous F I G U R E 3 Trends in various measures of thermal heterogeneity over the day (06:00-14:30 hr) for fine-scale temperature data collected using a thermal camera in primary (blue) and logged forests (orange). From left to right and top to bottom, the metrics are: median temperature (a); thermal Shannon Diversity Index (SHDI) (b); 95th percentile minus median temperature (c); median temperature minus 5th percentile (d); the average area (cm 2 ) per hot spot (e); the average area (cm 2 ) per cold spot (f); the number of hot spots per unit area (g); the number of cold spots per unit area (h); the Shape Index (SI) of hot spots (i); the Shape Index (SI) of cold spots (j); the Aggregation Index (AI) of hot spots (k); and the Aggregation Index (AI) of cold spots (l). Solid lines are model-predicted values with 95% confidence intervals. Semitransparent background points represent the raw data. Statistically significant differences are indicated by asterisks: .01 < p < .05 (*); .001 < p < .01 (**) and thermal images were taken in four orthogonal directions, with the camera held at breast height and pointing 45° downwards (relative to the ground). In total we collected 2,972 photos across 144 points. For full details see Scheffers et al. (2017) and Senior et al. (2018).
Data were extracted using batch _ extract and converted using batch _ convert (see Text S1 for a full worked example in the package vignette). We focused on the following summary statistics from Table 1: median temperature; Shannon Diversity Index; upper temperature range (95th percentile-median); and lower temperature range (median-5th percentile). For hot and cold spots separately, we calculated the following spatial statistics (Table 2): average area per patch (total area divided by number of patches, to correct for points with missing photos); average number of patches per unit area (density); Shape Index and Aggregation Index. We used stats _ by _ group to calculate each metric across all four photos taken each time a point was sampled (one photo in each direction).
To model the various thermal heterogeneity metrics against time of day and forest type (categorical: primary or logged) we fit Generalized Additive Mixed Effects Models (GAMMs), using the gamm4 package (Wood & Scheipl, 2017) in R. All measures of thermal heterogeneity were comparable between primary and unlogged forest (p > .05; Figure 3), but varied from the coolest time of day (~06:00 hr) to the hottest (~12:30 hr). Of particular note is the peak in thermal diversity (Shannon Diversity Index) at noon, and the different temporal patterns of hot spots (median temperature: 27°C) compared to cold spots (median temperature: 25°C). Our results suggest that logged forests provided comparable opportunities for thermoregulation as nearby primary forest, and that at noon, cold spots provided the greatest amount of temperature buffering.

| Averaging over repeated images
Averaging across repeated images of the same area could help to reduce unwanted variation, such as that introduced by the F I G U R E 4 Violin plots of temperature distribution across the raster of each distinct replicate in time. Each datapoint represents one pixel of the raster, with 100 pixels sampled at random. Each raster is the average of multiple thermal images collected at the same time and location, averaged within temporal replicates using the function average _ by _ group F I G U R E 5 Panel (a) shows the study area in the Stability of Altered Forest Ecosystems (SAFE) project, overlaid with the boundary of oil palm plantations (Ewers et al., 2011;Kahle & Wickham, 2013). Variation in modelled mean annual temperature is shown in panel (b) (Jucker et al., 2018). Panel (c) shows the temperature distribution of hot spots (orange) and cold spots (green), inside and outside plantations camera or photographer. In this example, thermal images were collected at nine points ('spatial replicate') spaced along ex- ideal for visualizing these data (Figure 4), as well as functions in the package rasterVis (Perpiñán & Hijmans, 2018).

| Using other gridded temperature data
Analytical functions in ThermStats are designed for any gridded temperature data. We illustrate this using a mean annual tempera- quantiles and the median, which together characterize the temperature distribution and can be used to construct boxplots ( Figure 5).
Qualitatively, we demonstrate that both hot and cold spots were hotter inside the oil palm plantation than the neighbouring forest.

| CON CLUDING REMARK S
Our R package presents users with a simple protocol for using thermography within ecology. The extraction of data from FLIR thermal cameras is made faster and easier, and for any gridded temperature data the package facilitates the calculation of different, biologically relevant metrics of thermal heterogeneity (Faye et al., 2016;Senior et al., 2018;Shi et al., 2016).
It is important to consider the strengths and weaknesses of thermography when deciding on the most appropriate methodology to answer research questions of interest. Thermal cameras cannot directly measure subsurface temperatures and are not as well suited for capturing temporal variation, compared with dataloggers (Senior et al., 2018). Although affordable smartphone attachments are now available, thermal cameras may still be more expensive than dataloggers (depending on the quantity of dataloggers required), and can be sensitive to extreme weather conditions common to regions such as the tropics and Arctic (FLIR, 2016). Bramer et al.
(2018) is an excellent resource for ecologists seeking best practice for using dataloggers; we intend that our study and the references herein offer a similar resource for ecologists using thermography.
As infrared technology develops, we anticipate that other applications and brands of thermal camera will be relevant to ecologists, and we therefore invite readers to suggest package enhancements via GitHub (https ://github.com/rasen ior/Therm Stats/ issues).
Fine-scale temperature variation across space and time has a huge influence on species' ecology, which will become increasingly pertinent as average temperatures rise under global climate warming. We showcase how our R package and framework can be used to quantify thermal heterogeneity using data at a fine spatial scale, collected using a FLIR thermal camera or modelled using remotely sensed data.
By simplifying and streamlining the processing of increasingly available thermal data, our approach enables researchers to more readily address important issues in ecology and conservation.

ACK N OWLED G EM ENTS
We thank Tommaso Jucker for providing the modelled mean annual temperature layer, and to Manoela S. Machado Mollinari for helpful discussions. We are very grateful to the reviewers for their insightful comments and suggestions. R.A.S. was funded by a NERC studentship through the ACCE (Adapting to the Challenges of a Changing Environment) Doctoral Training Partnership (Grant No. NE/L002450/1).

DATA AVA I L A B I L I T Y S TAT E M E N T
The development version of ThermStats can be downloaded from GitHub: https ://github.com/rasen ior/Therm Stats . Bug reports and suggested enhancements can be submitted to: https ://github.com/ rasen ior/Therm Stats/ issues. This manuscript used ThermStats version 0.9.6, which is deposited on Zenodo (https ://doi.org/10.5281/zenodo.3264016; Senior, 2019).