Parse lastiling line to bbox coordinates

[Marcel]

Hi, I have a tiled dataset, where:

1. I generate TIFs of the tiles + buffers,

2. process the TIFs outside of lastools

3. finally, combine the TIFs keeping only the original bounding box, so something like "-use_tile_bb".

I want to know:

a) the bounding box of tile+buffer and 

b) the bounding box of the tile (no buffer).

for each tile.

I believe the answer lies in this line from the output of lasinfo:

"LAStiling (idx A, lvl B, sub C, bbox D E F G, buffer) (size H x I, buffer J)"

From what I can tell, the line has the same structure for every tile and I replaced the actual numbers with the letters A..J.

What formula f(A..J) can I apply to get the info that I need?

I think this formula must have something to do with a quadtree, where "idx A, lvl B" is how we get the delta in x and y from the lower left corner, then just multiply these by tile size ("H x I") and add the corner coordinates  ("bbox D E F G"). Pleas help me put it into an equation,

Thank you!

Cheers,

Marcel

[Jochen Rapidlasso]

Hi Marcel,
let's take a small file and do some samples:

We take zurich.laz out of the LAStools demo folder:  

lasinfo -i zurich.laz -nc
  number of point records:    656837
  scale factor x y z:         0.01 0.01 0.01
  offset x y z:               0 0 0
  min x y z:                  676750.00 246000.00 547.29
  max x y z:                  676849.99 246099.99 573.87

first step always is indexing:

lasindex64 -i c:\lastools\data\zurich.laz -o tmp.laz

lastile64 -i tmp.laz -tile_size 80 -buffer 7 -odir tmp -olaz

creates 4 files:
  scale factor x y z:         0.01 0.01 0.01
  offset x y z:               0 0 0

based on our min/max values the total bounding box is calculated:
  bb-min = floor(min/tilesize)
  bb-max = ceil(min/tilesize)
  >> bbox 676720 246000 - 676880 246160

then, we see the bounding boxes of each file ONLY in the min/max values of the header:

676720_246000.laz
  number of point records:    304068
  min x y z:                  676750.00 246000.00 523.42
  max x y z:                  676806.99 246086.99 573.79
  LAStiling (idx 0, lvl 1, sub 0, bbox 676720 246000 676880 246160, buffer) (size 80 x 80, buffer 6.99)

676720_246080.laz
  number of point records:    80837
  min x y z:                  676750.00 246073.00 547.29
  max x y z:                  676806.99 246099.99 567.73
  LAStiling (idx 2, lvl 1, sub 0, bbox 676720 246000 676880 246160, buffer) (size 80 x 80, buffer 7)

676800_246000.laz
  number of point records:    373209
  min x y z:                  676793.00 246000.00 506.23
  max x y z:                  676849.99 246086.99 573.90
  LAStiling (idx 1, lvl 1, sub 0, bbox 676720 246000 676880 246160, buffer) (size 80 x 80, buffer 7)

676800_246080.laz
  number of point records:    114815
  min x y z:                  676793.00 246073.00 548.07
  max x y z:                  676849.99 246099.99 567.73
  LAStiling (idx 3, lvl 1, sub 0, bbox 676720 246000 676880 246160, buffer) (size 80 x 80, buffer 7)

So you see, the LAStiling entry is primary metadata about the WHOLE area.
idx: index of tile within a tile-set
lvl: level of index, increased, if a tile is sub-tiles (-refine_tiles)
bbox: bounding box of the whole tile set (original file or file set)
size: size of this tile, maybe not complete covered with data
buffer: buffer size

You can use floor/ceil functions of the min/max values to get the whole bounding box of the tile in full tile dimensions.
You can use the buffer value in the meta data to add/subtract this value from the min/max values to the the bounding box
without buffer. Please take care that the boundary is just added BETWEEN tiles, not at the outer bounding box.

Good luck!

Jochen @rapidlasso

[Marcel]

Hi, between posting the question and now I figured out that lastile seems to use a standard quadtree, so a possible solution to the above problem in python is presented below:

def get_quadtree_coords(level, cell_idx):

    """Standard quadtree, with sides of length 1"""

    # Initialize coordinates and cell size

    # Start in the bottom-left corner

    # x grows to the right, y grows upwards

    x, y = 0.0, 0.0

    cell_size = 1.0

    # Decode the ID

    for ctr in range(level):

        # halve each side of the cell to search

        cell_size /= 2

        # get the base-4 digit for the current level

        quadrant = (cell_idx // (4 ** (level - ctr - 1))) & 0b11

        # Update coordinates based on the quadrant

        if quadrant == 0b01:  # bottom-right

            x += cell_size

        elif quadrant == 0b10:  # top-left

            y += cell_size

        elif quadrant == 0b11:  # top-right

            x += cell_size

            y += cell_size

        # else stay in the bottom-left corner

    return x, y

def main():

    # use the numbers from an example "lastiling" line

    # LAStiling (idx 181, lvl 5, sub 0, bbox 407340 5578920 407980 5579560 buffer) (size 20 x 20, buffer 1)

    level = 5

    idx = 181

    sub = 0  # what is this used for?

    # coordinates encompassing all tiles

    left, bottom, right, top = 407340, 5578920, 407980, 5579560

    quadtree_size_x = right - left

    quadtree_size_y = top - bottom

    # size of this tile:

    tile_size_x, tile_size_y = 20, 20

    # buffer size:

    buffer = 1

    raw_x, raw_y = get_quadtree_coords(level, idx)

    # scale the coordinates

    raw_x *= quadtree_size_x

    raw_y *= quadtree_size_y

    # shift coordinates to actual location

    raw_x += left

    raw_y += bottom

    # tile bounds

    tile_left, tile_right = raw_x, raw_x + tile_size_x

    tile_bottom, tile_top = raw_y, raw_y + tile_size_y

    # tile+buffer bounds

    with_buffer_left = tile_left - buffer

    with_buffer_right = tile_right + buffer

    with_buffer_bottom = tile_bottom - buffer

    with_buffer_top = tile_top + buffer

    print(f"Tile bbox:\t{tile_left}\t{tile_bottom}\t{tile_right}\t{tile_top}")

    print(

        f"Tile+buff bbox:\t{with_buffer_left}\t{with_buffer_bottom}\t{with_buffer_right}\t{with_buffer_top}"

    )

if __name__ == "__main__":

    main()

For my dataset the "sub 0" part did not seem to play a role. but if it happens to not be 0, how should it be treated?

On Tuesday, July 9, 2024 at 8:09:31 PM UTC+2 Marcel wrote: