3D printing of digital elevation models with QGIS

Eiger area in Switzerland. A bit too much stringing

Introduction

Using QGIS, the GDAL library and the QGIS DEMto3D extension is an attractive toolset for 3D printing of digital elevation models. QGIS is a free GIS tool and its DEMto3D plugin allows to produce printable STL files. GDAL is a raster manipulation library available through the "Raster" menu in QGIS. It allows format translation, merging, cropping, etc.

We found this toolset to be both flexible and relatively easy to use. Recommended if you plan do create elevation models from a variety of data formats. In this article, we describe how to use QGIS as well as command line tools to:

• merge DEM tiles obtained from SRTM (or other) data servers
• to clip tiles
• and to translate DEM tiles to STL models

Read also:

• 3D printing of digital elevation models
• Digital elevation model

Updated: Feb 2019.

QGIS is a GIS program to "Create, edit, visualise, analyse and publish geospatial information on Windows, Mac, Linux, BSD (Android coming soon)". According to Wikipedia, “Similar to other software GIS systems, QGIS allows users to create maps with many layers using different map projections. Maps can be assembled in different formats and for different uses. QGIS allows maps to be composed of raster or vector layers. Typical for this kind of software, the vector data is stored as either point, line, or polygon-feature. Different kinds of raster images are supported, and the software can georeference images.”

GDAL is a “translator library for raster and vector geospatial data formats that is released under an X/MIT style Open Source license by the Open Source Geospatial Foundation. As a library, it presents a single raster abstract data model and single vector abstract data model to the calling application for all supported formats. It also comes with a variety of useful command line utilities for data translation and processing”. Besides QGIS and the command line there are other options for using it, see the list.

Using QGIS to manipulate GIS files

Before creating a printable STL terrain model, you may have to go through a series of file transformations.

In addition, you will need to retrieve DEM files as explained in the main 3D printing of digital elevation models article. A good start is using CGIAR-CSI SRTM files.

Getting and opening a DEM file

From most repositories you will have to download a "square" of the Earth and then extract the part you want to model. Sometimes you will have to merge two squares. If you are not familiar with various formats and repositories, we suggest getting a geoTIFF file from CGiar

• Click on one or more squares in the world map and then click "Search"
• In the next page displayed you will see a selection
• Select the one you want, download the zip file.
• There is no need to unpack the zip file

Now in QGIS:

• In the browser panel (to the left), browse in your directory structure until you find the *.zip
• click on the + to open it
• Select the *.tif file

Some canary islands (tile srtm_33_7)

Translate files

Sometimes you may get a file in a format that is not suitable. QGIS, i.e. the incorporated GDAL library, can read and translate several formats

In QGIS:

• Explore the raster menu

If you installed GDAL, you also can use a command line like this:

gdal_translate can translate many formats, e.g. all common raster and grid formats

gdal_translate -of output input

Merge images

Since terrain elevation data are available as squares, you may have to merge them in order to extract an area that sits "between tiles". The gdal_merge.py library can do that. Let's also mention that some software also allows working with virtual rasters, i.e. can process several tiles at the same time if they were a real picture. Merging images in QGIS is done with the raster menu:

QGIS: Menu Raster -> Miscellaneous -> Merge.....

• Select all the files you plan to merge
• Define the output file
• Set No data value = 0 for the output (better than -32676, see below)
• Make sure that "Load into Canvas when finished" is ticked.

Merging *.bil tiles with QGIS

As you can see, there seem to be some "holes" in the result, but that is another problem. If you see such holes, it can mean that you have "no data" values. In our case the number -32767 used for these turn out to be black.

• To see the nodata value of an input file, right-click on the file in the browser panel and look at properties.

Using the command line tool (btw. QGIS will always show command lines and arguments that you use):

gdal_merge.py -of file_format -o output_file input, input, input, .....

Example:

gdal_merge.py -n -32767 -a_nodata -32767 -of GTiff -o geneva_area.bil n46_e007_1arc_v3.bil n45_e007_1arc_v3.bil n46_e006_1arc_v3.bil n45_e006_1arc_v3.bil n46_e005_1arc_v3.bil n45_e005_1arc_v3.bil

Dealing with no data values (missing data)

Missing data (called "no data") is already marked as such. They may appear as totally black or white "holes" and you cannot and need not do anything about these if they are very small. Larger areas though need to converted to real values if you plan to 3D print.

In QGIS, load the file into a raster layer if not already done so (see merging above).

• Menu: Layer -> Add Raster Layer
• Select a supported file format, e.g. *.asc or *.bil

