Affichage de la corrélation spatiale et temporelle sur les cartes

J'ai des données pour un réseau de stations météorologiques à travers les États-Unis. Cela me donne un bloc de données qui contient la date, la latitude, la longitude et une certaine valeur mesurée. Supposons que les données soient collectées une fois par jour et dictées par la météo à l'échelle régionale (non, nous n'allons pas entrer dans cette discussion).

Je voudrais montrer graphiquement comment les valeurs mesurées simultanément sont corrélées dans le temps et l'espace. Mon objectif est de montrer l'homogénéité régionale (ou son absence) de la valeur étudiée.

Base de données

Pour commencer, j'ai pris un groupe de stations dans la région du Massachusetts et du Maine. J'ai sélectionné les sites par latitude et longitude à partir d'un fichier d'index qui est disponible sur le site FTP de la NOAA.

Vous voyez tout de suite un problème: il y a beaucoup de sites qui ont des identifiants similaires ou qui sont très proches. FWIW, je les identifie à l'aide des codes USAF et WBAN. En regardant plus profondément dans les métadonnées, j'ai vu qu'elles ont des coordonnées et des élévations différentes, et les données s'arrêtent sur un site puis commencent sur un autre. Donc, parce que je ne sais pas mieux, je dois les traiter comme des stations distinctes. Cela signifie que les données contiennent des paires de stations très proches les unes des autres.

Analyse préliminaire

J'ai essayé de regrouper les données par mois civil, puis de calculer la régression des moindres carrés ordinaires entre différentes paires de données. Je trace ensuite la corrélation entre toutes les paires comme une ligne reliant les stations (ci-dessous). La couleur de la ligne montre la valeur de R2 de l'ajustement OLS. La figure montre ensuite comment les 30+ points de données de janvier, février, etc. sont corrélés entre différentes stations dans la zone d'intérêt.

J'ai écrit les codes sous-jacents afin que la moyenne quotidienne ne soit calculée que s'il y a des points de données toutes les 6 heures, les données doivent donc être comparables d'un site à l'autre.

Problèmes

Malheureusement, il y a tout simplement trop de données pour donner un sens à une parcelle. Cela ne peut pas être résolu en réduisant la taille des lignes.

$k$

Le réseau semble être trop complexe, donc je pense que je dois trouver un moyen de réduire la complexité ou d'appliquer une sorte de noyau spatial.

Je ne sais pas non plus quelle est la mesure la plus appropriée pour montrer la corrélation, mais pour le public visé (non technique), le coefficient de corrélation d'OLS pourrait être le plus simple à expliquer. Il se peut que je doive également présenter d'autres informations comme le gradient ou l'erreur standard.

Des questions

J'apprends mon chemin dans ce domaine et R en même temps, et j'apprécierais des suggestions sur:

1. Quel est le nom le plus formel de ce que j'essaie de faire? Existe-t-il des termes utiles qui me permettraient de trouver plus de documentation? Mes recherches dessinent des blancs pour ce qui doit être une application courante.
2. Existe-t-il des méthodes plus appropriées pour montrer la corrélation entre plusieurs ensembles de données séparés dans l'espace?
3. ... en particulier, des méthodes dont il est facile de montrer visuellement les résultats?
4. Y en a-t-il dans R?
5. L'une de ces approches se prête-t-elle à l'automatisation?

Je pense qu'il y a quelques options pour montrer ce type de données:

La première option consisterait à effectuer une «analyse des fonctions orthogonales empiriques» (EOF) (également appelée «analyse en composantes principales» (ACP) dans les cercles non climatiques). Pour votre cas, cela doit être effectué sur une matrice de corrélation de vos emplacements de données. Par exemple, votre matrice de données datserait vos emplacements spatiaux dans la dimension de colonne et le paramètre mesuré dans les lignes; Ainsi, votre matrice de données sera constituée de séries chronologiques pour chaque emplacement. La prcomp()fonction vous permettra d'obtenir les principales composantes, ou modes de corrélation dominants, relatifs à ce domaine:

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.

La deuxième option serait de créer des cartes qui montrent une corrélation par rapport à un emplacement individuel d'intérêt:

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

EDIT: exemple supplémentaire

Bien que l'exemple suivant n'utilise pas de données gappy, vous pouvez appliquer la même analyse à un champ de données après une interpolation avec DINEOF ( http://menugget.blogspot.de/2012/10/dineof-data-interpolating-empirical.html ) . L'exemple ci-dessous utilise un sous-ensemble de données mensuelles de pression au niveau de la mer des anomalies provenant de l'ensemble de données suivant ( 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)

Mappez le principal mode 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

Créer une carte de corrélation

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

Je ne vois pas clairement derrière les lignes mais il me semble qu'il y a trop de points de données.

Puisque vous voulez montrer l'homogénéité régionale et pas exactement les stations, je vous suggère tout d'abord de les regrouper spatialement. Par exemple, superposer un "filet de pêche" et calculer la valeur mesurée moyenne dans chaque cellule (à chaque instant). Si vous placez ces valeurs moyennes dans les centres des cellules de cette façon, vous pixellisez les données (ou vous pouvez également calculer la latitude et la longitude moyennes dans chaque cellule si vous ne voulez pas superposer des lignes). Ou de faire la moyenne à l'intérieur des unités administratives, peu importe. Ensuite, pour ces nouvelles "stations" moyennes, vous pouvez calculer les corrélations et tracer une carte avec un plus petit nombre de lignes.

Cela peut également supprimer ces lignes de haute corrélation unique aléatoires traversant toute la zone.