problem with las2dem in parallel on Linux

[Qing XU]

Dear Martin,

Since the blast2dem-cli doesn't work on our Linux platform, we implement the CHM generation (proposed in one of your posts) using the las2dem-cli. The general rule is that when the tile has more than 5 million points, we retile it using lastile (the sub-tiles has less than 5 million points), and process each sub-tile sequentially. After the CHM for all sub-tiles are produced, we mosaic them together.  But tiles are processed in parallel on the Linux.

We use a lidar acquisition which has 20 tiles to test the procedure. When we serialize the tiles, it works fine, and give the correct results for each tile.

Here comes the problem when we parallel the tiles: for eight rasters (.bil) at different heights produced by las2dem, one or several of them (random too) have stripping, and those stripping is filled with abnormal high or low values. Sometimes, there is a huge black stripe that go through the whole tile, vertically or horizontally. Those striping or huge stripes are totally random, each time we run the 20 tiles in parallel, we've got stripping in different locations.

One week before, we just cleared the problem brought by GDAL, which we use for mosaicing those sub-tiles into an intact tile. Since the mosaicing are in parallel for the 20 tiles, black stripes are produced because gdal doesn't support writing in parallel. Please see the discussion below:

http://comments.gmane.org/gmane.comp.gis.gdal.devel/33736

We simply serialized the gdal. But unfortunately, we found that there are still random stripping in the CHM rasters, which led us to find out that the las2dem doesn't give correct outputs mentioned above.

Do you have any idea what is going on, and what could be the feasible solution? Does las2dem have the same "global block cashe" problem the same as gdal?

Expecting your help soon! Anyone who has any idea about it, please write back!

Best regards,

Qing

[Martin Isenburg]

Hello Qing,

I have a suspicion. The tile you are creating are just above the threshold where LAStools switches into "on-disk" buffering.

D:\LAStools\bin>more 41121D8410_tile_593000_4585500_chm_00.hdr

# created by LAStools (c) 2014 martin....@rapidlasso.com

nrows 2322

ncols 2322

nbands 1

nbits 32

layout bil

nodata -9999

byteorder I

ulxmap 592959.875000

ulymap 4586040.125000

xdim 0.25

ydim 0.25

pixeltype float

D:\LAStools\bin>more 41121D8410_tile_593000_4585500_chm_ground.hdr

# created by LAStools (c) 2014 martin....@rapidlasso.com

nrows 2320

ncols 2320

nbands 1

nbits 32

layout bil

nodata -9999

byteorder I

ulxmap 592960.125000

ulymap 4586039.875000

xdim 0.25

ydim 0.25

pixeltype float

Here is the corresponding piece in my code where the switch is made:

      // for large rasters use a disk buffer

      if (nrows * ncols > 2048*2048)

      {

        srwriter->set_buffer(new SRbufferRowBands()); // GIS 2006 paper

      }

      else

      {

        srwriter->set_buffer(new SRbufferInMemory()); // for small images

      }

When on-disk buffering is used a temporary file name is needed. All my temporry filenames include the process ID that is formed like this:

sprintf(base, "temp.%d", GetCurrentProcessId());

I assume that under Linux the Window's construct GetCurrentProcessId() does not work properly and always 0 (or so) is returned. Hence the temp files of the different processes over-write another. As a quick fix you can try the following: Add the command '-buffer inmemory' to your command string. Then no on-disk buffering is used.

By the way. I noticed that your input seems to be a tile but that you do not use any buffers around the tile edges (different kind of buffering) as described here.

http://rapidlasso.com/2015/08/07/use-buffers-when-processing-lidar-in-tiles/

I first noticed that when I saw that your BIL rasters are slightly different in size.

Finally ... here is the pit-free CHM LAStools command sequence that Qing is using:

las2dem.exe  -i input.laz -o chm_ground.bil -drop_z_above 0.1 -step 0.25 -v

lasthin.exe  -i input.laz -o temp_25.laz -subcircle 0.1 -step 0.125 -highest -olaz -v

las2dem.exe  -i temp_25.laz -o chm_00.bil -step 0.25 -kill 1.0 -v

las2dem.exe  -i temp_25.laz -o chm_02.bil -drop_z_below 2 -step 0.25 -kill 1.0 -v

las2dem.exe  -i temp_25.laz -o chm_05.bil -drop_z_below 5 -step 0.25 -kill 1.0 -v

las2dem.exe  -i temp_25.laz -o chm_10.bil -drop_z_below 10 -step 0.25 -kill 1.0 -v

las2dem.exe  -i temp_25.laz -o chm_15.bil -drop_z_below 15 -step 0.25 -kill 1.0 -v

las2dem.exe  -i temp_25.laz -o chm_20.bil -drop_z_below 20 -step 0.25 -kill 1.0 -v

las2dem.exe  -i temp_25.laz -o chm_25.bil -drop_z_below 25 -step 0.25 -kill 1.0 -v

I am curious. Do you like the result that the pit-free CHM algorithm currently produces? What were you using before to create CHMs? Does the absence of pits improve your subsequent analysis? I am asking in part because we have something new for all pit-free CHM fans this Silvilaser 2015 that happens at the end of September in France. I hope you are coming ... (-:

Martin

--
Download LAStools at
http://lastools.org
http://rapidlasso.com
Be social with LAStools at
http://facebook.com/LAStools
http://twitter.com/LAStools
http://linkedin.com/groups/LAStools-4408378
Manage your settings at
http://groups.google.com/group/lastools/subscribe

[Qing XU]

Dear Martin,

Suppressing the on-disk buffering works!!!

This is a product of lastile already. Reason for the different size of the two BIL rasters is that they are both wrong outputs of las2dem.

My feeling toward the pit-free CHM algorithm is quite ambiguous. The reason is "-kill 1.0". Do correct me if I am wrong. With "-kill 1", we are able to minimize the false canopy in the area without any lidar returns. But it also produces holes (NA) (quite big holes) in the canopy when the lidar returns are far away from each other. Point density varies across even a single tile, while the canopies looks really good in some parts of the tile, canopies in other parts just have holes. If the hole is in the middle of the canopy, it's not a difficult task to smooth it. But if the hole is in the boundary of the canopy (imaging the apple's logo), how to make it look like an intact tree crown? I tried to increase the triangle side length, and use "-kill 2", there won't be holes in the canopy any more, but the false canopy between true canopies is a real issue. It looks like a bridge between two true canopies, with similar height value as the boundary of the two canopies. Or is it related to the pixel size that I used, 0.25 m? My lidar data point density vary from 3 to 51 points/m2, the majority has about 10 points/m2. Do you have a better solution to this problem? I am testing a new strategy, in which I use "-kill 2" for the points higher than 10 meters. Hope it work out!

Previously, I apply the Gaussian height-based filtering (Pitkanen et al, 2004) on the raw CHM generated from the maximum height of the lidar returns in a pixel area. fill the holes is the first step to deal with before smoothing. Unfortunately, I can't go to Silvilaser 2015. Expecting your new things there!

Best regards,

Qing