Comments (4)
rogiersbart
commented on July 22, 2024
Well, a reprex would be a great start I guess. That would require the replacement to be made in a separate feature branch I think. It is probably best we replace it anyway, and then just try to figure out what is going wrong in those specific cases I suppose?
When I look at the function documentation, it seems the jittering is part of akima, while it is not of interp. I would however not immediately expect overlapping points in the point datasets from (R)MODFLOW ... But there might be cases where this does happen. For instance when you would plot HUF head contours for a vertical cross-section with multiple HUF layers that can be zero thickness? Seems not very likely that's what you were doing, but maybe it was? :) Or inactive cells with overlapping geometries? Not sure if MODFLOW complains about that ... In any case, if the missing jitter would be the problem, we could just apply it before using the interp function.
from rmodflow.
cneyens
commented on July 22, 2024
A reprex for a simple case with an irregular input grid:
xy <- expand.grid(c(50, 150, 350 , 450, 550, 650, 750, 850, 950), seq(950, 50, -100))
names(xy) <- c('x','y')
xy$z <- seq_len(nrow(xy))
akima::interp(xy$x, xy$y, xy$z, xo = seq(50, 950, 100), yo = seq(50, 950, 100))
#> $x
#> [1] 50 150 250 350 450 550 650 750 850 950
#>
#> $y
#> [1] 50 150 250 350 450 550 650 750 850 950
#>
#> $z
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 82.0 73.0 64.0 55.0 46.0 37.0 28.0 19.0 10.0 1.0
#> [2,] 83.0 74.0 65.0 56.0 47.0 38.0 29.0 20.0 11.0 2.0
#> [3,] 83.5 74.5 65.5 56.5 47.5 38.5 29.5 20.5 11.5 2.5
#> [4,] 84.0 75.0 66.0 57.0 48.0 39.0 30.0 21.0 12.0 3.0
#> [5,] 85.0 76.0 67.0 58.0 49.0 40.0 31.0 22.0 13.0 4.0
#> [6,] 86.0 77.0 68.0 59.0 50.0 41.0 32.0 23.0 14.0 5.0
#> [7,] 87.0 78.0 69.0 60.0 51.0 42.0 33.0 24.0 15.0 6.0
#> [8,] 88.0 79.0 70.0 61.0 52.0 43.0 34.0 25.0 16.0 7.0
#> [9,] 89.0 80.0 71.0 62.0 53.0 44.0 35.0 26.0 17.0 8.0
#> [10,] 90.0 81.0 72.0 63.0 54.0 45.0 36.0 27.0 18.0 9.0
interp::interp(xy$x, xy$y, xy$z, xo = seq(50, 950, 100), yo = seq(50, 950, 100))
#> $x
#> [1] 50 150 250 350 450 550 650 750 850 950
#>
#> $y
#> [1] 50 150 250 350 450 550 650 750 850 950
#>
#> $z
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 82.0 73.0 64.0 55.0 46.0 37.0 28.0 19.0 10.0 NA
#> [2,] 83.0 74.0 65.0 56.0 47.0 38.0 29.0 20.0 11.0 NA
#> [3,] 83.5 74.5 65.5 56.5 47.5 38.5 29.5 20.5 11.5 NA
#> [4,] 84.0 75.0 66.0 57.0 48.0 39.0 30.0 21.0 12.0 NA
#> [5,] 85.0 76.0 67.0 58.0 49.0 40.0 31.0 22.0 13.0 NA
#> [6,] 86.0 77.0 68.0 59.0 50.0 41.0 32.0 23.0 14.0 NA
#> [7,] 87.0 78.0 69.0 60.0 51.0 42.0 33.0 24.0 15.0 NA
#> [8,] 88.0 79.0 70.0 61.0 52.0 43.0 34.0 25.0 16.0 NA
#> [9,] 89.0 80.0 71.0 62.0 53.0 44.0 35.0 26.0 17.0 NA
#> [10,] NA NA NA NA NA NA NA NA NA NACreated on 2020-03-19 by the reprex package (v0.3.0)
interp generates NA's where akima doesn't.
According to the akima docs, the algorithm for interpolation is not suited for collinear points as is the case here:
"The points of x and y should not be collinear, i.e, they should not fall on the same line..."
Looking at the source of akima::interp, jitter is not added when linear = TRUE (default; what interp::interp tries to replicate).
from rmodflow.
rogiersbart
commented on July 22, 2024
Any idea if this would be a fix? Or does interp::interp() fail to mimic akima::interp() in other ways too?
xy <- expand.grid(c(50, 150, 350 , 450, 550, 650, 750, 850, 950),
seq(950, 50, -100))
names(xy) <- c('x','y')
xy$z <- seq_len(nrow(xy))
akima::interp(xy$x, xy$y, xy$z, xo = seq(50, 950, 100), yo = seq(50, 950, 100))
#> $x
#> [1] 50 150 250 350 450 550 650 750 850 950
#>
#> $y
#> [1] 50 150 250 350 450 550 650 750 850 950
#>
#> $z
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 82.0 73.0 64.0 55.0 46.0 37.0 28.0 19.0 10.0 1.0
#> [2,] 83.0 74.0 65.0 56.0 47.0 38.0 29.0 20.0 11.0 2.0
#> [3,] 83.5 74.5 65.5 56.5 47.5 38.5 29.5 20.5 11.5 2.5
#> [4,] 84.0 75.0 66.0 57.0 48.0 39.0 30.0 21.0 12.0 3.0
#> [5,] 85.0 76.0 67.0 58.0 49.0 40.0 31.0 22.0 13.0 4.0
#> [6,] 86.0 77.0 68.0 59.0 50.0 41.0 32.0 23.0 14.0 5.0
#> [7,] 87.0 78.0 69.0 60.0 51.0 42.0 33.0 24.0 15.0 6.0
#> [8,] 88.0 79.0 70.0 61.0 52.0 43.0 34.0 25.0 16.0 7.0
#> [9,] 89.0 80.0 71.0 62.0 53.0 44.0 35.0 26.0 17.0 8.0
#> [10,] 90.0 81.0 72.0 63.0 54.0 45.0 36.0 27.0 18.0 9.0
rmfi_interpolate <- function(x, y, z, xo, yo) {
result <- interp::interp(x, y, z,
xo = c(xo, NA),
yo = c(yo, NA))
result$z <- result$z[-nrow(result$z), -ncol(result$z)]
result
}
rmfi_interpolate(xy$x, xy$y, xy$z,
xo = seq(50, 950, 100),
yo = seq(50, 950, 100))
#> $x
#> [1] 50 150 250 350 450 550 650 750 850 950 NA
#>
#> $y
#> [1] 50 150 250 350 450 550 650 750 850 950 NA
#>
#> $z
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 82.0 73.0 64.0 55.0 46.0 37.0 28.0 19.0 10.0 1.0
#> [2,] 83.0 74.0 65.0 56.0 47.0 38.0 29.0 20.0 11.0 2.0
#> [3,] 83.5 74.5 65.5 56.5 47.5 38.5 29.5 20.5 11.5 2.5
#> [4,] 84.0 75.0 66.0 57.0 48.0 39.0 30.0 21.0 12.0 3.0
#> [5,] 85.0 76.0 67.0 58.0 49.0 40.0 31.0 22.0 13.0 4.0
#> [6,] 86.0 77.0 68.0 59.0 50.0 41.0 32.0 23.0 14.0 5.0
#> [7,] 87.0 78.0 69.0 60.0 51.0 42.0 33.0 24.0 15.0 6.0
#> [8,] 88.0 79.0 70.0 61.0 52.0 43.0 34.0 25.0 16.0 7.0
#> [9,] 89.0 80.0 71.0 62.0 53.0 44.0 35.0 26.0 17.0 8.0
#> [10,] 90.0 81.0 72.0 63.0 54.0 45.0 36.0 27.0 18.0 9.0Created on 2020-03-20 by the reprex package (v0.3.0)
from rmodflow.
cneyens
commented on July 22, 2024
This seems to fail in following test case where y represent the vertical dimension (i.e. cross-section).
xy <- expand.grid(c(50, 150, 350 , 450, 550, 650, 750, 850, 950),
seq(-5, -25, -10))
names(xy) <- c('x','y')
xy$z <- seq_len(nrow(xy))
rmfi_interpolate <- function(x, y, z, xo, yo) {
result <- interp::interp(x, y, z,
xo = c(xo, NA),
yo = c(yo, NA))
result$z <- result$z[-nrow(result$z), -ncol(result$z)]
result$x <- result$x[-length(result$x)] # added
result$y <- result$y[-length(result$y)] # added
return(result)
}
l <- 10
akima::interp(xy$x, xy$y, xy$z, xo = seq(50, 950, 100), yo = seq(-5, -25, length = l))
#> $x
#> [1] 50 150 250 350 450 550 650 750 850 950
#>
#> $y
#> [1] -5.000000 -7.222222 -9.444444 -11.666667 -13.888889 -16.111111
#> [7] -18.333333 -20.555556 -22.777778 -25.000000
#>
#> $z
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 1.0 3.0 5.0 7.0 9.0 11.0 13.0 15.0 17.0 19.0
#> [2,] 2.0 4.0 6.0 8.0 10.0 12.0 14.0 16.0 18.0 20.0
#> [3,] 2.5 4.5 6.5 8.5 10.5 12.5 14.5 16.5 18.5 20.5
#> [4,] 3.0 5.0 7.0 9.0 11.0 13.0 15.0 17.0 19.0 21.0
#> [5,] 4.0 6.0 8.0 10.0 12.0 14.0 16.0 18.0 20.0 22.0
#> [6,] 5.0 7.0 9.0 11.0 13.0 15.0 17.0 19.0 21.0 23.0
#> [7,] 6.0 8.0 10.0 12.0 14.0 16.0 18.0 20.0 22.0 24.0
#> [8,] 7.0 9.0 11.0 13.0 15.0 17.0 19.0 21.0 23.0 25.0
#> [9,] 8.0 10.0 12.0 14.0 16.0 18.0 20.0 22.0 24.0 26.0
#> [10,] 9.0 11.0 13.0 15.0 17.0 19.0 21.0 23.0 25.0 27.0
rmfi_interpolate(xy$x, xy$y, xy$z,
xo = seq(50, 950, 100),
yo = seq(-5, -25, length = l))
#> $x
#> [1] 50 150 250 350 450 550 650 750 850 950
#>
#> $y
#> [1] -5.000000 -7.222222 -9.444444 -11.666667 -13.888889 -16.111111
#> [7] -18.333333 -20.555556 -22.777778 -25.000000
#>
#> $z
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] NA NA NA 7.0 9.0 NA NA NA NA NA
#> [2,] NA 4.0 6.0 8.0 10.0 NA NA NA NA NA
#> [3,] NA 4.5 6.5 8.5 10.5 NA NA NA NA NA
#> [4,] NA 5.0 7.0 9.0 11.0 NA NA NA NA NA
#> [5,] NA 6.0 8.0 10.0 12.0 NA NA NA NA NA
#> [6,] NA 7.0 9.0 11.0 13.0 NA NA NA NA NA
#> [7,] NA 8.0 10.0 12.0 14.0 NA NA NA NA NA
#> [8,] NA 9.0 11.0 13.0 15.0 NA NA NA NA NA
#> [9,] NA 10.0 12.0 14.0 16.0 NA NA NA NA NA
#> [10,] NA 11.0 13.0 15.0 17.0 NA NA NA NA NAIt seems to be worse when a finer output mesh is used:
l <- 20
akima::interp(xy$x, xy$y, xy$z, xo = seq(50, 950, 100), yo = seq(-5, -25, length = l))
#> $x
#> [1] 50 150 250 350 450 550 650 750 850 950
#>
#> $y
#> [1] -5.000000 -6.052632 -7.105263 -8.157895 -9.210526 -10.263158
#> [7] -11.315789 -12.368421 -13.421053 -14.473684 -15.526316 -16.578947
#> [13] -17.631579 -18.684211 -19.736842 -20.789474 -21.842105 -22.894737
#> [19] -23.947368 -25.000000
#>
#> $z
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] 1.0 1.947368 2.894737 3.842105 4.789474 5.736842 6.684211 7.631579
#> [2,] 2.0 2.947368 3.894737 4.842105 5.789474 6.736842 7.684211 8.631579
#> [3,] 2.5 3.447368 4.394737 5.342105 6.289474 7.236842 8.184211 9.131579
#> [4,] 3.0 3.947368 4.894737 5.842105 6.789474 7.736842 8.684211 9.631579
#> [5,] 4.0 4.947368 5.894737 6.842105 7.789474 8.736842 9.684211 10.631579
#> [6,] 5.0 5.947368 6.894737 7.842105 8.789474 9.736842 10.684211 11.631579
#> [7,] 6.0 6.947368 7.894737 8.842105 9.789474 10.736842 11.684211 12.631579
#> [8,] 7.0 7.947368 8.894737 9.842105 10.789474 11.736842 12.684211 13.631579
#> [9,] 8.0 8.947368 9.894737 10.842105 11.789474 12.736842 13.684211 14.631579
#> [10,] 9.0 9.947368 10.894737 11.842105 12.789474 13.736842 14.684211 15.631579
#> [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
#> [1,] 8.578947 9.526316 10.47368 11.42105 12.36842 13.31579 14.26316 15.21053
#> [2,] 9.578947 10.526316 11.47368 12.42105 13.36842 14.31579 15.26316 16.21053
#> [3,] 10.078947 11.026316 11.97368 12.92105 13.86842 14.81579 15.76316 16.71053
#> [4,] 10.578947 11.526316 12.47368 13.42105 14.36842 15.31579 16.26316 17.21053
#> [5,] 11.578947 12.526316 13.47368 14.42105 15.36842 16.31579 17.26316 18.21053
#> [6,] 12.578947 13.526316 14.47368 15.42105 16.36842 17.31579 18.26316 19.21053
#> [7,] 13.578947 14.526316 15.47368 16.42105 17.36842 18.31579 19.26316 20.21053
#> [8,] 14.578947 15.526316 16.47368 17.42105 18.36842 19.31579 20.26316 21.21053
#> [9,] 15.578947 16.526316 17.47368 18.42105 19.36842 20.31579 21.26316 22.21053
#> [10,] 16.578947 17.526316 18.47368 19.42105 20.36842 21.31579 22.26316 23.21053
#> [,17] [,18] [,19] [,20]
#> [1,] 16.15789 17.10526 18.05263 19.0
#> [2,] 17.15789 18.10526 19.05263 20.0
#> [3,] 17.65789 18.60526 19.55263 20.5
#> [4,] 18.15789 19.10526 20.05263 21.0
#> [5,] 19.15789 20.10526 21.05263 22.0
#> [6,] 20.15789 21.10526 22.05263 23.0
#> [7,] 21.15789 22.10526 23.05263 24.0
#> [8,] 22.15789 23.10526 24.05263 25.0
#> [9,] 23.15789 24.10526 25.05263 26.0
#> [10,] 24.15789 25.10526 26.05263 27.0
rmfi_interpolate(xy$x, xy$y, xy$z,
xo = seq(50, 950, 100),
yo = seq(-5, -25, length = l))
#> $x
#> [1] 50 150 250 350 450 550 650 750 850 950
#>
#> $y
#> [1] -5.000000 -6.052632 -7.105263 -8.157895 -9.210526 -10.263158
#> [7] -11.315789 -12.368421 -13.421053 -14.473684 -15.526316 -16.578947
#> [13] -17.631579 -18.684211 -19.736842 -20.789474 -21.842105 -22.894737
#> [19] -23.947368 -25.000000
#>
#> $z
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
#> [1,] NA 1.947368 2.894737 NA NA NA NA NA NA NA NA NA
#> [2,] NA 2.947368 3.894737 NA NA NA NA NA NA NA NA NA
#> [3,] NA 3.447368 4.394737 NA NA NA NA NA NA NA NA NA
#> [4,] NA 3.947368 4.894737 NA NA NA NA NA NA NA NA NA
#> [5,] NA 4.947368 5.894737 NA NA NA NA NA NA NA NA NA
#> [6,] NA 5.947368 6.894737 NA NA NA NA NA NA NA NA NA
#> [7,] NA 6.947368 7.894737 NA NA NA NA NA NA NA NA NA
#> [8,] NA 7.947368 8.894737 NA NA NA NA NA NA NA NA NA
#> [9,] NA 8.947368 9.894737 NA NA NA NA NA NA NA NA NA
#> [10,] NA 9.947368 10.894737 NA NA NA NA NA NA NA NA NA
#> [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
#> [1,] NA NA NA NA NA NA NA NA
#> [2,] NA NA NA NA NA NA NA NA
#> [3,] NA NA NA NA NA NA NA NA
#> [4,] NA NA NA NA NA NA NA NA
#> [5,] NA NA NA NA NA NA NA NA
#> [6,] NA NA NA NA NA NA NA NA
#> [7,] NA NA NA NA NA NA NA NA
#> [8,] NA NA NA NA NA NA NA NA
#> [9,] NA NA NA NA NA NA NA NA
#> [10,] NA NA NA NA NA NA NA NACreated on 2020-03-25 by the reprex package (v0.3.0)
from rmodflow.
Related Issues (17)
- Error reading .dis file HOT 2
- RMODFLOW extension packages
- Consistency with ModelMuse HOT 2
- Add examples in function documentation
- Speed-up reading/writing arrays
- Reading recharge package fails when time-varying parameters are present HOT 3
- Weighted goodness-of-fit measures for use in optimization
- Automatic installation of MODFLOW variants, and finding paths to executables
- Package object structure: Dimensions, and maybe other things HOT 5
- Consistent S3 class names HOT 2
- Bug in rmf_create_hob (and rmf_write_hob) when dealing with transient observations HOT 1
- Metadata HOT 3
- Basic MODFLOW package/file object functions HOT 5
- Install codes and example models HOT 1
- Execute, optimize, analyze HOT 2
- rmf_convert_grid_to_xyz works not correctly HOT 3
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from rmodflow.