Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Confusion about the unit for maxdist in gstat? #113

Open
MatthieuStigler opened this issue Sep 7, 2022 · 2 comments
Open

Confusion about the unit for maxdist in gstat? #113

MatthieuStigler opened this issue Sep 7, 2022 · 2 comments

Comments

@MatthieuStigler
Copy link

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)

@edzer
Copy link
Member

edzer commented Sep 7, 2022

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?

@MatthieuStigler
Copy link
Author

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!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants