在地图上显示空间和时间相关性

机器算法验证 r 回归 数据可视化 主成分分析 空间的
2022-02-13 05:02:28

我有美国各地气象站网络的数据。这给了我一个包含日期、纬度、经度和一些测量值的数据框。假设每天收集一次数据,并由区域尺度的天气驱动(不,我们不打算进入那个讨论)。

我想以图形方式显示同时测量的值如何在时间和空间上相互关联。我的目标是展示正在调查的价值的区域同质性(或缺乏同质性)。

数据集

首先,我在马萨诸塞州和缅因州地区进行了一组站点。我从 NOAA 的 FTP 站点上提供的索引文件中按纬度和经度选择站点。

在此处输入图像描述

您会立即看到一个问题:有很多站点具有相似的标识符或非常接近。FWIW,我使用 USAF 和 WBAN 代码识别它们。深入研究元数据,我发现它们具有不同的坐标和高度,数据在一个站点停止,然后在另一个站点开始。所以,因为我不知道更好,我不得不将它们视为单独的站。这意味着数据包含彼此非常接近的成对站点。

初步分析

我尝试按日历月对数据进行分组,然后计算不同数据对之间的普通最小二乘回归。然后,我将所有对之间的相关性绘制为连接站点的线(下图)。线条颜色显示 OLS 拟合的 R2 值。然后,该图显示了 1 月、2 月等的 30 多个数据点如何在感兴趣区域的不同站点之间相互关联。

每个日历月的每日数据之间的相关性

我已经编写了基础代码,以便仅在每 6 小时有一个数据点时才计算每日平均值,因此数据应该在各个站点之间具有可比性。

问题

不幸的是,一张图上的数据太多了。这不能通过减小线条的大小来解决。

我尝试绘制该地区最近邻居之间的相关性,但这很快就变成了一团糟。下面的分面显示了没有相关值的网络,使用个最近邻居。这个数字只是为了测试这个概念。 k在此处输入图像描述

网络似乎太复杂了,所以我想我需要想办法降低复杂性,或者应用某种空间内核。

我也不确定什么是最适合显示相关性的指标,但对于目标(非技术)受众来说,OLS 的相关系数可能是最容易解释的。我可能还需要提供一些其他信息,例如梯度或标准误差。

问题

我正在同时学习进入这个领域和 R,并希望得到以下方面的建议:

  1. 我正在尝试做的更正式的名称是什么?是否有一些有用的术语可以让我找到更多的文献?我的搜索是为必须是一个常见的应用程序绘制空白。
  2. 是否有更合适的方法来显示空间分隔的多个数据集之间的相关性?
  3. ...特别是易于从视觉上显示结果的方法?
  4. 这些是否在 R 中实现?
  5. 这些方法中的任何一种都适合自动化吗?
2个回答

我认为有几个选项可以显示这种类型的数据:

第一种选择是进行“经验正交函数分析”(EOF)(在非气候圈中也称为“主成分分析”(PCA))。对于您的情况,这应该在您的数据位置的相关矩阵上进行。例如,您的数据矩阵dat将是您在列维度中的空间位置,以及在行中的测量参数;因此,您的数据矩阵将包含每个位置的时间序列。prcomp()函数将允许您获得与该字段相关的主成分或相关的主要模式:

res <- prcomp(dat, retx = TRUE, center = TRUE, scale = TRUE) # center and scale should be "TRUE" for an analysis of dominant correlation modes)
#res$x and res$rotation will contain the PC modes in the temporal and spatial dimension, respectively.

第二种选择是创建显示与单个感兴趣位置相关的地图:

C <- cor(dat)
#C[,n] would be the correlation values between the nth location (e.g. dat[,n]) and all other locations. 

编辑:附加示例

