# Examples

## Filtering

This example shows how you can extract points from a file and write them into a new one. We use the classification field to filter points, but this can work with the other fields.

    import laspy

    las = laspy.read('tests/data/simple.las')

    new_file = laspy.create(point_format=las.header.point_format, file_version=las.header.version)
    new_file.points = las.points[las.classification == 1]

    new_file.write('extracted_points.las')

## Creating from scratch

There are multiple ways to create new las files.

### Creating a new LasData

    import laspy
    import numpy as np

    # 0. Creating some dummy data
    my_data_xx, my_data_yy = np.meshgrid(np.linspace(-20, 20, 15), np.linspace(-20, 20, 15))
    my_data_zz = my_data_xx ** 2 + 0.25 * my_data_yy ** 2
    my_data = np.hstack((my_data_xx.reshape((-1, 1)), my_data_yy.reshape((-1, 1)), my_data_zz.reshape((-1, 1))))

    # 1. Create a new header
    header = laspy.LasHeader(point_format=3, version="1.2")
    header.add_extra_dim(laspy.ExtraBytesParams(name="random", type=np.int32))
    header.offsets = np.min(my_data, axis=0)
    header.scales = np.array([0.1, 0.1, 0.1])

    # 2. Create a Las
    las = laspy.LasData(header)

    las.x = my_data[:, 0]
    las.y = my_data[:, 1]
    las.z = my_data[:, 2]
    las.random = np.random.randint(-1503, 6546, len(las.points), np.int32)

    las.write("new_file.las")

### Using LasWriter

    import laspy
    import numpy as np

    # 0. Creating some dummy data
    my_data_xx, my_data_yy = np.meshgrid(np.linspace(-20, 20, 15), np.linspace(-20, 20, 15))
    my_data_zz = my_data_xx ** 2 + 0.25 * my_data_yy ** 2
    my_data = np.hstack((my_data_xx.reshape((-1, 1)), my_data_yy.reshape((-1, 1)), my_data_zz.reshape((-1, 1))))

    # 1. Create a new header
    header = laspy.LasHeader(point_format=3, version="1.2")
    header.offsets = np.min(my_data, axis=0)
    header.scales = np.array([0.1, 0.1, 0.1])

    # 3. Create a LasWriter and a point record, then write it
    with laspy.open("somepath.las", mode="w", header=header) as writer:
        point_record = laspy.ScaleAwarePointRecord.zeros(my_data.shape[0], header=header)
        point_record.x = my_data[:, 0]
        point_record.y = my_data[:, 1]
        point_record.z = my_data[:, 2]

        writer.write_points(point_record)

## Using chunked reading & writing

This example shows how to use the ‘chunked’ reading and writing feature to split potentially large LAS/LAZ file into multiple smaller file.

    import argparse
    import sys
    from pathlib import Path
    from typing import List, Optional

    import numpy as np

    import laspy

    def recursive_split(x_min, y_min, x_max, y_max, max_x_size, max_y_size):
        x_size = x_max - x_min
        y_size = y_max - y_min

        if x_size > max_x_size:
            left = recursive_split(
                x_min, y_min, x_min + (x_size // 2), y_max, max_x_size, max_y_size
            )
            right = recursive_split(
                x_min + (x_size // 2), y_min, x_max, y_max, max_x_size, max_y_size
            )
            return left + right
        elif y_size > max_y_size:
            up = recursive_split(
                x_min, y_min, x_max, y_min + (y_size // 2), max_x_size, max_y_size
            )
            down = recursive_split(
                x_min, y_min + (y_size // 2), x_max, y_max, max_x_size, max_y_size
            )
            return up + down
        else:
            return [(x_min, y_min, x_max, y_max)]

    def tuple_size(string):
        try:
            return tuple(map(float, string.split("x")))
        except:
            raise ValueError("Size must be in the form of numberxnumber eg: 50.0x65.14")

    def main():
        parser = argparse.ArgumentParser(
            "LAS recursive splitter", description="Splits a las file bounds recursively"
        )
        parser.add_argument("input_file")
        parser.add_argument("output_dir")
        parser.add_argument("size", type=tuple_size, help="eg: 50x64.17")
        parser.add_argument("--points-per-iter", default=10**6, type=int)

        args = parser.parse_args()

        with laspy.open(sys.argv[1]) as file:
            sub_bounds = recursive_split(
                file.header.x_min,
                file.header.y_min,
                file.header.x_max,
                file.header.y_max,
                args.size[0],
                args.size[1],
            )

            writers: List[Optional[laspy.LasWriter]] = [None] * len(sub_bounds)
            try:
                count = 0
                for points in file.chunk_iterator(args.points_per_iter):
                    print(f"{count / file.header.point_count * 100}%")

                    # For performance we need to use copy
                    # so that the underlying arrays are contiguous
                    x, y = points.x.copy(), points.y.copy()

                    point_piped = 0

                    for i, (x_min, y_min, x_max, y_max) in enumerate(sub_bounds):
                        mask = (x >= x_min) & (x <= x_max) & (y >= y_min) & (y <= y_max)

                        if np.any(mask):
                            if writers[i] is None:
                                output_path = Path(sys.argv[2]) / f"output_{i}.laz"
                                writers[i] = laspy.open(
                                    output_path, mode="w", header=file.header
                                )
                            sub_points = points[mask]
                            writers[i].write_points(sub_points)

                        point_piped += np.sum(mask)
                        if point_piped == len(points):
                            break
                    count += len(points)
                print(f"{count / file.header.point_count * 100}%")
            finally:
                for writer in writers:
                    if writer is not None:
                        writer.close()

    if __name__ == "__main__":
        main()

## COPC

    import time

    import laspy
    from laspy import Bounds, CopcReader

    def create_query(header):
        querys = []

        sizes = header.maxs - header.mins

        # Bottom left
        query_bounds = Bounds(mins=header.mins, maxs=header.mins + sizes / 2)
        query_bounds.maxs[2] = header.maxs[2]
        querys.append(query_bounds)

        # Top Right
        # Bounds can also be 2D (then all Z are considered)
        query_bounds = Bounds(mins=(header.mins + sizes / 2)[:2], maxs=header.maxs[:2])
        querys.append(query_bounds)

        return querys

    def main():
        # path = "http://localhost:8000/autzen-classified.copc(1).laz"
        # path = "http://localhost:8000/sofi.copc.laz"
        # path = "https://s3.amazonaws.com/hobu-lidar/autzen-classified.copc.laz"
        # path = "https://s3.amazonaws.com/hobu-lidar/sofi.copc.laz"
        # path = "https://s3.amazonaws.com/hobu-lidar/montreal-2015.copc.laz"
        path = "https://s3.amazonaws.com/hobu-lidar/melbourne-2018.copc.laz"

        with CopcReader.open(path) as crdr:
            print("copc_reader_ready")
            # querys = create_query(crdr.header)
            # for i, query_bounds in enumerate(querys):
            #     resolution = None
            #     points = crdr.query(query_bounds, resolution=resolution)
            #     print("Spatial Query gave:", points)
            #     print(len(points) / crdr.header.point_count * 100);

            # new_header = laspy.LasHeader(
            #     version=crdr.header.version,
            #     point_format=crdr.header.point_format
            # )
            # new_header.offsets = crdr.header.offsets
            # new_header.scales = crdr.header.scales
            # with laspy.open(f"output_{i}.las", mode="w", header=new_header) as f:
            #     f.write_points(points)

    if __name__ == "__main__":
        t0 = time.time()
        main()
        t1 = time.time()
        print("It took:", t1 - t0, " seconds")