Examples

Reading a LAZ file

We recommend you use laspy (pip install 'laspy[lazrs]' to install) to read a LAS/LAZ into a NumPy array, and then pass that array to startinpy:

import laspy
import numpy as np
import startinpy

las = laspy.read("../data/small.laz")
pts = np.vstack((las.x, las.y, las.z)).transpose()
pts = pts[::1]  # -- thinning to speed up, put ::10 to keep 1/10 of the points
dt = startinpy.DT()
dt.insert(pts)
print("number vertices:", dt.number_of_vertices())
# -- number vertices: 39231

Exporting the DT to GeoJSON

import numpy as np
import startinpy

# -- generate 100 points randomly in the plane
rng = np.random.default_rng(seed=42)
pts = rng.random((100, 3))
dt = startinpy.DT()
dt.insert(pts, insertionstrategy="AsIs")
dt.write_geojson("myfile.geojson")

Exporting the DT to several mesh formats with meshio

import startinpy
import meshio
import laspy
import numpy as np

las = laspy.read("../data/small.laz")
pts = np.vstack((las.x, las.y, las.z)).transpose()
dt = startinpy.DT()
dt.insert(pts)
vs = dt.points
vs[0] = vs[1] #-- to ensure that infinite vertex is not blocking the viz
cells = [("triangle", dt.triangles)]
meshio.write_points_cells("mydt.vtu", vs, cells)

Reading a GeoTIFF file with rasterio

We can use rasterio to read a GeoTIFF and triangulate the centre of the pixels/cells directly. Notice that retrieving the (x,y)-coordinates of the centres with the xy() function of rasterio is super slow and it’s better to use the code below.

Notice that we use the insertion strategy “BBox” because it is several orders of magnitude faster for gridded datasets. The code also randomly selects 1% of the points.

The no_data values are not inserted in the triangulation.

This code saves the resulting triangulation to a PLY file that can be opened directly in QGIS (with the newish MDAL mesh).

import startinpy
import rasterio
import random

d = rasterio.open('../data/dem_01.tif')
band1 = d.read(1)
t = d.transform
pts = []
for i in range(band1.shape[0]):
    for j in range(band1.shape[1]):
         x = t[2] + (j * t[0]) + (t[0] / 2)
         y = t[5] + (i * t[4]) + (t[4] / 2)
         z = band1[i][j]
         if (z != d.nodatavals) and (random.randint(0, 100) == 5):
             pts.append([x, y, z])
dt = startinpy.DT()
dt.insert(pts, insertionstrategy="BBox")
#-- exaggerate the elevation by a factor 2.0
dt.vertical_exaggeration(2.0)
dt.write_ply("mydt.ply")
_images/mdal.jpg

3D visualisation with Polyscope

You need to install Polyscope (basically pip install polyscope).

import laspy
import numpy as np
import polyscope as ps
import startinpy

las = laspy.read("../data/small.laz")
pts = np.vstack((las.x, las.y, las.z)).transpose()
pts = pts[::10]  # -- thinning to speed up, put ::10 to keep 1/10 of the points
dt = startinpy.DT()
dt.insert(pts)

pts = dt.points
pts[0] = pts[1]  # -- first vertex has inf and could mess things
trs = dt.triangles

ps.init()
ps.set_program_name("mydt")
ps.set_up_dir("z_up")
ps.set_ground_plane_mode("shadow_only")
ps.set_ground_plane_height_factor(0.01, is_relative=True)
ps.set_autocenter_structures(True)
ps.set_autoscale_structures(True)
pc = ps.register_point_cloud(
    "mypoints", pts[1:], radius=0.0015, point_render_mode="sphere"
)
ps_mesh = ps.register_surface_mesh("mysurface", pts, trs)
ps_mesh.reset_transform()
pc.reset_transform()
ps.show()
_images/polyscope_gui.jpg

Plotting the DT with matplotlib

import numpy as np
import startinpy

# -- generate 100 points randomly in the plane
rng = np.random.default_rng(seed=42)
pts = rng.random((100, 3))
# -- scale to [0, 100]
pts = pts * 100
t = startinpy.DT()
t.insert(pts)
pts = t.points
trs = t.triangles
# -- plot
import matplotlib.pyplot as plt

plt.triplot(pts[:, 0], pts[:, 1], trs)
# -- the vertex "0" shouldn't be plotted, so start at 1
plt.plot(pts[1:, 0], pts[1:, 1], "o")
plt.show()
_images/matplotlib.png

Gridding the dataset with spatial interpolation

import json
import math
import laspy
import numpy as np
import rasterio
import startinpy
from tqdm import tqdm


def main():
    las = laspy.read("../data/small.laz")
    dt = startinpy.DT()
    dt.duplicates_handling = "Highest"
    d = las.xyz
    print("Constructing the TIN with {} points".format(len(d)))
    for each in tqdm(d):
        dt.insert_one_pt(each)

    # -- grid with 50cm resolution the bbox
    bbox = dt.get_bbox()
    cellsize = 0.5
    deltax = math.ceil((bbox[2] - bbox[0]) / cellsize)
    deltay = math.ceil((bbox[3] - bbox[1]) / cellsize)
    centres = []
    i = 0
    for row in range((deltay - 1), -1, -1):
        j = 0
        y = bbox[1] + (row * cellsize) + (cellsize / 2)
        for col in range(deltax):
            x = bbox[0] + (col * cellsize) + (cellsize / 2)
            centres.append([x, y])
            j += 1
        i += 1
    centres = np.asarray(centres)
    print("Interpolating at {} locations".format(centres.shape[0]))
    zhat = dt.interpolate({"method": "TIN"}, centres)

    # -- save to a GeoTIFF with rasterio
    write_rasterio(
        "grid.tiff", zhat.reshape((deltay, deltax)), (bbox[0], bbox[3]), cellsize
    )


def write_rasterio(output_file, a, bbox, cellsize):
    with rasterio.open(
        output_file,
        "w",
        driver="GTiff",
        height=a.shape[0],
        width=a.shape[1],
        count=1,
        dtype=np.float32,
        crs=rasterio.crs.CRS.from_string("EPSG:28992"),
        nodata=np.nan,
        transform=(cellsize, 0.0, bbox[0], 0.0, -cellsize, bbox[1]),
    ) as dst:
        dst.write(a, 1)
    print("File written to '%s'" % output_file)


if __name__ == "__main__":
    main()
_images/grid.png

Creating a DSM

To create the DSM from an aerial lidar dataset of a city, one wants to preserve the highest z-value for each location.

With startinpy, this can be performed easily, as the example below shows. The LAZ file is read using the Python library laspy, a rather large tolerance for xy-duplicates is set (0.10m), and the highest z-value for each xy-location is kept. Furthermore, the intensity property of the input LAZ points is preserved.

The resulting file is exported to the PLY format, which can be read by several software, the open-source QGIS being one of them.

import laspy
import numpy as np
import startinpy

las = laspy.read("../data/small.laz")
pts = np.vstack((las.x, las.y, las.z, las.intensity)).transpose()
dt = startinpy.DT(np.dtype([("intensity", float)]))
dt.snap_tolerance = 0.10
dt.duplicates_handling = "Highest"
for pt in pts:
    dt.insert_one_pt([pt[0], pt[1], pt[2]], intensity=pt[3])
dt.write_ply("mydt.ply")
_images/qgis.png