如何計算在R中繪製矢量場的矩陣的梯度？

Question

我環顧四周，但我能找到的是如何使用diff(d)找到矩陣的導數，其中d是矩陣。這並不能爲我提供矢量，只是一堆標量。我不確定該怎麼做。

我想找到一種方法來計算矩陣表示的整個表面上的幾個點的梯度。該漸變可以顯示爲矢量場。有關於在R中創建矢量場的question here，但我不知道如何計算梯度。

編輯：我會盡力詳細說明我在找什麼。比方說，我有這樣一個矩陣：

現在，我想很簡單：：在x的一定的時間間隔和

 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 

當您繪製它，它看起來是這樣的Ÿ，我希望能夠找到表面的斜度。因此，例如，從x = 0開始，y = 0，我希望能夠以稍後可以使用的矢量形式查找斜率。然後，找出x = 0，y = 1處的斜率，以此類推y的所有值。然後找出x = 1的所有y值，依此類推。

這樣做的目的是爲了在矢量字段like this中繪製一組向量。

這可以在R中完成嗎？

Answer 1

這是一些開始的東西。

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

你必須決定是否在NA或0填充在開始或結束，或在diff(x)複製第一個或最後一個值，或...

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

梯度X （行）方向：

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

在Y（列）方向：

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

對於你的榜樣，東西大約是這樣工作的：

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

這裏的raster一攬子方針。開始用相同的基質本的答案：

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) 

轉換爲柵格（注意轉置與一般小提琴演奏是因爲矩陣從左上角開始，但是從底部協調工作的權利）：

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) 

現在按照Ben的說法，你需要給它一個座標系做。目前，它只有1到10和1到15的數字行和列座標。如果這是一個真實世界地圖，那麼raster需要知道這是經緯度，米還是英尺，以及X座標和Y座標是相同的比例。這一點很重要，即使對於我懷疑數據是真實世界的數據而言，即使對於不是也是如此。

如果您的X和Y座標不在相同單位中，則漸變無意義。如果X是以歐姆爲單位的電阻，Y是以安培爲單位的電流，而Z是以伏特爲單位的測量電位，那麼斜率是多少？那麼，它可能是X軸每歐姆2V，Y方向每安培-3V。所以呢？你不能說，因爲你不能結合歐姆和安培來獲得方向。

所以我會假設無論X和Y的單位在你的例子中，它們是相同的單位（也許它們是電阻器A上的歐姆和電阻器B上的歐姆）並且它們從1到10和1到15.

現在我認爲有一個投影代碼只是說「這些是x和y座標，沒有真正的地理意義」，但我不記得它是什麼或找到它。所以我只會說謊，並使用任何我知道的舊座標參考系統是一個普通的笛卡爾網格。在這種情況下，GB國家電網。如果你想在地圖上繪製該光柵這將是一個小方休英格蘭西南部的海岸，因爲那裏的座標原點是，你的數據在本系統10米由15米：

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

我們繪製它，以確保我們沒有亂了嗎：

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

注意，X和Y座標在同一方向的樣地指向，所以我知道我到目前爲止，這一切都很順利。

現在我可以從rasterVis做levelplot，但我必須做一個輕微的縮放。這是因爲真實地圖上的梯度是根據具有相同單位（米或英尺）的高度和距離來計算的，但您的數據只是數字。因此，在自然整數座標系中，梯度實際上非常小。所以：

vectorplot(r, scaleSlope=.1) 

爲您提供：

注意坡度一般是向下的，因爲這是你的X軸和Y軸是在你的榜樣情節的方式（因此在我的光柵）。還要注意，單元格是正方形的，因爲我們保留了數據的高寬比（因爲我們將X和Y座標視爲相等）。 Ben的答案顯示了一個一般的L-R流程，這意味着他的X和Y座標不是按照常規順序。

此外，在vectorplot梯度搜索算法做一些平滑化程度，所以很少間斷右上方看起來並不極端，如本的區別算法：

但你必須決定你是否真的想繪製平滑的梯度或有限差異...