虽然以下示例不使用 gappy 数据,但您可以在使用 DINEOF ( http://menugget.blogspot.de/2012/10/dineof-data-interpolating-empirical.html )插值后对数据字段应用相同的分析. 下面的示例使用来自以下数据集 ( http://www.esrl.noaa.gov/psd/gcos_wgsp/Gridded/data.hadslp2.html )的月度异常海平面压力数据的子集:

library(sinkr) # https://github.com/marchtaylor/sinkr

# load data
data(slp)

grd <- slp$grid
time <- slp$date
field <- slp$field

# make anomaly dataset
slp.anom <- fieldAnomaly(field, time)

# EOF/PCA of SLP anom
P <- prcomp(slp.anom, center = TRUE, scale. = TRUE)

expl.var <- P$sdev^2 / sum(P$sdev^2) # explained variance
cum.expl.var <- cumsum(expl.var) # cumulative explained variance
plot(cum.expl.var)

映射领先的 EOF 模式

# make interpolation
require(akima)
require(maps)

eof.num <- 1
F1 <- interp(x=grd$lon, y=grd$lat, z=P$rotation[,eof.num]) # interpolated spatial EOF mode


png(paste0("EOF_mode", eof.num, ".png"), width=7, height=6, units="in", res=400)
op <- par(ps=10) #settings before layout
layout(matrix(c(1,2), nrow=2, ncol=1, byrow=TRUE), heights=c(4,2), widths=7)
#layout.show(2) # run to see layout; comment out to prevent plotting during .pdf
par(cex=1) # layout has the tendency change par()$cex, so this step is important for control

par(mar=c(4,4,1,1)) # I usually set my margins before each plot
pal <- jetPal
image(F1, col=pal(100))
map("world", add=TRUE, lwd=2)
contour(F1, add=TRUE, col="white")
box()

par(mar=c(4,4,1,1)) # I usually set my margins before each plot
plot(time, P$x[,eof.num], t="l", lwd=1, ylab="", xlab="")
plotRegionCol()
abline(h=0, lwd=2, col=8)
abline(h=seq(par()$yaxp[1], par()$yaxp[2], len=par()$yaxp[3]+1), col="white", lty=3)
abline(v=seq.Date(as.Date("1800-01-01"), as.Date("2100-01-01"), by="10 years"), col="white", lty=3)
box()
lines(time, P$x[,eof.num])
mtext(paste0("EOF ", eof.num, " [expl.var = ", round(expl.var[eof.num]*100), "%]"), side=3, line=1) 

par(op)
dev.off() # closes device

在此处输入图像描述

创建相关图

loc <- c(-90, 0)
target <- which(grd$lon==loc[1] & grd$lat==loc[2])
COR <- cor(slp.anom)
F1 <- interp(x=grd$lon, y=grd$lat, z=COR[,target]) # interpolated spatial EOF mode


png(paste0("Correlation_map", "_lon", loc[1], "_lat", loc[2], ".png"), width=7, height=5, units="in", res=400)

op <- par(ps=10) #settings before layout
layout(matrix(c(1,2), nrow=2, ncol=1, byrow=TRUE), heights=c(4,1), widths=7)
#layout.show(2) # run to see layout; comment out to prevent plotting during .pdf
par(cex=1) # layout has the tendency change par()$cex, so this step is important for control

par(mar=c(4,4,1,1)) # I usually set my margins before each plot
pal <- colorRampPalette(c("blue", "cyan", "yellow", "red", "yellow", "cyan", "blue"))
ncolors <- 100
breaks <- seq(-1,1,,ncolors+1)
image(F1, col=pal(ncolors), breaks=breaks)
map("world", add=TRUE, lwd=2)
contour(F1, add=TRUE, col="white")
box()

par(mar=c(4,4,0,1)) # I usually set my margins before each plot
imageScale(F1, col=pal(ncolors), breaks=breaks, axis.pos = 1)
mtext("Correlation [R]", side=1, line=2.5)
box()

par(op)

dev.off() # closes device

在此处输入图像描述

我看不清楚线条后面,但在我看来,数据点太多了。

由于您想显示区域同质性而不是确切的站点,我建议您首先将它们按空间分组。例如,用“渔网”覆盖并计算每个单元格中的平均测量值(在每个时刻)。如果您以这种方式将这些平均值放置在单元格中心,则可以对数据进行栅格化(或者,如果您不想覆盖线,也可以计算每个单元格中的平均纬度和经度)。或者平均内部行政单位,无论如何。然后对于这些新的平均“站”,您可以计算相关性并绘制具有较少线数的地图。

在此处输入图像描述

这也可以去除那些穿过所有区域的随机单条高相关线。