Darstellung der räumlichen und zeitlichen Korrelation auf Karten

16

Ich habe Daten für ein Netzwerk von Wetterstationen in den USA. Dies gibt mir einen Datenrahmen, der Datum, Breite, Länge und einige Messwerte enthält. Angenommen, die Daten werden einmal pro Tag erfasst und sind abhängig vom regionalen Wetter (nein, wir werden nicht auf diese Diskussion eingehen).

Ich möchte grafisch zeigen, wie zeitlich und räumlich gleichzeitig gemessene Werte korrelieren. Mein Ziel ist es, die regionale Homogenität (oder das Fehlen derselben) des untersuchten Wertes aufzuzeigen.

Datensatz

Zunächst nahm ich eine Gruppe von Stationen in der Region Massachusetts und Maine. Ich habe Websites nach Breiten- und Längengrad aus einer Indexdatei ausgewählt, die auf der FTP-Site von NOAA verfügbar ist.

Bildbeschreibung hier eingeben

Auf Anhieb sehen Sie ein Problem: Es gibt viele Websites, die ähnliche Bezeichner aufweisen oder sehr nahe beieinander liegen. FWIW, ich identifiziere sie mit den USAF- und WBAN-Codes. Bei genauerem Hinsehen der Metadaten stellte ich fest, dass sie unterschiedliche Koordinaten und Höhen haben und die Daten an einer Stelle anhalten und dann an einer anderen beginnen. Da ich es nicht besser weiß, muss ich sie als separate Stationen behandeln. Dies bedeutet, dass die Daten Paare von Stationen enthalten, die sehr nahe beieinander liegen.

Voruntersuchung

Ich habe versucht, die Daten nach Kalendermonaten zu gruppieren und dann die gewöhnliche Regression der kleinsten Quadrate zwischen verschiedenen Datenpaaren zu berechnen. Ich zeichne dann die Korrelation zwischen allen Paaren als eine Linie, die die Stationen verbindet (unten). Die Linienfarbe zeigt den Wert von R2 aus der OLS-Anpassung. Die Abbildung zeigt dann, wie die über 30 Datenpunkte von Januar, Februar usw. zwischen verschiedenen Stationen im interessierenden Bereich korreliert sind.

Korrelation zwischen Tagesdaten während jedes Kalendermonats

Ich habe die zugrunde liegenden Codes so geschrieben, dass der Tagesmittelwert nur berechnet wird, wenn alle 6 Stunden Datenpunkte vorliegen. Daher sollten die Daten standortübergreifend vergleichbar sein.

Probleme

Leider gibt es einfach zu viele Daten, um auf einem Plot einen Sinn zu ergeben. Dies kann nicht durch Verringern der Zeilengröße behoben werden.

kBildbeschreibung hier eingeben

Das Netzwerk scheint zu komplex zu sein, daher denke ich, dass ich einen Weg finden muss, um die Komplexität zu reduzieren, oder eine Art räumlichen Kernel anzuwenden.

Ich bin mir auch nicht sicher, welche Metrik am besten geeignet ist, um die Korrelation darzustellen, aber für das beabsichtigte (nicht technische) Publikum ist der Korrelationskoeffizient von OLS möglicherweise am einfachsten zu erklären. Möglicherweise muss ich auch andere Informationen wie den Gradienten oder den Standardfehler angeben.

Fragen

Ich lerne mich gleichzeitig in dieses Feld und in R ein und würde mich über Vorschläge freuen zu:

  1. Was ist der formalere Name für das, was ich versuche zu tun? Gibt es einige hilfreiche Begriffe, mit denen ich mehr Literatur finden könnte? Meine Suche ist das Zeichnen von Leerzeichen für eine übliche Anwendung.
  2. Gibt es geeignetere Methoden, um die Korrelation zwischen mehreren im Raum getrennten Datensätzen darzustellen?
  3. ... insbesondere Methoden, die sich visuell leicht darstellen lassen?
  4. Sind einige davon in R implementiert?
  5. Bietet sich einer dieser Ansätze für die Automatisierung an?
