SAGA-based bicubic spline interpolation

[mosscoder]

Greetings, I've been working on a pipeline to make gridded terrain models from classified point cloud based upon the catalog engine and SAGA. I was interested in a bicubic spline interpolation method found here.

Here is the function (and some others) in its present state:

bicubicSplineInterp <- function(cluster, res, buffer) {
  las <- readLAS(cluster)
  if (is.empty(las))
    return(NULL)

  bbox  <- raster::extent(cluster)
  bufferedExtent <- raster::extent(las)

  dir <- paste0(tempdir(), '/', bbox[1], '_', bbox[3])
  dir.create(dir)

  xyz <-lasfilter(las,
                  Classification == 2L)@data[ , c('X', 'Y', 'Z')] 
  xyz <- st_as_sf(as.data.frame(xyz), 
                  coords = c('X', 'Y'), 
                  crs = crs(las))
  point_dir <- paste0(dir, "/irr_pts.shp")
  st_write(xyz, point_dir, driver = "ESRI Shapefile", overwrite = TRUE)

  ras_dir <- paste0(dir, '/interp.sgrd')

  system(
    paste0(
      'saga_cmd grid_spline 4', # http://www.saga-gis.org/saga_tool_doc/2.1.3/grid_spline_4.html
      ' -SHAPES ', point_dir, # Location of point cloud in vector form
      ' -FIELD Z', # Layer containing Z
      ' -METHOD 1', # Default 'refinement', seems faster than alternative
      ' -LEVEL_MAX 14', # Highest level ~looks best, doesn't have huge performance disadvantage
      ' -TARGET_USER_XMIN ', bufferedExtent[1],
      ' -TARGET_USER_XMAX ', bufferedExtent[2],
      ' -TARGET_USER_YMIN ', bufferedExtent[3],
      ' -TARGET_USER_YMAX ', bufferedExtent[4],
      ' -TARGET_USER_SIZE ', res, 
      ' -TARGET_OUT_GRID ', ras_dir
    )
  )

  interp <- raster(paste0(dir, '/interp.sdat'))
  interp <- raster::crop(interp, bbox)
  unlink(dir, recursive=TRUE)

  return(interp)

}

Here is a reproducible example (assuming you have SAGA running) of its use in a catalog pipeline:

library(lidR)
source('https://raw.githubusercontent.com/mosscoder/lidR_supplements/master/las_helper_functions.R')

setwd('~/Downloads')

set_lidr_threads(1L)
data.table::setDTthreads(1L)
plan('multisession', workers = 4L)

t <- tempfile()

download.file(url = 'https://storage.googleapis.com/mpg-files/classified_example.las.zip',
              destfile = t)

ctg <- catalog(unzip(t)[1])

buffer <- 30
res <- 0.2

opt_chunk_buffer(ctg) <- buffer
opt_chunk_size(ctg) <- 200
opt_stop_early(ctg) <- TRUE
opt_select(ctg) <- 'xyzc'
opt_output_files(ctg) <- './bicubic/{XLEFT}_{YBOTTOM}_{ID}'

opt <- list(raster_alignment = res,
            automerge = TRUE)

dtm <- catalog_apply(ctg, 
              FUN = bicubicSplineInterp, 
              buffer = buffer,
              res = res, 
              .options = opt)

hillPlot(readLAS(ctg@data$filename), main = 'TIN interpolation') 
hillPlot(dtm, main = 'Bicubic spline interpolation') 

It largely performs as I'd hope, and goes quickly enough for my needs. I do suffer from tile edge effects, however. Evident in the second plot that the above code chunk produces. Particularly, evident in the interpolated regions. I'm not sure if this is because of my misuse of the tiling functionality. The function does crop by the extent of the cluster (bbox) at the end. Perhaps it is an alignment issue?

[Jean-Romain]

That is definitively a question and this should be asked on https://gis.stackexchange.com/ with the tag lidr. Quickly (if you need more details please ask on the forum):

The function does crop by the extent of the cluster (bbox) at the end.

You did is yourself

interp <- raster::crop(interp, bbox)

[Jean-Romain]

Sorry, I think I did not answer the question actually. But what is the question? And what is the problem? This look like a pretty advanced but correct usage of catalog_apply(). I don't have SAGA so I can't reproduce.

[mosscoder]

First, I thought I might prompt a feature request for bicubic spline interpolation methods and offer my (currently flawed) solution. It is perhaps doubly flawed because it adds dependency outside of R. I have been looking for a smooth in appearance interpolation algorithm that is also fast. Exploring options within R in the Akima and Fields packages reveal that those can offer bicubic interpolation methods, but can't handle data of point cloud size. Kriging also takes quite some time at the resolutions I'm interested in. Thus, I went down the external dependency path of SAGA GIS, which is also freely available.

The question I had relates to edge effects, and I can ask it on stack exchange. But cropping by the bbox extent is not eliminating my edge artifacts, and no amount of additional buffering mitigates them either. I suspect it is an alignment issue, and am digging deeper.

[Jean-Romain]

There is no alignment issue according to your code.You properly used the option for ensuring correct alignment. Please show the output to see what do you mean by edge artifacts.

[mosscoder]

Thanks for taking the time.

[Jean-Romain]

To me it looks good. There are indeed edge artifacts but only in regions with a very low density of ground points. The bicubic spline interpolation seems weak at these locations. My advices are

