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