Getting Edge Artifacts in Processing DEM in HTC system

[pma...@g.clemson.edu]

Hi, 

I am having issues in processing LAZ files to create a DEM. I am getting edge artifacts and white spots for the area where there is no data. We have a high throughput computing system (HTCondor) here which means I can split my big dataset into pieces and send it to individual core of each computer to process it. We have 34 computers here with 8 cores each and the way my process works is through the following steps:

1. To split this big data set, I zip each laz file into a separate folder. The reason I do this is to split my job in a way that each core of each computer gets one laz file. 

2. Then on each zipped lazfile, I unzip it first then do the following:

- Las2Las (Proj) to project my laz file

- lastile Pro to tile my laz file into smaller pieces of 1000 m size and 25 m buffer ( I have tried smaller buffers too it does the same)

- las ground with ultrafine granularity and nature or forest type.

- lasheight with ground(2) class

- lasclassify everything default including gutters and no tiny buildings

- las2dem pro - using tile bb with ground points only

then I mosaic all the pieces on each core using ArcGIS Mosaic to new raster as tif, and once I have all the .tif files back to the computer where I submitted my job from. I do Mosaic to New Raster again to mosaic all the final pieces. 

Can anyone please help me in solving the edge problem, I attached my final results and my scripts. I have a big deadline next week. I am not sure if the problem is with splitting or if the tiles need to see their neighborhoods. 

Any help is much appreciated. 

Best, 

Palak

 

[Martin Isenburg]

Hello Palak,

I am having issues in processing LAZ files to create a DEM. I am getting edge artifacts and white spots for the area where there is no data.

The PNG file you sent does not show (obvious) edge artifacts but I can see many "empty" areas. It may be that there are sometime large gaps between the ground points (e.g. lack in coverage, water bodies, clouds, ...). If there is a gap larger than 100 meters then you will see empty pixels in the generated raster as the default '-kill' is 100. Set the '-kill 1000' and those areas will be closed ... but the interpolation across those areas will obviously be a poor approximation of the terrain.

We have a high throughput computing system (HTCondor) here which means I can split my big dataset into pieces and send it to individual core of each computer to process it.

Along those "splits" of the data (into 34  pieces? one per computer?) you will get edge artifacts if you do not use buffers. An easy solution is to make "overlapping splits" and then remove the overlapping parts.

 

We have 34 computers here with 8 cores each and the way my process works is through the following steps:

1. To split this big data set, I zip each laz file into a separate folder. The reason I do this is to split my job in a way that each core of each computer gets one laz file. 

You got to give me more detail here.  What is the input? Can you send a lasinfo report of two three typical input files?

Is it LiDAR flightlines in LAZ format? Your 'process_laz_files.py' scripts suggests that as you use the '-files_are_flightlines' option in the lastile step.

What (original) projection is the data in? Your 'process_laz_files.py' scripts suggests that the tiles must have a correctly populated source projection (stateplane?) as you use the '-target_utm 17N' option in the las2las step without specifying the source projection.

 

2. Then on each zipped lazfile, I unzip it first then do the following:

- Las2Las (Proj) to project my laz file

- lastile Pro to tile my laz file into smaller pieces of 1000 m size and 25 m buffer ( I have tried smaller buffers too it does the same)

- las ground with ultrafine granularity and nature or forest type.

- lasheight with ground(2) class

- lasclassify everything default including gutters and no tiny buildings

- las2dem pro - using tile bb with ground points only

All that looks good except that your script uses the term DEM for the ground model and DTM for the raster generated from the first returns. The latter is usually called a Digital *Surface* Model or DSM. I suggest you drop the term DEM and use DTM instead and rename what you currently call the DTM as the DSM.