Terrain data into SketchUp

Recently, a client wanted some UK 3D terrain models - in Google sketchup format.

I didn't suspect this would be a problem as I know Sketchup well having used it for many years. I knew it had the ability to display some terrain via the Tools > Google Earth menu > Toggle Terrain (as long as you've set the location). But this wasn't detailed enough for my client, the data they wanted to use was Ordnance Survey Panorama.

The second option was to use the DEM import function - but that would have involved converting Panorama data into the USGS DEM format, this would probably work, but sounded overly complicated.

So I had a dig around for SketchUp plug-ins that would import meshes of 3D elevation data, and of course, there was something that did exactly that job - Cloud v8.

So I Downloaded and dropped it into the plugins folder and gave the demo a try.

All I needed now was to get the Panorama data into that format.

Once again, gdal v1.8 to the rescue...

I already had the Panorama data in geotiff format (using gdal to convert). So then it was just a matter of using the gdal XYZ ascii export format and specifying a window (in OS national grid coordinates in this case).

gdal_translate panorama.tif snowdon.txt -projwin 260100 355500 262100 353500 -of XYZ

Then, opening up Sketchup and the cloud plug-in

File > Points Cloud > Import > (pick a point on screen for origin) then select file snowdon.txt, and use the space delimiter.

The triangulation can take a few minutes...

Eh voila...

Digital terrain mapping using free tools and data

There are a variety of public domain terrain data sets and tools available which can be a great resource to digital cartographers, but when I started out looking I struggled to find any information on how to process this data simply. This post is aimed at cartographers who would like to start using this data, but don't really know where to begin. I hope people find it useful.

Background

But first, why? What are we trying to achieve by using digital terrain data in our maps? Well, there is a long and respected tradition of representing elevation in maps in various ways. Contours are the most familiar, but there are also hyspometric tints and hill-shading effects that can be very effectively employed to indicate to the reader the elevation.

See www.reliefshading.com for more examples of how these cartographic techniques.
See wikipedia for the history of relief shading.

These techniques, pioneered before the advent of digital cartography, required laborious manual processes to create the images. Today we can achieve the same using various open-source tools and data - and it is this process that I want to describe. Primarily with reference to creating online maps compatible with google maps etc.

Resources

The approach I describe uses open-source tools which I run on Mac OSX. The tools will run under other operating systems of course, and there shouldn't be any differences based on operating system used.

GDAL
The primary data processing tool is gdal which is a fantastic collection of raster processing routines. It can look a little intimidating to start with, but it's worth getting to grips with because of it's power and speed. Installation for Mac users is best undertaken by going to kyngchaos and getting the gdal complete package. Other OS users should go to gdal.org to download and install, but insure that you get version 1.7 or above.

Mapnik
For the map production I use Mapnik, again it can be quite daunting to master, but pays huge dividends in the long term. Again for Mac users there is a lovely installer available from dbsgeo.com.

Data
There are two primary dataset that are of most use for global map-makers:

• ETOPO1
• SRTM data

ETOPO1 is a global database of elevation and bathimetry, in a single file, suitable for maps up to a scale of about 1:15m. It comes in a couple of flavours, I use the cell-registered ice-surface geotiff version.

SRTM (Shuttle Radar Topography Mission) data isn't global - but close to it, just missing data at the poles. It's an incredible data set, but it does have data holes, on which more later. It comes as over 1400 1 degree tiles, each 1201 x 1201 pixels in dimension.

Both these datasets are digital terrain models of the earth where each pixel contains an elevation in meters rather a color as would be the case with an image. Each pixel has a geographic extent, usually measured in digital degrees. In the case of ETOPO1 each pixel is 1 arc-minute square, the SRTM data is 3 arc seconds square (approximately 90m).

Using GDAL...

So you've got gdal installed, got some SRTM tiles and you're not afraid of the linux command line... let's get started...

In these examples I am going to use SRTM tile N54W004, but of course any should work.

