How to perform an anisotropic Elliptical Inverse Distance Weighting (EIDW) interpolation in R?

River bathymetry is generally more variable in the transverse than the streamwise direction, so I need to use anisotropic interpolation technique like Elliptical Inverse Distance Weighting (EIDW) (see Anisotropic considerations while interpolating river channel bathymetry) to construct a digital elevation model (DEM) of the surveyed channel. I only found how to perform an Inverse Distance Weighted using the gstat package in R:

library(sp)
data(meuse)
coordinates(meuse) = ~x+y
data(meuse.grid)
gridded(meuse.grid) = ~x+y

# inverse distance weighted:
x <- krige(log(zinc)~1, meuse, meuse.grid)

Answer

I could apply an Elliptical Inverse Distance Weighting (EIDW) in R using the function gdal_grid of the gdalUtils package. I modified a little bit the gdal_grid help example in R documentation as follows:

# Load packages
library('sp')
library('gdalUtils')

# Load data
data('meuse')
coordinates(meuse) = ~x+y

# Create a properly formatted CSV:
temporary_dir <- tempdir()

tempfname_base <- file.path(temporary_dir, "var")
tempfname_csv <- paste(tempfname_base, ".csv", sep = "")

pts <- data.frame('Easting' = meuse@coords[,1], 'Northing' = meuse@coords[,2], 'Variable' = meuse@data$cadmium)

write.csv(pts, file = tempfname_csv, row.names = FALSE)

# Now make a matching VRT file
tempfname_vrt <- paste(tempfname_base, ".vrt", sep = "")


vrt_header <- c(
  '',
  '\t',
  paste("\t", tempfname_csv, "", sep = ""),
  '\twkbPoint',
  '\t',
  '\t',
  '\t'
)


vrt_filecon <- file(tempfname_vrt, "w")
writeLines(vrt_header, con = vrt_filecon)
close(vrt_filecon)

tempfname_tif <- paste(tempfname_base, ".tiff", sep = "")

# Interpolation algorithm and parameters
parameters <- "invdist:power=2.0:smoothing=1.0:radius1=120:radius2=360:angle=45:max_points=0"

# Now run gdal_grid:

meuseEIDW <- gdal_grid(src_datasource = tempfname_vrt, dst_filename = tempfname_tif, of = "GTiff", ot = "Float64", a = parameters, l = "var", outsize = c(40, 40), output_Raster = TRUE)

# Plot
spplot(obj = meuseEIDW)