3D Printing Topographical Maps from Space Shuttle Data

In February 2000 the Space Shuttle Endeavor embarked on an 11 day mission known as the Shuttle Radar Tomography Mission (SRTM). Its purpose was to collect RADAR interferometric data to generate a high-resolution digital topographic database of Earth. The resulting data is a collection of files called Digital Elevation Models (DEMs). These files contain elevations in meters above mean sea level below 60 degrees of latitude. The resolution of the data is 1 arc second over North America and 3 arc seconds over North America and (most of) the rest of the Earth’s land masses.

NOTE: As of the fall of 2014 the agency responsible for SRTM data has begun releasing SRTM data for the rest of the world at 1 arc second resolution. The first continent released was Africa with the other regions of the world released in the months following.

The goal of this project is to convert SRTM elevation 3 arc second data into a 3D surface that can be used as the basis of creating Raised Relief Maps using a 3D Printer. The software created for this purpose is called SRTM2STL. The software reads a 3 arc second binary SRTM file containing elevation data and converts it into triangles and writes an output file in Standard Tessellation Language (STL), a common 3D file format, in either binary or text format. Once an STL file is generated it can easily be imported into 3D printer vendor-provided slicing software. The slicing software converts the data to G-Code, a file format used by 3D printers to control printing.

SRTM Data

The SRTM 3 arc second data is encoded as signed 16 bit 2’s compliment binary data. At 3 arc seconds a one degree by one degree square of the earth is broken into 1200 by 1200 data points. SRTM data is stored as 1201 by 1201 points. Each data point is stored as a pair of bytes. These byte pairs form a 2,884,802 byte file (1201 X 1201 * 2 bytes per point). Each top and right edge of the data contains the same data as the bottom and left edges of the adjacent square (1200+1). This can be used as registration information when files are concatenated.

Software

The software to create STL files, SRTM2STL, is written in plain (very plain) vanilla C. The software can be compiled into an executable program on virtually any platform i.e any platform with a C compiler, for example, gcc.

Algorithm

The process of converting the SRTM data to an STL file begins by converting every group of 4 points into two triangles. Each triangle requires 3 triple points x, y, and z to define the triangle’s location in space. Additionally, each triangle needs a normal vector x, y, and z to define its orientation in space. The normal vector should point out of the object being defined according to the Right Hand Rule (described later). The resulting data forms the surface of the relief map.

Besides the relief map surface the program also has to create the four sides and a bottom making the 3D model closed or ‘watertight.’ This completes the 3D model which can then be imported into an STL file viewer for inspection or into slicing software to create the data necessary to 3D print the model.

For SRTM data containing 1201 X 1201 data points, the surface of the area being created (along with the sides and bottom) will contain 5759996 triangles (or facets/polygons).

Filling Holes

For a variety of reasons some of the data in the original SRTM files is missing. Missing points are marked by making the point equal to the largest negative value possible for a 16 bits or -32768. The marked data points can be easily located within the data and, if the option is selected, the holes can be filled in by the software. The software implements a simple linear algorithm for this purpose.

The linear fill algorithm works by finding a hole in the data on a row in the array containing the data. The algorithm then looks for the next data point on the row that isn’t a hole. The slope is then computed between the two end points. The elevation is computed for each hole based on the slope and, the holes are filled in with the computed value.

The linear fill algorithm has limitations. For example, when a lot of data is missing, the program creates long linear fill lines (that can look bad). Usually the data is missing because of the scattering of the RADAR signal along rough irregular surface areas. Long linear filled lines don’t usually match the surface very well in these instances and therefore often look bad. A better algorithm to deal with holes is one place where improvements can be made to this software.

There is another problem that you can encounter with the linear fill algorithm (as implemented in the initial release of the software). This is when the hole starts at the beginning or the end of a row (an edge). In order to correct the problems encountered in this situation, the user should extract the data away from holes that appear along the edge. This is yet another place where there is room for improvement in the original algorithm.

Another way to avoid the problem with holes is to download corrected files. Much work has been done by the SRTM team to fill in holes using a variety of techniques. A link to the corrected files can be found in the references.

Command-line Options

The program uses command line options to control its behavior. The syntax of the command line is as follows:

srtm2stl input output [Options]

Options
/A – Amplification – Factor to multiply the heights by. This is applied before the bias is added/subtracted.
/B – Bias – Must be an integer. Negative to lower the map (thin).
/C – Coordinate subset – Latitude, Longitude, Latitude, Longitude (decimal minutes).
/F – Fix holes (Linear Fill).
/N – Change the Solid Name in the STL ASCII files.
/H+ means print the array coordinates of hole (missing data). S|s – Subset – Start_Row, Start_Column, Stop_Row, Stop_Column
/T+ means output ASCII format
/T- means output binary format.
V – Turn on verbose output

Creating Polygons