Andy Clifton
quelle
[Räumliche Beschreibung der zeitlichen Korrelation in einer Visual Analytics-Umgebung, "Abish Malik et al.] [1] [1]: google.com/…
pat
2
yWy
Was ist, wenn Sie versuchen, den Plot-Schwellenwert (0,5) zu erhöhen und mehr als 4 Farbschritte zu verwenden? Oder dünnere Linien anstelle von Farben verwenden.
Nadya
norder((n2)/2)
1
Daraus ist mir klar geworden, dass ich viel Arbeit in der Vorverarbeitung der Daten habe, bevor ich mit der Analyse beginne, die ich hier skizziert habe. Wenn ich die Antwort von @nadya lese, denke ich, dass es klar ist, dass ich eine Art räumliche Aggregation betrachten muss, aber das wird eine Herausforderung sein, da es falsch ist, Land- und Ozeandaten zu aggregieren. Dann muss ich mir Strategien zum Füllen von Lücken ansehen. Dann (und nur dann) kann ich anfangen, mir die Mapping- / Visualisierungsarbeit anzuschauen.
Andy Clifton

Antworten:

10

Ich denke, es gibt ein paar Möglichkeiten, um diese Art von Daten anzuzeigen:

Die erste Option wäre die Durchführung einer "Empirical Orthogonal Functions Analysis" (EOF) (auch als "Principal Component Analysis" (PCA) in Nicht-Klimakreisen bezeichnet). Für Ihren Fall sollte dies auf einer Korrelationsmatrix Ihrer Datenpositionen durchgeführt werden. Beispielsweise datwäre Ihre Datenmatrix Ihre räumlichen Positionen in der Spaltendimension und der gemessene Parameter in den Zeilen. Ihre Datenmatrix enthält also Zeitreihen für jeden Standort. Mit dieser prcomp()Funktion können Sie die Hauptkomponenten oder dominanten Korrelationsmodi für dieses Feld ermitteln:

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.

Die zweite Möglichkeit wäre, Karten zu erstellen, die die Korrelation zu einem bestimmten Ort von Interesse anzeigen:

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

EDIT: zusätzliches Beispiel

Im folgenden Beispiel werden keine Gappy-Daten verwendet, Sie können jedoch nach der Interpolation mit DINEOF dieselbe Analyse auf ein Datenfeld anwenden ( http://menugget.blogspot.de/2012/10/dineof-data-interpolating-empirical.html ). . Das folgende Beispiel verwendet eine Teilmenge der monatlichen Anomaliedruckdaten aus dem folgenden Datensatz ( 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)

Ordnen Sie den führenden EOF-Modus zu

# 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

Bildbeschreibung hier eingeben

Erstellen Sie eine Korrelationskarte

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

Bildbeschreibung hier eingeben

Marc in der Kiste
quelle
Wie gut gehen diese Funktionen mit fehlenden Daten um? Ich habe ziemlich oft Lücken in der Zeitreihe.
Andy Clifton
2
Es gibt EOF-Methoden, die für den von Ihnen beschriebenen Spezialfall "Gappy Data" entwickelt wurden. Hier ist ein Link zu einem Artikel , der diese Methoden behandelt: dx.doi.org/10.6084/m9.figshare.732650 . Sie werden sehen, dass die RSEOF- und DINEOF-Methoden am genauesten zum Ableiten von EOFs aus Gappy-Datensätzen geeignet sind. Der DINEOF-Interpolationsalgorithmus ist hier zu finden: menugget.blogspot.de/2012/10/…
Marc in der Box
1
Ich denke, dies ist die beste Antwort auf eine schreckliche Frage (im Nachhinein).
Andy Clifton
3

Ich sehe nicht klar hinter den Linien, aber es scheint mir, dass es zu viele Datenpunkte gibt.

Da Sie die regionale Homogenität und nicht genau die Stationen zeigen möchten, empfehle ich Ihnen, diese zunächst räumlich zu gruppieren. Beispiel: Überlagerung mit einem "Fischnetz" und Berechnung des durchschnittlichen Messwerts in jeder Zelle (zu jedem Zeitpunkt). Wenn Sie diese Durchschnittswerte auf diese Weise in die Zellmitten platzieren, rastern Sie die Daten (oder Sie können auch den mittleren Breiten- und Längengrad in jeder Zelle berechnen, wenn Sie keine überlagernden Linien möchten). Oder um in Verwaltungseinheiten zu arbeiten, was auch immer. Dann können Sie für diese neuen gemittelten "Stationen" Korrelationen berechnen und eine Karte mit einer geringeren Anzahl von Linien zeichnen.

Bildbeschreibung hier eingeben

Dies kann auch jene zufälligen einzelnen Linien mit hoher Korrelation entfernen, die durch den gesamten Bereich verlaufen.

Nadya
quelle
x×xx
Ja, die Koordinaten zu projizieren ist eine gute Idee. Viel Glück!
Nadya