• Bigger chunks. 200 m is not relevant. Why you don't use 1 km?
• Bigger buffer. It seems you have regions with very low densities of ground points. Use a buffer of 100 m. Maybe 200.

Also optimize you code to compute faster and save memory

• Use opt_filter(ctg) <- "-keep_class 2" and remove lasfilter(las, Classification == 2L)@data[ , c('X', 'Y', 'Z')]
• Use as.spatial instead of st_as_sf(...). It is optimized to recycle the memory
• Index your point cloud

[mosscoder]

The buffer suggestion worked for me thanks. I tried 100m buffers earlier, and thought that was crazy high, but it looks like the 300m is needed. I do like how it looks, even if interpolated understory is just fantasy anyway. Thanks for the memory optimization tips, I'll need them.

[Jean-Romain]

This is indeed crazy especially if you still use a chunk size of 200m. But it seems that your methods is weak when there are to few points. Look at the crater on the top. I'm sure it is an artifact of the algorithm. Maybe you should also chnage the parameters of the SAGA algorithm.

Next time you have a question please ask on gis.stackexchange. Usually I do not answer to questions on github. But you triggered my curiosity.

[mosscoder]

Regarding the crater, this is a structure from motion point cloud, not lidar in origin. So, canopy penetration is pretty bad.

To add some additional closure and for future experimentation I found an R package that implements this interpolation method called MBA that eliminates the SAGA dependency and revised the interpolation function as follows:

bicubicSplineInterp <- function(cluster, res, h = 12, zDecimals = 2) {
  las <- readLAS(cluster)
  if (is.empty(las))
    return(NULL)

  bbox  <- raster::extent(cluster)
  bufferedExtent <- raster::extent(las)

  template <- raster(bufferedExtent, crs = crs(las), res = res)

  interp <- mba.surf(xyz = las@data[,c('X', 'Y', 'Z')], 
           no.X = nrow(template),
           no.Y = ncol(template),
           h = h,
           extend = TRUE
           )$xyz.est$z

  interp <- raster(apply(interp, 1, rev))

  extent(interp) <- extent(template)
  crs(interp) <- crs(template)

  interp <- raster::crop(interp, bbox)

  return(interp)

}

Note there is a new tuning parameter h, which controls the complexity of the interpolation method. I see diminishing returns beyond 11, and memory demands increase as well. More information here.

The following can now be used to reproduce my results without SAGA:

library(lidR)
library(MBA)
library(future)

source('https://raw.githubusercontent.com/mosscoder/lidR_supplements/master/las_helper_functions.R')

set_lidr_threads(1L)
plan('multisession', workers = 4L)

t <- tempfile()

fUrl <- 'https://github.com/mosscoder/lidR_supplements/blob/master/classified_example.las.zip?raw=true'

download.file(url = fUrl,
              destfile = t)

ctg <- catalog(unzip(t)[1])

res <- 0.2

opt_chunk_buffer(ctg) <- 300 #absurd buffer needed to avoid edge artifcats
opt_chunk_size(ctg) <- 300
opt_stop_early(ctg) <- TRUE
opt_select(ctg) <- 'xyz'
opt_filter(ctg) <- '-keep_class 2'

opt <- list(raster_alignment = res,
            automerge = TRUE)

dtm <- catalog_apply(
  ctg,
  FUN = bicubicSplineInterp,
  res = res,
  .options = opt
)

hillPlot(readLAS(ctg@data$filename), main = 'TIN interpolation') 
hillPlot(dtm, main = 'Bicubic spline interpolation') 

[Jean-Romain]

Let me show you something undocumented

mba = function(n = 1, m = 1, h = 8, extend = TRUE)
{
  n <- lazyeval::uq(n)
  m <- lazyeval::uq(m)
  h <- lazyeval::uq(h)
  extend <- lazyeval::uq(extend)

  f = function(what, where, scales = c(0,0), offsets = c(0,0))  {
    return(MBA::mba.points(what, where, n, m , h, extend)$xyz.est[,3])
  }

  class(f) <- c("SpatialInterpolation", "Algorithm", "lidR", "function")
  return(f)
}

LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las = readLAS(LASfile)
dtm = grid_terrain(las, algorithm = mba())

You no longer need to think about catalog_apply() and its options. You no longer need to take care of cropping the raster or raster alignment. You no longer need to take care of memory optimization. grid_terrain() does all of that for you. Just read the documentation of grid_terrain().

Enjoy

---

ps: the following will work as well

las = lasnormalize(las, mba())
las = las - mba()
# and so on

[mosscoder]

Wow, yes thanks for taking the time to share that, it is much faster than what I'd written.

Generally, I'm curious about your opinions of this method as compared to other interpolation techniques. The results from this paper piqued my curiosity in spline-based approaches.

Looking at figures 3 and 5, there appears to be a marginal benefit to splines at higher resolutions in steep slopes and low point densities, though they are more outlier prone according to the authors. Maybe a followup outlier screening is a good idea post-gridding?

[Jean-Romain]

No opinion. I'm providing tools but I don't tell people what they should or should not do. I'm only telling people how to do what they want to do.

[jp136769]

Hello @Jean-Romain ,

With lidR 4.0.2, using rasterize_terrain, I get the following error :

rasterize_terrain(clipped, res=0.5, algorithm=mba())
Error: Invalid function provided as algorithm.

Any clues to make it works ? Thanks,