如何计算矩阵的梯度以在R中绘制矢量场？

Question

我环顾四周,但我能找到的是如何找到矩阵的导数diff(d),其中d是矩阵.这不会给我矢量,只是一堆标量.我不太确定如何处理这些问题.

我想找到一种方法来计算整个表面中由矩阵表示的几个点处的梯度.该梯度可以显示为矢量场.这里有一个关于在R中制作矢量场的问题,但我不知道如何计算梯度.

编辑:我将尝试详细说明我正在寻找的东西.假设我有一个这样的矩阵:

     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和y的某些间隔,我希望能够找到表面的斜率.因此,例如,从x = 0,y = 0开始,我想找到我可以用于稍后绘制的矢量形式的斜率.然后,找到x = 0,y = 1处的斜率,依此类推y的所有值.然后找到y的所有值,x = 1,依此类推.

这样做的目的是为了有一堆可在矢量场被绘制矢量的这样.

这可以在R中完成吗？

Answer 1

这是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)

现在,作为本暗示,你就需要给它一个坐标系.目前它只有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坐标指向与样本图相同的方向,所以我知道到目前为止我已经完成了这一切.

现在,我可以做levelplot的rasterVis,但我必须做轻微的脱屑.这是因为真实地图上的渐变是根据具有相同单位(可能是米或英尺)的高度和距离计算的,但您的数据只是数字.因此,在自然整数坐标系中,梯度实际上非常小.所以:

vectorplot(r, scaleSlope=.1)

给你:

请注意,斜率通常是向下的,因为这是X轴和Y轴在示例图中的方式(因此在我的栅格中).另请注意,单元格是正方形的,因为我们保留了数据的纵横比(因为我们将X和Y坐标视为度量相等).Ben的答案显示了一般的LR流量,这意味着他的X和Y坐标不是传统的顺序.

此外,梯度发现算法在vectorplot某种程度上平滑,因此右上角的小不连续性看起来不像Ben的差分算法那样极端:

但你必须决定是否真的要绘制平滑的渐变或有限的差异......

Answer 2

这里有一些东西可以开始。

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)