How to: Extract building heights from LiDAR data and make 3D buildings

Screen Shot 2015-09-25 at 11.19.00

The Environment Agency recently released their LiDAR as Open Data meaning it is now free to use and without restrictions. You can read a bit more about that here.

For those unfamiliar, LiDAR, which stands for Light Detection and Ranging, is a technology which involves firing a laser at a feature (usually the ground from a plane) and analysing the reflected light. The timing and pattern of the reflected light returned can be used to determine distance and type of feature scanned. A typical LiDAR system can scan hundreds of thousands of points a second producing incredibly detailed models of terrain, buildings, cars or whatever else they are pointed at.

The LiDAR data released by EA comes in 2 flavours:

  • DSM – Digital Surface Model
  • DTM – Digital Terrain Model

The DSM is exactly what you see on the surface of the earth with your own eyes; so depending on the resolution, this will include all man made features such as buildings, cars, infrastructure and so on as well as all the natural features such as vegetation.
The DTM is derived from the DSM where all the surface features (buildings, trees etc) have been stripped out leaving just the bare terrain.

Each product has it’s own uses, but one of the benefits of having access to both is the ability to extract the surface features by comparing the two datasets.

This is what I set out to do using some data I downloaded for the Bath area – it was surprisingly straightforward.

This is the end result:


Full screen version

For the step by step details…

You will need:

1) Download LiDAR data (both DSM and DTM) from here.  You will notice that depending on the area you click, different resolution products are available, e.g 25cm, 50cm, 1m, 2m. This is the distance between sampled points so the smaller the distance, the more detailed and the larger (file size) of the dataset. I went with 1m, the highest resolution for my area.

2) Extract the DTM and DSM zip files into separate folders. The contents will be a bunch of .asc files. These represent smaller tiles which make up the larger tile you downloaded. If you open one of these files in a text editor you will notice a bit of metadata at the top explaining where the data is located, resolution etc and then followed by a long list of values which are essentially just spot heights in meters.

3) In QGIS go to Raster > Misc > Build Virtual Raster (Catalog) and select all your DSM asc files and set an output file name. This will produce a .vrt which essentially creates a a single large virtual dataset from your individual .asc files.

 

Once complete, the VRT will be loaded onto the map and look something like this. (If prompted for a Coordinate Reference System (CRS), you will need to select EPSG27700, aka British National Grid.) You can see buildings as well as natural terrain features. Depending on your area, you may have gaps in LiDAR coverage like I do in the south east corner. Note that the bright white areas are water features – water absorbs LiDAR so no light is reflected and therefore no values can be determined.

Screen Shot 2015-09-23 at 20.22.22

 

4) Repeat the above process to create a DTM .vrt file.

5) The next step is to subtract the DTM from the DSM which will produce a new dataset containing the elevation of the surface features. To do this goto Raster > Raster Calculator and enter the following expression:

“combined_DSM_1m@1″ – “combined_DTM_1m@1″

Obviously substitute with your vrt layer names and also set the Output CRS to EPSG27700/British National Grid

Screen Shot 2015-09-23 at 20.28.41

This will produce a new layer which really highlights the surface features – even things like individual trees and field/garden boundaries are visible in my 1m dataset. If you are using 50 or 25cm data you should be able to pick out even finer details.

In this example we can see the parapet walls on the roof of The Circus and vegetation in the gardens

Screen Shot 2015-09-25 at 10.11.54

 

6) This step is optional –  we can view these surface features in 3D using the Qgis2threejs plugin (You can install it by going to Plugins > Manage and Install Plugins and searching for Qgis2threejs – you may need to restart QGIS after installing)

Here we have Bath Abbey – not bad from 1m LiDAR.

Screen Shot 2015-09-23 at 20.56.04

 

We can also drape aerial imagery on top of this model to produce a poor man’s version of Google and Apple’s 3D buildings layers.

Click the images below to see the full 3D interactive versions (Each one is around 8mb so may take a moment to load and requires a WebGL compatible browser)

Screen Shot 2015-09-23 at 21.30.24

Screen Shot 2015-09-25 at 11.19.26

Screen Shot 2015-09-25 at 11.19.00

