gstat idw: is the point to polygon interpolation with maxdist restricted to centroid-covered polygons?

[MatthieuStigler]

I am trying to understand how the point to polygon interpolation with gstat is working wiht maxdist. My intuition was that it would be equivalent to compute the distance say with st_dist, thresholding to NA distances d>maxdist, then doing the weighted mean. However, I notice that polygons that have a d<maxdist might still get NA values. I suspect there is an inclusion criterion, whereby a point only gives a positive weights if d<maxdist AND d(point, poly centroid)<maxdist, is that correct?

If yes, what is the difference between doing to to point and point to polygon? It still seems the results are different, so it is not just that gstat is using a point to polygon centroid concept?

Thanks!!

In the example below, I expected the two left cells in the second row to receive an interpolating value, since they cover buffer of points?

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(gstat)
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.0; sf_use_s2() is TRUE

## prepare points and grid
nc = st_read(system.file("shape/nc.shp", package="sf"), quiet=TRUE) %>% 
  st_transform("ESRI:102008")
set.seed(1234)
nc_points = st_sf(x=runif(6), geometry=st_sample(nc[1:3, ], 6) )
nc_grid <- st_as_sf(st_make_grid(nc, n = 3)) %>% 
  rename(geometry=x)
nc_grid_pts <- st_centroid(nc_grid)

## visual buffedr
thresh <- 80000
thresh <- 100000
nc_points_buf <- nc_points %>% st_buffer(thresh) %>% st_geometry() %>% 
  st_union()



## run gstat on sf objects
gs_out <- gstat(formula = x ~ 1, data = nc_points, 
                maxdist = thresh,
                set = list(idp = 2))
z <- suppressWarnings(predict(gs_out, nc_grid))
#> [inverse distance weighted interpolation]

plot(z %>% st_geometry())
plot(z %>% select(var1.pred), add=TRUE)
plot(nc_points, add=TRUE)
plot(nc_points_buf, add=TRUE)
plot(nc_grid_pts, add=TRUE)

^{Created on 2022-09-13 with reprex v2.0.2}

[edzer]

From the top of my head the distances are evaluated from the point to the centroid of the target polygon, not as the distance to the polygon (boundary). Would that answer your question?

[MatthieuStigler]

Thanks, this would answer the first part, yet leave open the second question: if gstat takes centroids of polygons, why does one see (arguably small) differences compared to when taking explicitly the centroid?

library(dplyr, warn.conflicts=FALSE)
library(gstat)
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.0; sf_use_s2() is TRUE

## prepare points and grid
nc = st_read(system.file("shape/nc.shp", package="sf"), quiet=TRUE) %>% 
  st_transform("ESRI:102008")
set.seed(1234)
nc_points = st_sf(x=runif(6), geometry=st_sample(nc[1:3, ], 6) )
nc_grid <- st_as_sf(st_make_grid(nc, n = 3)) %>% 
  rename(geometry=x)
nc_grid_pts <- st_centroid(nc_grid)


## run gstat on sf objects
thresh <- 150000
gs_out <- gstat(formula = x ~ 1, data = nc_points, 
                maxdist = thresh,
                set = list(idp = 2))
z <- suppressWarnings(predict(gs_out, nc_grid)) |> pull(var1.pred)
#> [inverse distance weighted interpolation]
z2 <- suppressWarnings(predict(gs_out, nc_grid_pts))|> pull(var1.pred)
#> [inverse distance weighted interpolation]

z
#> [1]        NA        NA        NA 0.4275601 0.7097799        NA 0.4910593
#> [8] 0.6872785        NA
z2
#> [1]        NA        NA        NA 0.4443680 0.7113301        NA 0.4918421
#> [8] 0.6921974        NA

## run gstat on sf objects
gs_out <- gstat(formula = x ~ 1, data = nc_points, 
                maxdist = 200000,
                set = list(idp = 2))
z <- suppressWarnings(predict(gs_out, nc_grid)) |> pull(var1.pred)
#> [inverse distance weighted interpolation]
z2 <- suppressWarnings(predict(gs_out, nc_grid_pts))|> pull(var1.pred)
#> [inverse distance weighted interpolation]

z
#> [1]        NA        NA        NA 0.5221320 0.6160211        NA 0.5491632
#> [8] 0.6284615        NA
z2
#> [1]        NA        NA        NA 0.5267554 0.6262449        NA 0.5469904
#> [8] 0.6352545        NA

^{Created on 2022-09-13 with reprex v2.0.2}