How to obtain leaf temperature estimates world-wide using the R packages, tealeaves and nasapower.

The incidence and severity of plant diseases are dependant on environmental effects which is emphasised by the disease triangle (Francl 2001). Arguably the effect of temperature on disease can be the most important of environmental influences. One problem with evaluating the effect of temperatures on disease incidence and severity is the difference between air temperatures measured at a weather station and the micro-environment of the plant canopy.

When modelling plant disease outbreaks, being able to accurately estimate the conditions which lead to localised plant disease epidemics can be very useful. This led me to the R package *tealeaves* (Muir 2019), which provides functions for calculating leaf temperatures from ambient environmental inputs.

If you have not already, install the packages *tealeaves* and *nasapower* with `install.packages(c("tealeaves", "nasapower"))`

.

```
library(tealeaves)
```

Before I begin, it is worth noting *tealeaves* uses the package, *units*, to format the inputs and outputs (Pebesma et al. 2016). Detail about the package *tealeaves* and it’s functions can be found in Muir’s publication (Muir 2020), and in the vignettes which can be obtained by running the following commands in R:

`vignette("tealeaves-introduction")`

,`vignette("tealeaves-intermediate")`

For brevity, I will not go in depth into *tealeaves* but provide a simple example using default parameters, then how to change default environmental parameters with ones obtained from the *nasapower* package.

```
lf_par <- make_leafpar()
env_par <- make_enviropar()
const <- make_constants()
Tm_leaf <-
tleaf(leaf_par = lf_par,
enviro_par = env_par,
constants = const)
Tm_leaf[, c("T_leaf", "E")]
```

```
T_leaf E
1 301.4181 [K] 0.00794791 [mol/m^2/s]
```

The function `tleaf()`

calculates leaf temperature and evapotranspiration by solving an energy budget equation from incoming solar radiation, air temperatures, reflectance, and heat loss to wind and evapotranspiration for a leaf with certain characteristics (Muir 2019). This equation requires these parameter inputs which describe this energy budget for calculating leaf temperature. The parameters are created or altered through the package functions:

`make_leafpar()`

defines the leaf parameters and characteristics,`make_enviropar()`

sets environmental parameters,`make_constants()`

defines the constants used in the`tleaf()`

equation.

The `tleaf()`

function returns an output containing details of the thermal energy absorbed and emitted, including an estimate of evapotranspiration (`"E"`

), which assists in defining the estimated leaf temperature (`"T_leaf"`

), given in Kelvin.

We can replace the default environmental parameters with values parsed from a weather station or remote sensing satellite data retrieved by using the R package *nasapower* (Sparks 2018; Sparks 2022).

Let’s download some hourly data using the *nasapower* (Sparks 2018; Sparks 2022) package. More information on using the *nasapower* package can be found by calling the vignette `vignette(nasapower)`

.

```
library(nasapower)
np <- get_power(
community = "ag",
pars = c(
"RH2M", # Relative humidity at 2 meters
"T2M", # Air temperature at 2 meters
"WS2M", # wind speed at 2 meters
"ALLSKY_SFC_SW_DWN", # short wave downward solar radiation at the earth surface
"ALLSKY_SFC_LW_DWN", # long wave downward solar radiation at the earth surface
"PS", # surface pressure
"ALLSKY_SRF_ALB" # Surface albedo
),
temporal_api = "hourly",
lonlat = c(150.88, -26.83),
dates = c("2021-02-01", "2021-02-08")
)
```

In the `get_power()`

function above, I have requested data from the NASA “POWER” project from the agricultural weather data collection (`"ag"`

). These data contain all the energy balance inputs needed from a 0.5° grid cell at the specified latitude and longitude (`latlong`

). The coordinates I have used are just north of Dalby in Queensland, Australia. Additionally, I have requested hourly data for the first week in February 2021 with the `"temporal_api"`

argument.

```
np
```

