I feel like this should be a simple fix but can't figure it out. I have two rasters, one with temperature data and one with humidity data. I need to calculate the heat index using a set of equations and am trying to get a raster output. I put the heat index equations into a function, which is as follows (feel free to ignore the equations themselves, I know it's a lot of numbers, this is not where my issue is):
fun_heat_index <- function(RH, T){
# https://www.wpc.ncep.noaa.gov/html/heatindex_equation.shtml
#RH <- relative humidity in percent
#T <- maximum temperature in degrees F
# if HI is below 80 degrees, apply only this simple formula
iheat_index <- 0.5 * (T + 61 + ((T-68) * 1.2) + (RH * 0.094))
if (iheat_index > 80){
iheat_index <- -42.379 + 2.04901523 * T + 10.14333127 * RH - .22475541 * T *
RH - .00683783 * T * T - .05481717 * RH * RH + .00122874 * T * T * RH +
.00085282 * T * RH * RH - .00000199 * T * T * RH * RH
# ADJUSTMENT 1:
#if RH less than 13% and temp is between 80 and 112
if (RH < 13 && between(T, 80, 112)){
adj_1 <- ((13 - RH) / 4) * sqrt((17-abs(T-95))/17)
iheat_index <- iheat_index + adj_1
} # if statement Adjustment 1
# ADJUSTMENT 2:
#if RH is greater than 85% and temp is between 80 and 87
if (RH > 85 && between(T, 80, 87)){
adj_2 <- ((RH - 85) / 10) * ((87 - T) / 5)
iheat_index <- iheat_index + adj_2
} # if statement Adjustment 2
} # if statement simple equation over 80% HI
# Output the result after going through adjustments & equations
return(iheat_index)
} # fun_heat_index
This is what I was trying to run using the function above which was just coming up with an error and not running through.
fun_heat_index(relative_humidity_data, tmax_data)
My rasters are the same extent and each have the same number of grid cells. I'm trying to apply the function to calculate from the two rasters to get a third output raster with the heat index values. I found the calc() function that makes it go through each grid cell individually, which I thought that would fix the problem. I tried this next:
calc(tmax_data, fun = fun_heat_index(relative_humidity_data, tmax_data))
It started working a little after I added that, but it seems to be tripping up after the initial iheat_index variable creation still, which is the first step in the function. It will output a raster from that line, but once it goes to the next line into the if statement if (iheat_index > 80){ then it says argument is not interpretable as logical, which I'm assuming is because it's expecting a single value not a raster. Is there some way to apply the calc() function again in the if statements? Or is there another way to get R to just look at some grid cells and apply additional equations to only grid cells that meet a certain criteria?
In looking around for how to do this, I found overlay() as well and tried that as follows:
overlay(iloop_tmax, iloop_rh, fun = fun_heat_index(iloop_rh, iloop_tmax))
That didn't work either. I also tried going through each cell and removing the values using getValues() and applying the function to any cells that met the criteria for the if statements, but I couldn't figure out how to get it back into raster form rather than tabular form after that.
Any help is much appreciated! :)
I try to get simple things working in a step by step manner as it makes my many mistakes easier to detect and debug. Taking your above dropped files ending with .gri as T2 & RH2, then creating a range of values via sampling in each:
set.seed(42) # for reproducibility as we're sampling
values(RH2) <- sample(5:95, ncell(RH2), replace= TRUE)/100 # % relative humidity
values(T2) <- sample(70:115, ncell(T2), replace = TRUE)
# data T2 and RH2 at bottom
Moving to the heatindex_equiation, we see that traditionally, the simplest equation is calculated first and this result is averaged with T2:
HI_tst1 = 0.5 * (T2 + 61.0 +((T2-68.0) * 1.2) + (RH2 * 0.094))
HI_tst1_avg = 0.5 * (T2 + HI_tst1)
So it appears that HI_tst1_avg is our base data set, and further calculations and adjustments are applied to specific cells in this data set, depending on the stated criteria. And it is probably more accurate to say that HI_tst1_avg is an index to cells where further calcs will be applied (conceptually). So now, perhaps we do some fishing within our avg to get a sense of what calc and adjustments should likely be done:
# ? How many HI_tst1_avg >= 80
length(which(HI_tst1_avg@data@values >= 80))
[1] 402
# ? which cells in T2 and RH2 are these, here just indexed
gt80_idx <- which(HI_tst1_avg@data@values >= 80)
# ? How many HI_tst1_avg < 80 and subject to simple calc
length(which(HI_tst1_avg@data@values < 80))
[1] 138
# which cells are those indexed
HI_simple_idx <- which(HI_tst1_avg@data@values < 80)
# ? How many subject to first adjustment (RH <= 0.13 cells also T btw 80:87)
length(which(HI_tst1_avg@data@values >= 80 & HI_tst1_avg@data@values <=87 & RH2@data@values <= 0.13))
[1] 7
# ? Which cells in T2 & RH2 do these mean
intersect((which(HI_tst1_avg@data@values >= 80 & HI_tst1_avg@data@values <=87)),(which(RH2@data@values <= 0.13)))
[1] 20 45 72 177 184 422 441
# so, `adj1_idx = intersect((which...` save as index
# ? how many will be subject to adjustment 2
length(which(HI_tst1_avg@data@values >= 80 & HI_tst1_avg@data@values <=87 & RH2@data@values >= 0.85))
[1] 11
# ? which cells to index for adj2
intersect((which(HI_tst1_avg@data@values >= 80 & HI_tst1_avg@data@values <=87)),(which(RH2@data@values >= 0.85)))
[1] 56 221 230 300 307 346 353 393 412 470 485
And now I think I have the right indexes to start plugging into creating a heatindex raster. So, create it with NA values.
heatindex <- T2
values(heatindex) <- NA
names(heatindex) <- 'summer_day'
heatindex
class : RasterLayer
dimensions : 18, 30, 540 (nrow, ncol, ncell)
resolution : 0.0625, 0.0625 (x, y)
extent : -78.125, -76.25, 38.375, 39.5 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : memory
names : summer_day
values : NA, NA (min, max)
Working in the suggested order, 1)simple 2)regression 3) adj1 4) adj2
#simple
values(heatindex)[HI_simple_idx] = 0.5 * (values(T2)[HI_simple_idx] + 61.0 + ((values(T2)[HI_simple_idx] - 68.0) * 1.2) + (values(RH2)[HI_simple_idx] * 0.094))
# regression coefficients - Yikes, trusting to precedence
values(heatindex)[gt80_idx] = -42.379 + 2.04901523 * values(T2)[gt80_idx] + 10.14333127 * values(RH2)[gt80_idx] - .22475541 * values(T2)[gt80_idx] *
values(RH2)[gt80_idx] - .00683783 * values(T2)[gt80_idx] * values(T2)[gt80_idx] - .05481717 * values(RH2)[gt80_idx] * values(RH2)[gt80_idx] + .00122874 * values(T2)[gt80_idx] * values(T2)[gt80_idx] * values(RH2)[gt80_idx] +
.00085282 * values(T2)[gt80_idx] * values(RH2)[gt80_idx] * values(RH2)[gt80_idx] - .00000199 * values(T2)[gt80_idx] * values(T2)[gt80_idx] * values(RH2)[gt80_idx] * values(RH2)[gt80_idx]
# god help me, appears to work
# adjustment 1
values(heatindex)[adj1_idx] = values(heatindex)[adj1_idx] - ((13 - values(RH2)[adj1_idx]) /4) * sqrt((17 - abs(values(T2)[adj1_idx]- 95)/17))
#adjustment 2
values(heatindex)[adj2_idx] = values(heatindex)[adj2_idx] + ((values(RH2)[adj2_idx] - 85) / 10) * ((87 - values(T2)[adj2_idx]) / 5)
# are we nearly done?
which(is.na(values(heatindex)) == TRUE)
[1] 20 45 72 177 184 422 441
length(which(is.na(values(heatindex)) == TRUE))/ncell(heatindex)
[1] 0.01296296
all.equal(which(is.na(values(heatindex)) == TRUE), adj1_idx)
[1] TRUE
As they say, 'Heat Kills'. And there remains noodling around to figure out why adjustment 1 is not working, but these are the steps I'd walk through prior to making a function that wraps these steps. I also find these base approaches easier to pick up and reason about when revisited many months later.
Debugging the error meant first reversing the introduction of NAN(s) in the adj_idx cells, fiddling with ) placement for reasonable values
# redo gt80_idx
values(heatindex)[gt80_idx] = -42.379 + 2.04901523 * values(T2)[gt80_idx] + 10.14333127 * values(RH2)[gt80_idx] - .22475541 * values(T2)[gt80_idx] *
values(RH2)[gt80_idx] - .00683783 * values(T2)[gt80_idx] * values(T2)[gt80_idx] - .05481717 * values(RH2)[gt80_idx] * values(RH2)[gt80_idx] + .00122874 * values(T2)[gt80_idx] * values(T2)[gt80_idx] * values(RH2)[gt80_idx] +
.00085282 * values(T2)[gt80_idx] * values(RH2)[gt80_idx] * values(RH2)[gt80_idx] - .00000199 * values(T2)[gt80_idx] * values(T2)[gt80_idx] * values(RH2)[gt80_idx] * values(RH2)[gt80_idx]
# start values
(values(heatindex)[adj1_idx])
[1] 83.25591 82.37403 82.37816 80.57744 81.48331 84.12430 80.57744
# seem reasonable?
(values(heatindex)[adj1_idx]) - ((13 - values(RH2)[adj1_idx]) / 4) * sqrt((17 - abs(values(T2)[adj1_idx] - 95.0) / 17))
[1] 70.14729 69.32935 69.28284 67.58968 68.45192 70.96179 67.58968
# assign
values(heatindex)[adj1_idx] = (values(heatindex)[adj1_idx]) - ((13 - values(RH2)[adj1_idx]) / 4) * sqrt((17 - abs(values(T2)[adj1_idx] - 95.0) / 17))
> which(is.na(values(heatindex)) == TRUE)
integer(0)
Data:
# RH2
new("RasterLayer", file = new(".RasterFile", name = "", datanotation = "INT4U",
byteorder = c(value = "little"), nodatavalue = 2147483647,
NAchanged = FALSE, nbands = 1L, bandorder = c(value = "BIL"),
offset = 0L, toptobottom = TRUE, blockrows = 0L, blockcols = 0L,
driver = "", open = FALSE), data = new(".SingleLayerData",
values = c(0.38, 0.21, 0.65, 0.87, 0.47, 0.83, 0.19, 0.09,
0.5, 0.93, 0.78, 0.51, 0.57, 0.15, 0.76, 0.45, 0.32, 0.12,
0.88, 0.08, 0.86, 0.84, 0.84, 0.66, 0.67, 0.43, 0.67, 0.28,
0.17, 0.6, 0.06, 0.6, 0.66, 0.63, 0.32, 0.46, 0.75, 0.71,
0.86, 0.49, 0.32, 0.61, 0.89, 0.11, 0.12, 0.44, 0.49, 0.61,
0.48, 0.82, 0.84, 0.39, 0.3, 0.15, 0.49, 0.91, 0.55, 0.41,
0.18, 0.16, 0.18, 0.61, 0.48, 0.8, 0.52, 0.49, 0.72, 0.9,
0.24, 0.23, 0.05, 0.07, 0.37, 0.29, 0.21, 0.21, 0.89, 0.69,
0.42, 0.08, 0.87, 0.72, 0.12, 0.6, 0.67, 0.44, 0.54, 0.9,
0.93, 0.43, 0.65, 0.6, 0.49, 0.63, 0.18, 0.75, 0.78, 0.35,
0.53, 0.4, 0.62, 0.82, 0.7, 0.09, 0.81, 0.82, 0.08, 0.47,
0.67, 0.89, 0.67, 0.91, 0.1, 0.77, 0.91, 0.25, 0.58, 0.24,
0.34, 0.4, 0.52, 0.73, 0.31, 0.3, 0.39, 0.68, 0.21, 0.07,
0.82, 0.05, 0.62, 0.39, 0.29, 0.44, 0.15, 0.11, 0.27, 0.31,
0.61, 0.59, 0.28, 0.38, 0.93, 0.91, 0.74, 0.61, 0.44, 0.29,
0.35, 0.94, 0.11, 0.54, 0.15, 0.39, 0.14, 0.18, 0.68, 0.95,
0.65, 0.73, 0.48, 0.5, 0.52, 0.23, 0.24, 0.54, 0.66, 0.72,
0.63, 0.91, 0.24, 0.14, 0.33, 0.53, 0.61, 0.9, 0.13, 0.85,
0.73, 0.49, 0.78, 0.83, 0.08, 0.11, 0.65, 0.06, 0.52, 0.55,
0.77, 0.31, 0.44, 0.11, 0.42, 0.91, 0.71, 0.88, 0.58, 0.83,
0.12, 0.37, 0.12, 0.78, 0.88, 0.76, 0.66, 0.12, 0.28, 0.5,
0.48, 0.68, 0.53, 0.29, 0.19, 0.12, 0.42, 0.86, 0.7, 0.36,
0.34, 0.44, 0.87, 0.45, 0.2, 0.71, 0.77, 0.39, 0.23, 0.28,
0.41, 0.95, 0.57, 0.13, 0.44, 0.95, 0.24, 0.17, 0.21, 0.6,
0.79, 0.1, 0.05, 0.17, 0.42, 0.83, 0.81, 0.5, 0.26, 0.74,
0.47, 0.95, 0.81, 0.35, 0.49, 0.76, 0.95, 0.94, 0.24, 0.85,
0.42, 0.31, 0.46, 0.58, 0.85, 0.52, 0.41, 0.09, 0.39, 0.65,
0.54, 0.28, 0.18, 0.35, 0.11, 0.87, 0.18, 0.37, 0.36, 0.43,
0.66, 0.71, 0.83, 0.23, 0.45, 0.36, 0.93, 0.85, 0.73, 0.32,
0.65, 0.4, 0.34, 0.63, 0.15, 0.85, 0.21, 0.78, 0.68, 0.54,
0.28, 0.91, 0.73, 0.95, 0.23, 0.63, 0.74, 0.86, 0.93, 0.84,
0.55, 0.8, 0.63, 0.6, 0.88, 0.51, 0.47, 0.88, 0.31, 0.27,
0.73, 0.14, 0.88, 0.4, 0.95, 0.77, 0.32, 0.44, 0.5, 0.86,
0.74, 0.07, 0.88, 0.61, 0.79, 0.41, 0.07, 0.83, 0.55, 0.31,
0.37, 0.34, 0.66, 0.51, 0.78, 0.45, 0.35, 0.87, 0.51, 0.34,
0.26, 0.89, 0.76, 0.33, 0.9, 0.8, 0.76, 0.87, 0.53, 0.29,
0.9, 0.69, 0.54, 0.46, 0.14, 0.39, 0.85, 0.2, 0.47, 0.65,
0.41, 0.49, 0.59, 0.28, 0.52, 0.55, 0.76, 0.3, 0.56, 0.64,
0.05, 0.25, 0.86, 0.69, 0.64, 0.46, 0.18, 0.6, 0.88, 0.43,
0.48, 0.08, 0.7, 0.38, 0.95, 0.48, 0.55, 0.6, 0.44, 0.61,
0.05, 0.12, 0.44, 0.95, 0.23, 0.78, 0.92, 0.59, 0.82, 0.94,
0.58, 0.38, 0.5, 0.87, 0.27, 0.24, 0.53, 0.16, 0.76, 0.41,
0.91, 0.61, 0.32, 0.05, 0.18, 0.8, 0.15, 0.13, 0.1, 0.09,
0.56, 0.53, 0.33, 0.24, 0.21, 0.59, 0.95, 0.11, 0.79, 0.57,
0.91, 0.8, 0.13, 0.58, 0.8, 0.24, 0.08, 0.92, 0.24, 0.14,
0.22, 0.46, 0.82, 0.48, 0.64, 0.66, 0.91, 0.69, 0.75, 0.16,
0.63, 0.5, 0.85, 0.75, 0.58, 0.51, 0.74, 0.79, 0.61, 0.31,
0.42, 0.86, 0.58, 0.55, 0.43, 0.83, 0.49, 0.53, 0.69, 0.32,
0.66, 0.63, 0.09, 0.17, 0.56, 0.79, 0.95, 0.72, 0.79, 0.23,
0.8, 0.13, 0.16, 0.89, 0.91, 0.2, 0.35, 0.79, 0.45, 0.45,
0.9, 0.28, 0.2, 0.16, 0.28, 0.87, 0.55, 0.46, 0.08, 0.11,
0.22, 0.37, 0.81, 0.44, 0.24, 0.2, 0.23, 0.38, 0.63, 0.64,
0.75, 0.45, 0.34, 0.82, 0.94, 0.64, 0.19, 0.08, 0.77, 0.25,
0.18, 0.76, 0.08, 0.81, 0.62, 0.31, 0.28, 0.36, 0.39, 0.22,
0.13, 0.39), offset = 0, gain = 1, inmemory = TRUE, fromdisk = FALSE,
isfactor = FALSE, attributes = list(), haveminmax = TRUE,
min = 0.05, max = 0.95, band = 1L, unit = "", names = "X73147"),
legend = new(".RasterLegend", type = character(0), values = logical(0),
color = logical(0), names = logical(0), colortable = logical(0)),
title = character(0), extent = new("Extent", xmin = -78.125,
xmax = -76.25, ymin = 38.375, ymax = 39.5), rotated = FALSE,
rotation = new(".Rotation", geotrans = numeric(0), transfun = function ()
NULL), ncols = 30L, nrows = 18L, crs = new("CRS", projargs = "+proj=longlat +datum=WGS84 +no_defs"),
history = list(), z = list())
And T2:
new("RasterLayer", file = new(".RasterFile", name = "", datanotation = "INT4U",
byteorder = c(value = "little"), nodatavalue = 2147483647,
NAchanged = FALSE, nbands = 1L, bandorder = c(value = "BIL"),
offset = 0L, toptobottom = TRUE, blockrows = 0L, blockcols = 0L,
driver = "", open = FALSE), data = new(".SingleLayerData",
values = c(104L, 87L, 90L, 89L, 78L, 97L, 105L, 98L, 79L,
94L, 86L, 111L, 72L, 85L, 101L, 114L, 77L, 75L, 72L, 86L,
105L, 70L, 86L, 84L, 114L, 93L, 94L, 70L, 110L, 75L, 108L,
101L, 108L, 102L, 93L, 76L, 83L, 111L, 109L, 91L, 111L, 108L,
106L, 72L, 85L, 115L, 93L, 85L, 107L, 83L, 115L, 93L, 76L,
99L, 90L, 87L, 90L, 113L, 73L, 113L, 106L, 78L, 86L, 112L,
94L, 102L, 79L, 100L, 114L, 90L, 88L, 85L, 84L, 97L, 71L,
98L, 72L, 109L, 98L, 97L, 100L, 102L, 70L, 105L, 96L, 94L,
84L, 98L, 99L, 110L, 72L, 102L, 74L, 83L, 94L, 72L, 93L,
73L, 102L, 95L, 95L, 76L, 83L, 91L, 72L, 92L, 115L, 72L,
110L, 101L, 107L, 94L, 88L, 113L, 102L, 83L, 90L, 111L, 113L,
102L, 100L, 88L, 112L, 81L, 83L, 77L, 104L, 100L, 100L, 93L,
111L, 93L, 93L, 80L, 81L, 72L, 74L, 88L, 71L, 105L, 108L,
110L, 108L, 100L, 102L, 73L, 89L, 79L, 95L, 73L, 107L, 79L,
85L, 103L, 100L, 107L, 105L, 102L, 78L, 102L, 73L, 72L, 77L,
70L, 114L, 86L, 100L, 93L, 109L, 113L, 91L, 102L, 92L, 85L,
105L, 76L, 83L, 76L, 100L, 103L, 70L, 80L, 110L, 84L, 73L,
70L, 74L, 102L, 107L, 97L, 85L, 101L, 89L, 78L, 108L, 74L,
106L, 73L, 98L, 99L, 104L, 95L, 94L, 74L, 91L, 91L, 96L,
71L, 84L, 84L, 101L, 83L, 74L, 108L, 94L, 101L, 103L, 107L,
106L, 85L, 84L, 100L, 98L, 88L, 104L, 105L, 111L, 85L, 82L,
84L, 115L, 97L, 86L, 78L, 89L, 99L, 73L, 112L, 99L, 75L,
81L, 104L, 80L, 73L, 82L, 96L, 84L, 96L, 108L, 106L, 82L,
86L, 78L, 92L, 101L, 88L, 87L, 75L, 100L, 92L, 100L, 83L,
78L, 76L, 80L, 106L, 101L, 88L, 93L, 82L, 113L, 91L, 90L,
102L, 93L, 94L, 102L, 99L, 100L, 101L, 107L, 87L, 93L, 114L,
73L, 78L, 97L, 107L, 79L, 84L, 94L, 77L, 112L, 75L, 74L,
83L, 94L, 92L, 101L, 87L, 77L, 103L, 82L, 107L, 89L, 114L,
84L, 112L, 115L, 76L, 89L, 94L, 92L, 98L, 111L, 101L, 79L,
81L, 92L, 90L, 100L, 74L, 113L, 115L, 70L, 105L, 102L, 110L,
107L, 97L, 113L, 80L, 81L, 93L, 76L, 96L, 94L, 105L, 98L,
96L, 90L, 88L, 75L, 78L, 104L, 82L, 79L, 77L, 105L, 73L,
80L, 78L, 85L, 99L, 107L, 90L, 91L, 110L, 98L, 100L, 78L,
93L, 71L, 73L, 101L, 76L, 110L, 96L, 96L, 115L, 70L, 115L,
71L, 82L, 74L, 90L, 97L, 84L, 111L, 75L, 101L, 103L, 113L,
103L, 113L, 87L, 76L, 77L, 95L, 71L, 89L, 72L, 83L, 100L,
84L, 90L, 83L, 110L, 77L, 103L, 93L, 80L, 70L, 89L, 100L,
85L, 84L, 96L, 82L, 106L, 74L, 85L, 102L, 109L, 77L, 79L,
113L, 108L, 89L, 101L, 85L, 87L, 78L, 80L, 104L, 90L, 104L,
79L, 89L, 103L, 110L, 102L, 91L, 112L, 105L, 77L, 115L, 113L,
114L, 96L, 83L, 113L, 71L, 103L, 110L, 111L, 101L, 102L,
81L, 98L, 74L, 100L, 82L, 110L, 73L, 109L, 71L, 77L, 91L,
75L, 72L, 115L, 71L, 75L, 103L, 107L, 96L, 84L, 113L, 86L,
73L, 107L, 88L, 79L, 113L, 110L, 110L, 102L, 95L, 93L, 79L,
95L, 79L, 108L, 87L, 84L, 115L, 111L, 94L, 105L, 99L, 98L,
113L, 90L, 92L, 82L, 106L, 73L, 96L, 70L, 81L, 91L, 102L,
114L, 78L, 90L, 70L, 78L, 106L, 94L, 102L, 108L, 97L, 83L,
83L, 78L, 73L, 90L, 78L, 101L, 73L, 112L, 77L, 115L, 93L,
71L, 74L, 97L, 103L, 109L, 99L, 114L, 102L, 105L, 94L, 73L,
111L, 109L, 81L, 115L), offset = 0, gain = 1, inmemory = TRUE,
fromdisk = FALSE, isfactor = FALSE, attributes = list(),
haveminmax = TRUE, min = 70L, max = 115L, band = 1L, unit = "",
names = "X2000.04.09"), legend = new(".RasterLegend", type = character(0),
values = logical(0), color = logical(0), names = logical(0),
colortable = logical(0)), title = character(0), extent = new("Extent",
xmin = -78.125, xmax = -76.25, ymin = 38.375, ymax = 39.5),
rotated = FALSE, rotation = new(".Rotation", geotrans = numeric(0),
transfun = function ()
NULL), ncols = 30L, nrows = 18L, crs = new("CRS", projargs = "+proj=longlat +datum=WGS84 +no_defs"),
history = list(), z = list())