Estimating leaf temperatures for spatial epidemiology

R4PlantPath Reproducible Research

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

Paul Melloy


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")).



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:

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.

First, a brief tealeaves example

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:

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).

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.

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  
 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>

Estimate leaf temperatures

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]

Leaf temperature gradients

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)

[1] 64 25
head(Tm_leaves[, c("T_leaf", "E")])
# 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.

Leaf temperature time-series

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.

  mean(np$ALLSKY_SRF_ALB, na.rm = TRUE)
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")
      leaf_par = lf_par,
      enviro_par = env_par,
      constants = const,
      quiet = TRUE

# bind the list into a data.frame
Leaf_temperatures <-"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 <-

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")

Plotting leaf temperature

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


       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() +
    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:
Muir, C. D. 2020. tealeaves: An R package for modelling leaf temperature using energy budgets. Available at:
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:
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, 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 ...".