Interpolating DEM using Kriging in R [closed]

Question

I have several DEMs as geotiff files which I load into R using raster::raster. Now some of these DEMs contain patches of missing values, represented by NA. I would like to interpolate these missing values from the surrounding cells using Kriging.

One approach I tried was the raster::interpolate function. However, the examples in the documentation use spatial points data frames as predictors (from the "meuse" data set) and I can't figure out how it works for my raster.

I have also seen this, but the approach suffers from the same problem (interpolation of a raster from points). A similar example is this.

Several other questions on https://gis.stackexchange.com address similar issues, however none of them provided a working solution for me.

What I also tried is the R package filling. After converting my raster to a matrix, several functions therein work for me, yet I am not sure if these methods are appropriate for the problem at hand.

Answer 1

Try this. It works by converting the raster out to a data frame of x,y,value columns, making a spatial points data frame, extracting the missing and non-missing points, fitting on the non-missing points with prediction on the missing points, and then filling in the missing points with the predictions at the missing points in a new raster:

library(raster)
library(sp)
library(automap) # for autoKrige - or adapt to use gstat

fillRaster <- function(r){

    xyV = as.data.frame(r,xy=TRUE)
    sp::coordinates(xyV)=~x+y
    miss = is.na(xyV$layer)

    m = automap::autoKrige(
                     layer~1,
                     input_data = xyV[!miss,],
                     new_data=xyV[miss,])

    rfill = raster(r)
    rfill[] = r[]
    rfill[miss] = m$krige_output$var1.pred

    return(rfill)
}

To test, a data generator function:

testRaster  <- function(seed,n=40,m=50,pmiss=0.02){
    set.seed(seed)
    rmiss = raster(matrix(runif(n*m),n,m))
    rfull = rmiss
    rmiss[runif(n*m)<pmiss]=NA
    list(full = rfull,miss= rmiss)
}

Then:

Make a raster with 10% missing values and fill them via kriging:

> rr = testRaster(123,pmiss=0.1)
> rrfilled = fillRaster(rr$miss)
[using ordinary kriging]
> plot(rr$full-rrfilled)

that should show zero where you had data, and your prediction error where there was missing data. r$full is the raster before the data was made missing.

For the test rasters generated by that function the plot looks like noise because it is all noise, there's no spatial structure. You're lucky if the kriging doesn't complain with that sort of data. Feed it something sensible and test.