2

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.

I have the impression this should be a pretty straightforward issue, yet I could not solve it so far despite considerable effort. Any help would be appreciated.

3

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.

6
  • Thanks a lot for your help! I made one small adjustment: after xyV = as.data.frame(r,xy=TRUE), I added: colnames(xyV) <- c("x", "y", "layer"). With this, your code runs for me. It seems to be quite slow on large rasters though, my testing suggests because of the automap::autoKrige function (i.e. the Kriging). – Abdirizak Apr 25 '20 at 12:13
  • If you have a large raster and your spatial correlation is a lot bigger than the raster cell size then you can use a sample of the known points and you shouldn't get a much worse answer. Random thinning or maybe keeping more points near your unknown points...Hmmm – Spacedman Apr 25 '20 at 12:49
  • Ok. But what do you mean by the correlation being larger than the raster cell size? Correlation ranges from -1 to +1, whereas the cell size can be anything in meters. I don't think the two are readily comparable. Or is there some formula that sets them into relation? – Abdirizak Apr 28 '20 at 9:51
  • I mean the spatial range of the correlation. If you have a 1000x1000 1meter grid and your spatial correlation range is about 100m then the amount of information got from two adjacent grid cells is going to be not much different from one cell. Hence if computing with the 1000000 cells is a problem you can sample a subset of locations and not do much worse with your predictions. – Spacedman Apr 28 '20 at 10:32
  • You can estimate the spatial range from the variogram, and you can produce that from subsets of the locations if you have too many grid points to do the whole thing. I used autokrige in my answer because its an all-in-one solution and scaling to large grids wasn't mentioned as an issue. – Spacedman Apr 28 '20 at 10:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.