lastile
Gordon Frazer
Hi Martin,
I have a lidar dataset that is composed of 4 billion classified (ground vs. non-ground) points and currently partitioned into 105 separate 1000 x 1000 m tiles. I would like to try using LAStools (i.e., lastile vs. lasclip vs. las2las) to retile the entire dataset using a 500 x 500 m tile size and user-defined grid origin. Ideally, I would like the tiling procedure to follow a consistent boundary logic where points are included in a tile: IF (pt_x >= x_min & pt_x < x_max & pt_y >= y_min & pt_y < y_max), so that no lidar points will be duplicated or excluded from the dataset. The tiles should also be temporally buffered by some small margin (10 to 25 m) to avoid edge artifacts when normalizing laser elevations to heights above the ground surface using lasheight. Finally, I’ll want to remove the buffer from the tiles containing laser heights so that I can compute a broad range of gridded vegetation metrics using existing software. Could you please let me know the most efficient way to do this with LAStools.
Best regards and many thanks,
Gordon
Martin Isenburg
excellent question. lastile uses exactly the tiling logic you describe
below as I like the tiling you describe. i have come across many
tilings (generated off shapefile descriptions?) that duplicate points
that fall exactly onto a tile edge (or quadruple points falling
exactly onto a corner). I dislike such tilings a lot as this
"duplication" of points is completely arbitrary and unnecessary as it
happens completely "by chance".
You can implement the workflow you describe as follows:
(1) run lastile to retile your existing tiles. afterward each point
belongs to exactly one tile and only buffer points are appear
duplicated that can easily be removed later.
lastile -i old_tiles*.las -merged -o tiles -olaz -tile_size 500 -
buffer 25
(1a) optional: if your original tiles had duplicated points along the
boundaries you may remove them now with
mkdir duplicate
lasduplicate -i tiles_*.laz -unique_xyz -odir duplicate -olaz
(2) run lasheight on all the tiles and putting the output into a the
height directory.
mkdir height
lasheight -i tiles_*.laz -replace_z -odir height -olaz
if you have a multi-core machine then this is an excellent time to put
all of them to use by adding, for example, the '-cores 3' or '-cores
7' option:
lasheight -i tiles_*.laz -replace_z -odir height -olaz -cores 3
(3) remove the buffer
cd height
mkdir no_buffer
lastile -i tiles_*.laz -remove_buffer -odir no_buffer -olaz
If you use lasgrid or las2dem to process the tiles you do not need to
remove the buffer. You can simply instruct the tools to respect the
tile bounding box with the '-use_tile_bb'. For example, you can create
a false color tree canopy rasters:
cd height
lasgrid -i tiles_*.laz -use_tile_bb -highest -false -opng -set_min_max
0 40 -step 2
Or raster a hill-shaded model of the normalized heights (DSM - DTM)
directly from the points
cd height
las2dem -i tiles_*.laz -clip_z_below 0 -use_tile_bb -hillshade -opng -
step 1
We are currently running similar work-flows to create canopy height
and biomass models for the forested areas on a Canary island. However,
I did not try this exact work-flow as described, so please report back
to me if something does not work correctly.
Cheers,
Martin @lastools
PS: Please note that for any sort of production use, be it commercial
or government, you would need to license LAStools for this and many
similar kind of work-flows.
Gordon Frazer
Hi Martin,
Thank you for the helpful response. The one detail that did not appear was how to invoke a user-defined origin. It is often the case that gridded lidar data will be spatially integrated with other types of remotely sensed data (digital colour photography, satellite imagery, radar, etc.). Therefore, it is helpful that for a lidar analyst to generate a grid that shares the same origin point and cell size as the other spatial datasets (i.e., grid cells from disparate dataset will overlay perfectly upon one another). Is it possible to input a specific grid origin in lastile or is it arbitrarily defined by the bounding box (range rectangle) the entire dataset?
I’ll happily report any success or failure with the workflow you kindly supplied. I would also license LAStools for any type of production work. This effort is purely to see how easily I can accomplish the task with LAStools. Your efforts are greatly appreciated.
Thanks,
Gordon
Gordon Frazer
Hi Martin,
A user-defined grid origin (ll_x, ll_y) would simply provide a predetermined starting point to build a tiling scheme that would cover the (x,y)-dimensions of the entire lidar dataset. The benefit here (assuming that lidar points will be used to generate gridded surfaces; e.g., vegetation metrics) is that all lidar-derived grid cells (pixels), once mosaicked, will sit directly upon other existing GIS raster layers or remotely sensed imagery without having to resample any of the datasets. Starting the tiling scheme on the lower left-hand corner of the point cloud bounding box does not allow the flexibility required to spatially integrate lidar-derived gridded data with other raster datasets unless grid-cell resampling takes place. I wrote a R program last evening that builds a tiling scheme with a user-defined origin (see below). The lidar point cloud would still fall entirely within the bounds of this tiling scheme; however, the resultant grid cells (pixels) would now assume the correct position according to the other datasets with which they will be spatially integrated with. Please let me know if this still does not make sense to you.
Thanks,
Gordon
# This R program builds a lidar tiling scheme in shape file format.
# load libraries.
library(shapefiles)
# set input parameters.
# range rectangle of tiling grid.
x_min = 1063323.805
y_min = 490409.75
x_max = 1086323.805
y_max = 508909.75
# tile size
tile_size = 500
buffer = 10
# number of columns and rows in tiling scheme.
num_col = ceiling(x_max - x_min)/tile_size
num_row = ceiling(y_max - y_min)/tile_size
# create data table to store coordinates of polygon vertices.
dd = data.frame(matrix(0, nrow = (5 * num_col * num_row), ncol= 3))
names(dd) = c("Id", "X", "Y")
# create data table to store table attributes.
ddTable = data.frame(matrix(0, nrow = (num_col * num_row), ncol = 11))
names(ddTable) = c("Id", "Block_id","buffer", "ll_x", "ll_y", "ul_x", "ul_y", "ur_x", "ur_y", "lr_x", "lr_y")
# loop through grid and assign xy coordinates for vertices, block_id, and tile_id.
tile_num = 0
row_id = 0
for (i in 1:num_col){
for (j in 1:num_row){
# compute table values.
# sequential tile number.
tile_num = tile_num + 1
# concatentate tile id.
tile_id = paste("chm_", i - 1, "_", j - 1, sep="")
# compute vertices of key points.
# lower left corner.
ll_x = x_min + ((i - 1) * tile_size) - buffer
ll_y = y_min + ((j - 1) * tile_size) - buffer
# upper left corner.
ul_x = ll_x
ul_y = ll_y + tile_size + 2 * buffer
# upper right corner.
ur_x = ul_x + tile_size + 2 * buffer
ur_y = ul_y
# lower right corner.
lr_x = ur_x
lr_y = ll_y
# load data in matrices for points.
for (k in 1:5){
row_id = row_id + 1
if (k == 1){
dd$Id[row_id] = tile_num
dd$X[row_id] = ll_x
dd$Y[row_id] = ll_y
} else if (k == 2){
dd$Id[row_id] = tile_num
dd$X[row_id] = ul_x
dd$Y[row_id] = ul_y
} else if (k == 3){
dd$Id[row_id] = tile_num
dd$X[row_id] = ur_x
dd$Y[row_id] = ur_y
} else if (k == 4){
dd$Id[row_id] = tile_num
dd$X[row_id] = lr_x
dd$Y[row_id] = lr_y
} else if (k == 5){
dd$Id[row_id] = tile_num
dd$X[row_id] = ll_x
dd$Y[row_id] = ll_y
}
}
# load data matrices for attribute table.
ddTable$Id[tile_num] = tile_num
ddTable$Block_id[tile_num] = tile_id
ddTable$buffer[tile_num] = buffer
# range rectangle for each tile
ddTable$ll_x[tile_num] = ll_x
ddTable$ll_y[tile_num] = ll_y
ddTable$ul_x[tile_num] = ul_x
ddTable$ul_y[tile_num] = ul_y
ddTable$ur_x[tile_num] = ur_x
ddTable$ur_y[tile_num] = ur_y
ddTable$lr_x[tile_num] = lr_x
ddTable$lr_y[tile_num] = lr_y
}
}
# create shape file.
ddShapefile = convert.to.shapefile(dd, ddTable, "Id", 5)
write.shapefile(ddShapefile, "e:\\temp\\bb_analysis_grid_10m_buffer", arcgis=F)
Terje Mathisen
processing pipeline, starting with lastile to split the input into
250x250 m tiles, with 35 m boundaries.
After ground point have been determined, the next step is lasheight, and
here I've found a small glitch:
If lasheight gets to a tile where all the points are classified as
ground, you get a warning message
all points have classification 2. skipping
'tiles_raw\tile_455250_6602750.laz'
and (as it says), that tile is skipped.
The problem is of course that for all subsequent processing, those
skipped tiles will turn up as holes with zero data. :-(
My current workaround is to add three lines that CD into the source
directory, then copies any tile file that doesn't exist in the
tiles_classified directory, and then CDs back up:
rem run from the source dir:
cd tiles_ground
rem The next is all on one line!
for %%f in (tile_*.laz) do if not exist ..\tiles_classified\%%f copy %%f
..\tiles_classified\%%f
rem get back to parent dir
cd ..
If the number of ground points is too low, you also get an error message.
I suspect though that most users of lastools will expect such a pipeline
to pass through all the tiles with any data in them!
Terje
--
- <Terje.M...@tmsw.no>
"almost all programming can be viewed as an exercise in caching"
Martin Isenburg
yes. If you look in the source code (lasdefinitions.hpp) you find that for LAStools points that fall exactly onto the lower and the left edge of a tile or a rectangle are considered part of the area of interest whereas points that fall exactly onto the upper or right edge are not.
See the code samples below for more details on the filter logic.
Regards,
Martin @rapidlasso
---
BOOL inside_rectangle(const F64 r_min_x, const F64 r_min_y, const F64 r_max_x, const F64 r_max_y) const
{
F64 xy;
xy = get_x();
if (xy < r_min_x || xy >= r_max_x) return FALSE;
xy = get_y();
if (xy < r_min_y || xy >= r_max_y) return FALSE;
return TRUE;
}
BOOL inside_tile(const F32 ll_x, const F32 ll_y, const F32 ur_x, const F32 ur_y) const
{
F64 xy;
xy = get_x();
if (xy < ll_x || xy >= ur_x) return FALSE;
xy = get_y();
if (xy < ll_y || xy >= ur_y) return FALSE;
return TRUE;
}
BOOL inside_box(const F64 min_x, const F64 min_y, const F64 min_z, const F64 max_x, const F64 max_y, const F64 max_z) const
{
F64 xyz;
xyz = get_x();
if (xyz < min_x || xyz >= max_x) return FALSE;
xyz = get_y();
if (xyz < min_y || xyz >= max_y) return FALSE;
xyz = get_z();
if (xyz < min_z || xyz >= max_z) return FALSE;
return TRUE;
}
Hi Martin,Does this same logic apply for the "-inside_rectangle" option for las2las and las2txt?ie: IF (pt_x >= x_min & pt_x < x_max & pt_y >= y_min & pt_y < y_max)RegardsAndrew
http://rapidlasso.comhttp://facebook.com/LAStools
Be social with LAStools at
http://twitter.com/LAStools
http://linkedin.com/groups/LAStools-4408378
Manage your settings at
http://groups.google.com/group/lastools/subscribe
Ty Kennedy-Bowdoin
--
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
Ty Kennedy-Bowdoin
Carnegie Institution for Science
260 Panama St.
Stanford, CA 94305
USA
+1 831 227-3885
Skype: Ty.bowdoin
Reginald Argamosa
Martin Isenburg
--
Download LAStools at
http://lastools.org
Terje Mathisen
> Hello Reginald, > > as a rule-of-thumb I use a 25 meter buffer for a 1000km by 1000km
> tile. If you have only forest and few waterbodies (=> areas without > any returns) then a smaller buffer may do such as 10 or 15 meters.
If > you have large buildings (e.g. IKEA's) or lots of waterbodies in
your > survey I would go up to 50 meter buffers for a 1000km by 1000km
tile.
:-)
> I have not personally investigated what buffer size guarantees > edge-artifact free DTMs / CHMs on different types of terrains. Might
> be an interesting study for a student to publish a paper on. The >
buffer should generally a lot smaller than the tile size. You cannot >
set it bigger as lastile won't let you. Bigger buffers mean more data >
needs to be processed and that will slow you down on larger areas so >
don't be too generous with your buffer sizes ... (-:
(320x320 including the buffers), which is about 60% more than the
internal size of 62500 sq m.
Intuitively the buffers need to be wide enough that we'll always find
ground points inside the buffer zone, so that all border TINs have an
anchor location, this means that it needs to be a bit larger than the
worst case ground point distances.
With even wider (relatively) buffers, the overhead becomes pretty huge.