Now you can decide which values are missing values. In the picture above (merge images), we have single "band" of color: grey. Darker means lower regions, lighter means higher regions. Black is likely "freaky" data, i.e. impossibly high or low. Data points that point deep into the earth, i.e. Min = -32767 are "missing" data and should be ignored:

Min=-32767.000 Max=4267.000 Mean=410.805, StdDev=4753.033

• Right click on the file name in the layers panel
• Select properties
• See if there are "no data value". In our case, we had -32767 which is correct.

no data value STL

Converting missing data to another value

If your missing data represent the bottom of an model that was clipped in the z direction, e.g. you got an island but no sea, you could transform the missing data values to 0. "Reclassifying" no data values is neccessary for using the Demto3D plugin, just described below.

This situation is very easy to spot before you even attempt to produce the STL:

• The background of your raster is white (while the highest points also look white)
• In addition, in the DEMto3D popup, you also would see the highest point as 32767

A good solution is therefore to transform the 32767 value (or whatever else) into the lowest value of your model, e.g. 0. To do so use the Reclassify values tool that is part of the SAGA tools (a bit difficult to find):

Processing menu -> Toolbox -> Select SAGA in the panel to the right -> Raster tools -> Reclassify values

Now, in the popup window, scroll down and uncheck replace other values. If you want another value than 0, you can change that. The resulting new raster now should have a black background.

No value data (missing data) being translated to 0

Otherwise, e.g. if you got missing areas within a mountain, you may have to use an algorithm that will interpolate values (I did not investigate how to do this yet, simple substitution will not work, unless you later can smooth them out with Meshlab.).

Crop/extract tiles

We suggest cropping a tile somewhat before using the DEMtoSTL extension. Since you also can clip in the DEMto3D extension, you can clip to a larger size than needed at this point.

Identify the area

Firstly we suggest to identify precisely or roughly the coordinates. In google, typing for example "coordinates el capitan" will give this:

Coordinates of El Capitan

You then could enter the following values in the popup window. We chose to extract a 0.4 x 0.4 sized square, i.e. roughly a 44km by 44km area. One degree is about 111km (if we are right!) and here we got 0.4.

-119.8,-119.4,37.5,37.9 [EPSG:4326]

xmin = -119.8 (longitude W, left, more west)

xmax = -119.4 (longitude W, right, less west)

ymin = 37.5 (lattitude N, bottom, closer to equator !)

ymax = 37.9 (lattitude N, top, closer to north pole)

Later, in the Demto3D tool we can clip smaller areas, e.g. -119.7, -119.55, 37.65, 37.8

Clipping with QGIS

In the QGIS GUI, you can crop (extract an square area) from a tile.

• In QGIS V.2x: Menu Raster -> Extraction -> Clipper
• QGIS V.3.4: Menu Raster -> Extraction -> Clip Raster by Extent

Procedure:

• In the Popup window, click on the button to the right of Clipping extent and either select "Select extent on Canvas" or enter some values, e.g. something like the values above. You also can combine the two methods. First select with the mouse, then fix values. Make sure to understand in which direction the coordinates go...

Make sure to get the order right if you enter your own values and do not forget the "-" if you are in the western or southern hemisphere. To get a feel for right type of coordinates, you can manually select some area, then modify. It usually takes 4-5 trials to get the right extent. If you select with the mouse, zoom into the area and make it as big as you can, the select the rectangle

• Press OK

QGIS popup clip raster by extent, typing values for an el capitan square

Saving a clipped area (not mandatory)

• Right-click on the new created layer (bottom left panel)

export -> save as

• Define the name of an output file. Make sure that the file path is correct (by default the program attempts to save in a wrong place). File extension will define the output type. By default the "tif" format is used.

QGIS Save clipped layer to file (optional)

Alternatively, use the command line: Using the command line tool you could provide coordinates, e.g.

gdal_translate -projwin 5.68197537628 46.6941320621 7.28048142324 45.6929800494 geneva_region.bil geneva_small.bil

Polygonize raster files

In default installations, QGIS and GDAL cannot directly produce usable 3D formats for printing. However, one can produce some polygone formats that one could further translate to STL. We do not suggest using this feature, since we found the resulting files fairly unusable. See the next item.

Create STL from DEM files with the QGIS DEMto3D extension

Prerequisites

• DEMto3D must be installed, see the instructions at the end of this article.
• You must have a DEM file, see the instructions at the start of this article.

DEMto3D plugin with missing height base

Procedure for using DEMto3D

Tif to STL for Gran Canaria, correctly clipped (DEMto3D display problem on a 4K screen)

DEMto3D plugin. Clipping El Capitan area from a larger square

The plugin is available in the Raster menu (once installed):

