在科学或实验中，探究两个变量之间的关系是一项常见任务。通过相关性测试，可以获得相关系数和p值。相关系数表明关系的强度，而p值则告诉相关性测试是否显著。基于这些结果，可以得出一个因变量是否影响目标变量的结论。本文将讨论如何对表格数据进行相关性测试，特别是空间栅格数据。栅格数据是具有空间属性的图像，对空间栅格数据进行相关性测试与表格数据类似。
假设有两组栅格数据：(1) NDVI（归一化植被指数）和(2) 亮度温度（BT）。想要测试这两者之间的关系，以了解植被是否影响温度。NDVI代表植被状况，NDVI值高于0.6通常表示高密度植被，如森林；NDVI值在0.4到0.6之间通常表示农作物；NDVI值在0.2到0.4之间通常表示土壤、裸地或房屋；NDVI值为负数通常表示水体。
首先，需要将这两个栅格数据集叠加在一起。接下来，将它们从栅格转换为表格数据。每一对单元格或像素将并排排列在表格中。因此，表格有两列，分别对应NDVI和BT。表格的行数与单元格数量相同。之后，就可以运行相关性测试了。
但是，并不总是能得到具有相同分辨率或单元格大小、覆盖区域的栅格数据集。在这个例子中，NDVI是使用Landsat 8 OLI的Band 5和Band 4处理的，单元格大小为30米。亮度温度的单元格大小为90米。这两个栅格数据集具有不同的单元格大小，并且覆盖区域也不同。将它们叠加在一起后，它们不能被转换为数据框，因为单元格之间不匹配。图1展示了两组栅格图像。黑色栅格的分辨率高于红色栅格。单元格不匹配。
这是空间分析的挑战。在运行相关性测试之前，必须处理这两个栅格图像，使它们具有相同的单元格大小和覆盖区域。NDVI栅格的每个单元格必须与BT栅格的每个单元格分别对应。为了做到这一点，在叠加栅格图像之前，需要根据另一个栅格的单元格大小对其中一个栅格数据进行重采样。现在的选择是，是按照BT栅格对NDVI栅格进行重采样，还是反之。
如果根据BT栅格对NDVI栅格进行重采样，那么30米的单元格大小将被泛化到90米。其中将丢失一些信息。如果将分辨率为90米的BT栅格重采样到30米分辨率，将会产生插值的合成信息。在空间数据中，更倾向于将复杂的数据泛化为更简单的数据，而不是反之。因此，在这种情况下，将对NDVI栅格进行重采样，使单元格分辨率为90米。
在“raster”包中有两种重采样方法。第一种方法“ngb”用于根据“最近邻”对分类值进行重采样。另一种方法是“双线性”插值，用于对连续值进行重采样。在这种情况下，使用双线性插值。
“ngb”方法创建一个新的栅格，其单元格值根据旧栅格的最近值确定。双线性插值创建一个新的栅格，其单元格值基于旧栅格的4个最近值。请观察图3。创建了一个新的单元格。新栅格的值在蓝点上。有四个单元格围绕着蓝点。它们是7、9、1和11。蓝点的值是4个周围值的加权平均值，根据距离而定。
重采样NDVI图像后，现在可以将两个图像转换为具有配对单元格的数据框。下一个挑战是解决无值单元格的问题。有些区域有NDVI值，但没有BT值。也有一些区域有BT值，但没有NDVI值。不希望这些区域被包含在相关性测试中。通常，无值区域用特殊数字填充，如0、9999或-9999。这将破坏相关性测试。
为了过滤掉这些区域，需要用彼此的掩膜处理两个栅格数据集。这意味着，对于NDVI或BT无值的单元格将被移除。这与交集相同。只返回两个栅格数据集都有可用值的单元格。以下是执行栅格掩膜的代码。
# 掩膜栅格
ndvi_m <- mask(ndvi, bt)
bt_m <- mask(bt, ndvi)
# 绘图
plot(ndvi_m)
plot(bt_m)
现在重采样和掩膜已经完成，两个栅格数据集可以转换为数据框。如前所述，数据框有一列用于NDVI值，另一列用于BT值。每一行是配对的单元格。现在可以进行相关性测试以测试关系。在这种情况下，使用皮尔逊相关性，因为两个变量都是连续的。以下是获取相关系数和p值的代码。
# 叠加两个栅格并创建数据框
overlay <- stack(ndvi_m, bt_m)
overlay <- data.frame(na.omit(values(overlay)))
names(overlay2) <- c("ndvi", "bt")
# 相关性测试
cor.test(overlay[,1], overlay[,2])
输出结果如下：
Pearson's product-moment correlation
data: overlay[, 1] and overlay[, 2]
t = 6.1095, df = 31557, p-value = 1.011e-09
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.02334775 0.04538748
sample estimates:
cor
0.0343718
以下是数据的前7行：
ndvi
bt
1
0.16652
31.49438
2
0.1741
31.42719
3
0.223416
31.60288
4
0.273133
32.10332
5
0.28525
32.57689
6
0.250479
32.85951
7
0.23371
33.02891
可以看到，相关系数是0.0343718，p值是1.011e-09。这意味着NDVI不影响亮度温度。
还可以进行其他类型的分析。例如，让以NDVI作为自变量，BT作为目标变量进行线性回归。以下是代码。
# 线性回归（备选）
linear <- lm(overlay[,1] ~ overlay[,2])
summary(linear)
Call:
lm(formula = overlay[, 1] ~ overlay[, 2])
Residuals:
Min 1Q Median 3Q Max
-0.48114 -0.14726 0.01716 0.16513 0.28196
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.4453772 0.0091248 48.81 < 2e-16 ***
overlay[, 2] 0.0017075 0.0002795 6.11 1.01e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1733 on 31557 degrees of freedom
Multiple R-squared: 0.001181, Adjusted R-squared: 0.00115
F-statistic: 37.33 on 1 and 31557 DF, p-value: 1.011e-09