How to do spatial interpolation that extends beyond location of sample points?

I think this is a common question and, having just hit upon it, thought it worth asking the issue here. Reproducible example:

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
library(gstat)
library(stars)
#> Loading required package: abind
u = "https://github.com/saferactive/saferactive/releases/download/0.1.1/rnet_lnd_1pcnt.Rds"
rnet_lnd = readRDS(url(u))
rnet_lnd
#> Simple feature collection with 698 features and 9 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 506201.1 ymin: 157489.3 xmax: 554024 ymax: 197327
#> projected CRS:  OSGB 1936 / British National Grid
#> First 10 features:
#>       local_id bicycle govtarget_slc govnearmkt_slc gendereq_slc dutch_slc
#> 36469    36469       2             8              4            4        38
#> 72154    72154     296           464            509          415       992
#> 13692    13692      64            75             85           99       170
#> 51194    51194      10            10             15           20        38
#> 45023    45023       5            20             14           10        78
#> 14664    14664     140           166            189          198       270
#> 62204    62204      31            41             48           49       118
#> 48161    48161       7            18             16           14        83
#> 37306    37306       2            17             12            3        88
#> 8329      8329      10            34             32           19       144
#>       ebike_slc    length km_cycled_yr                  geometry
#> 36469        40  75.14911     76.05090 POINT (540554.4 184261.5)
#> 72154      1191  21.59051   3233.73959 POINT (532707.7 182641.9)
#> 13692       224  50.93146   1649.36427 POINT (526445.2 177201.8)
#> 51194        50 236.08443   1194.58721 POINT (533507.4 186557.7)
#> 45023        91 211.86509    536.01868 POINT (541408.3 180908.9)
#> 14664       342  14.72725   1043.27843 POINT (532361.8 186576.6)
#> 62204       193 550.77091   8639.39256 POINT (536944.7 176095.4)
#> 48161       104 194.18121    687.78986   POINT (534208 190667.8)
#> 37306        98  52.39301     53.02172 POINT (534251.2 195203.9)
#> 8329        249  81.10461    410.38932 POINT (540147.3 192165.4)
grd = st_bbox(rnet_lnd) %>%
  st_as_stars(dx = 500, dy = 500) %>%
  st_set_crs(27700) %>%
  st_crop(rnet_lnd) 
grd
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>     values     
#>  Min.   :0     
#>  1st Qu.:0     
#>  Median :0     
#>  Mean   :0     
#>  3rd Qu.:0     
#>  Max.   :0     
#>  NA's   :7051  
#> dimension(s):
#>   from to offset delta                       refsys point values    
#> x    1 96 506201   500 OSGB 1936 / British Natio...    NA   NULL [x]
#> y    1 80 197327  -500 OSGB 1936 / British Natio...    NA   NULL [y]
v = variogram(bicycle~1, rnet_lnd, cutoff = 10000)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
#> = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
#> definition
#> Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum OSGB 1936 in CRS definition
plot(v)

vm = fit.variogram(v, vgm(psill = "Sph", model = "Exp", range = 10000, nugget = 1))
plot(vm, cutoff = 10000)

rnet_krige = gstat::krige(bicycle~1, rnet_lnd, grd, vm, nmax = 100)
#> [using ordinary kriging]
plot(rnet_lnd$geometry)
plot(rnet_krige, add = TRUE)

^{Created on 2020-10-26 by the reprex package (v0.3.0)}

Any ideas @Nowosad ? I see you have done some work on gstat. FYI I've read the documentation here https://keen-swartz-3146c4.netlify.app/interpolation.html

The problem is that you have simply copy-pasted the code, and you replaced de which is a polygon layer by your point layer rnet_lnd:

grd = st_bbox(rnet_lnd) %>%
  st_as_stars(dx = 500, dy = 500) %>%
  st_set_crs(27700) %>%
  st_crop(rnet_lnd)

Just get rid off the st_crop(rnet_lnd) line and you will get a raster which covers the bounding box of your bicycle points. Please note that ordinary kriging is quite a simple approach, basically only taking distance into account, maybe you would like to use universal kriging, however, for this you would need a predictor which you also have spatially available, i.e., as a raster or stars objects.