• If you already clipped the tile, then select full extent (click on the little globe), else give coordinates
• Now make sure to untick all other models from the layers panel, we also suggest killing the ones that you won't need. On the screen you only should see a single layer. Make sure that the red dotted lines surround the object. Else, close QGIS and just load this file.
• Define print extent first, else there will be error messages !!
• Define width or length, e.g. 150mm (15cm) for a good build plate or 32mm for a 4x4 lego block. Alternatively you could enter the scale and then decide how large you can/want print and adjust the width now.
• Define spacing (we suggest 0.25 or 0.3 mm for larger models)
• Exaggeration factor: 2.0 (or less or more) for large models and 1.5 for smaller models
• Model height: Enter the lowest point in meters, i.e. 0 in most cases. If you model the K2, you could start at 3500 meters though. Entering a negative number will have the plugin crash.

We tested this plugin with several data sources and several file sizes. It can produce large and detailed STL files.

Exporting a larger mode will take a few minutes and much more on a slow computer. Resulting STL files can have twice this size. We extracted Gran Canaria from a CGiar tile (srtm_33_7) that included also Tenerifa and other "middle" islands. The tif was 625KB (rather small) and the resulting STL was 775 MB, i.e. very large. I used both Meshlab and Meshmixer to reduce the triangles as explained in the main article.

Below is the area around Geneva made from 30m SRTM data. The raw result as seen in Meshlab is rather pretty. Of course, it may be pointless to create a model of such a large area with 30m data. Creating model of small areas with more precise data, however, is an option one must consider.

DEMto3D plugin. View of the Montblanc from somwhere over the Jura. Geneva is to the right and Lac d'Annecy top right

Below is a screenshot of the almost ready-to-print 10x10 cm model of the Salève, the Geneva house mountain. It was done with the following parameters:

• Full extent (I clipped the area before)
• Size width = 100
• Lowest point = 100
• Exageration factor =1

I did not make any changes to the output file. Of course it would be wise to do a preventive Meshrepair. Once also could consider augmenting the exaggeration factor to 1.2 or 1.5 maybe.

DEMto3D plugin. View of the Salève (Netfabb). As you can see model fits on the platform as expected and is well positioned.

DEMto3D plugin. View of the Salève (Meshlab)

Post processing of the generated STL

The STL output is ascii (not binary), i.e. it could be compressed quite a lot to binary STL. If you started with low resolution data, you have too many triangles which is a not problem per se, but can drastically augment slicing time. Also you cannot email 500MB to 1TB files even if you compress these to 50 to 100MB files . Therefore, read in the overview article how to reduce meshes with the free Meshlab program. Alternatively, use another tool, e.g. Meshmixer or Microsoft 3D builder. In addition, I recommend repairing the resulting file with a tool like Netfabb.

Also, I do have the feeling that some DEM models should be emptied with a tool like Meshmixer, using its "make solid" functions. I had more than one print abort because the second hot end got stuck because too much plastic got deposited in models that were not optimally sliced. But I did not have time to look into the problem, i.e. analyse g-code and test different slicers. In any case, volano printing is much harder than printing nice elegant Voronoi structures.

For starters, I printed two models of the famous Vulcano island using both 2m and 15m resolution *.tif files. Before printing I did reduce a large STL file (made from 2m data) using Meshlab. Before reduction the ascii STL file size was 590MB, the reduced binary STL was 16MB, i.e. more manageable. After reducing I had to fix the STL in Netfabb and it was good for printing. I also printed a model of just the volcano made from 15m data and in addition made some small custom Legos to distribute to the children of the island. Later, I also produced smaller models of Geneva area, Canari islands (e.g. Gran Canaria) and the Vercors area in France.

DEMto3D plugin. View of Vulcano island (2m data) with a 0.1 reduced mesh (Meshlab)

I do recommend this tool chain, if you work with high resolution tiles and/or if you have data that need some processing (e.g. tiles that need to be assembled and clipped, or no data value files)..

Issues

(as of March 2017)

Missing data

• If your STL model is deeply hidden at the bottom of a huge tunnel block then you have a no data problem. E.g. if you have a little island and the sea around the island has no data coded with "32767" then the island will sit at the bottom of a 32km high tunnel. See above how to deal with no data values. IMHO this plugin should not compute "no data values", but it does :)

Messed up display with high resolution screens

An older version of this plugin did not adapt to high resolution screens under Windows (QGIS 2.18.4 to 2.18.14). This has been fixed in summer 2019, therefore upgrade'. If you still must use an old version or still got problem in your configuration:

