Come si calcola il gradiente di una matrice per disegnare un campo vettoriale in R?

Question

Mi sono guardato intorno, ma tutto quello che ho potuto trovare era come trovare la derivata di una matrice usando diff(d), dove d è una matrice. Questo non mi dà i vettori, solo un mucchio di scalari. Non sono davvero sicuro di cosa fare con quelli.

Mi piacerebbe trovare un modo per calcolare il gradiente in pochi punti su tutta la superficie rappresentata da una matrice. Questo gradiente può essere visualizzato come un campo vettoriale. C'è un question here sulla creazione di campi vettoriali in R, ma non so come calcolare il gradiente.

Modifica: cercherò di elaborare quello che sto cercando. Diciamo che ho una matrice simile a questo:

 X0 X1.5 X3.1 X4.3 X5.9 X7.3 X8.6 X9.8 X11 X12.3 X13.6 X14.9 X16.4 X17.9 X20 
[1,] 0 1.4 3.0 4.5 6.0 7.3 8.6 9.7 10.9 12.2 13.4 14.9 16.4 18.1 20 
[2,] 0 1.6 3.2 4.9 6.4 7.6 8.7 9.6 10.6 11.8 13.2 14.7 16.4 18.1 20 
[3,] 0 1.7 3.5 5.2 7.0 8.3 9.0 9.4 9.9 11.1 12.7 14.6 16.3 18.2 20 
[4,] 0 1.8 3.7 5.8 8.0 9.3 9.3 9.3 9.4 10.2 12.1 14.1 16.2 18.0 20 
[5,] 0 1.7 3.9 6.0 8.8 9.3 9.3 9.4 9.6 9.9 11.8 14.0 16.2 18.1 20 
[6,] 0 1.8 3.8 5.7 8.1 9.3 9.3 9.4 9.6 10.1 12.3 14.4 16.3 18.0 20 
[7,] 0 1.6 3.5 5.2 7.0 8.4 9.1 9.5 10.1 11.3 13.0 14.6 16.4 18.2 20 
[8,] 0 1.5 3.2 4.9 6.4 7.7 8.7 9.7 10.7 11.9 13.3 14.9 16.5 18.3 20 
[9,] 0 1.5 3.1 4.6 6.0 7.4 8.6 9.7 10.9 12.1 13.5 15.1 16.6 18.3 20 
[10,] 0 1.5 3.0 4.6 6.0 7.3 8.5 9.7 10.9 12.4 13.6 13.1 16.6 18.2 20 

Sembra qualcosa di simile quando si traccia è:

Ora, quello che voglio è semplicemente questo: a determinati intervalli di xe y, mi piacerebbe essere in grado di trovare la pendenza della superficie. Quindi, per esempio, iniziando con x = 0, y = 0, vorrei trovare la pendenza sotto forma di un vettore che posso usare per tracciare in seguito. Quindi, trova la pendenza in x = 0, y = 1 e così via per tutti i valori di y. Quindi trova tutti i valori di y per x = 1 e così via.

Lo scopo di questo è di avere un gruppo di vettori che possono essere tracciati in un campo vettoriale like this.

Questo può essere fatto in R?

Answer 1

Ecco alcune cose per iniziare.

m <- matrix(1:9,nrow=3) 

Dovete decidere se compilare NA o 0 all'inizio o alla fine, o replicare il primo o l'ultimo valore diff(x), o ...

bdiff <- function(x) c(NA,diff(x)) 

gradienti nel x (riga) direzione:

t(apply(m,1,bdiff)) 
##  [,1] [,2] [,3] 
## [1,] NA 3 3 
## [2,] NA 3 3 
## [3,] NA 3 3 

Nella direzione y (colonna):

apply(m,2,bdiff) 
##  [,1] [,2] [,3] 
## [1,] NA NA NA 
## [2,] 1 1 1 
## [3,] 1 1 1 

Per esempio, qualcosa circa come questo funziona:

m2 <- matrix(c(
0,1.4,3.0,4.5,6.0,7.3,8.6,9.7,10.9,12.2,13.4,14.9,16.4,18.1,20, 
0,1.6,3.2,4.9,6.4,7.6,8.7,9.6,10.6,11.8,13.2,14.7,16.4,18.1,20, 
0,1.7,3.5,5.2,7.0,8.3,9.0,9.4,9.9,11.1,12.7,14.6,16.3,18.2,20, 
0,1.8,3.7,5.8,8.0,9.3,9.3,9.3,9.4,10.2,12.1,14.1,16.2,18.0,20, 
0,1.7,3.9,6.0,8.8,9.3,9.3,9.4,9.6,9.9,11.8,14.0,16.2,18.1,20, 
0,1.8,3.8,5.7,8.1,9.3,9.3,9.4,9.6,10.1,12.3,14.4,16.3,18.0,20, 
0,1.6,3.5,5.2,7.0,8.4,9.1,9.5,10.1,11.3,13.0,14.6,16.4,18.2,20, 
0,1.5,3.2,4.9,6.4,7.7,8.7,9.7,10.7,11.9,13.3,14.9,16.5,18.3,20, 
0,1.5,3.1,4.6,6.0,7.4,8.6,9.7,10.9,12.1,13.5,15.1,16.6,18.3,20, 
0,1.5,3.0,4.6,6.0,7.3,8.5,9.7,10.9,12.4,13.6,13.1,16.6,18.2,20), 
byrow=TRUE,nrow=10) 

