interpIDW offset issue.

[h-a-graham]

Hey,

I recently stumbled across this issue where it appears that increasing the radius argument in interpIDW results in an increased offset between the original and interpolated raster.

I've put together a minimal example below (terra v ‘1.7.35’):

library(terra)
#> terra 1.7.29
pfunc <- function(.x, .idw) {
    par(mfrow = c(1, 2))
    plot(.x)
    plot(.idw)
    par(mfrow = c(1, 1))
}
x <- rast(system.file("ex/logo.tif", package="terra"))$red
x[x>200] <- NA

idw5 <- interpIDW(x, as.points(x), field="red", radius=5, power=3, smooth=0, maxPoints=5)
pfunc(x, idw5)

idw25 <- interpIDW(x, as.points(x), field="red", radius=25, power=3, smooth=0, maxPoints=5)
pfunc(x, idw25)

^{Created on 2023-06-15 with reprex v2.0.2}

As always many thanks for this excellent package!

[rhijmans]

Thanks. The problem seems to stem from setting the maxpoints argument with near=FALSE. The results look good when omitting the argument (using the default maxpoints=Inf).

idw25 <- interpIDW(x, as.points(x), field="red", radius=25, power=3, smooth=0)
pfunc(x, idw25)

I believe this is a bug in GDALGRID

[h-a-graham]

That's great, many thanks!

[rhijmans]

As suggested by Even Roualt, this is caused by taking the first five points that are within the circle; in order of appearance. For that reason you do not see this effect if you first randomize the order of the points

p  <- as.points(x)  |> sample()
idw5 <- interpIDW(x, p, field="red", radius=5, power=3, smooth=0, maxPoints=5)

[h-a-graham]

That is super interesting! That solution solves things for me. I'd be happy to submit a PR with some suggestions for docs if it's useful but I suppose it is a fairly specific issue that occurs where the points to interpolate are derived from a grid? perhaps something in details along the lines of:

"Where maxPoints is finite and y contains point locations that are orderd by spatial location or derived from a grid, the point ordering should be randomised to avoid offset values in the retutned SpatRaster."

Cheers!

[rhijmans]

I think the solution is to use near=TRUE and I have made that the default. If you read the manual carefully you will see that with near=FALSE you will get up to maxPoints observations, but there is no guarantee that these will be near to the target grid cell. That is not intuitive.
However, there was bug when using near=TRUE that I have now fixed.

I now get:

y <- interpIDW(x, p, field="red", radius=10, power=3, maxPoints=5, near=T)
pfunc(x, y)

Thank you for reporting this and for the nice repex.

[h-a-graham]

Perfect thanks so much for resolving this so quickly! Have a good weekend.