How to calculate the matrix gradient for drawing a vector field in R?

Question

I looked around, but all I could find was to find the derivative of the matrix using diff(d) , where d is the matrix. This does not give me vectors, just a bunch of scalars. I'm not quite sure what to do with them.

I would like to find a way to calculate the gradient at several points across the surface, which is represented by a matrix. This gradient can be displayed as a vector field. Here's the question of creating vector fields in R, but I don't know how to calculate the gradient.

Edit: I will try to elaborate on what I'm looking for. Say I have a matrix like this:

  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 

It looks something like this when you close it up:

Now, what I want is just this: at certain intervals x and y, I would like to find the slope of the surface. So, for example, starting with x = 0, y = 0, I would like to find the slope as a vector, which I can use to build later. Then find the slope at x = 0, y = 1, etc. For all y values. Then find all y values ​​for x = 1, etc.

The goal is to have a bunch of vectors that can be built in a vector field like this .

Can this be done in R?

Answer 1

Here are some things to get you started.

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

You need to decide whether to fill in NA or 0 at the beginning or end, or copy the first or last value to diff(x) or ...

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

Gradients in the x direction (row):

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

In y direction (column):

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

In your example, something like this works:

 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

This uses the raster batch approach. Start with the same matrix as 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) 

Convert to raster (note that transposition and general driving is because matrices start from the upper left corner, and the coordinates work from the lower right):

 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) 

Now that Ben has hinted, you need to give him a coordinate system. At the moment, it just has the numeric coordinates of rows and columns from 1 to 10 and from 1 to 15. If it was a real-world map, then raster should know if it is long-long, or meters or feet, and X coordinates and coordinates Y are on the same scale. This is important even for data that is not displayed in the real world, as I suspect, your data.

A gradient does not make sense if your X and Y coordinates are not in the same units. If X is the resistance in ohms, and Y is the current in amperes, and Z is your measured potential in volts, then what is the bias? Well, it could be 2 V per ohm on the X axis and -3 V per amplifier in the Y direction. So what? You cannot say because you cannot combine ohms and amplifiers to get direction.

So, I proceed from the assumption that no matter what units X and Y are in your example, they are the same units (maybe they are ohms on resistor A and ohms on resistor B), and they go from 1 to 10 and 1 to 15.

Now I think that there is a projection code that simply says: “These are the x and y coordinates without a real geographical value,” but I cannot remember what it is or find. So I just lie and use any old coordinate system that I know is a normal Cartesian grid. In this case, the National GB Network. If you tried to build this bitmap on a map, it will be a tiny square off the coast of southwest England, because where the grid takes place, and your data is 10 m by 15 m in this system:

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

Let me talk to make sure we haven't messed up yet:

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

Note that the X and Y coordinates are in the same direction as your sample, so I know that everything is still fine.

Now I can just do a levelplot from rasterVis , but I need to do a little scaling. This is because the gradient on a real map is calculated from heights and distances that have the same units (e.g. meters or feet), but your data is just numbers. Consequently, the gradients are actually quite small in a natural integer coordinate system. So:

 vectorplot(r, scaleSlope=.1) 

gives you:

Note that the slope is usually downward, because the way your X and Y axes are in your example (and therefore in my raster). Also note that the cells are square because we keep the aspect ratio of the data (because we treat the X and Y coordinates as equal in measure). Ben's answer shows the general LR flow, which means that his X and Y coordinates are not in the usual order.

In addition, the gradient search algorithm in vectorplot performs some degree of smoothing, so a small gap in the upper right corner does not look as extreme as in the Ben differential algorithm:

but you need to decide if you really want to build a smooth gradient or final differences ...