gstat icon indicating copy to clipboard operation
gstat copied to clipboard

Confusion about the unit for maxdist in gstat?

Open MatthieuStigler opened this issue 3 years ago • 2 comments

I am confused on the unit in which maxdist is to be specified? It seems the answer depends on whether the data is projected or not?

In the example below, I find heuristically that I need to give in [m] for projected data, and possibly in [km] for unprojected? But I might as well be confused! The related question would be how to specify maxdist in meter whenever one uses un-projected data?

Thanks a lot!

library(gstat)
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.0; sf_use_s2() is TRUE

nc = st_read(system.file("shape/nc.shp", package="sf"), quiet=TRUE)# %>% 

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)) 

nc_points_prj <- st_transform(nc_points, "ESRI:102008")
nc_grid_prj <- st_transform(nc_grid, "ESRI:102008")



## predict
gs_out <- gstat(formula = x ~ 1, data = nc_points, maxdist = 80)
gs_out_prj <- gstat(formula = x ~ 1, data = nc_points_prj, maxdist = 80000)
suppressWarnings(predict(gs_out, nc_grid)$var1.pred)
#> [inverse distance weighted interpolation]
#> [1]        NA        NA        NA        NA        NA        NA        NA
#> [8] 0.8609154        NA
suppressWarnings(predict(gs_out_prj, nc_grid_prj)$var1.pred)
#> [inverse distance weighted interpolation]
#> [1]        NA        NA        NA        NA        NA        NA        NA
#> [8] 0.8609154        NA

Created on 2022-09-07 by the reprex package (v2.0.1)

MatthieuStigler avatar Sep 07 '22 15:09 MatthieuStigler

in [m] for projected data, and possibly in [km] for unprojected?

Yes, that is right, and I agree it is confusing. Today I would do this differently; suporting units from the units package might be a way out?

edzer avatar Sep 07 '22 15:09 edzer

thanks for the fast answer!

I can imagine using units might be more elegant indeed. In the meanwhile, maybe providing a discussion in the help file might be useful?

PS: Also, is the warning message in In proj4string(obj) : CRS object has comment, which is lost in output; predict(gs_out, nc_grid) necessary in this context? Thanks!

MatthieuStigler avatar Sep 07 '22 16:09 MatthieuStigler