rr <- row(m2) 
cc <- col(m2) 
dx <- t(apply(m2,1,bdiff)) 
dy <- apply(m2,2,bdiff) 
sc <- 0.25 
off <- -0.5 ## I *think* this is right since we NA'd row=col=1 
plot(rr,cc,col="gray",pch=16) 
arrows(rr+off,cc+off,rr+off+sc*dx,cc+off+sc*dy,length=0.05) 

Answer 2

Ecco l'approccio raster pacchetto. Inizia con la stessa matrice di risposta di Ben:

m2 <- matrix(c(
0,1.4,3.0,4.5,6.0,7.3,8.6,9.7,10.9,12.2,13.4,14.9,16.4,18.1,20, 
0,1.6,3.2,4.9,6.4,7.6,8.7,9.6,10.6,11.8,13.2,14.7,16.4,18.1,20, 
0,1.7,3.5,5.2,7.0,8.3,9.0,9.4,9.9,11.1,12.7,14.6,16.3,18.2,20, 
0,1.8,3.7,5.8,8.0,9.3,9.3,9.3,9.4,10.2,12.1,14.1,16.2,18.0,20, 
0,1.7,3.9,6.0,8.8,9.3,9.3,9.4,9.6,9.9,11.8,14.0,16.2,18.1,20, 
0,1.8,3.8,5.7,8.1,9.3,9.3,9.4,9.6,10.1,12.3,14.4,16.3,18.0,20, 
0,1.6,3.5,5.2,7.0,8.4,9.1,9.5,10.1,11.3,13.0,14.6,16.4,18.2,20, 
0,1.5,3.2,4.9,6.4,7.7,8.7,9.7,10.7,11.9,13.3,14.9,16.5,18.3,20, 
0,1.5,3.1,4.6,6.0,7.4,8.6,9.7,10.9,12.1,13.5,15.1,16.6,18.3,20, 
0,1.5,3.0,4.6,6.0,7.3,8.5,9.7,10.9,12.4,13.6,13.1,16.6,18.2,20), 
byrow=TRUE,nrow=10) 

convertire in raster (nota la trasposizione e giocherellare generale è perché matrici Inizia dalla parte superiore sinistra, ma le coordinate di lavoro dal basso a destra):

require(raster) 
require(rasterVis) 
r=raster(t(m2[,ncol(m2):1]), xmn=0.5,xmx=nrow(m2)+.5, ymn=0.5,ymx=ncol(m2)+0.5) 

Ora come Ben ha suggerito, è do necessario dargli un sistema di coordinate. Al momento ha solo coordinate numeriche di righe e colonne da 1 a 10 e da 1 a 15. Se questa era una mappa del mondo reale, allora raster ha bisogno di sapere se questo è lat-long, o metri, o piedi, e se la X le coordinate e le coordinate Y sono sulla stessa scala.Questo è importante, anche per i dati che non è mappato al mondo reale come sospetto che i tuoi dati siano.

Il gradiente non ha significato se le coordinate X e Y non si trovano nelle stesse unità. Se X è una resistenza in ohm e Y è corrente in ampere e Z è il tuo potenziale misurato in Volt allora qual è la pendenza? Bene, potrebbe essere 2 V per ohm sull'asse X e -3 V per amplificatore nella direzione Y. Quindi del tutto? Non si può dire, perché non è possibile combinare ohm e amplificatori per ottenere una direzione.

Quindi assumerò che tutte le unità di X e Y siano nel tuo esempio, sono le stesse unità (forse sono ohm sul resistore A e ohhm sul resistore B) e vanno da 1 a 10 e da 1 a 15.

Ora penso che ci sia un codice di proiezione che dice semplicemente "Queste sono le coordinate xey senza un reale significato geografico" ma non riesco a ricordare cosa sia o trovarlo. Quindi mi limiterò a mentire e ad usare qualsiasi vecchio sistema di riferimento di coordinate che conosco sia una normale griglia cartesiana. In questo caso, la griglia nazionale GB. Se si è tentato di tracciare il raster su una mappa sarebbe una piccola piazza al largo della costa sud-ovest dell'Inghilterra, perché è lì che l'origine della griglia è, ei dati sono 10m da 15m in questo sistema:

projection(r)=CRS("+init=epsg:27700") 

Diciamo tracciarla per assicurarsi che non abbiamo ancora incasinato:

persp(r,theta=-50,phi=20, shade=0.23,col="red") 

nota che le coordinate X e Y siano rivolte nella stessa direzione come trama del campione, quindi so io' ce l'ho fatta fin qui.

Ora posso solo fare levelplot da rasterVis, ma devo fare un leggero ridimensionamento. Questo perché il gradiente su una mappa reale viene calcolato da altezze e distanze che hanno le stesse unità (forse metri o piedi) ma i tuoi dati sono solo numeri. Quindi i gradienti sono in realtà piuttosto piccoli nel sistema di coordinate intero naturale. Quindi:

vectorplot(r, scaleSlope=.1) 

ti dà:

Nota la pendenza è generalmente verso il basso, perché questo è il modo in cui l'assi X e Y sono nel tuo esempio plot (e quindi nel mio raster). Nota anche che le celle sono quadrate perché stiamo preservando le proporzioni dei dati (perché stiamo trattando le coordinate X e Y come uguali nella misura). La risposta di Ben mostra un flusso L-R generale, il che significa che le sue coordinate X e Y non sono nell'ordine convenzionale.

Inoltre, l'algoritmo gradiente trovando in vectorplot fa un certo grado di lisciatura, in modo che il piccolo discontinuità in alto a destra non sembra così estreme come in algoritmo di differenziazione di Ben:

, ma si deve decidi se vuoi veramente tracciare il gradiente levigato o le differenze finite ...