Each polygon that forms the facets of the 3D surface created by SRTM2STL is composed of 3 three-dimensional coordinates, X, Y, and Z. The X and Y coordinates represent longitude and latitude while the Z coordinate represents height of the point above the geoid or roughly, mean sea level. To create the top surface of the map we used 4 adjacent points. For example, if we are starting in the upper left hand coordinates of the map we will use the first two points on the first row and the first two points on the second row. The first polygon is created from the first point on row 1 and the two points on row two. The second polygon is created from the second point on the first row and the two points on the second row. This implies that the two polygons share a common edge, and they do.

Here is some sample data produced by the program which makes up two polygons. The data representing a single polygon is comprised of three sets of points. Each point is made up of three coordinates X, Y, and Z. The first polygon in this example is made up of these three points:

Point 1: 0.000000, 0.000000, 223.000000
Point 2: 0.000000, 92.408531, 219.000000
Point 3: 68.984360, 0.000000, 182.000000

The second polygon is made up of the following three points (X, Y, and Z):

Point 1: 68.984360, 92.408531, 180.000000
Point 2: 0.000000, 92.408531, 219.000000
Point 3: 68.984360, 0.000000, 182.000000

So, the top surface is created from all of the polygons generated by the program based on the options selected. The elevation data is extracted from the SRTM data and stored in a 2D array in the program’s memory. The distances and therefore the coordinates, of the latitude and longitude data are created by multiplying the integer indexes of the array rows and columns by the distance per point in each direction. The distances for longitude (north and south) are different from latitude (east and west) because the latitude distance gets smaller as you go from the equator to the poles. The distance for the Latitude (east and west) is calculated for the bottom of the 1 degree by 1 degree square. It would be more accurate to re-calculate this value as you traverse the array so that the ‘taper’ that exists naturally is represented accurately. Practically speaking, this distance difference is mostly unnoticeable. It could be added in the future as a program option.

After the top surface is created the program closes the sides and the bottom by creating three-point polygon pairs as can be seen in the figure below.

Calculations

One of the things required in the output STL file is a normal to the surface of each polygon. This vector points out of the object being created, i.e. the vector indicates the direction of the outside of the object.

The normal is computed using the right-hand rule. Graphically this looks like this:

So, given that as you create the polygons you know where outside and inside is, choose a vertex of the polygon. According to the picture above, the vectors a and b will correspond to the index and middle finger respectively and the vertex is the point where they begin. So, the cross-product will compute the normal vector which is represented by the thumb in the picture above. The thumb, after this calculation, will point ‘out’ of the object.

The calculation of the Cross Product assumes that the vertex coordinates are (0,0,0). In order to make the calculation from three sets of coordinates the vectors a and b are calculated by subtracting the component coordinates of the vertex from the component coordinates of the end points i.e. the coordinates of the tips of the fingers in the picture above. We will use these normalized coordinates for the calculation of the Cross Product.

Given three sets of coordinates that define a polygon as r, s and t where r is the vertex and, given that the right-hand rule is followed in that r->s represents vector a and r->t represents the vector b, then the calculation of the vectors a and b is as follows:

ax = sx – rx ay = sy – ry az = sz – rz

bx = tx – rx by = ty – ry bz = tz – rz

The Cross Product, c, is calculated from a and b as:

cx = (ay * bz) – (az * by) cy = (ax * bz) – (az * bx) cz = (ax * by) – (ay * bx)

The Cross Product is written to the output STL file along with the coordinates of the ends of the polygon.

File Format

There are two output file formats for STL files: text and binary. The text file format represents the data as an ASCII file of text. Here is an example of what an STL data point looks like as text:

solid some_name

facet normal 61.678 -65.695 -675.324

outer loop

vertex 65708.836 30.839 999.000
vertex 65708.836 61.678 1002.000
vertex 65730.734 30.839 1001.000

endloop

endfacet

endsolid

In the example above there would of course be more than one facet (polygon). The three numbers following facet normal are the coordinates of the normal to the facet x, y, and z. Each vertex is a coordinate of one of three points that represents a corner of the facet and therefore each has an x, y, and a z component.

Binary STL file format is much more compact than text format. It is also harder to understand if you don’t ‘do’ binary numbers. In binary format each of the numbers is stored as a single-precision 32 bit floating point number. The file begins with a fixed 80 byte header which can contain anything but usually contains a short description of the file. Next is a single 32 bit integer that holds the value of the number of triangles contained in the file. Following this are four sets of three floating point numbers. The first set of floats is the normal to the surface x, y and z. After that comes the three points that describe the corners of the triangle and each point has three 32 bit floating point numbers x, y and z.

Conclusion

The source code for SRTM2STL is available on GitHub. It has been released under the GNU license. My hope is that others will modify and extend the software. I also hope that the raised relief maps created with it will be used to educate as well as fascinate.

Below are some screencaptures from the 3D viewing software alongside the 3D printed pieces. The better the resolution of the printer, the better the output matches the model.