我觉得这应该是一个简单的修复程序，但无法解决。我有两个光栅，一个有温度数据，另一个有湿度数据。我需要使用一组方程来计算热指数，并试图获得光栅输出。我把热指数方程放入一个函数中，如下所示(可以忽略方程本身，我知道这是很多数字，这不是我的问题所在(：
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
这就是我试图使用上面的函数运行的内容，该函数只是出现了一个错误，没有运行完。fun_heat_index(relative_humidity_data, tmax_data)
我的光栅是相同的范围，每个都有相同数量的网格单元。我正在尝试应用该函数从两个光栅进行计算，以获得具有热指数值的第三个输出光栅。我找到了calc()函数，它使它单独通过每个网格单元，我认为这会解决问题。我尝试了下一个：
calc(tmax_data, fun = fun_heat_index(relative_humidity_data, tmax_data))
在我添加之后，它开始工作了一点，但在最初创建iheat_index变量之后，它似乎仍然会出错，这是函数的第一步。它将从该行输出一个光栅，但一旦它进入if语句if (iheat_index > 80){的下一行，它就会说参数不能解释为逻辑，我认为这是因为它期望的是一个值而不是光栅。有没有办法在if语句中再次应用calc()函数？或者，有没有另一种方法可以让R只看一些网格单元，并将额外的方程仅应用于满足特定标准的网格单元？
在四处寻找如何做到这一点时，我也找到了overlay()，并尝试如下：
overlay(iloop_tmax, iloop_rh, fun = fun_heat_index(iloop_rh, iloop_tmax))
这也没用。我还尝试过遍历每个单元格，使用getValues()删除值，并将该函数应用于符合if语句条件的任何单元格，但我不知道如何将其恢复为光栅形式而不是表格形式。
非常感谢您的帮助！：(
我试图让简单的事情一步一步地进行，因为这会让我的许多错误更容易检测和调试。将您上面删除的以.gri结尾的文件作为T2&RH2，然后通过在每个中采样创建一个值范围
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
转到heatindex_equiation，我们可以看到，传统上，最简单的方程首先计算，这个结果用T2:进行平均
HI_tst1 = 0.5 * (T2 + 61.0 +((T2-68.0) * 1.2) + (RH2 * 0.094))
HI_tst1_avg = 0.5 * (T2 + HI_tst1)
因此，HI_tst1_avg似乎是我们的基本数据集，根据所述标准，将对该数据集中的特定单元格进行进一步的计算和调整。更准确的说法可能是，HI_tst1_avg是将应用进一步计算的单元格的索引(概念上(。所以现在，也许我们在我们的平均值范围内进行一些研究，以了解可能应该进行的计算和调整：
# ? 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
现在，我想我有了合适的索引，可以开始插入创建热索引光栅。所以，用NA值创建它。
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)
按照建议的顺序工作，1(简单2(回归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
正如他们所说，"热会致命"。还有很多问题要弄清楚为什么调整1不起作用，但在制作一个包含这些步骤的函数之前，我会经历这些步骤。我还发现，当几个月后再次访问时，这些base方法更容易理解和推理。
调试错误意味着首先扭转在adj_idx单元格中引入NAN的趋势，篡改)的位置以获得合理的值
# 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)
数据：
# 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())
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())