terra:interpNear interpolates in a square window, not in a circle, beyond 100 points

[jldupouey]

I get a result I don't understand, with interpNear. As soon as the number of starting points exceeds 100, it interpolates in square windows, not in circles (and therefore does not respect the radius provided). Using the example given in the documentation:

library(terra)
r <- rast(ncol=100, nrow=100, crs="local", xmin=0, xmax=50, ymin=0, ymax=50)
set.seed(100)
x <- runif(101, 10, 40)
y <- runif(101, 10, 40)
z <- sample(101)
xyz <- cbind(x,y,z)

x <- interpNear(r, xyz, radius=5)
plot(x)

It works well below 100 points:

xyz100 <- xyz[-1, ]

x <- interpNear(r, xyz100, radius=5)
plot(x)

Perhaps I've misunderstood something? Thanks,

Jean-Luc

[rhijmans]

Thank you for reporting this. This calls GDALGridCreate and I see that in the source code it has a parameter nPointCountThreshold created like this:

const unsigned int nPointCountThreshold =
        atoi(CPLGetConfigOption("GDAL_GRID_POINT_COUNT_THRESHOLD", "100"));

If the number of points is above this threshold a Quadtree based algorithm is used. That explains the variation you see, and leads to this work-around in R:

setGDALconfig("GDAL_GRID_POINT_COUNT_THRESHOLD", nrow(xyz))
x <- interpNear(r, xyz, radius=5)

I could of course implement this but at this point, I am not sure if this is a bug in GDAL; or an undocumented option to speed up the computation (with a, perhaps undesired, side-effect).

[jldupouey]

Thanks, that's clear. The terra::interpNear documentation states that the radius parameter is “The radius of the circle”. Perhaps it should be indicated that it is no longer a circle beyond 100 points or, better still, indicate in the documentation this problem and give the trick to get around it.