Author: David Zelený & Tsung-Yi Lin

The function uses aspect, slope and latitude of the plot and calculates heat load or radiation according to equations published by McCune & Keon (2002).

Note: prior to 4/16/2019, the definition of this function had a bug and returned incorrect values; thanks to Tsung-Yi Lin (林宗儀) for noticing the problem and fixing it!

heatload.r

# Function heatload
# Calculates heatload or potential annual direct incident radiation, using the formulas published in 
# McCune & Keon (2002) based on aspect, slope and latitude.
# Arguments: 
# aspect - aspect of the plot, vector (either in degrees, or radians)
# slope - inclination of the plot, vector (in degrees)
# latitude - latitude of the plot, either vector of the same length as aspect or slope, or single value (in case that all plots are from the same relatively small region)
# method - default is "heatload", alternative is "radiation".
# units - default id "degrees", alternative is "radians".
# equation - the number of equation (1, 2 or 3); default is 3.
# Author: David Zeleny & Tsung-Yi Lin
heatload <- function (aspect, slope, latitude, method = 'heatload', units = 'degrees', equation = 3)
{
  if (units == 'degrees')   # convert degrees to radians
    {
    aspect <- aspect/180*pi
    slope <- slope/180*pi
    aspect[slope == 0] <- 0
    latitude <- latitude/180*pi
    }  
  A <- if (method == 'heatload') abs (pi - abs (aspect - (5*pi/4))) else pi - abs (aspect-pi)
  S <- slope
  L <- if (length (latitude) == 1) rep (latitude, length (A)) else latitude
  if (equation == 1) res <- exp (-1.467 +1.582*cos(L)*cos(S) -1.500*cos(A)*sin(S)*sin(L) -0.262*sin(L)*sin(S) +0.607*sin(A)*sin(S))
  if (equation == 2) res <- exp (-1.236 +1.350*cos(L)*cos(S) -1.376*cos(A)*sin(S)*sin(L) -0.331*sin(L)*sin(S) +0.375*sin(A)*sin(S))
  if (equation == 3) res <-      +0.339 +0.808*cos(L)*cos(S)                             -0.196*sin(L)*sin(S)                       - 0.482*cos(A)*sin(S)
  return (res)
}

Example of the use on the dataset grasslands.env (acidophilous grasslands in Trebic region, Czech Republic). The grasslands.env contains variables aspect, slope and latitude, all in degrees. Dataset represents small plots (16-25m2) located on convex outcrops in the agricultural landscape, covered by seminatural grassland vegetation.

source ('http://anadat-r.davidzeleny.net/doku.php/en:customized_functions:heatload?do=export_code&codeblock=0') # reads the function definition from above
grasslands.env <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/grasslands-env.txt') # loads the file
# Calculate heat load and add it to the grasslands.env data.frame:
grasslands.env$heatload <- heatload (aspect = grasslands.env$aspect, slope = grasslands.env$slope, latitude = grasslands.env$latitude)