lazer.pdal

The Point Data Abstraction Library PDAL is a new software system for reading and writing LiDAR data, including LAS/LAZ files.

PDAL aims to emulate GDAL and become the premier open source standard library for working with LiDAR data. It has been successful insofar as it is a complete nightmare to understand how it works, and nothing makes sense until you dig (and dig, and dig) into the code.

But we've wrapped a few commands to make it easier to work with PDAL command line tools that retrieve las file information, like projection and classification status.

---

Convenience functions for working with PDAL.

info(input_file, opts=['--stats', '--metadata', '--enumerate Classification'], raw=False)

Retrieves file info from a las/laz dataset.

Useful for getting information like the number of points, the spatial reference system, the point classes, etc.

Parameters:

Name	Type	Description	Default
input_file	str	path to the file to get info from.	required
opts	list	the command line options to pass to pdal info.	['--stats', '--metadata', '--enumerate Classification']
raw	bool	return the raw pdal info output, which is a shitshow.	False

Returns:

Type	Description
dict	a dictionary of the file metadata.

Source code in lazer/pdal.py

def info(
    input_file: str,
    opts: list = ["--stats", "--metadata", "--enumerate Classification"],
    raw: bool = False,
) -> dict:
    """Retrieves file info from a las/laz dataset.

    Useful for getting information like the number of points, the spatial
        reference system, the point classes, etc.

    Args:
        input_file: path to the file to get info from.
        opts: the command line options to pass to `pdal info`.
        raw: return the raw `pdal info` output, which is a shitshow.

    Returns:
        a dictionary of the file metadata.
    """
    command = f"pdal info {input_file} {' '.join(opts)}"
    rawinfo = run_command_line(command)
    info = json.loads(rawinfo)
    if raw:
        return info
    else:
        # store key spatial reference info
        spatial = {}
        spatial["proj4"] = info["metadata"]["srs"]["proj4"]
        spatial["wkt"] = info["metadata"]["srs"]["wkt"]
        spatial["units"] = info["metadata"]["srs"]["units"]
        spatial["bounds"] = [
            round(info["metadata"]["minx"] + info["metadata"]["offset_x"], 2),
            round(info["metadata"]["miny"] + info["metadata"]["offset_y"], 2),
            round(info["metadata"]["maxx"] + info["metadata"]["offset_x"], 2),
            round(info["metadata"]["maxy"] + info["metadata"]["offset_y"], 2),
        ]
        spatial["xsize"] = round(info["metadata"]["maxx"] - info["metadata"]["minx"], 2)
        spatial["ysize"] = round(info["metadata"]["maxy"] - info["metadata"]["miny"], 2)

        # store key point return info
        points = {}
        points["count"] = info["metadata"]["count"]
        points["compressed"] = info["metadata"]["compressed"]
        stats = info["stats"]["statistic"]
        for stat in stats:
            if stat["name"] == "Classification":
                classes = stat["values"]
                class_labels = []
                known_classes = list(LASClasses.keys())
                for number in classes:
                    if number in known_classes:
                        label = f"{number} ({LASClasses[number]})"
                    else:
                        label = f"{number} (user-defined)"
                    class_labels.append(label)
                points["classes"] = class_labels
            elif stat["name"] == "Z":
                points["z_min"] = stat["minimum"]
                points["z_max"] = stat["maximum"]
                points["z_avg"] = stat["average"]

        metadata = {
            "points": points,
            "spatial": spatial,
        }

        return metadata

reproject(input_file, output_file, epsg)

Translates xyz locations to a new projection, writing a new output file.

Parameters:

Name	Type	Description	Default
input_file	str	path to the file to reproject.	required
output_file	str	path to the reprojected file.	required
epsg	int	the output file EPSG authority spatial reference code.	required

Source code in lazer/pdal.py

def reproject(
    input_file: str,
    output_file: str,
    epsg: int,
) -> None:
    """Translates xyz locations to a new projection, writing a new output file.

    Args:
        input_file: path to the file to reproject.
        output_file: path to the reprojected file.
        epsg: the output file EPSG authority spatial reference code.
    """
    command_list = [
        "pdal translate",
        f"-i {input_file}",
        f"-o {output_file}",
        "-f reprojection",
        f"--filters.reprojection.out_srs='EPSG:{epsg}'",
        "--writers.las.compression='laszip'",
        "--writers.las.scale_x=0.01",
        "--writers.las.scale_y=0.01",
        "--writers.las.scale_z=0.01",
    ]
    command = " ".join(command_list)
    run_command_line(command)