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.
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.
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.
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.
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 