The raw SRTM data is in zipped HGT format which is unique to SRTM data and not very widely supported. Markus Neteler kindly wrote a bash script to convert these files into geotiff. (Windows users, I am sorry but I don't know if there is anything to unpack and convert the tiles as simply for you). However, I found a small problem with this script not dealing correctly with negative elevations, so here is a modified version that corrects that issue. Get the script, use chmod u+rwx to make executable, then enter the following at the terminal:

./srtm_hgt_to_gtiff.sh N54W004.hgt.zip

and N54W004.tif should be created. This is the geotiff data file that we'll be using in the following examples.

Hillshading
Now, at the terminal enter:

gdaldem hillshade

And you'll see all the options listed.

The basic hillshading command can be run very simply...

gdaldem hillshade N54W004.tif hillshaded.tif

and you'll get something like this...



While this obviously works, the image is very contrasty. This is a result of the data being in raw lon/lat coordinates, effectively exaggerating the vertical scale. A quick solution to get something more appealing is to add a scaling parameter to compensate -

gdaldem hillshade N54W004.tif hillshaded.tif -s 111120

(The scaling value is an approximate value to convert degrees to meters, to match the elevation values)

Try that and you should get this:



That may be good enough for you, if not experiment with changing the rendering options - in particular the vertical exaggeration (z factor) and the sun angle (-az), for example...

gdaldem hillshade -z 3 -az 90 N54W004.tif hillshaded.tif

Hypsometric Tinting
To tint the data you will need a color palette that associates each elevation with a color. You can make your own based on the gdaldem format, or grab this one which uses solid bands of color to represent different elevation levels.

Then run..

gdaldem color-relief N54W004.tif colorpalette.txt colortinted.tif

And you should generate this...



a nice color tinted elevation map.

Compositing
Once you have a hillshaded and color tinted image you need to composite them together to create a single shaded and tinted image. This can be done using a variety of image manipulation tools, one approach is to use Photoshop.

The method used to composite the two images together can have a big impact on the final map quality, so it is worth experimenting with different merging techniques. But the method described in the link above using Photoshop 'multiply layers' works well. Imagemagick can also be used if you prefer the command line approach. I will describe how to use Mapnik to composite images in a future post. Until then see http://trac.mapnik.org/wiki/RasterSymbolizer

Either way, you should end up with something like this:



A single hillshaded and tinted relief map


Mosaicing data tiles together
If you are using SRTM data it is unlikely that one tile will cover your whole area of interest, so prior to creating your hillshaded imagery, you will need to mosaic several tiles together to create one large elevation file covering your map area. Fortunately there is a python script installed as part of gdal to do that for you. First, you will need to identify and download the data tiles that you require and put them in a directory. Uncompress and convert them all to geotiff using the script above. Then, in that directory enter:

gdal_merge.py -o mymap.tif *.tif

This combines all the tif files in that directory into a single file (mymap.tif). Clearly, if you have a lot of tiles that you need to mosaic then the output image can be very large. Fortunately, there are options to clip out the area needed based on bounding box coordinates.

Once you have a merged data tile you can then hillshade and tint this merged file as described above.

Projecting
The raw terrain data is provided in a geographic projection, which is unlikely to be suitable for many mapping needs - so you are going to need to reproject the data into something more useful. This time it is gdalwarp to the rescue.

Gdalwarp can accept projection information using the EPSG codes which is probably the simplest way to capture all the projection parameters. This example reprojects the data tile from the raw geographic projection (EPSG:4326) the SRTM data is provided in, into the Google Spherical Mercator projection (EPSG:900913), but you could use any of the supported projections.

A simple reprojection can be executed like this:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:900913 N54W004.tif N54W004_projected.tif

s_srs is the source projection, t_srs the output projection.

Important
If you are projecting a tiled dataset such as the SRTM data, tile by tile, then you need to be aware that using non-cylindrical projections will result in missing data areas when Gdal needs data from adjacent tiles to interpolate from. I am not aware of any solution to this issue other than merging the data first, projecting, and then re-tiling. Cylindrical projections do not suffer from this problem because the output tile is rectangular and tiles can therefore be mosaiced together, after re-projecting, without overlaps.

Resampling
One of the easiest things to get wrong when generating hillshaded images is to use the data at its native resolution and then resize/resample any derived images to fit. However, resampling a hillshaded image doesn't always produce the best visual result. Fortunately this is easily avoided by resampling the terrain data prior to rendering to an hill-shaded image and considering the resolution of the media you're intend to deliver your map on.

Gdalwarp is the tool to use to resample the data to fit your need - usually the SRTM data has too much information for most mapping needs, so we resample down. There are two ways of doing that using the -ts or -tr parameter.

The -ts parameter sets the output data file dimensions so would be used in the form -ts 600 600 if a 600x600 pixel image was required. The -tr parameter sets the resolution in meters of each pixels and would be used in the form -tr 100 100 if 100m square pixels were required. A sampling method needs to be specified in conjuction with these parameters using the -r switch. For height data I use 'cubic' resampling which delivers high quality resampling but it isn't the fastest method.

How to calculate the correct resampling values...
But how do you work out the correct values for these parameters? The best approach I have found is to match the data resolution to the image resolution you require. Image resolution is measured in DPI (dots per inch) and will depend on whether you are outputting to print or screen media. DPI for print media can vary greatly, but I tend to work with a value of 150, but experiment with your own values if you think that is too low or are creating a particularly high resolution image. For screen the DPI is 96 (which is what I will be using in the following examples).

To calculate the tr parameter (meters in a pixel) you use the formula map scale * 0.0254 / dpi (0.0254 is the number of meters in an inch.)

So to get the resampling parameter for a world map at 1:100 million to be used on screen:
100,000,000 * (0.0254/96) = 26,458.33 m/pix

When used in gdal with ETOPO1 data this gives us...
gdalwarp -tr 26458.33 26458.33 -r cubic -s_srs EPSG:4326 -t_srs EPSG:900913 ETOPO1_Ice_c.tif ETOPO1_Ice_c900913_L4.tif

If you are generating maps for use as tiles in online mapping applications you will need to render separate images for each zoom level. Fortunately, the meters/pixel for each zoom levels for the Google maps and Microsoft Bing maps tile structure is documented online.

Terrain data can be rendered successfully at scales greater than the native level, but as a rule of thumb I tend not to go higher than double the resolution (SRTM data will render nice imagery down to zoom level 11)

Combining Projecting and Resampling
Projection and resampling is usually undertaken simultaneously. A complete projecting/resampling command looks like this:
gdalwarp -srcnodata -32768 -dstnodata -32768 -tr 305.7481 305.7481 -r cubic -s_srs EPSG:4326 -t_srs EPSG:900913 N54W004.tif N54W004_zoom9.tif

(the -srcnodata/dstnodata switch tells gdal what value represents null data).

This will create a single data tile, in spherical mercator projection, sampled for use at zoom level 9. This tile can now be shaded and tinted as required, ready for use in your maps.

Issues

SRTM voids and how to fill them
Unfortunately, the SRTM data isn't glitch free. Because of the way the Shuttle collected the data, there are holes. There have been various attempts to fill these holes using a range of techniques and one of the best sources of void-filled data is CGIAR. However, unless you want to pay license fees, you are better to download from source and void-fill yourself.

Gdal has a simple python route to fill the voids - gdal_fillnodata.py.

This program does little more than patch the holes by interpoloating from adjacent existing pixel values, so it's not particularly intelligent, and other more sophisticated solutions exist where you can patch in data from other sources.

There is also a source of improved SRTM data at http://www.viewfinderpanoramas.org/dem3.html

However for most cartographic purposes gdal_fillnodata.py works OK. You will need to batch process all the tiles through this routine like this:

gdal_fillnodata.py -md 100 infile.tif outfile.tif


Black edge pixels after hillshading
There is one frustrating issue with the hillshading routine in gdaldem 1.7, the output shaded image has a 1 pixel black border, presumably because the program can't the calculate the grey value without knowing adjacent pixel values. If you are creating a single map image from a merged set of tiles then this isn't really a problem. But if you want to process tiles and keep the output in a tiled structure (as you would if you are creating a large scale global map set for online use) then you need a solution to this issue. My approach involved writing a simple script that created a virtual database of all the data tiles (using gdalbuildvrt) then grabbing the NSEW coordinates of each tile, extending the coordinates by a small amount in each direction, extracting this buffered tile from the database, then shading the buffered tile, then clipping the image back to original coordinates. Not ideal, but it works.

Update...
This issue has been addressed in Gdal 1.8 with the switch -compute_edges


Other related web resources

Here are some other great resources for processing and mapping with digital elevation data

• http://wiki.openstreetmap.org/wiki/Contours (info on SRTM data processing)
• http://www.gravitystorm.co.uk/shine/archives/2010/01/05/hill-shading-on-opencyclemap-org/