What is causing the error "Z matrix not specified" when making a graph of vectors

Question

I am making a graph for some vectors, and I have the following code (r script):

calculating the derivatives

#calculating the derivatives
#x direction`#x direction
gradientx2d_dd=function(x,dx){}
dd=uwind
xdim=dd[1]
ydim=dd[2]
gradx=0.0*uwind
for (i in 1:ydim){
gradx[2:xdim,i]=diff(uwind[1:uwind,i])/dx}
return(gradx)

#y direction
gradienty2d_dd=function(x,dy){}
dd=uwind
xdim=dd[1]
ydim=dd[2]
grady=0.0*uwind
for (i in 1:xdim){
grady[i,2:ydim]=diff(uwind[i, 1:ydim])/dy}
return(grady)

Divergence=gradientx2d_dd(uwind,dx)+gradienty2d_dd(uwind,dy)

#plot divergence
filled.contour(lon,lat,Divergence,zlim=c=
(-10,10),nlevels=124,colour=rainbow, main="Mean wind Divergence",
plot.axes={axis(1);axis(2);map(add=T)})

I know the error is within the y derivative section, which I have taken from a previous code where it worked. I have double, triple checked it, however it is saying that it cannot plot the graph, as the "no proper 'z' matrix specified". However, I think I have specified it as divergence. I cannot figure out where I went wrong here.

Answer 1

there are a couple of things that might cause problems in your code.

1. The body of your function should be placed between the curly braces
2. The functions take an argument x, but within the function, the data is hard coded as uwind (dd = uwind instead of dd = x)
3. The dimensions of your matrix should be specified using xdim = nrow(dd) and ydim = ncol(dd) (not xdim = dd[1])
4. In the calculation of the gradx, the index should be uwind[1:xdim,i] (or better yet uwind[,i]), not uwind[1:uwind,i]

Your function for the derivative in x-direction should therefore look something like this:

gradientx2d_dd = function(x, dx){ 
    dd = x
    xdim = nrow(dd) 
    ydim = ncol(dd)
    gradx = 0.0*uwind 
    for (i in 1:ydim){ 
        gradx[2:xdim, i] = diff(uwind[1:xdim, i])/dx} 
    return(gradx) 
}

The same function can, however, be implemented in a simpler way using apply:

gradientx2d_dd <- function(x, dx){ 
    gradx <- apply(x, 2, function(xx) diff(xx)/dx)
} 

The same goes for the derivative in y-direction:

gradientxyd_dd = function(x, dy){ 
    dd = x
    xdim = nrow(dd) 
    ydim = ncol(dd)
    grady = 0.0*uwind 
    for (i in 1:xdim){ 
        grady[i, 2:ydim]=diff(uwind[i, 1:ydim])/dy}
    return(grady) 
}

Or simpler:

gradienty2d_dd <-  function(z, dy){ 
    grady <- t(apply(z, 1, function(xx) diff(xx)/dy))
}

Using these functions you can compute the derivatives but you will find that the dimensions of the return values of the functions do not match. The gradx has one fewer row than your uwind and the grady has one fewer column, as your code evaluates the derivatives in the midpoints of the x-grid (for gradx) or the midpoints of the y-grid (for grady). Therefore the line

Divergence = gradientx2d_dd(uwind,dx) + gradienty2d_dd(uwind,dy)

will likely throw an error.

In order to fix this you could do the following (my maths is a bit rusty, so this might not be 100 % correct):

1. assume some boundary conditions for your function uwind (e.g. value = 0 on the boundary [Dirichlet] or derivative is 0 on the boundary [Neumann])
2. move the whole evaluation of the gradients to the mid points of the grid using e.g. linear interpolation (then you will also have to adapt the grid for the plot)

An example for approach (2) could look like this:

Sample data:

## Set up your grid:
dx <- 1/100
dy <- 1/100
x <-  seq(0, 1, by = dx)
y <-  seq(0, 0.5, by = dy)
## Evaluate some differentiable function:
uwind <- outer(x,y,function(x, y) 6*x*y^2)

Code:

## Compute the gradients in x- and y-direction (at the mid-points of the grid)
gradientx2d_dd <- function(x, dx){ 
    gradx <- apply(x, 2, function(xx) diff(xx)/dx)
    gradx <- (gradx[, -1] + gradx[, -ncol(gradx)])/2
} 
gradienty2d_dd <-  function(x, dy){ 
    grady <- t(apply(x, 1, function(xx) diff(xx)/dy))
    grady <- (grady[-1, ] + grady[-nrow(grady), ])/2
}

## Compute divergence:
Divergence <- gradientx2d_dd(uwind, dx) + gradienty2d_dd(uwind, dy) 

## Plot divergence 
xMidpoints <- (x[-1] + x[-length(x)])/2
yMidpoints <- (y[-1] + y[-length(y)])/2
filled.contour(xMidpoints, yMidpoints, Divergence, 
        nlevels = 124, main = "Mean wind Divergence")