Confirmed, that fixes it. Many thanks for the quick reply. Feeling :facepalm: now but happy it's a simple solution!

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
library(gstat)
library(stars)
#> Loading required package: abind
u = "https://github.com/saferactive/saferactive/releases/download/0.1.1/rnet_lnd_1pcnt.Rds"
rnet_lnd = readRDS(url(u))
rnet_lnd
#> Simple feature collection with 698 features and 9 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 506201.1 ymin: 157489.3 xmax: 554024 ymax: 197327
#> projected CRS:  OSGB 1936 / British National Grid
#> First 10 features:
#>       local_id bicycle govtarget_slc govnearmkt_slc gendereq_slc dutch_slc
#> 36469    36469       2             8              4            4        38
#> 72154    72154     296           464            509          415       992
#> 13692    13692      64            75             85           99       170
#> 51194    51194      10            10             15           20        38
#> 45023    45023       5            20             14           10        78
#> 14664    14664     140           166            189          198       270
#> 62204    62204      31            41             48           49       118
#> 48161    48161       7            18             16           14        83
#> 37306    37306       2            17             12            3        88
#> 8329      8329      10            34             32           19       144
#>       ebike_slc    length km_cycled_yr                  geometry
#> 36469        40  75.14911     76.05090 POINT (540554.4 184261.5)
#> 72154      1191  21.59051   3233.73959 POINT (532707.7 182641.9)
#> 13692       224  50.93146   1649.36427 POINT (526445.2 177201.8)
#> 51194        50 236.08443   1194.58721 POINT (533507.4 186557.7)
#> 45023        91 211.86509    536.01868 POINT (541408.3 180908.9)
#> 14664       342  14.72725   1043.27843 POINT (532361.8 186576.6)
#> 62204       193 550.77091   8639.39256 POINT (536944.7 176095.4)
#> 48161       104 194.18121    687.78986   POINT (534208 190667.8)
#> 37306        98  52.39301     53.02172 POINT (534251.2 195203.9)
#> 8329        249  81.10461    410.38932 POINT (540147.3 192165.4)
grd = st_bbox(rnet_lnd) %>%
  st_as_stars(dx = 500, dy = 500) %>%
  st_set_crs(27700) 
grd
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>     values  
#>  Min.   :0  
#>  1st Qu.:0  
#>  Median :0  
#>  Mean   :0  
#>  3rd Qu.:0  
#>  Max.   :0  
#> dimension(s):
#>   from to offset delta                       refsys point values    
#> x    1 96 506201   500 OSGB 1936 / British Natio...    NA   NULL [x]
#> y    1 80 197327  -500 OSGB 1936 / British Natio...    NA   NULL [y]
v = variogram(bicycle~1, rnet_lnd, cutoff = 10000)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
#> = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
#> definition
#> Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum OSGB 1936 in CRS definition
plot(v)

vm = fit.variogram(v, vgm(psill = "Sph", model = "Exp", range = 10000, nugget = 1))
plot(vm, cutoff = 10000)

rnet_krige = gstat::krige(bicycle~1, rnet_lnd, grd, vm, nmax = 100)
#> [using ordinary kriging]
plot(rnet_lnd$geometry)
plot(rnet_krige, add = TRUE)

^{Created on 2020-10-26 by the reprex package (v0.3.0)}

No worries. I think it was a good question, and there is really no need to feel :facepalm:.

@jannes-m made a lot of great points. I also added some tiny changes to the code: (1) I cropped the map to the data area, (2) plotted model together with the points, (3) updated the map:

library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.0.4, PROJ 6.3.2
library(gstat)
library(stars)
#> Loading required package: abind
library(tmap)
library(concaveman)
u = "https://github.com/saferactive/saferactive/releases/download/0.1.1/rnet_lnd_1pcnt.Rds"
rnet_lnd = readRDS(url(u))
# rnet_lnd

grd = st_bbox(rnet_lnd) %>%
  st_as_stars(dx = 250, dy = 250) %>%
  st_set_crs(27700) %>%
  st_crop(concaveman(rnet_lnd))
# grd

v = variogram(bicycle~1, rnet_lnd, cutoff = 10000)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
#> = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
#> definition
#> Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum OSGB 1936 in CRS definition
plot(v)

vm = fit.variogram(v, vgm(model = "Sph", range = 5000, nugget = 1))
plot(v, model = vm)

rnet_krige = krige(bicycle~1, rnet_lnd, grd, vm, nmax = 30)
#> [using ordinary kriging]

tm_shape(rnet_krige[1, ]) +
        tm_grid() +
        tm_raster(style = "cont",
                  breaks = c(0, 10, 100, 200, 300, 500, 1000, 1800),
                  palette = rev(hcl.colors(n = 99, palette = "Ag_Sunset")),
                  title = "Bicycles:") +
        tm_shape(rnet_lnd) +
        tm_dots() +
        tm_layout(legend.outside = TRUE)
#> Warning: Values have found that are less than the lowest break

^{Created on 2020-10-27 by the reprex package (v0.3.0)}

Very nice. Good basis for a section in Geocomputation with R version 2?

geocompx / geocompr

How to do spatial interpolation that extends beyond location of sample points? #585