Background

For test purposes, I got LIDAR + FKB-data (buildings, roads, rivers, lakes etc.) for an area north of Bergen in Norway. The LIDAR data is on LAS-format and the FKB data is on a Norway specific SOSI format. The making an orienteering map of this is described here - now I want to show this in 3D in the browser.

Simple test: 10m data from Kartverket

First up is to repeat the work of the ThematicMapping blog - but with data from the area of interest. For the gdal stuff I use the FWTools windows binary which was linked to in the Attackpoing discussion (note, only the mirror download site worked when I tried to download it).

I first downloaded the 10m data from here - the data is free (you need to register to download, though)
Then I did the conversion for the grid described in this blog post for my data:
```
gdalwarp holsnoy.vrt holsnoy.tif
gdal_translate -scale 0 500 0 255 -outsize 1000 1000 -of PNG holsnoy.tif holsnoy.png
gdal_translate -scale 0 500 0 65535 -ot UInt16 -outsize 200 200 -of ENVI holsnoy.tif holsnoy.bin
```
Note that compared to the blog post, I changed the second parameter of "scale" to a smaller number. This is the min/max meters. My area had a maximum height of around 500 meters.

Then we also want to overlay the map image. I did the following to get a nice graphical overview:

gdal_translate -outsize 1000 1000 -of GTiff holsnoy.tif holsnoy_small.tif
gdaldem  color-relief holsnoy_small.tif color-relief-holsnoy.txt -of png holsnoy-relief.png

Here is the color-relief file I used (saved as "color-relief-holsnoy.txt"):

0 80 80 255
5 110 220 110
100 240 250 160 
200 230 220 170 
300 220 220 220
400 250 250 250

Finally I downloaded the sample files described on the ThematicMapping blog from Github, and did a few small modifications in the file to fit my data,

[..]
  terrainLoader.load('holsnoy.bin' ... // my bin-file from above
[..]
  geometry.vertices[i].z = data[i] / 65535 * 5; // Adjust scale of height
[..]
  map: THREE.ImageUtils.loadTexture('holsnoy-relief.png') // my relief
[..]

The resulting 3D mode in the browser is available here

Karttapullautin: Single tile

The process for a single Karttapullautin tile is not much more difficult. I use this laz-file which I renamed to "laztest.laz". I decided to take the easy way and just convert from the LAZ-format to the DEM format which was used for the first sample. There might be other / easier ways to do it, but I decided to make it easy for myself first.

So here is how I did it:

Convert the laz-file to a DEM file using laz2dem (from lastools).
```
las2dem laztest.laz -o laztest.asc
```
This however gives me a black diagonal in the resulting file as I need to have a license to work with this many las points in laz2dem. So insteaed of finding an alternative way to convert to DEM, I reduce the number of points in my las-file by thinning out my las-data using las2las. I don't need that many points for the remaining process anyway:
```
las2las -o laztest_thin.laz laztest.laz -keep_class 2 -keep_random_fraction 0.1
```
Note that I also use "-keep-class 2" to only keep the ground points.
Now I can make my DEM-file (asc-format) using the following command:
```
las2dem laztest_thin.laz -o laztest.asc
```

Proceeding with the asc-file, I can do exactly as I did above for the DEM-file:

gdalbuildvrt laztest.vrt laztest.asc
gdalwarp laztest.vrt laztest.tif
gdal_translate -scale 0 200 0 65535 -ot UInt16 -outsize 200 200 -of ENVI laztest.tif laztest.bin

Instead of using the color relief, I this time simply use my "pullautin.png" from processing the same laz-tile as we are using now.
Putting it all together, we get this 3D model of the Pullautin tile. Note that the Pullautin sample in the actual 3D file is with roads etc. while the one pictures below is the older version without roads.

Adding hillshading

Adding hillshading should be fairly easy. See this 3D model of the Pullautin tile with hillshading Some comments about hillshading:

Note that for generating nice looking hillshaded DTMs it is often it is better to first create an elevation raster in BIL format and then hillshade this raster. Why? Because shading the (often small and irregular) triangles is lightly to generate jitter and aliasing during the lighting computation of the (often instable) surface normals of the (often tiny or slivery) triangles. After rasterization onto a BIL the LiDAR gets essentially resampled and using these regular spaced points avoids all these shading irregularities. Hence instead of las2dem -i lidar.laz -keep_class 2 -hillshade -o dtm.png try running las2dem -i lidar.laz -keep_class 2 -o dtm.bil las2dem -i dtm.bil -hillshade -o dtm.png

Complete map and/or part of map

The final step is to take a complete map / several tiles and/or a normal orienteering map. As long as my las-file (and therefore the resulting bin-file) is exactly the same size as the full orienteering map/Pullautin map, everything should be 100% OK and no change is needed for the procedure. It should be possible to do this for any Pullautin map. For a normal orienteering map, I need to be a bit more careful, but that should be possible as well.

Edit: Turns out it should be possible to work directly with xyz-files without converting from laz to dem using e.g. this or this. I won't dive into this now, but maybe later some time... Maybe one could even get the elevation data from a WMS-server like they do here. This would only give the 10m data and not the accurate elevation data used in the Pullautin generation, but would still be quite elegant. Note that it also looks like I can get the 10m data for Norway directly using vrt-files in GDAL.

Edit 2: Looking a bit further, I also found this nice service for plotting a GPS-track in Norway on a 3D map - some nice elements there also.

Following a track / rotation