Dutch AHN3 point cloud visualization

The current Dutch Elevation (Actueel Hoogtebestand Nederland, AHN) map is a digital elevation map of the whole of the Netherlands. The map is created using massive laseraltimetry pointclouds.

There are multiple versions of AHN, the most recent and detailed one is AHN3.  AHN3 will have a complete coverage next year (2019) with aimed point density of 8 points/m2.

Image

In AHN3 each point has an x,y and z value and is classified into: ground surface (2), water (9), buildings (6), artificial objects (26), vegetation, and unclassified (1). In addition each point contains scanning information like number of returns, intensity and gps time.

The nice thing is that all the raw data is now available as open data. As these are huge files it takes some care to process. In this blog a method is described to get from a raw pointcloud file to a 3D visualization for a specific area (Johan Cruijff Arena  in Amsterdam) with open source tooling.

The following steps are performed:

1] Download raw pointcloud LAZ file

2] Get summary information about the file using PDAL

3] Crop the file to an area of interest using PDAL and bboxfinder.com

4] get information about 1 sample point using PDAL

5] Visualize in 3D viewer (plas.io)

Note: The PDAL steps are executed using Docker. As an alternative you can install PDAL tooling yourself instead of using Docker.

0] Preparation

Create a directory like d:\gisdata\ahn3:

$ cd d:\gisdata\ahn3

1] download a file

Go to  https://github.com/bertt/ahn3/blob/master/ahn3.geojson and select your favorite area.

Image(1)

We’ll download file C_25GZ1.LAZ as this file contains the Johan Cruijff Arena:

$ wget https://geodata.nationaalgeoregister.nl/ahn3/extract/ahn3_laz/C_25GZ1.LAZ

warning: its a 2GB file!

2] get summary info with PDAL

After downloading, run the following docker command to get  info about the LAZ file.

$ docker run -v d:/gisdata/ahn3:/data pdal/pdal pdal info /data/C_25GZ1.LAZ –summary

returns:

{ "filename": "\/data\/C_25GZ1.LAZ", "pdal_version": "1.7.2 (git-version: 07d19a)", "summary": { "bounds": { "X": { "max": 124999.999, "min": 120000 }, "Y": { "max": 481249.999, "min": 475000 }, "Z": { "max": 108.495, "min": -10.416 } }, "dimensions": "X, Y, Z, Intensity, ReturnNumber, NumberOfReturns, ScanDirectionFlag, EdgeOfFlightLine, Classification, ScanAngleRank, UserData, PointSourceId, GpsTime", "num_points": 470536792 } }

So it contains 470536792 (470 million) points with dimensions X, Y, Z, Intensity, ReturnNumber, NumberOfReturns, ScanDirectionFlag, EdgeOfFlightLine, Classification, ScanAngleRank, UserData, PointSourceId, GpsTime.

3] Find AOI

As the file is quite big, we’ll crop it to our area of interest. First we have to define the area of our interest in Dutch projection (epsg code 28992).

Go to http://bboxfinder.com/ and select your area of interest (Johan Cruijff Arena) in the map.

Switch to 28992 in the ‘epsg’ box.

Image(2)

Our area of interest is: 123606,480216,125816,481280 (xmin, ymin, xmax., ymax)

4] Crop area

Now with the raw pointcloud file and area of interest defined we can crop the file using PDAL pipeline command.

First create a crop.json file containing the following code:

{ "pipeline":[ "/data/C_25GZ1.LAZ", { "type":"filters.crop", "bounds":"([123606, 125816],[ 480216, 481280])" }, "/data/arena.laz" ] }

NB: the bounds are defined as [xmin, xmax], [ymin, ymax]

run with the PDAL pipeline with:

$ docker run -v d:/gisdata/ahn3:/data pdal/pdal pdal pipeline /data/crop.json

A new file ‘arena.laz’ with cropped area will be created.

5] Inspect a point

Now we can inspect a point in the arena.laz file to see its attributes.

$ docker run -v d:/gisdata/ahn3:/data pdal/pdal pdal info /data/arena.laz -p 0

returns:

{ “filename”: “\/data\/arena.laz”, “pdal_version”: “1.7.2 (git-version: 07d19a)”, “points”: { “point”: { “Blue”: 0, “Classification”: 1, “EdgeOfFlightLine”: 0, “GpsTime”: 328710.2203, “Green”: 0, “Intensity”: 20, “NumberOfReturns”: 1, “PointId”: 0, “PointSourceId”: 30811, “Red”: 0, “ReturnNumber”: 1, “ScanAngleRank”: 11, “ScanDirectionFlag”: 0, “UserData”: 2, “X”: 123985.86, “Y”: 480216.33, “Z”: -4.21 } } }

So the first point height is -4.21 (below sealevel) and it’s classification is 1.

6] Visualize in plas.io

Using a 3D viewer (like http://plas.io/) we can visualize the cropped point cloud in 3D.

Open the ‘arena.laz’ file in plas.io with ‘choose data to display’ option.

Image(3)

To get the colors fiddle around with the parameters on the right side:

set intensity source -> heightmap greyscale

set colorization -> classification

set intensity blending -> slide halfway

Another frequently used 3D desktop tool for pointclouds is CloudCompare (https://www.danielgm.net/cc/).

The same method can be applied to another area or another set of pointcloud data so have fun with that Smile

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s