• Workaround 1: use the TAB to move forward and type blindly. Look at the screenshot above. You probably still can fill in "model size" paramters. It is important to fill in at least the "width", e.g. 150 (mm). You will have to fill in Exaggeration factor (type something like "2.0" and hit TAB, then type e.g. "0" for the lowest point.
• Workaround 2: Adjust screen resolution to a low one (e.g. HD) and reduce font size to 100% (typically you would use between 150% and 250% on a high resolution screen and this is the problem). This will require login out and back in.
• Workaround 3 (my preferred solution for Windows):
  • Read in the Inkscape article how to lower DPI Scaling for just some applications.
  • In short, right click on QGIS shortcut in the file manager, and select Properties
  • Select Compatibility TAB and click on Change high DPI settings. Then, play with the settings. I ticked both boxes and selected System (enhanced.

DEMto3D problem on a 4K screen

The author of the plugin can probably not see this. Since it took Adobe several years to figure that people actually do use screens with lots of pixels, we won't complain here.

Crash under Ubuntu

• The application may crash under Ubuntu 16/QGIS 2.14.1. Workaround: Use Windows (or maybe install a newer version by hand). Also, avoid using network drives.

Empty STL

• If you use a Network drive (hapened both on Ubuntu and Windows), try saving locally. A file with 347053 MB only worked that way.

Installation tips

Installation of QGIS under Windows 64bit

Take the latest 64 bit version installer from here.

The program is rather large: download is 395MB and installation is 1.6GB.

Installation of QGIS under Ubuntu 16x

To install the huge QGIS Geographic Information Systems Package, you may need to configure a GIS repository

sudo apt-get install python-software-properties
sudo add-apt-repository ppa:ubuntugis/ppa
sudo apt-get update

Now you can select dozens of packages described in Packages in “ubuntugis-stable”. To install QGIS:

sudo apt-get install qgis

Installation of QGIS under Ubuntu 18x

The required "Palma" version is included in the standard distribution

sudo apt install qgis

If you want a better version, i.e. QGIS 3.4.x Madeira LTR

• edit /etc/apt/sources.list and add:

deb https://qgis.org/debian-ltr bionic main

• Add the signature

• Add a key. do something like this (please, search for information, I forgot)

wget -O - https://qgis.org/downloads/qgis-2017.gpg.key | gpg --import
gpg --fingerprint CAEB3DC3BDF7FB45
gpg --export --armor CAEB3DC3BDF7FB45 | sudo apt-key add -

• update

sudo apt update

sudo apt install qgis

Installation of the DEMto3D extension

To produce directly STL from raster files, install the DEMto3D plugin. From the plugin installation dialog: DEMto3D is the first tool that links GIS (Geographic Information System) and 3D printing. DEMto3D allows export DEM to STL format ready to 3D printing. [..] Author: Francisco Javier Venceslá Simón

As of April 2018, this plugin only worked with QGIS 2x "Palma" (which is the official stable long term version). As of Jan 2019, DEMto3D version 3.1 works with QGIS 3.x

Installation steps:

• Menu Plugins
• In the popup, select Not installed in the left panel
• Select DEMto3D

Alternative installation steps you can try if the above does not work, e.g. for QGIS 3x on April 2018

• Get the zip file: https://plugins.qgis.org/plugins/DEMto3D/
• Menu Plugins
• Install from Zip

Trouble

As of July 2018, this plugin did not work with the most recent QGIS version. E.g. on April 5 2018 we got:

Couldn't load plugin 'DEMto3D' due to an error when calling its classFactory() method 
ModuleNotFoundError: No module named 'PyQt4' 

The easiest solution was to install the older stable "Palma" QGIS version. That is what I did :)

As of Jan 2019 DEMto3D works just fine with QGIS 3.4.5 Madeira (tested with windows 10 and Ubuntu 18). There is still a display problem though...

Installation of GDAL only under Ubuntu 16x

We rather recommend installing just QGIS as describe above, but the following is more practical for people who just need a few command line options for manipulating terrain data files.

To install a complete GDAL package, follow these instructions as summarized below:

Requirements, some python extension and the GIS PPA

• See above for GIS PPA

The following will install GDAL, including Python bindings

 sudo apt-get install gdal-bin
 sudo apt-get install libgdal-dev
 sudo apt-get install python-gdal
 sudo easy_install gdal

Test if something works. Get a DEM file and see if you can extract information. Examples, using files I extracted from various archives. Each should work and provide about 20-30 lines of information (file content, coordinates, etc.)

gdalinfo n45_e006_1arc_v3.bil
gdalinfo n46_e006_1arc_v3.tif
gdalinfo srtm_38_03.asc