Interpolation extent issue

[jfprieur]

Hello,

I am trying to understand the behaviour of the interpolation when I am creating surfaces such as a DTM. As mentionned in my email discussion, our point clouds are already classified so I filter out the ground returns then use the following to create the DTM raster.

chm = tile.chm(cell_size=0.25, interp_method = "linear")

We are processing 1600 buffered 1km^2 las tiles (40m buffer) to create DTM, DSM and CHM which then gets used by our ITC segmentation script. The rasters are at 25cm resolution so the output raster tiles should be 4160*4160. Our previous method using the lidR package in R produced this result.

When I use pyfor (and another python package called whitebox), the rasters have uneven dimensions such as 4079*4159 for example. I understand what is happening, the ground las points do not reach the corner of the declared extent of the las file and the interpolation stops a bit short. Here are some screenshots that explain the issue

LAS ground points (1000m * 1000m)

lidR output (4160*4160)

pyfor output (4079*4079)

whitebox output for comparison (4080*4080)

Your package does a better job as it uses the points at the top to interpolate whereas wbt does not

Arcmap TIN display

I believe what may be happening is the scipy interpolation is not getting the extent from the LAS file but calculating it from the actual points which makes it fall short of the 'official' extent. Do not know what R does different in this respect but this would be the desired output.

I will be digging into the code eventually to see if I can force this myself.

Thank you!

[brycefrank]

I see you are using the linear interpolation method. My experience with this is that it requires points on the "inside" of the cloud to produce a prediction. Using the "nearest" interpolator should fill the output raster completely.

I can't speak to what JR is doing, he may be using a nearest interpolator or some other method. To my knowledge he implements his own methods in C++, where I am more or less implementing a wrapper for scipy.interpolate as we have discussed.

One trick that I used to use (for better or worse, I will let you decide) is to create synthetic points on the outer rim of the tile to get the linear interpolator to work. Even then, it introduces uncertainty at the edge of the tiles.

One other option is to buffer a tile using its neighbor tiles and do the linear interpolation on the buffered tile. This would relegate the edge effect to only tiles at the edge of the study region, rather than the edges of all tiles.

Automated buffering tasks like this are a work in progress, with some preliminary code up (god willing) some time today, either on 0.3.4 or 0.3.4a (a staging branch for a very messy update I am working on)

[brycefrank]

I have also ran into some other issues in Grid that will be fixed today - for example the cell size of the output raster may not always equal the input cell_size argument.

[jfprieur]

Understood about the linear option. I tried with nearest and it mitigates the extent problem (4159*4159 raster produced) but the quality is not optimal. These tiles are already buffered for the segmentation process so we could always increase the buffer ;)

JRR is a colleague so I will reach out to him to find out. Will let you know if I find anything interesting.

[brycefrank]

Sounds good.

In any case I think I need to add some avenues for using already classified points. I will add it to do the to-do list.

[brycefrank]

JF,

I will close for now, and open a separate issue for adding the use of already-classified ground points for creating bare earth models.