Author: David Zelený & Tsung-Yi Lin

The function uses aspect, slope and latitude of the plot and calculates potential annual direct incident radiation and heat load, 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!

Arguments:

• aspect - aspect of the plot, vector (either in degrees, or radians)
• slope - the inclination of the plot, vector (either in degrees, or radians)
• latitude - latitude of the plot, either a vector of the same length as aspect or slope, or a single value (in case that all plots are from the same relatively small region)
• method - default is “heatload”, alternative is “radiation”
• units - default is “degrees”, alternative is “radians”
• equation - the equation number (1, 2 or 3); default is 1 (the most general one). The three equations have slightly different uses. Eq. 1 (default in the function below) is broadest in the application, covering all slopes <= 90° in steepness at latitudes 0-60°N, but has the lowest precision. Eq. 2 increases the precision by excluding slopes steeper than 60°, an inconsequential omission for almost all data sets. Eq. 3 uses only three parameters to produce a slightly stronger model but is applicable only to latitudes 30-60°N (ie not for the tropical and subtropical region).

heatload.r

# Function "heatload"
# Calculates heat load or potential annual direct incident radiation, using the formulas published in 
# McCune & Keon (2002) based on aspect, slope and latitude.
 
heatload <- function (aspect, slope, latitude, method = 'heatload', units = 'degrees', equation = 1)
{
  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 Třebíč 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, equation = 3)