Menampilkan korelasi spasial dan temporal pada peta

16

Saya memiliki data untuk jaringan stasiun cuaca di seluruh Amerika Serikat. Ini memberi saya bingkai data yang berisi tanggal, lintang, bujur, dan beberapa nilai yang diukur. Asumsikan bahwa data dikumpulkan sekali sehari dan didorong oleh cuaca skala regional (tidak, kami tidak akan masuk ke dalam diskusi itu).

Saya ingin menunjukkan secara grafis bagaimana nilai-nilai yang diukur secara bersamaan berkorelasi lintas waktu dan ruang. Tujuan saya adalah untuk menunjukkan homogenitas regional (atau ketiadaan nilai) dari nilai yang sedang diselidiki.

Himpunan data

Untuk memulainya, saya mengambil sekelompok stasiun di wilayah Massachusetts dan Maine. Saya memilih situs berdasarkan garis lintang dan bujur dari file indeks yang tersedia di situs FTP NOAA.

masukkan deskripsi gambar di sini

Langsung Anda melihat satu masalah: ada banyak situs yang memiliki pengidentifikasi serupa atau sangat dekat. FWIW, saya mengidentifikasi mereka menggunakan kode USAF dan WBAN. Melihat lebih dalam ke metadata saya melihat bahwa mereka memiliki koordinat dan ketinggian yang berbeda, dan data berhenti di satu situs kemudian mulai dari yang lain. Jadi, karena saya tidak tahu yang lebih baik, saya harus memperlakukan mereka sebagai stasiun terpisah. Ini berarti data berisi pasangan stasiun yang sangat dekat satu sama lain.

Analisis awal

Saya mencoba mengelompokkan data berdasarkan bulan kalender dan kemudian menghitung regresi kuadrat terkecil biasa antara pasangan data yang berbeda. Saya kemudian memplot korelasi antara semua pasangan sebagai garis yang menghubungkan stasiun (di bawah). Warna garis menunjukkan nilai R2 dari kecocokan OLS. Gambar tersebut kemudian menunjukkan bagaimana 30+ titik data dari Januari, Februari, dll. Berkorelasi antara stasiun yang berbeda di bidang yang diminati.

korelasi antara data harian selama setiap bulan kalender

Saya telah menulis kode yang mendasarinya sehingga rata-rata harian hanya dihitung jika ada poin data setiap periode 6 jam, sehingga data harus dapat dibandingkan di seluruh situs.

Masalah

Sayangnya, ada terlalu banyak data yang masuk akal pada satu plot. Itu tidak bisa diperbaiki dengan mengurangi ukuran garis.

kmasukkan deskripsi gambar di sini

Jaringan tampaknya terlalu rumit, jadi saya pikir saya perlu mencari cara untuk mengurangi kompleksitas, atau menerapkan semacam kernel spasial.

Saya juga tidak yakin apa metrik yang paling tepat untuk menunjukkan korelasi, tetapi untuk audiens yang dimaksud (non-teknis), koefisien korelasi dari OLS mungkin hanya yang paling sederhana untuk dijelaskan. Saya mungkin perlu menyajikan beberapa informasi lain seperti kesalahan gradien atau standar juga.

Pertanyaan

Saya belajar cara saya ke bidang ini dan R pada saat yang sama, dan akan sangat menghargai saran tentang:

  1. Apa nama yang lebih formal untuk apa yang saya coba lakukan? Adakah beberapa istilah yang membantu yang memungkinkan saya menemukan lebih banyak lektur? Pencarian saya kosong untuk apa yang harus menjadi aplikasi umum.
  2. Apakah ada metode yang lebih tepat untuk menunjukkan korelasi antara beberapa set data yang dipisahkan dalam ruang?
  3. ... khususnya, metode yang mudah untuk menunjukkan hasil dari secara visual?
  4. Apakah ini diimplementasikan dalam R?
  5. Apakah ada dari pendekatan ini yang mengarah pada otomatisasi?
