Skip to content

Confusion about the unit for maxdist in gstat? #113

@MatthieuStigler

Description

@MatthieuStigler

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)

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions