Fix heatload function

[shuysman]

I noticed that the way heatload is calculated in get_OudinPET doesn't seem correct, particularly for steep, south slopes. The math looks like it should be correct, but I think there is an issue with parenthesis and order of operations. As written the function is difficult to read, so I refactored it to calculate heat load in a separate step. I think having a separate function for heat load would be useful as it is an output I'd like to include with my spatial analyses. I also added helper functions for degrees to radian conversions (deg2rad) and folded aspect calculation (fold_aspect). get_OudinPET is still a mess, in terms of readability, and could use a few more changes.

I also noted that the heatload function we use only accepts inputs from 1-60 degrees.

Old and new functions can be compared with this script:

library(dplyr)
library(plotly)
library(WaterBalance)

slope = 0:60
aspect = 0:180
lat = 43

vars <- expand.grid(slope = slope, aspect = aspect, lat = lat)

vars$hl <- get_heatload(vars$slope, vars$aspect, vars$lat)
vars$old_hl <- (0.339 + 0.808 * cos(deg2rad(vars$lat)) * cos(deg2rad(vars$slope))) -
    (0.196 * sin(deg2rad(vars$lat)) * sin(deg2rad(vars$slope))) -
    (0.482 * cos(deg2rad(fold_aspect(vars$aspect))) * sin(deg2rad(slope)))

plot_ly(data = vars, x = ~slope, y = ~aspect, z = ~hl, type = "mesh3d")
plot_ly(data = vars, x = ~slope, y = ~aspect, z = ~old_hl, type = "mesh3d")

Compare new output:

To old output:

Steep, south slopes were getting a HL < 1, which doesn't make sense. This new output roughly matches what I get with the python WB script for a similar latitude, but I have not verified formally.

[shuysman]

So I noticed I screwed up some parens in my new function. Now the output seems to match what I was getting with the old method. I don't think this output is correct, because using the python script I get HL ranging from 0.3555403553162483 - 1.2846587281750943 for an example site around 43 degrees latitude. Also the McCune & Keon 2002 paper specifies HL should range from 0.03 - 1.11, so there are a few things going on here. Might have to think more about this. I think my commits here are an improvement in readability, but I think the output is still incorrect.

[shuysman]

I had been setting up the python function incorrectly (giving it Lat in degrees instead of radians), which explains the mismatched results. When the right units are used both Python and R seem to match. I do still think that heat load should be it's own function as it would be useful to save the results from for diagnostic purposes.