Andy Clifton
sumber
[Menggambarkan Korelasi Temporal secara Spasial di Lingkungan Visual Analytics, "Abish Malik et al.] [1] [1]: google.com/…
pat
2
yWy
Bagaimana jika Anda mencoba meningkatkan ambang ploting (0,5) dan menggunakan lebih dari 4 langkah warna? Atau untuk menggunakan garis yang lebih tipis-lebih tebal daripada warna.
nadya
nmemesan((n2)/2)
1
Saya menyadari dari sini bahwa saya memiliki banyak pekerjaan yang harus dilakukan pada pra-pemrosesan data sebelum saya memulai analisis yang telah saya uraikan di sini. Membaca jawaban dari @nadya saya pikir jelas saya perlu melihat beberapa jenis agregasi spasial, tetapi itu akan sulit karena salah untuk menggabungkan data darat dan laut. Maka saya perlu melihat strategi mengisi celah. Kemudian (dan hanya kemudian) saya bisa mulai melihat pekerjaan pemetaan / visualisasi.
Andy Clifton

Jawaban:

10

Saya pikir ada beberapa opsi untuk menampilkan tipe data ini:

Opsi pertama adalah melakukan "Analisis Fungsi Orthogonal Empiris" (EOF) (juga disebut sebagai "Analisis Komponen Utama" (PCA) dalam lingkaran non-iklim). Untuk kasus Anda, ini harus dilakukan pada matriks korelasi lokasi data Anda. Misalnya, matriks data Anda datakan menjadi lokasi spasial Anda di dimensi kolom, dan parameter yang diukur di baris; Jadi, matriks data Anda akan terdiri dari deret waktu untuk setiap lokasi. The prcomp()fungsi akan memungkinkan Anda untuk mendapatkan komponen utama, atau mode dominan korelasi, yang berkaitan dengan bidang ini:

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.

Opsi kedua adalah membuat peta yang menunjukkan korelasi relatif terhadap lokasi tertentu yang diminati:

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

Sunting: contoh tambahan

Meskipun contoh berikut tidak menggunakan data gappy, Anda bisa menerapkan analisis yang sama ke bidang data setelah interpolasi dengan DINEOF ( http://menugget.blogspot.de/2012/10/dineof-data-interpolating-empirical.html ) . Contoh di bawah ini menggunakan subset data tekanan permukaan laut anomali bulanan dari kumpulan data berikut ( 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)

Petakan mode EOF terkemuka

# 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

masukkan deskripsi gambar di sini

Buat peta korelasi

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

masukkan deskripsi gambar di sini

Marc di dalam kotak
sumber
Seberapa baik fungsi-fungsi ini menangani data yang hilang? Saya cukup sering memiliki celah dalam deret waktu.
Andy Clifton
2
Ada metode EOF yang dirancang untuk kasus khusus "data gappy" yang Anda gambarkan. Berikut ini tautan ke makalah yang mengulas metode ini: dx.doi.org/10.6084/m9.figshare.732650 . Anda akan melihat bahwa metode RSEOF dan DINEOF adalah yang paling akurat untuk memperoleh EOF dari set data gappy. Algoritma interpolasi DINEOF dapat ditemukan di sini: menugget.blogspot.de/2012/10/…
Marc di dalam kotak
1
Saya pikir ini adalah jawaban terbaik untuk apa pertanyaan yang mengerikan (di belakang).
Andy Clifton
3

Saya tidak melihat dengan jelas di balik garis tetapi bagi saya tampaknya ada terlalu banyak titik data.

Karena Anda ingin menunjukkan homogenitas regional dan bukan stasiun, saya sarankan Anda terlebih dahulu untuk mengelompokkannya secara spasial. Misalnya, overlay oleh "jala" dan menghitung nilai rata-rata yang diukur di setiap sel (setiap saat). Jika Anda menempatkan nilai rata-rata ini di pusat sel dengan cara ini Anda meraster data (atau Anda dapat menghitung juga berarti garis lintang dan bujur di setiap sel jika Anda tidak ingin overlay garis). Atau rata-rata di dalam unit administrasi, apa pun. Kemudian untuk "stasiun" rata-rata baru ini Anda dapat menghitung korelasi dan plot peta dengan jumlah garis yang lebih kecil.

masukkan deskripsi gambar di sini

Ini juga dapat menghapus garis korelasi tunggal acak tinggi yang melewati semua area.

nadya
sumber
x×xx
Ya, memproyeksikan koordinat adalah ide yang bagus. Semoga berhasil!
nadya