r-spatial / gstat

Spatial and spatio-temporal geostatistical modelling, prediction and simulation
http://r-spatial.github.io/gstat/
GNU General Public License v2.0
195 stars 49 forks source link

Can't krige with 3 variables, only with 2 #119

Closed CGenie closed 1 year ago

CGenie commented 1 year ago

Hello,

I have this function:

krige <- function(data, bbox, xlength = 256, ylength = 256, nmax=15, maxdist=200) {
    data_sf <- st_as_sf(data,
                        coords=c("lon", "lat"),
                        crs=4326,
                        agr="constant") %>%
        as_Spatial(.)

    dt.fit <- automap::autofitVariogram(value~coords.x1+coords.x2,
                                        data_sf)

    grd <- tidyr::expand_grid(
                      lon = seq(bbox$xmin, bbox$xmax, length.out = xlength),
                      lat = seq(bbox$ymin, bbox$ymax, length.out = ylength)
                  )
    grd_sf <- grd %>%
        st_as_sf(coords=c("lon", "lat"),
                 crs=4326,
                 agr="constant") %>%
        as_Spatial(.)

    result <- gstat::krige(value~coords.x1+coords.x2+height,  # formula
                           data_sf,  # locations
                           grd_sf,  # data
                           model=dt.fit$var_model,
                           nmax=nmax,
                           maxdist=maxdist)

    list(xlength = xlength,
         ylength = ylength,
         result = result,
         sf = data_sf,
         grd = grd)
}

and I call it with this simple data:

krige(data.frame(lat=c(50.4400978, 50.18796, 50.50484, 50.3132324, 53.33001), lon=c(17.9800358, 17.81987, 18.2546, 17.5486336, 16.308094), height=c(171.57905308564648, 284.1249314275102, 247.21124382499872, 270.71395841781765, 122.6001442065739), value=c(4.12, 2.46, 4.43, 4.09, 3.56)), list(xmin=50.4, xmax=50.6, ymin=17.8, ymax=18.0))

However, I'm getting:

Error in eval(predvars, data, env) : object 'height' not found

This seems to be caused by the formula used in gstat::krige. If I change it to value~coords.x1+coords.x2 then things work, but I don't really like the result, as the height does seem to be an important component here.

Could anyone tell me what's wrong with my code can't it find the height element even if it's there? Here's the output of summary(data_sf) for my data:

Object of class SpatialPointsDataFrame
Coordinates:
               min      max
coords.x1 16.30809 18.25460
coords.x2 50.44010 53.33001
Is projected: FALSE 
proj4string : [+proj=longlat +datum=WGS84 +no_defs]
Number of points: 5
Data attributes:
     height          value      
 Min.   :122.6   Min.   :2.460  
 1st Qu.:171.6   1st Qu.:3.560  
 Median :247.2   Median :4.090  
 Mean   :219.2   Mean   :3.732  
 3rd Qu.:270.7   3rd Qu.:4.120  
 Max.   :284.1   Max.   :4.430  
edzer commented 1 year ago

This is not a question about this package, but a question about writing R functions.