eliocamp / metR

Tools for Easier Analysis of Meteorological Fields
https://eliocamp.github.io/metR/
139 stars 22 forks source link

Is Laplacian function working as expected? #170

Closed pascaloettli closed 1 year ago

pascaloettli commented 1 year ago

There is an error in the Derivate/Laplacian/Divergence/Vorticity help page. As far as I understand, we should get the same figures at the end of the example, which is not the case.

metR version 0.13.0

library(metR)

theta <- seq(0, 360, length.out = 20)*pi/180
theta <- theta[-1]
x <- cos(theta)
dx_analytical <- -sin(theta)
dx_finitediff <- Derivate(x ~ theta, cyclical = TRUE)[[1]]

# Curvature (Laplacian)
# Note the different boundary conditions for each dimension
variable <- expand.grid(lon = seq(0, 360, by = 3)[-1],
                        lat = seq(-90, 90, by = 3))
variable$z <- with(variable, cos(lat*pi/180*3) + sin(lon*pi/180*2))
variable <- cbind(
  variable,
  as.data.frame(Derivate(z ~ lon + lat, data = variable,
                         cyclical = c(TRUE, FALSE), order = 2)))
library(ggplot2)
ggplot(variable, aes(lon, lat)) +
  geom_contour(aes(z = z)) +
  geom_contour(aes(z = z.ddlon + z.ddlat), color = "red")
#> Warning: Removed 480 rows containing non-finite values (`stat_contour()`).


# The same as
ggplot(variable, aes(lon, lat)) +
  geom_contour(aes(z = z)) +
  geom_contour(aes(z = Laplacian(z ~ lon + lat, cyclical = c(TRUE, FALSE))),
               color = "red")

Created on 2023-01-12 with reprex v2.0.2

eliocamp commented 1 year ago

Well, it seems that it is not working as expected xD.

There might be an issue with the underlying derivate code. Right now I'm a bit busy with real life stuff, so I cannot promise a quick fix (:disappointed:) but I wanted to acknowledge that I've seen your issue.

eliocamp commented 1 year ago

This is now fixed in the development version. Thanks!

library(metR)
# Curvature (Laplacian)
# Note the different boundary conditions for each dimension
variable <- expand.grid(lon = seq(0, 360, by = 3)[-1],
                        lat = seq(-90, 90, by = 3))
variable$z <- with(variable, cos(lat*pi/180*3) + sin(lon*pi/180*2))
variable <- cbind(
    variable,
    as.data.frame(Derivate(z ~ lon + lat, data = variable,
                           cyclical = c(TRUE, FALSE), order = 2)))
variable$lap <- Laplacian(z ~ lon + lat, cyclical = c(TRUE, FALSE), data = variable)

library(ggplot2)
ggplot(variable, aes(lon, lat)) +
    geom_contour(aes(z = z)) +
    geom_contour(aes(z = z.ddlon + z.ddlat), color = "red") + 
    geom_contour(aes(z = lap), color = "magenta") 
#> Warning: Removed 480 rows containing non-finite values (`stat_contour()`).
#> Removed 480 rows containing non-finite values (`stat_contour()`).


with(variable, plot(lap, z.ddlon + z.ddlat))

Created on 2023-03-25 with reprex v2.0.2