Terra rasterize interpolation output produces horizontal and vertical lines in output raster

Question

I'm trying to produce a continuous interpolated surface from point data. I have a data frame of longitude, latitude and response values (mud). I've tried both Inverse Distance Weighting and Ordinary Kriging but get the same problem. Here's some reproducible code below:

library(terra)
library(sp)
library(gstat)

extent <- c(-7.67151099570083, -2.34424614123549, 50.159278977608, 56.2717492369751)
blankRaster <- rast(ext(extent))
cell_size <- c(0.00117, 0.002075) * 10
res(blankRaster) <- cell_size

x_range <- c(ext(blankRaster)[1], ext(blankRaster)[2])
y_range <- c(ext(blankRaster)[3], ext(blankRaster)[4])

grd <- expand.grid(x = seq(from = x_range[1], to = x_range[2], by = cell_size[1]), y = seq(from = y_range[1], to = y_range[2], by = cell_size[2]))
coordinates(grd) <- ~ x+y
gridded(grd) <- T

lon <- c(-6.054578, -6.030297, -6.148603, -6.001423, -6.029923, -6.041965, -6.026100, -6.013715, -6.027235, -6.026288, -6.004180, -6.020842, -5.997938, -6.030913, -6.025397, -6.107727, -6.089913, -6.054938, -6.042965, -6.072297)

lat <- c(53.93057, 53.93222, 53.83308, 53.57506, 53.57339, 53.57612, 53.56305, 53.54610, 53.60224, 53.61145, 53.52641, 53.51799, 53.51425, 53.72778, 53.76036, 53.90055, 53.89335, 53.90751, 53.90092, 53.90038)

mud <- c(0.032, 0.039, 0.126, 0.146, 0.536, 0.225, 0.222, 0.408, 0.145, 0.112, 0.241, 0.031, 0.186, 0.340, 0.074, 0.162, 0.379, 0.147, 0.482, 0.220)

df <- data.frame(lon, lat, mud)

obs_rast <- rasterize(x = vect(df), y = blankRaster, field = "mud", fun = median)
obs_rast <- focal(obs_rast, w = 9, fun = median, na.policy="only", na.rm = T)
obs_rast <- resample(obs_rast, blankRaster)
out_df <- as.data.frame(obs_rast, xy = T)
names(out_df) <- c("lon", "lat", "mud")
coordinates(out_df) <- c("lon", "lat")

#IDW interpolation for each interp_df   
idw <- idw(formula = mud ~ 1, locations = out_df, newdata = grd)
idw_output <- as.data.frame(idw)
names(idw_output)[1:3] <- c("lon", "lat", "mud")
interp_output <- idw_output[, 1:3]

interp_rast <- terra::rasterize(x=vect(interp_output), y=blankRaster, field='mud', fun=median)
interp_rast <- resample(interp_rast, blankRaster)
plot(interp_rast)

The problem is that the final output raster object has regularly spaced horizontal and vertical white lines, but I need it to be continuous.

I've tried resample() on the final interp_rast object but no change. The only other thing I can think of is that the grd and out_df objects don't quite match up. If I plot them on top of each other and zoom in they seem ok.

plot(interp_rast, xlim = c(-6.4, -5.8), ylim = c(53.2, 54.1))
plot(grd, add = T, pch = 20)
plot(out_df, add = T, pch = 20, col = "red")

Answer 1

Here is a simplified workflow that does not produce these artifacts.

library(terra)
library(gstat)

df <- data.frame(
    lon = c(-6.054578, -6.030297, -6.148603, -6.001423, -6.029923, -6.041965, -6.026100, -6.013715, -6.027235, -6.026288, -6.004180, -6.020842, -5.997938, -6.030913, -6.025397, -6.107727, -6.089913, -6.054938, -6.042965, -6.072297),
    lat = c(53.93057, 53.93222, 53.83308, 53.57506, 53.57339, 53.57612, 53.56305, 53.54610, 53.60224, 53.61145, 53.52641, 53.51799, 53.51425, 53.72778, 53.76036, 53.90055, 53.89335, 53.90751, 53.90092, 53.90038),
    mud = c(0.032, 0.039, 0.126, 0.146, 0.536, 0.225, 0.222, 0.408, 0.145, 0.112, 0.241, 0.031, 0.186, 0.340, 0.074, 0.162, 0.379, 0.147, 0.482, 0.220)
)

extent <- ext(c(range(df$lon), range(df$lat))) + .1
blankRaster <- rast(extent, res=c(0.0117, 0.02075))

model <- gstat(id = "mud", formula = mud~1, locations = ~lon+lat, data=df)
out <- interpolate(blankRaster, model, xyNames=c("lon", "lat"))

Look at the output

plot(out[[1]])
points(df[, c("lon", "lat")])