Pokazywanie korelacji przestrzennej i czasowej na mapach

Mam dane dla sieci stacji pogodowych w Stanach Zjednoczonych. To daje mi ramkę danych, która zawiera datę, szerokość, długość i pewną zmierzoną wartość. Załóżmy, że dane są gromadzone raz dziennie i zależą od pogody w skali regionalnej (nie, nie będziemy wchodzić w tę dyskusję).

Chciałbym pokazać graficznie, jak jednocześnie mierzone wartości są skorelowane w czasie i przestrzeni. Moim celem jest pokazanie regionalnej jednorodności (lub jej braku) badanej wartości.

Zbiór danych

Na początek wziąłem grupę stacji w regionie Massachusetts i Maine. Wybrałem witryny według szerokości i długości geograficznej z pliku indeksu, który jest dostępny na stronie FTP NOAA.

Od razu widzisz jeden problem: istnieje wiele witryn, które mają podobne identyfikatory lub są bardzo blisko. FWIW, identyfikuję je za pomocą kodów USAF i WBAN. Zagłębiając się w metadane, zauważyłem, że mają one różne współrzędne i wysokości, a dane zatrzymują się w jednym miejscu, a następnie zaczynają w innym. Ponieważ nie wiem nic lepszego, muszę traktować je jako osobne stacje. Oznacza to, że dane zawierają pary stacji, które są bardzo blisko siebie.

Wstępna analiza

Próbowałem pogrupować dane według miesiąca kalendarzowego, a następnie obliczyć regresję metodą najmniejszych kwadratów między różnymi parami danych. Następnie rysuję korelację między wszystkimi parami jako linię łączącą stacje (poniżej). Kolor linii pokazuje wartość R2 z dopasowania OLS. Na rysunku pokazano następnie, w jaki sposób ponad 30 punktów danych ze stycznia, lutego itp. Jest skorelowanych między różnymi stacjami w obszarze zainteresowania.

Napisałem podstawowe kody, aby średnia dzienna była obliczana tylko wtedy, gdy istnieją punkty danych co 6 godzin, więc dane powinny być porównywalne w różnych witrynach.

Problemy

Niestety, jest po prostu zbyt wiele danych, aby można je było zrozumieć na jednym wykresie. Nie można tego naprawić, zmniejszając rozmiar linii.

Próbowałem wykreślić korelacje między najbliższymi sąsiadami w regionie, ale to bardzo szybko zmienia się w bałagan. Poniższe aspekty pokazują sieć bez wartości korelacji, używając najbliższych sąsiadów z podzbioru stacji. Ta liczba miała tylko przetestować koncepcję. $k$

Sieć wydaje się być zbyt złożona, więc myślę, że muszę znaleźć sposób na zmniejszenie złożoności lub zastosowanie jakiegoś przestrzennego jądra.

Nie jestem również pewien, która metoda jest najbardziej odpowiednia do pokazania korelacji, ale dla zamierzonej (nietechnicznej) grupy odbiorców współczynnik korelacji z OLS może być najprostszy do wyjaśnienia. Może być konieczne przedstawienie innych informacji, takich jak gradient lub błąd standardowy.

pytania

W tym samym czasie uczę się tej dziedziny i R. Docenię sugestie dotyczące:

1. Jaka jest bardziej formalna nazwa tego, co próbuję zrobić? Czy są jakieś przydatne terminy, które pozwoliłyby mi znaleźć więcej literatury? Moje wyszukiwania rysują puste miejsca dla tego, co musi być powszechną aplikacją.
2. Czy istnieją bardziej odpowiednie metody pokazania korelacji między wieloma zestawami danych oddzielonymi w przestrzeni?
3. ... w szczególności metody, które łatwo pokazać wizualnie?
4. Czy którekolwiek z nich są zaimplementowane w języku R?
5. Czy któreś z tych podejść nadaje się do automatyzacji?

Myślę, że istnieje kilka opcji wyświetlania tego typu danych:

Pierwszą opcją byłoby przeprowadzenie „Empirycznej analizy funkcji ortogonalnych” (EOF) (zwanej również „główną analizą składową” (PCA) w kręgach innych niż klimatyczne). W twoim przypadku należy to przeprowadzić na macierzy korelacji lokalizacji danych. Na przykład macierzą danych datbędą lokalizacje przestrzenne w wymiarze kolumny i zmierzony parametr w wierszach; Zatem macierz danych będzie się składać z szeregów czasowych dla każdej lokalizacji. Ta prcomp()funkcja umożliwia uzyskanie głównych składników lub dominujących trybów korelacji dotyczących tego pola:

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.

Drugą opcją byłoby stworzenie map pokazujących korelację w stosunku do indywidualnej interesującej lokalizacji:

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

EDYCJA: dodatkowy przykład

Podczas gdy w poniższym przykładzie nie używa się nieciekawych danych, można zastosować tę samą analizę do pola danych po interpolacji za pomocą DINEOF ( http://menugget.blogspot.de/2012/10/dineof-data-interpolating-empirical.html ) . W poniższym przykładzie wykorzystano podzbiór miesięcznych danych dotyczących ciśnienia na poziomie morza z anomalii z następującego zestawu danych ( 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)

Odwzoruj wiodący tryb 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

Utwórz mapę korelacji

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

Nie widzę wyraźnie za liniami, ale wydaje mi się, że jest zbyt wiele punktów danych.

Ponieważ chcesz pokazać regionalną jednorodność, a nie dokładnie stacje, sugeruję najpierw pogrupować je przestrzennie. Na przykład nałóż „kabaretkę” i oblicz średnią zmierzoną wartość w każdej komórce (w każdej chwili). Jeśli umieścisz te średnie wartości w centrach komórek w ten sposób, zrasteryzujesz dane (lub możesz obliczyć średnią szerokość i długość geograficzną w każdej komórce, jeśli nie chcesz nakładać linii). Lub uśredniać wewnątrz jednostek administracyjnych, cokolwiek. Następnie dla tych nowych uśrednionych „stacji” możesz obliczyć korelacje i wykreślić mapę z mniejszą liczbą linii.

Może to również usunąć losowe pojedyncze linie o wysokiej korelacji przechodzące przez cały obszar.