r-lib / isoband

isoband: An R package to generate contour lines and polygons.
https://isoband.r-lib.org
Other
130 stars 14 forks source link

Inconsistencies with `isolines` results #18

Closed gueyenono closed 3 years ago

gueyenono commented 3 years ago

I create a grid of values with a simple linear function and I attempt to get the coordinates of the contour at a specific z level (i.e. z = 150). I do so with 3 different packages and the results obtained by the isoband::isolines() function very unexpected. Note that I installed the development version in the master branch.

Of the 3 methods in the code below, the results obtained by contoureR::getContourLines() are the best in the sense that they are the most accurate.

# Simple linear function
myfun <- function(x, y){
  5*x + 15*y
}

# Compute grid values
x <- y <- 0:100
z1 <- outer(X = x, Y = y,  myfun) # grid format for isoband::isolines() and grDevices::contourLines()
df <- expand.grid(x = x, y = y)
z2 <- with(df, 5*x + 15*y) # grid format for contoureR::getContourLines()

# Method 1: grDevices::contourLines()
q1 <- contourLines(x = x, y = y, z = z1, levels = 150)
myfun(q1[[1]]$x, q1[[1]]$y)

 [1] 148.2353 150.0000 150.0000 148.5714 148.2353 150.0000 150.0000 148.5714 148.2353 150.0000 150.0000 148.5714 148.2353 150.0000
[15] 150.0000 148.5714 148.2353 150.0000 150.0000 148.5714 148.2353 150.0000 150.0000 148.5714 148.2353 150.0000 150.0000 148.5714
[29] 148.2353 150.0000 150.0000 148.5714 148.2353 150.0000 150.0000 148.5714 148.2353 150.0000 150.0000 148.5714

# Method 2: contoureR::getContourLines()
q2 <- contoureR::getContourLines(df$x, df$y, z2, levels = 150)
myfun(q2$x, q2$y)

[1] 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150
[34] 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150
[67] 150 150 150 150 150 150 150

# Method 3: isoband::isolines()
q3 <- isoband::isolines(x = x, y = y, z = z1, levels = 150)
myfun(q3[[1]]$x, q3[[1]]$y)

[1]  50.00000  63.33333  76.66667  90.00000  90.00000 103.33333 116.66667 130.00000 130.00000 143.33333 156.66667 170.00000 170.00000
[14] 183.33333 196.66667 210.00000 210.00000 223.33333 236.66667 250.00000 250.00000 263.33333 276.66667 290.00000 290.00000 303.33333
[27] 316.66667 330.00000 330.00000 343.33333 356.66667 370.00000 370.00000 383.33333 396.66667 410.00000 410.00000 423.33333 436.66667
[40] 450.00000
clauswilke commented 3 years ago

Without looking in detail, I'm quite certain this is caused by confusion over the order in which x and/or y are given (probably y). I think isoband assumes the coordinate origin is in the bottom left, since that's what all grid graphics functions assume. If you plot the isobands you get you'll certainly see that they're just flipped relative to your input.