0

I have two vectors representing x and y-coordinates in a scatter plot, and a thrid variable (z) for each (x,y)-coordinate representing the variable from which to draw contour lines. Example data are given as follows:

df<-data.frame(x=runif(n=30,min=-6,max=6),
               y=runif(n=30,min=-6,max=10),
               z=seq(1,100,length.out=30))

I use the R-package akima to generate the z-matrix for the contour plot

library(akima)
M1 <- interp(x=df$x,y=df$y,z=df$z)
contour(x=M1$x,y=M1$y,z=M1$z)

I now want to draw arrows perpendicular to the contourlines, preferably using something like the function "quiver" in the R-package pracma, with the origin of an arrow at every (x,y)-coordinate and with the arrow pointing in the direction of the gradient of the contourlines. Is there a way to do this?

My best idea so far is to somehow extract (x,y)-gradients of the contourlines and use these as velocities in the quiver function.

Grateful for any assistance.

user2554330
  • 37,248
  • 4
  • 43
  • 90
user09034
  • 15
  • 6
  • 2
    Could you please clarify: do you want the arrows to be attached to the contour lines, or at every `(M1$x, M1$y)` pair? – user2554330 Oct 23 '18 at 13:11
  • Good point. I want the arrows attached at every (M1$x,M1$y) pair – user09034 Oct 23 '18 at 13:21
  • 1
    You can use the raster package "terrain" function with the "aspect" option to get the angle from a point on an elevation raster. e.g. `plot(terrain(raster(M1),"aspect"))` - then sample the raster at the points, draw arrows with sines and cosines... – Spacedman Oct 23 '18 at 13:53
  • @Spacedman Thank you for the suggestion, I will try it out – user09034 Oct 23 '18 at 13:58

1 Answers1

2

The pracma package has a gradient function that can do this for you using the original M1$z values. For example, using your code to get M1 after set.seed(123):

contour(x=M1$x,y=M1$y,z=M1$z, asp = 1) # asp = 1 needed so things look perpendicular

library(pracma)
g <- gradient(M1$z, M1$x, M1$y)

x <- outer(M1$x, M1$y, function(x, y) x)
y <- outer(M1$x, M1$y, function(x, y) y)

quiver(x, y, g$Y, g$X, scale = 0.02, col = "blue")

Note that the gradient labels in the quiver plot have been swapped. Maybe I set up the x and y values transposed from the way the package expects. Here's what you get:

enter image description here

user2554330
  • 37,248
  • 4
  • 43
  • 90
  • Follow-up question: Do you know of a way to attach the arrows to the original (x,y)-coordinates? (df$x,df$y) – user09034 Oct 25 '18 at 05:45
  • 1
    You might be able to use `pracma::interp2` to interpolate gradient values to those points, but in a quick try just now, I found it fails on any points on the convex hull of your dataset. With 30 points, that's half of the points. – user2554330 Oct 25 '18 at 09:54
  • OK will try to use interp2, grateful for your efforts – user09034 Oct 26 '18 at 10:04