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

Question

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 1

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(
  '<OGRVRTDataSource>',
  '\t<OGRVRTLayer name="var">',
  paste("\t<SrcDataSource>", tempfname_csv, "</SrcDataSource>", sep = ""),
  '\t<GeometryType>wkbPoint</GeometryType>',
  '\t<GeometryField encoding="PointFromColumns" x="Easting" y="Northing" z="Variable"/>',
  '\t</OGRVRTLayer>',
  '\t</OGRVRTDataSource>'
)

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)