7) As you can see from the above images, the buildings are not very well defined due to the 1m resolution of the LiDAR. We can get more discrete building footprints from another Open Data product, OS Open Map Local. Download the relevant tile in ESRI shape format from here.

8) Once complete, unzip and open the XX_Building.shp file in QGIS, setting the CRS to EPSG27700/British National Grid.

9) Since the building layer probably covers an area much larger than your LiDAR data we will clip the buildings to the LiDAR extents. Create a new vector layer by clicking this button on the left: Screen Shot 2015-09-25 at 10.22.35Then select Polygon, set CRS to EPSG27700/BNG and hit ok and save somewhere.

Screen Shot 2015-09-25 at 10.22.20

 

10) With your new layer selected, select the edit button (yellow pencil) then the add feature (green blob with yellow *) Screen Shot 2015-09-25 at 10.27.10Zoom out until you can see all of your LiDAR data and then draw a rough polygon around your LiDAR extents by clicking each corner. To close the shape, right click the last point. Once done, press the save button above (blue disk with red pencil).

11)  Now we can clip the buildings dataset to the polygon you just created by going to Vector > Geoprocessing Tools > Clip

Screen Shot 2015-09-25 at 10.33.06

12) We now need to calculate the height for each of the building footprints in your area. For this we need to install the Zonal Statistics plugin by going to Plugins > Manage and Install Plugins

13) Once installed, run the plugin by going to Raster > Zonal Statistics. Enter your relevant layers and calculate the mean.

Screen Shot 2015-09-25 at 10.37.05

14) Done!

We now have building footprints with an average height attribute.

From here you can use your building heights for whatever you want, I decided to render them into simple 3D buildings to display on top of OpenStreetMap.
While the building footprints from Ordnance Survey are pretty much up to date, the LiDAR data appears to be at least 6-7 years old (you can tell this by the lack of a few recent housing and shopping developments). This meant there are a few newer building footprints which don’t appear in the LiDAR so have no or very low heights. Therefore I applied a filter to my dataset to remove all building footprints with a height less than 1m.


Full screen version

The process to create these 3D buildings is actually pretty straightforward – I used Mapbox’s TileMill but you could also use Studio.

Import your building footprint shapefile into TileMill and use the ‘building-height’ cartoCSS style along with the column name containing your heights. Then apply a fill colour.

Here’s the full CartoCSS I used

#bathbuildings{
     building-height: '[mean]';
     building-fill:rgba(189,168,144,.9); 
}

Once you’re happy with how it looks, hit ‘Export’ and ‘upload’ to your Mapbox account.

Posted in GIS
  • Toneman1984

    Hey the link to the Lidar no longer works. Can you provide a working link?
    Thanks a lot!

    • brendan

      Thanks for letting me know, looks like they changed the URL slightly, I’ve updated the post accordingly. Cheers.

  • Phil Knight

    Great post. I’d love to give this a go.

    • brendan

      Thanks Phil, yeah you should definitely give it a go.

  • Ollie Brown

    Hi Nadnerb, thanks a lot for the post! I’ve created the DSM_combined, DTM_combined, and DSM_minus_DTM layers, and have installed the Qgis2threejs plugin. However, when I go to run the 3D viewer, the dialogue box says “The unit of current CRS is degrees, so terrain may not appear well.”, and the map comes out flat in my web browser. Any idea why this might be?

    • brendan

      Have you you got the CRS set to WGS84/EPSG4326? Try changing it to British National Grid / EPSG27700

      • Ollie Brown

        Yep that was it! I had the DSM_combined and DTM_combined files set to EPSG27700 but not the DSM_minus_DTM. Where did you get your aerial map from?

        • brendan

          I used the OpenLayers plugin if I remember correctly.

          • Ollie Brown

            Thanks!

  • Interested of Yeovil

    Hi Brendan … an excellent post. I’m particularly interested in getting the slope of the roof of a building. I tried adding the _max and _min attributes with the Zonal Statistics, but _min usually comes out as zero, whereas I would be looking for the height to the lowest part of the roof. The ultimate goal would be to calculate the volume of the roof. Cheers. Charles.