Inverse Distance Weighting with gstat

Question

I am attempting to execute an idw interpolation with gstat but am unable to.

text =
name,long,lat,waterLevel,elevation,depth

EM_01,18.553392,-34.07027,14.4,20.358,63.0

EM_27,18.574777,-34.068709,16.196,19.966,48.0

EM_29,18.613985,-34.053271,18.766,25.477,39.0

EM_20,18.654089,-34.045177,20.102,36.502,45.0

EM_23,18.643495,-34.037468,22.915,32.715,33.0

EM_13,18.637267,-34.029194,25.516,29.716,33.0

install.packages("sf")
install.packages("gstat")
install.packages("units")
install.packages("terra")

library(sf) 
library(gstat) # geostatistics
library(units) # units of measure
library(terra) 

file = 'text.txt'

#-- import
cfaq <- read.csv(file, header = 1, sep = ',', dec = '.') 

#- set as a spatial feature with xy coords an existing projection
cfaq.sf <- st_as_sf(cfaq, coords=c("long", "lat"), crs = 4326) 

#- transform to local crs
cfaq.sf <- st_transform(cfaq.sf, crs = 32734)

#-- construct empty grid
(n.col <- length(seq.e <- seq(min.x <- floor(min(st_coordinates(cfaq.sf)[,2])/1000)*1000,
                              max.x <- ceiling(max(st_coordinates(cfaq.sf)[,2])/1000)*1000, by=1000)))

(n.row <- length(seq.n <- seq(min.y <- floor(min(st_coordinates(cfaq.sf)[,1])/1000)*1000,
                              max.y <- ceiling(max(st_coordinates(cfaq.sf)[,1])/1000)*1000, by=1000)))


#- 2km grid
grid2km <- rast(nrows = n.row, ncols = n.col,
                xmin=min.x, xmax=max.x,
                ymin=min.y, ymax=max.y, crs = st_crs(cfaq.sf)$proj4string,
                resolution = 2000, names="waterLevel")

values(grid2km) <- NA_real_

#-- grid to df
grid2km.df <- as.data.frame(grid2km, xy = TRUE, na.rm = FALSE)
names(grid2km.df)[1:2] <- c("X", "Y") # match the names of the point dataset

#- idw
grid2km.df$idwp05 = idw(log(waterLevel)  ~ 1, cfaq.sf, 
                        grid2km.df, idp = 0.5)$var1.pred

The idw does not succeed.
Error in .local(formula, locations, ...) : sf::st_crs(locations) == sf::st_crs(newdata) is not TRUE

How can this be? The crs is defined. How do I execute the idw successfully please?

score 0 · Accepted Answer · answered Aug 10 '23 at 13:18

That error is not wrong; though note that it does not complain about CRS of cfaq.sf (locations param). You are using grid2km.df for newdata parameter, and that's the one with no CRS. Besides, newdata should should be Spatial, sf or stars, so perhaps try to use stars for grid:

library(sf) 
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(gstat) 
library(stars)

text_ <- "name,long,lat,waterLevel,elevation,depth
          EM_01,18.553392,-34.07027,14.4,20.358,63.0
          EM_27,18.574777,-34.068709,16.196,19.966,48.0
          EM_29,18.613985,-34.053271,18.766,25.477,39.0
          EM_20,18.654089,-34.045177,20.102,36.502,45.0
          EM_23,18.643495,-34.037468,22.915,32.715,33.0
          EM_13,18.637267,-34.029194,25.516,29.716,33.0"

# -- import
cfaq <- read.csv(text = text_, header = 1, sep = ',', dec = '.') 

# - set as a spatial feature with xy coords an existing projection
cfaq.sf <- st_as_sf(cfaq, coords=c("long", "lat"), crs = 4326) 

# - transform to local crs
cfaq.sf <- st_transform(cfaq.sf, crs = 32734)

cfaq.bb <- st_bbox(cfaq.sf)
cfaq.bb[c("xmin","ymin")] <-   floor(cfaq.bb[c("xmin","ymin")]/1000)*1000
cfaq.bb[c("xmax","ymax")] <- ceiling(cfaq.bb[c("xmax","ymax")]/1000)*1000

grid2km <- st_as_stars(cfaq.bb, dx = 2000)
grid2km
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>         Min. 1st Qu. Median Mean 3rd Qu. Max.
#> values     0       0      0    0       0    0
#> dimension(s):
#>   from to  offset delta                refsys x/y
#> x    1  5  274000  2000 WGS 84 / UTM zone 34S [x]
#> y    1  3 6233000 -2000 WGS 84 / UTM zone 34S [y]

i <- idw(log(waterLevel) ~ 1, cfaq.sf, grid2km, idp = 0.5)
#> [inverse distance weighted interpolation]
as.data.frame(i)
#>         x       y var1.pred var1.var
#> 1  275000 6232000  2.941439       NA
#> 2  277000 6232000  2.957383       NA
#> 3  279000 6232000  2.982699       NA
#> 4  281000 6232000  3.031934       NA
#> 5  283000 6232000  3.036951       NA
#> 6  275000 6230000  2.918279       NA
#> 7  277000 6230000  2.938086       NA
#> 8  279000 6230000  2.963870       NA
#> 9  281000 6230000  2.995139       NA
#> 10 283000 6230000  3.010777       NA
#> 11 275000 6228000  2.875184       NA
#> 12 277000 6228000  2.907489       NA
#> 13 279000 6228000  2.948413       NA
#> 14 281000 6228000  2.972109       NA
#> 15 283000 6228000  2.986109       NA

^{Created on 2023-08-10 with reprex v2.0.2}

Thank you @margusl. I would appreciate the ability to interpolate into the df ---like the code originally wants: grid2km$idw05, etc.--- how can I assign a `crs` to the `df` please. I am failing too with `grid2km.df <- st_set_crs(grid.df, st_crs(32734))` — arkriger, Aug 11 '23 at 13:37
Maybe a not the best choice of words from me, but you can't . Error you included in your question states that in your current settings, `grid2km.df` is evaluated as `newdata` argument, but that can't be a data.frame, it must be either `Spatial`, `sf` or `stars`. You can check `?gstat::idw` for reference. — margusl, Aug 11 '23 at 15:35

Inverse Distance Weighting with gstat

1 Answers1