```
NASA/POWER CERES/MERRA2 Native Resolution Hourly Data
Dates (month/day/year): 02/01/2021 through 02/08/2021
Location: Latitude -26.83 Longitude 150.88
Elevation from MERRA-2: Average for 0.5 x 0.625 degree lat/lon region = 336.69 meters
The value for missing source data that cannot be computed or is outside of the sources availability range: NA
Parameter(s):
Parameters:
RH2M MERRA-2 Relative Humidity at 2 Meters (%) ;
T2M MERRA-2 Temperature at 2 Meters (C) ;
WS2M MERRA-2 Wind Speed at 2 Meters (m/s) ;
ALLSKY_SFC_SW_DWN CERES SYN1deg All Sky Surface Shortwave Downward Irradiance (MJ/hr) ;
ALLSKY_SFC_LW_DWN CERES SYN1deg All Sky Surface Longwave Downward Irradiance (W/m^2) ;
PS MERRA-2 Surface Pressure (kPa) ;
ALLSKY_SRF_ALB CERES SYN1deg All Sky Surface Albedo (dimensionless)
# A tibble: 192 × 13
LON LAT YEAR MO DY HR RH2M T2M WS2M
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 151. -26.8 2021 2 1 10 39.8 30.9 3.72
2 151. -26.8 2021 2 1 11 34.1 32.6 3.51
3 151. -26.8 2021 2 1 12 30.5 33.8 3.37
4 151. -26.8 2021 2 1 13 28.3 34.5 3.39
5 151. -26.8 2021 2 1 14 27.3 34.6 3.6
6 151. -26.8 2021 2 1 15 27.7 33.8 3.9
7 151. -26.8 2021 2 1 16 29.1 32.6 4.13
8 151. -26.8 2021 2 1 17 31.9 30.9 4.08
9 151. -26.8 2021 2 1 18 38.5 28.3 3.38
10 151. -26.8 2021 2 1 19 43 26.3 4.13
# … with 182 more rows, and 4 more variables:
# ALLSKY_SFC_SW_DWN <dbl>, ALLSKY_SFC_LW_DWN <dbl>, PS <dbl>,
# ALLSKY_SRF_ALB <dbl>
```

Next, I will characterise the environmental parameters using the `np`

data at a single time point into the `tleaf()`

.

```
env_par <- make_enviropar(
replace = list(
P = set_units(np$PS[1], "kPa"), # Air pressure
r = set_units(np$ALLSKY_SRF_ALB[1]), # Albedo
RH = set_units(np$RH2M[1]/100), # Relative Humidity (as a proportion)
S_sw = set_units(np$ALLSKY_SFC_SW_DWN[1],
"W/m^2"), # Short wave radiation
T_air = set_units(
set_units(np$T2M[1], "DegC"),
"K"), # Air temperature, converted from Celsius to Kelvin
wind = set_units(np$WS2M[1], "m/s") # Wind speed
)
)
```

Now I can run the `tleaf()`

function again with the updated environmental parameters.

```
Tm_leaf <-
tleaf(leaf_par = lf_par,
enviro_par = env_par,
constants = const)
Tm_leaf[, c("T_leaf", "E")]
```

```
T_leaf E
1 299.5096 [K] 0.006604466 [mol/m^2/s]
```

Because the outputs are specified using the *units* package I can easily convert from Kelvin back to degrees Celsius using the `set_units()`

function again.

```
set_units(Tm_leaf$T_leaf, "DegC")
```

`26.35956 [°C]`

Imputing vectors as parameters returns leaf temperature gradients across the supplied parameter permutations and not a vector of times. Below is an example of the permutations `t_leaves()`

processes from a vector of environmental parameters. This is described further in the *tealeaves* vignette.

```
env_par <- make_enviropar(
replace = list(
P = set_units(np$PS[1:2], "kPa"),
r = set_units(np$ALLSKY_SRF_ALB[1:2]),
RH = set_units(np$RH2M[1:2] / 100), # needs to be specified as a proportion
S_sw = set_units(np$ALLSKY_SFC_SW_DWN[1:2], "W/m^2"),
T_air = set_units(set_units(np$T2M[1:2], "DegC"),
"K"), # convert from degrees to Kelvin
wind = set_units(np$WS2M[1:2], "m/s")
)
)
Tm_leaves <-
tleaves(leaf_par = lf_par,
enviro_par = env_par,
constants = const)
dim(Tm_leaves)
```

`[1] 64 25`

```
# A tibble: 6 × 2
T_leaf E
[K] [mol/m^2/s]
1 299. 0.00710
2 299. 0.00724
3 300. 0.00760
4 300. 0.00775
5 299. 0.00710
6 299. 0.00724
```

We can see that the `tleaves()`

function uses the multiple values specified in the environmental parameters to calculate all possible gradients. As such 64 values are returned (\(2^6\)), permutations of the input environmental parameters.

*However in this example we want leaf temperature estimates over time, not space.*

To obtain a time series I can use the `apply()`

function to pass the environmental parameters to `tleaf()`

one time point at a time.

First, I need to clean up the data a little by substituting `NA`

albedo values generated by night-time measurements. Here I have substituted it with the mean.

```
tm_leaf_list <-
apply(np, 1, function(nasa_p) {
env_par <- make_enviropar(
replace = list(
P = set_units(nasa_p["PS"], "kPa"),
r = set_units(nasa_p["ALLSKY_SRF_ALB"]),
RH = set_units(nasa_p["RH2M"] / 100), # needs to be specified as a proportion
S_sw = set_units(nasa_p["ALLSKY_SFC_SW_DWN"], "W/m^2"),
T_air = set_units(set_units(nasa_p["T2M"], "DegC"),
"K"), # convert from degrees to Kelvin
wind = set_units(nasa_p["WS2M"], "m/s")
)
)
return(tleaf(
leaf_par = lf_par,
enviro_par = env_par,
constants = const,
quiet = TRUE
))
})
# bind the list into a data.frame
Leaf_temperatures <-
do.call("rbind", c(tm_leaf_list, make.row.names = FALSE))
Leaf_temperatures[1:10, c("T_leaf", "E")]
```

```
T_leaf E
1 299.5096 [K] 0.006604466 [mol/m^2/s]
2 300.3629 [K] 0.007602218 [mol/m^2/s]
3 300.9245 [K] 0.008296177 [mol/m^2/s]
4 301.2553 [K] 0.008803558 [mol/m^2/s]
5 301.3194 [K] 0.009132864 [mol/m^2/s]
6 300.9297 [K] 0.009061479 [mol/m^2/s]
7 300.2532 [K] 0.008639131 [mol/m^2/s]
8 299.0727 [K] 0.007694088 [mol/m^2/s]
9 297.1887 [K] 0.005870064 [mol/m^2/s]
10 296.0459 [K] 0.005306279 [mol/m^2/s]
```

I can now pair up the time and air temperature data from the NASA power data with the leaf temperature data. Also, I will convert temperature in Kelvin to degrees Celsius.

```
Leaf_temperatures$air_temp <-
np$T2M
Leaf_temperatures$times <-
as.POSIXct(paste0(np$YEAR, "-",
np$MO, "-",
np$DY, " ",
np$HR, ":00:00"),
tz = "Australia/Brisbane")
Leaf_temperatures$T_leaf <-
set_units(Leaf_temperatures$T_leaf, "DegC")
```

Finally, I can visualise the change in air temperature and leaf temperatures over time.

```
library(ggplot2)
ggplot(Leaf_temperatures,
aes(x = times)) +
geom_line(aes(y = air_temp, colour = "darkblue"), size = 1) +
geom_line(aes(y = drop_units(T_leaf),
colour = "forestgreen"), size = 1) +
labs(y = "Temperature (Degrees Celcius)",
color = "Estimates") +
theme_minimal() +
scale_color_manual(
values = c("darkblue", "forestgreen"),
labels = c("Air Tm",
"Leaf Tm")
)
```

Hopefully this will help get you started for using these two packages to estimate leaf and canopy temperatures. However, further parameterisation would be necessary prior to any in-depth analyses. Leaf parameters, such as leaf size, reflectance, absorbance, stomatal conductance and ratio should be carefully considered in addition to vertical environmental gradients across the height of the plant. Shading in the canopy is known to effect pathogen infection, colonisation and sporulation in addition to leaf wetness.

Francl, L. J. 2001. The disease triangle: A plant pathological paradigm revisited. The Plant Health Instructor.

Muir, C. D. 2019. tealeaves: an R package for modelling leaf temperature using energy budgets. AoB PLANTS. 11 Available at: https://academic.oup.com/aobpla/article/11/6/plz054/5666155.

Muir, C. D. 2020. tealeaves: An R package for modelling leaf temperature using energy budgets. Available at: https://CRAN.R-project.org/package=tealeaves.

Pebesma, E., Mailund, T., and Hiebert, J. 2016. Measurement units in R. R Journal. 8:486–494.

Sparks, A. 2022. *nasapower: NASA-POWER Data from R*. Available at: https://CRAN.R-project.org/package=nasapower.

Sparks, A. H. 2018. nasapower: A NASA POWER Global Meteorology, Surface Solar Energy and Climatology Data Client for R. The Journal of Open Source Software. 3:1035.

If you see mistakes or want to suggest changes, please create an issue on the source repository.

Text and figures are licensed under Creative Commons Attribution CC BY 4.0. Source code is available at https://github.com/openplantpathology/OpenPlantPathology, unless otherwise noted. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".