comparison pracma::gradient() and metR::Derivate()

[chrisdane]

Hi

How do I need to interpret the different 2d-gradient results from pracma::gradient() and metR::Derivate()?

# 2d example data from ?pracma::gradient
v <- seq(-2, 2, by=0.2)
X <- meshgrid(v, v)$X
Y <- meshgrid(v, v)$Y
Z <- X * exp(-X^2 - Y^2)
# pracma 2d-gradient
grXY_pracma <- pracma::gradient(Z, v, v)
# metR 2d-gradient
XY_vec <- expand.grid(X=v, Y=v)
data <- data.frame(X=XY_vec$X, Y=XY_vec$Y, Z=as.vector(Z))
grXY_metR <- metR::Derivate(formula=Z ~ X + Y, data=data, order=1, cyclical=c(F, F), sphere=F)
png("pracma_gradient_vs_metr_derivate.png", width=2000, height=1000, res=200)
par(mfrow=c(1, 2))
image(v, v, t(Z), main="pracma::gradient()")
contour(v, v, t(Z), col="black", add = TRUE)
grid(col="white")
quiver(X, Y, grXY_pracma$X, grXY_pracma$Y, scale = 0.2, col="blue")
image(v, v, t(Z), main="metR::Derivate()")
contour(v, v, t(Z), col="black", add = TRUE)
grid(col="white")
quiver(XY_vec$X, XY_vec$Y, grXY_metR$Z.dX, grXY_metR$Z.dY, scale = 0.2, col="blue")
dev.off()

Thanks a lot! Chris

[eliocamp]

Hi! I've never used pracma but the issue you're having here with metR is something to to with ordering. Because your X and Y coordinates are in separate variables from your values, it gets rotated somehow.

Here's and example of using metR using tidy data which works as expected:

v <- seq(-2, 2, by=0.2)

XY_vec <- data.table::as.data.table(expand.grid(X=v, Y=v))
XY_vec[, Z := X * exp(-X^2 - Y^2)]
XY_vec[, c("dx", "dy") := metR::Derivate(Z ~ X + Y)]

library(ggplot2)

ggplot(XY_vec, aes(X, Y)) +
    geom_contour(aes(z = Z, color = ..level..)) + 
    metR::geom_arrow(aes(dx = dx, dy = dy))

^{Created on 2020-04-11 by the reprex package (v0.3.0)}

[eliocamp]

Here's an example closer to yours to illustrate what's going on.

v <- seq(-2, 2, by=0.2)
X <- pracma::meshgrid(v, v)$X
Y <- pracma::meshgrid(v, v)$Y
Z <- X * exp(-X^2 - Y^2)
# pracma 2d-gradient
grXY_pracma <- pracma::gradient(Z, v, v)
# metR 2d-gradient
XY_vec <- expand.grid(X=v, Y=v)
data <- data.frame(X=XY_vec$X, Y=XY_vec$Y, Z=as.vector(Z))
grXY_metR <- metR::Derivate(formula=Z ~ X + Y, data=data, order=1, cyclical=c(F, F), sphere=F)

par(mfrow=c(1, 2))
image(v, v, t(Z), main="pracma::gradient()")
contour(v, v, t(Z), col="black", add = TRUE)
grid(col="white")

pracma::quiver(X, Y, grXY_pracma$X, grXY_pracma$Y, scale = 0.2, col="blue")
#> Warning in arrows(x, y, x + scale * u, y + scale * v, angle = 10, length =
#> length, : zero-length arrow is of indeterminate angle and so skipped

#> Warning in arrows(x, y, x + scale * u, y + scale * v, angle = 10, length =
#> length, : zero-length arrow is of indeterminate angle and so skipped

#> Warning in arrows(x, y, x + scale * u, y + scale * v, angle = 10, length =
#> length, : zero-length arrow is of indeterminate angle and so skipped

#> Warning in arrows(x, y, x + scale * u, y + scale * v, angle = 10, length =
#> length, : zero-length arrow is of indeterminate angle and so skipped
image(v, v, Z, main="metR::Derivate()")     # Not transposed
contour(v, v, Z, col="black", add = TRUE)   # Not transposed
grid(col="white")
pracma::quiver(XY_vec$X, XY_vec$Y, grXY_metR$Z.dX, grXY_metR$Z.dY, scale = 0.2, col="blue")

^{Created on 2020-04-11 by the reprex package (v0.3.0)}

[chrisdane]

Great, thank you very much for figuring this out!

[chrisdane]

One follow-up with real data

# real data
f <- system.file("extdata", "temperature.nc", package="metR")
nc <- ncdf4::nc_open(f)
lon <- as.vector(nc$dim$lon$vals)
lat <- as.vector(nc$dim$lat$vals)
tair <- ncdf4::ncvar_get(nc, "air", start=c(1, 1, 1, 1), count=c(length(lon), length(lat), 1, 1))
lat <- rev(lat)
tair <- tair[,length(lat):1]

# pracma gradient
dxy_tair <- pracma::gradient(t(tair), lon, lat) 
png("pracma_gradient_vs_metr_derivate_real_data.png", width=2000, height=1000, res=200)
par(mfrow=c(1, 2))
image(lon, lat, tair, xlim=c(350, 357.5), ylim=c(35, 45), main="pracma::gradient()")
maps::map("world2", add=T)
lonlat_mat <- pracma::meshgrid(lon, lat)
pracma::quiver(lonlat_mat$X, lonlat_mat$Y, dxy_tair$X, dxy_tair$Y, scale = 0.2, col="blue")

# metR::Derivate() needs vector input, not matrix
data <- expand.grid(lon=lon, lat=lat)
data$z <- as.vector(tair)
dxy_tair <- metR::Derivate(formula=z ~ lon + lat, data=data,
                           order=1, cyclical=c(T, F), sphere=T)
image(lon, lat, tair, xlim=c(350, 357.5), ylim=c(35, 45), main="metR::Derivate()")
maps::map("world2", add=T)
pracma::quiver(lonlat_mat$X, lonlat_mat$Y, dxy_tair$z.dlon, dxy_tair$z.dlat, scale = 0.2, col="blue")
dev.off()

With many warnings:

There were 50 or more warnings (use warnings() to see the first 50)
1: In arrows(x, y, x + scale * u, y + scale * v, angle = 10,  ... :
  zero-length arrow is of indeterminate angle and so skipped              

Some factor missing? Thanks a lot and stay healthy! Chris

[eliocamp]

Probably also something to do with the methods you're using to plot them. Maybe quiver needs matrixes as the u and v inputs?

Using data.frames and ggplot2, it all lines up perfectly.

# real data
f <- system.file("extdata", "temperature.nc", package="metR")
nc <- ncdf4::nc_open(f)
lon <- as.vector(nc$dim$lon$vals)
lat <- as.vector(nc$dim$lat$vals)
tair <- ncdf4::ncvar_get(nc, "air", start=c(1, 1, 1, 1), count=c(length(lon), length(lat), 1, 1))
lat <- rev(lat)
tair <- tair[,length(lat):1]

library(data.table)  # Just to subset more easily and the %between% operator
library(ggplot2)
data <- expand.grid(lon=lon, lat=lat)
data$z <- as.vector(tair)
data[c("z.dlon", "z.dlat")] <- metR::Derivate(formula=z ~ lon + lat, data=data,
                           order=1, cyclical=c(T, F), sphere=T)
data <- as.data.table(data) 

ggplot(data[lon %between% c(330, 360) & lat %between% c(35, 45)], aes(lon, lat)) +
  geom_raster(aes(fill = z)) +
  metR::geom_arrow(aes(dx = z.dlon, dy = z.dlat)) +
  coord_cartesian(xlim = c(349, 358), ylim= c(35, 45))

^{Created on 2020-04-11 by the reprex package (v0.3.0)}

A more "metR-like" approach would be

library(ggplot2)
library(data.table)

f <- system.file("extdata", "temperature.nc", package="metR")
data <- metR::ReadNetCDF(f, vars = c(z = "air"), subset = list(level = 1000))
data[, c("z.dlon", "z.dlat") := metR::Derivate(z ~ lon + lat, cyclical = c(TRUE, FALSE), sphere = TRUE), 
     by = level]

ggplot(data[lon %between% c(330, 360) & lat %between% c(35, 45)], aes(lon, lat)) +
  geom_raster(aes(fill = z)) +
  metR::geom_arrow(aes(dx = z.dlon, dy = z.dlat)) +
  coord_cartesian(xlim = c(349, 358), ylim= c(35, 45))

^{Created on 2020-04-11 by the reprex package (v0.3.0)}