I experimented with extracting the heights values from the drone imagery and add them to the building polygons traced in OpenStreetMap. Read along if you want to understand the process or just go here to see it in action here!

3D buildings around Foundation University, Dumaguete

Day 3 of Pista ng Mapa last Aug 3, 2019 was a series of field mapping activities. A couple of people did on the ground survey with field papers and mapillary. Our resident drone Mapher Leigh deployed her DJI Phantom 4 to survey the event venue and its surrounding community. Leigh uploaded all the drone derived data into OpenAerialMap including the elevation models! Using QGIS 3.6, I decided to explore ways to extract the heights from the derived DSM/DTM and use it for visualizing building polygons from OpenSteetMap.

For this exploration I used QGIS 3.6 with the GRASS GIS, Profile tool and QGIS2threejs plugin.

Downloading the drone derived layers from OpenAerialMap

• The drone derived raster layers are available for download at OpenAerialMap. You will need the 3 rasters listed below:

Orthomosaic	Digital Surface Model (DSM)	Digital Terrain Model (DTM)

We will use the orthomosaic as a background base image for our visualization. The DSM and DTM layers will be used for extracting the height of buildings. I will not talk about the difference between DSM and DTM here, but here’s a nice explanation from gis.stackexchange.

Resampling the DTM and DSM to 1 meter pixel resolution

The DSM and DTM we downloaded have different pixel resolution, to simplify the process, we will interpolate both raster pixel values to a coarser resolution of 1 meter.

• To resample to a coarser grid, we will use the GRASS GIS modules in QGIS. In the Processing toolbox, search for the r.resamp.stats.

• Select the DTM raster layer to resample, choose average for the Aggregation method and 1.0 for GRASS GIS region cell size.

• Save the file in your machine. The parameters I used are shown below:

Once you finished resampling the DTM, do the same process for the DSM.

Extracting height difference from DTM and DSM

Nest step is to subtract the values from DSM to the DTM to get a raster layer that only has the difference between the two. This will be the height value we will use for the buildings.

• In the Processing dialog, search for r.mapcalc.simple
• In the r.mapcalc.simple dialog, choose your dsm layer for Raster Layer A and your dtm layer for Raster Layer B. Select any raster layer for C-F (we will not use any of these layers but r.mapcalc.simple dialog requires all dropdown values to have a valid layer).
• In the Formula field, type A-B. This formula is simple difference between the two rasters.
• In the GRASS GIS region cell size, use 1.0
• Finally, choose to either use a temporary file or to save to file in the Calculated field.

• Once the process completes, a new raster layer will be displayed in your map canvass.

To visualize the data, you can use the Profile Tool plugin like shown below:

Height in meters for each layer along the transect (red line in the map canvas). Red - calculated difference, blue - DSM, purple - DTM.

Download OSM buildings

Now that we have the raster layer for the height, we will add these values to OSM buildings.

• Extract the buildings around the area using overpass-tubo
• I used the following overpass code to get the buildings around Foundation University.

[out:xml]/*fixed by auto repair*/[timeout:25]; ( node["building"]({{bbox}}); way["building"]({{bbox}}); relation["building"]({{bbox}}); ); out meta;/*fixed by auto repair*/ >; out meta qt;/*fixed by auto repair*/

• Once the extraction is complete, click Export and choose download/copy as GeoJSON
• Add the data in QGIS.

Re-project the buildings to the same CRS of the raster layers

In order to get the height values from the calculated height difference raster, we need both layers in the same coordinate reference system (CRS).

• To get the CRS of your rasters, select any of the raster within the Layers panel. Right-click and choose Properties. In the Layer Properties dialog, click on the Information tab. Our rasters are in EPSG:32651 - WGS 84 / UTM zone 51N - Projected CRS. We will use the same CRS for our building vector layer.
• To export the building layer to the same CRS, select the building within the Layers panel. Right-click and choose Export > Save Feature As …
• Choose any file format and filename but make sure to choose *EPSG:32651 - WGS 84 / UTM zone 51N as the CRS. Then click OK.

Get height of buildings from the calculated difference raster

Now that we have all data in the same CRS, we can start adding the height data from the raster layer to the building polygon.

• In the Processing dialog, search for v.rast.stats.
• In the v.rast-stats dialog, choose your building polygon in the Name of vector polygon map and the calculated heigh difference in the Name of raster map to calculate statistics from.
• In the Column prefix for new attribute columns type height_.
• Choose a filename to save your output and click Run.

• Once the process completes, a new building vector layer will be added to your map canvas that includes the univariate statistics of each polygon in the attribute table.

• Symbolize your map to show different colors of building according to height. In the image below, I used the median values (height__median) to show height for each building. I chose median value based on previous evaluation I did when we imported building heights in San Francisco from LiDAR derived data.

Visualize with QGIS2threejs plugin

Our processing is finally complete! We can start visualizing the heights in QGIS using the QGIS2ThreeJS plugin.

• In the QGIS menu, click Web > QGIS2threejs > QGIS2threejs Exporter. A new window will be shown to configure the rendering of your 3D map.
• Play around with the various parameters of the plugin to show the best results for your purpose. The following settings are what I used:
  • DEM - I used the dtm raster layer;
  • Polygon - I used the building vector with the extracted height values (Rast stats).
• In the Rast stats Properties, I chose the following:
  • Object type - Extruded;
  • Z coordinate - Relative to the DTM layer;
  • Height - Using the expression dialog I added "height__median" * 1.2, this formula exaggerates the height of building by a factor of 1.2,

• Once you are happy with the results, you can either export the visualization as an image or a webpage in the File menu of the QGIS2threejs Exporter

And its done! Run a local webserver then open the save index.html. See it in action here!

Some notes

• Buildings partially covered by trees can have wrong height. This is because the DSM layer includes not only height of buildings but also of trees and other features for example, in this small building (median_height = 11.6 meters).

In such cases where there is a lot vegetation obscuring the building, median_height may not be be best metric to use.

• This is not a full 3D, instead, it’s a “simple block model representation” (LOD1). My process cannot model roof structures or detailed architectural building models.

• Can we use this process for mapping building:height in OpenStreetMap? Definitely, but this requires further calibration and ground validation. For example, I noticed height anomalies along the edge of the raster data. This is likely due to the height extraction used in processing the raw image files along the edges where there is not enough data.

• Please share DSM/DTM data in OpenAerialMap. Extracting DSM/DTM from drone survey are fairly straightforward with the common processing tools like pix4d or OpenDroneMap. Although it will take more processing time, the benefits of using the DSM/DTM products along with the orthomosaic are huge. My appeal to all drone mappers uploading data in OpenAerialMap, please include the DSM/DTM. ;)

Give it a try with your own data and comment in this diary if I miss something. Happy mapping!