Making shaded relief from digital elevation models (DEMs) in QGIS: a British Columbia perspective

Who would not agree with Tom Patterson, creator of the fabulous shadedrelief.com website, when he said, “There is no more important component of a map than the shaded relief.”

But the topic of creating your own shaded relief from a DEM is rather complex, so I’ve made a few assumptions in this how-to. I’m assuming:

  • you know your way around QGIS (I am using QGIS 2.18 for this post)
  • you’re familiar with the idea of projections, and re-projecting raster data
  • you know that each projection/datum combination has an EPSG number, which is a convenient way to refer to it. For example
    • Lat/Long/WGS84 = 4326
    • UTM Zone 9N/WGS84 = 32609
    • BC Albers = 3005

If that’s the case, there are really only three things you need to know in order to made your own hillshades:  where to get the data, how to transform it, and what pitfalls to avoid.

Why even make your own hillshade?

Usually when I want to include shaded relief in a map that is in British Columbia, the first place I will turn is the WMS service at openmaps.gov.bc.ca: http://openmaps.gov.bc.ca/imagex/ecw_wms.dll?

001_WMS hillshade example

The 315 black-and-white hillshade from openmaps.gov.bc.ca. The black-and-white hillshade with the sun at an azimuth of 315° is typically the most useful of the layers on this WMS server.

It’s got some limitations:

  • there are only two sun azimuths to select from: 315° (northwest) and 225° (southwest). More flexibility would be good, because each landscape seems to have a different ideal azimuth to bring out the landforms that you want to bring out. More about this down below.
  • there’s no height exaggeration possible
  • occasionally there are odd artifacts, like straight lines running across the hillshade

On the other hand, you can adjust its brightness and contrast, and use the Multiply blend mode, which means you can do some nice things with it.

002_WMS hillshade example with overlay

When you set the blend mode of a hillshade to MULTIPLY, it shows through upper layers.

But, if you want to adjust the azimuth or exaggerate height, you’ll need to find a DEM and make your own hillshade.  It’s well-established (though not well-explained) that the human eye needs to see light from above, preferably from above and to the left. If your map is, say,  south-up, you will need a hillshade where the light comes from the southeast (which will be in the upper left).

004_WMS hillshade example rotated 180

Now south is at the top (map is rotated 180°) and the shaded relief appears flat or inverted. For this map you would need a hillshade where the sun’s azimuth is southeast, or 135°.

005_DEM hillshade rotated 180

Here the map is still rotated 180°, but the hillshade has been manufactured with an azimuth of 135°. The mountains look like mountains again.

006_DEM hillshade rotated 180 3x

A 3x vertical exaggeration of the terrain emphasizes the detail in the valley bottoms.

The overall process

007a_flow chart

Let’s talk about the three principles.

  1. The hillshading algorithms require a DEM in a metric projection. That means that DEMs projected in degrees won’t work: you have to re-project them first. Unfortunately, just about all DEMS come projected in degrees.
  2. The scale of your final map determines what sort of cell size you want in your re-projected DEM. A DEM with 10m cells is far too detailed for a map at 1:500,000, and the file would be enormous. On the other hand, a cell size of 500m would make a very coarse hillshade at 1:500,000. As a general guideline, divide the denominator of your scale by 5000 to get roughly what cell size you want. So if your map is 1:100,000, you’d be looking for a DEM with (roughly) 20m cells.
  3. It is an option to style any DEM as “Hillshade” (other options are singleband grey, multiband colour, paletted, etc.) but the hillshading algorithm in the toolbox (Processing>Toolbox) produces a better result.

Acquiring DEM data

The first thing you need to do is pick your resolution. Because DEMs usually come projected in 4326, DEM resolution is typically expressed in arc-seconds, or seconds of latitude. Because degrees of latitude in BC are bigger than degrees of longitude, these cells are not square. They are upright rectangles.

What you want to know is how these non-square cells measured in arc-seconds will convert to square cells measured in metres. Here’s a handy chart.

resolution degrees resolution metres pixel size in degrees source limitations recommended scale Notes
1/3 “ 6 9.26E-05 ftp://rockyftp.cr.usgs.gov/vdelivery/Datasets/Staged/Elevation/13/ArcGrid/ USA only ~ 1:30,000 It comes as 1°x 1° tiles. Files have names like “USGS_NED_13_n19w068_ArcGrid.zip,” which would be the 1°x 1° tile northeast of 44°N, 110°W.
1” 17 0.000208333 SRTM 1 for all of North America at http://rmw.recordist.com/
Europe at https://www.eea.europa.eu/data-and-maps/data/eu-demre
North America and Europe ~ 1:100,000 From recordlist these come as 1°x1° tiles in HGT format, each about 25MB. The file N55W128.hgt would be north and east of 55°N, 128°W
3” 45 0.000833333 This coverage, based on SRTM data, is available for the world in two different versions.

worldwide~ 1:300,000v4 data is better, but it comes as 5°x5° tiles. The SRTM 3 comes as 1°x1°tiles in HGT format. The old CDED dems were this resolution.

15”4500.00416667http://www.viewfinderpanoramas.org/dem3.htmlworldwide~ 1:1,500,000  250m and 500m https://hc.box.net/shared/1yidaheouv (Password: ThanksCSI!) worldwide~ 1:1,000,000 and 1:2,000,000re-samplings of the finer resolution data30”1 km0.00833333https://lta.cr.usgs.gov/GTOPO30 Or https://earthexplorer.usgs.gov/worldwide~ 1:3,000,000at this point you should be considering the shaded relief at http://www.naturalearthdata.com/5’10 km https://www.eea.europa.eu/data-and-maps/data/world-digital-elevation-model-etopo5worldwide~ 1:30,000,000 

You’ll notice that I favour the DEMs created from the Shuttle Radar Topography Mission (SRTM). These have a few advantages over the DEMs produced by Natural Resources Canada or the provincial mapping agencies:

  • they are usually more recent (dating to about 2000)
  • they represent actual measurements, as opposed to a grid generated from contour lines
  • they go right across provincial and international boundaries

From here, I’ll demonstrate how this works with an actual example. In this case I want shaded relief for a series of maps that I’m making. All the maps are in the same area, and they range in scale from 1:45,000 to 1:130,000. Looking at the handy chart above, this scale suggests I’m going to want to use SRTM 1 DEM.

007c determining the tiles you want2

Displaying an area in a print composer with grid of lines in CRS 4326 is one way to figure out which tiles you’ll need.

So I go to Recordlist, enter the bounding box and select SRTM1.

008_recordlist selecting DEMs

Downloading DEM tiles from Recordlist

Once they are downloaded and unzipped, I’ll read them into QGIS to confirm that they cover the right area.

Merge and Clip

You want to merge all of the individual DEMs into one, using Raster>Miscellaneous>Merge. (Incidentally, they do not have to be read into QGIS to do this.)

016_merging DEMs

QGIS’s raster merge dialogue

I tend to name the resulting DEM with “_4326” on the end so that later I will know what projection it’s in.

It’s tempting to display as hillshade now, but don’t. The hillshade styling is not meant for DEMs projected in degrees.

040_displaying as hillshade when projected in degrees

Styling the DEM as “hillshade” produces ugly results when the DEM is projected in degrees

If you want to clip the merged DEM, now is the time to do it. Remember that with Raster>Extraction>Clipper you will need to change the QGIS projection to the DEM’s native projection (4326) before you draw the clipping box. Be sure to check, once you come back to your map’s projection, that the clip you made covers your whole print composer.

Re-project and re-sample

Reprojection is the process of giving raster cells coordinates in a new projection. Resampling, on the other hand, is the process of creating new cells based on old cells at a coarser or finer resolution. The two are essentially inseparable, since as you reproject from 4326 to, say, 32609, you will also want to go from the rectangular cells of the degree-projected DEM to nice square cells in the UTM projection.

The first thing to do is to right-click the DEM and choose Save As… so you can see the dimensions of the original DEM cells.

050_Save As dialogue

Note that once you set the CRS for the saved copy to be 32609, you get a suggested resolution of (roughly) 17 x 31m cells. That’s the native cell size of this SRTM DEM at this latitude. Note down the “17” (the smaller dimension)  somewhere, and close this dialogue.

You can reproject and resample using QGIS’s Save As… feature, but you don’t get control over the resampling algorithm used, and the results, once you get to the final hillshade, are ugly.

060_hillshade produced with QGIS Save As

When you reproject using Save As…, the final hillshade sometimes has these weird cross-hatching artifacts in it

Instead, you want to re-project and re-sample with the Warp tool in the toolbox. (Go Processing>Toolbox and search on “warp.”)

084_reprojecting using gdalwarp through QGIS toolbox

The toolbox’s Warp tool

In this dialogue…

  • Source SRS should be 4326
  • Destination SRS should be 32609 (or whatever metric projection you are making your final map in)
  • Output file resolution can be whatever you want, but I get good results with the smaller of the two cell dimensions I saw in the Save As dialogue — in this case, 17.
  • Trial and error with making hillshades has convinced me that the best resampling method to choose here is “bilinear.”
  • I name the resulting DEM with a “_32609” on the end so I will know its projection in the future.

The new 32609 DEM should look the same as the 4326 DEM when read into QGIS, but if you go to the Metadata tab in Layer Properties you’ll see it has a cell size of 17 metres, and quite different pixel dimensions.

For your final hillshade, the one you use in your map, you will want something better that what you get if you just style this DEM as “hillshade.”  But for now, go ahead and style it as “hillshade.” This enables you to play with the sun azimuth and elevation, and the vertical exaggeration, to see what is going to work best for your terrain. The human eye probably wants an azimuth around north-northwest (337.5) but there are a lot of azimuths on either side that will work.

Notice how changing azimuth changes what’s brought out in this piece of terrain:

Make a static hillshade

Now if this hillshade has no strange artifacts, you’re done. But if you want a smoother results, it’s time to take the azimuth, sun elevation and vertical exaggeration you’ve chosen, and go over to the Hillshade tool in the toolbox (go Processing>Toolbox and search on “hillshade”). This will produce a raster hillshade that you just display as singleband grey. It also tend to be about half the file size of the DEM itself.

086 running toolbox hillshade2

The toolbox’s Hillshade tool

And here’s the result

090_hillshade produced with gdalwarp gdaldem hillshade

Look, no artifacts!

Displaying hillshades

Typically you put the hillshade as the bottom layer in your map, and adjust brightness and contrast — because hillshades tend to be a bit darker than you want. Usually it’s brightness up, contrast down. Sometimes you will also make it semi-transparent.

More important, in terms of getting the rest of you map layers to appear to drape over your hillshade, is to set the blend mode of the hillshade to Multiply.

091_adjusting brightness etc

Multiply causes the values (lightness, darkness) in the hillshade to be adopted by layers above it.

So you can get something like this.

095_final result with multiply

Overlaying a hillshade whose blend mode is set to Multiply

That’s it. Using these techniques, you should be able to manufacture hillshades with any azimuth, pretty much all over the world. It gets the most challenging when you are mapping at large scales like 1:20,000.

Finally, if you master this, and you are really in love with shaded relief, you will want to experiment with making shaded relief in Blender.

Advertisements

Tracing the Great Southern and Western

The “Ring of Kerry” is a 180 km tourist route around one of the rugged peninsulas on the west coast of Ireland. Although its moniker probably arose in a 20th century tourism promotion, the  route itself is one that travellers in the neighbourhood of the Irish town of Killarney have been encouraged to go round since the 1800s.

Here is Samuel Carter Hall, in his 1858 traveller’s guide, A Week At Killarney, expanding in full Victorian eloquence on why one might make this journey:

We shall ask the reader to accompany us to the wild sea-coast of the South-west, and the Tourist to follow us into a district where the graceful beauties of Killarney may be contrasted with the wild grandeur of scenery certainly unsurpassed in Ireland. That district is now visited by a large number of those who visit Killarney; and one of our special objects in our latest tour–in 1858– was to describe the routes to it, with the facilities for travelling and accommodation; and at the same time to picture its peculiarities as well as our limited space and opportunities permit us to do.

The district Hall is referring to is the Iveragh peninsula, jutting west into the Atlantic from the area around Killarney. Hall’s tourist is to begin from Killarney and proceed south to the town of Kenmare — which anchors the southeast corner of the peninsula; then out along its southern coast, up around the end, and back along the north coast to Killarney. With horse and jaunting cart the tour took the Victorian traveller two days. (Hall recommended hiring a single set of horses for the entire journey.) Today most visitors drive it in one.

But not all visitors have explored the Iveragh peninsula by horse-drawn car or horseless carriage.  Some 35 years after Hall published his guide, a new option appeared: rail. The West Kerry Branch of the Great Southern and Western Railway began its service in 1893, running along the north coast of the Iveragh from Farranfore to Reenard Point. Farranfore was a connection with the main line running from Tralee to Killarney and on to Cork; Reenard Point was a ferry terminal on the outermost coast where passengers could board a boat and continue out to Valentia Island, a place sometimes promoted with the tantalizing (but geographically incorrect) fact of being the westernmost point in Europe. However it is fair to say that the GS&W was the westernmost railway in Europe.

1902 route map

The West Kerry Branch of the GS&W highlighted in yellow. (Map from Project Gutenberg: “Great_Southern_and_Western_Railway_-_1902_Ireland_routemap_-_Project_Gutenberg_eText_19329.jpg”)

Tourism was only one reason for this thirty-mile rail line out to the end of the world. Valentia Island hosted both famous slate quarries and an active fishing industry, and there was demand for the products of both in Ireland’s big neighbour to the east. No less a symbolic building than the Houses of Parliament in London was floored with Valentia slate.

This branch line ceased operations in 1960, and today most stations and track – which was the Irish Gauge, spaced 5′ 3” apart – have been removed. But, if you’re paying attention, you can still encounter bits and pieces of the railway, which visited the towns of Killorglin, Glenbeigh and Cahersiveen on its way.

 

DSCF0178

In Killorglin a plaque on the sidewalk, near the Aldi supermarket on Iveragh Road, tells the story of this vanished line.

For example, in the village of Glenbeigh, if you cross the River Behy on the old bridge and go just a short distance left, you find the road cutting through the old railgrade, still some ten feet high.  The Gleensk Viaduct looms above your head as the road does a tight turn to contour through a small drainage on the steep-sided hill.  At Cahersiveen the old rail bridge stands just upstream as you cross the river Ferta.

 

DSCF0042

On the way to Cahersiveen, the road snakes high above the Atlantic along the side of Drung Hill, and here a roadside pullout features an interpretive sign that details the history of the GS&W.

DSCF0043

At the same location: the mouths of two railway tunnels are above the road but below the electric line running across the slope

In many ways though the most remarkable remains of the railway are the mute raised railbeds that run through farmers’ fields.

 

DSCF0182

Do the cows know that a railway used to run along that strange, straight mound that cuts through their field? (Muingaphuca townland, near Killorglin)

As well, the railway is detectable in small kinks in present-day roads, like the one in the Caragh Lake Road 200m from where it meets the Ring of Kerry. Here the line crossed over the road, which jogs briefly to pass under a viaduct that is now completely gone.

There’s a remarkable web resource for tracing this old route. The Historic Environment Viewer of the Irish government’s Department of Arts, Heritage, Regional, Rural and Gaeltacht Affairs features a detailed basemap from Ordnance Survey Ireland (OSI) which shows the current status in superb detail at large scales.

 

Screenshot from 2017-08-23 15-54-25

OSI’s superbly detailed and up-to-date map of the area between the Caragh River (right) and the village of Glenbeigh (left). Portions of the old rail line are shown.

But here’s the magic: you can change the base map to the Cassini 6″ mapping. This is mapping at six inches to the mile (or 1:10,560) from the late nineteenth or early twentieth century, perfect for seeing where the railway ran. It turns out that, as one drives the Ring of Kerry from Glenbeigh to Cahersiveen, the old rail grade is never far away, sometimes running on the right, sometimes the left, side of the road.

 

Screenshot from 2017-08-23 10-39-19

The same region, but now seen against the Cassini 6″ basemap, c. 1900

As well, there are two other historic map layers here. The Historic 6″ dates from between 1829 and 1841; the Historic 25″ (twenty-five inches to the mile, or 1:2534!) from the end of the nineteenth century.

There was a hope at one time that trans-Atlantic vessels would depart from Valentia Island, passengers preferring to go as far towards North America as possible by rail before boarding ship. It never came to pass, but today the old rail line may perhaps have another life. Between Glenbeigh and Reenard Point the rail bed is presently the locus of debate about whether it will be turned into a cycleway. The conflict focuses on how local landowners will be compensated. If this can be worked out, yet a fourth mode of transport will be added to the choices travellers have had for exploring that peninsula west of Killarney.

 

 

Downloading OpenStreetMap data through Trimble

OpenStreetMap data is often the best large-scale data you can find in regions of the world where governments are not yet distributing free, open geodata.

There are several ways to download OSM data (QGIS plugin, Geofabrik, direct from OSM), but if the format you prefer is shapefile, and you have a specific area you’re interested in, the old Weogeo service  may be the easiest and the best way.

Weogeo seems to have been bought or otherwise taken over by the for-profit company Trimble. I appreciate Trimble hosting this service and keeping it alive, but since Trimble integrated it  into their pre-existing sales system, you have to go through a number of odd steps to get your free data. Since it’s free, open data, you hope for one of those good-feeling interactions that suggest the sharing-without-strings that OSM represents. Instead be prepared here for a less comfortable, more corporate, interaction where they’ll want to you to go to their “data marketplace,” make an account, “add to cart,” and so on. But it works really well.

Two important notes up front:

  1. If you use the Firefox plug-in “Privacy Badger,” it breaks this site. You have to disable it for trimbledata.com.
  2. You won’t get your data immediately. Depending on the order size you may have to wait up to 24 hours to receive the email saying that you can download it.

Here’s the procedure:

Go to http://market.weogeo.com/datasets/osm-openstreetmap-planet.html

Click Shop OpenStreetMap data now. You are now at the Trimble Data Marketplace.

Sign In (upper right corner). (Registration is free, and necessary.) Your name should now appear in the upper right hand corner, and below it, over the northwest portion of the map, you should see a summary of your current order, with headings for Region, Layers, Datum-Projection and File Format.

Initially the map shows the entire world (in a ghastly pseudo-mercator projection, but that does go hand-in-hand with OSM). Region is “Entire map,” Layers are “26/26” (26 of 26 available layers), Datum-Projection is “Lat/Long-WGS84 (Native),” File Format is “ESRI Shape,” Estimated size is “1.36 TB” and cost is “Free.Trimble1

Don’t be fooled by the 1.36 TB. That’s the size of the entire world database for OSM. Trimble will restrict you to a 5 GB limit in a single order.

On the map you can now pan and zoom to the region you’re interested in. Let’s say you go to the area around Budapest, Hungary.

Notice that even though you are zoomed in, the Region is still “Entire Map” (meaning the whole world) and Estimated Size is still 1.36 TB. To narrow down your data area you have to either upload a KML file with a polygon of your area of interest, or draw a polygon on the screen.Trimble2

To draw a polygon, click the pencil icon next to Region, click Draw, and begin placing points around your polygon. Clicking on the initial point closes the polygon. Now Region and Estimated Size change to something more reasonable.

You may not want all 26 layers, and you can click the pencil icon next to Layers to select the subset you want. The 26 layers are explained in depth at the OSM wiki page on Map Features, but for simplicity I’ll just list the 26 categories here, with links to the OSM wiki page:

Note that these categories are often subdivided into separate shapefiles for point, line and polygon features.

In this case let’s say I want only Waterway, Highway, Railway. I deselect all, select these three, and click the “X” to close the Layers list. Layers has now changed to “3/26” and Estimated Size has now dropped to 381 MB.Trimble3

Then click Order. Accept the Content License and click Add To Cart. Click View Cart.

Everything should still read “Free,” so click Checkout. Go through the [annoying] Address Information and Select Payment steps (“No Payment Required”) and then finally Place Order.

The next step is that you receive a series of emails acknowledging your order, telling you that your order is being prepared, and finally that your order is ready for download. It can be almost immediate for small orders, or up to 24 hours for large ones.

When you do download your data it will be in a ZIP file called “weogeo_<order number>.zip.”

Mapping Gottfried Merzbacher’s “The Central Tian-Shan Mountains 1902-1903”

In 1905 Gottfried Merzbacher, a geologist, published an account of two years exploring the Tian-Shan.  The Central Tian-Shan Mountains, 1902-1903 (John Murray, London) is an engaging read, partly because of Merzbacher’s insights about the landscape, and descriptions of the people he meets, but also because it is an account of his efforts to solve the puzzle of where the prominent mountain of Khan Tengri actually was. This 7000m peak was readily visible from lowlands outside the Tian-Shan, but previous mappers had merely been able to theorize about how one might get to its base.

Merzbacher’s book originally came with a map, but, as is so commonly the case, the map was not properly scanned for the digital copies of the book that one can download today from archive.org. You get a fragment like this:

merzbacher-map-fragment-2-colour

It’s a shame. Hopefully someone who owns the book will post a good scan of the complete map one day.

What fascinates me about these exploration accounts is the question How would we follow his route today? To that end I’ve begun mapping Merzbacher’s travels as described in The Central Tian-Shan Mountains. So far I’ve mapped Chapter 1.

merzbacher_chapter_01

(Download as PDF)

The key decision in this kind of mapping is what base map to use. Use a period base map (something from around 1900) and you may not be able to follow it today: place names change, roads move, and mountains get mapped more precisely. (See my previous post on the yet-moving location of Khan Tengri.) Use a modern map and you lose the sense of travel in an era of inaccurate mapping.

I tend more to the modern map. To me, the question Where did he go? means Where was his route relative to modern landmarks? But many modern base maps are unsuitable. Most online mapping services do not label rivers or mountain ranges. A modern road map (like Freytag+Berndt’s Central Asia) is too small-scale (1:1,750,000) to see what’s going on. So in this case I settled for a 1980s aeronautical chart (Tactical Pilotage Chart F-6C) which gives shaded relief at a reasonable scale (1:500,000) plus the added benefit of named rivers and mountain ranges. It ain’t that pretty though.

Representing Timberline

When I taught map reading in the U.S., there was a piece of folklore that we used to tell our students: that the green areas on the USGS topographic map indicated there was a sufficiently dense forest there that you could hide a platoon of soldiers (about 40 people) per acre. It was a good way to explain why small clumps of trees (or “tree islands”) didn’t appear on the map.

This story lives on on the Internet, but there’s no evidence that it’s really true. I like the implied subtext — cartographers producing the maps for military officers involved in some kind of domestic war, and needing to know where they could hide their men from aircraft –- but it doesn’t take much reflection to realize that the U.S. Geological Survey never could have visited all those places, looked at the tree cover, and decided where you could or couldn’t hide 40 guys. Or even how to divide it up into suitable one acre blocks.

In Canada, the green area on the topographic maps has a specific definition. According to Natural Resources Canada it’s “An area at least 35 per cent covered by perennial vegetation of a minimum height of 2 m.” And they probably estimate that 35% coverage from air photos.

NRCan topo_sm

So, when you’re ascending a mountain and the green ends on the map, is that treeline, or is it timberline? I’d always thought these were just variant terms for the same thing, but then I read Jim Pojar’s book Alpine Plants of British Columbia. In the Introduction to this photo-rich handbook of all the plants you’re likely to see up there, I learned that treeline and timberline are different:

The term treeline designates the upper limit of the occurrence of tree species, regardless of their stature, whereas timberline refers to the upper limit of forest, of continuous cover of upright trees 3 m or more in height.

So timberline, being where the solid forest ends, is the end of green on our classic NRCan topos. Treeline is the last little, twisted, stunted tree.

Neither, of course, is really a line. As map scale decreases and you zoom in, the timberline becomes impossibly complex, and has to be generalized somehow. And no map, I think it’s fair to say, tries to represent treeline, since this would somehow be defined by many isolated clumps of krummholtz (“twisted-wood”, the bonzai-like tree clumps also known affectionately as shintangle) that you see after ascending past timberline. And just to make things a bit more complicated, as the climate changes it has become easy (in northern BC at least) to find areas above treeline where dozens of tiny seedlings are coming up and now surviving. How big would they have to be before one moves treeline?

Timberline however is a very important landmark for hikers and skiers, and how to represent it is a question that comes up frequently in topographic mapping. Of course using using generalized green and white is not the only option. I first learned this when I was bushwhacking across a 1:50,000 scale “provisional” series, black-and-white topographic map in northern BC. On these there is actually a black line that snakes across the elevation contours. It has “W” on one side of it (for wooded) and “C” on the other (for clear). Little pairs of tick marks pointed into the forested side. It was hard enough to read that I took a pencil crayon and shaded in some green on the W side.

NRCan BW topo_sm

There’s also the solution that National Geographic used years ago in mapping northern Canada; the Northern Limit of Wooded Country is represented by a line of tiny tree symbols.

NatGeo_sm

Sometimes using green is just out of the question. If you want to use a range of colours to represent elevation (hyposgraphic tinting), having green forest is going to be quite confusing. In these cases I have I tried a technique of having a dashed green line at timberline, bordered by some fill on the downhill side, fill that quickly fades out. You can do this in QGIS by using shapeburst fill, shading to a set distance of a few millimetres.

McDonell_Lake_detail_sm

I have also played around with not marking timberline at all, and just putting a note beside the trail at the point when one would clear the trees.

OpalRidge_sm

Another option available to you if you have access to landcover data, is to give map-readers more information about how the forest makes that complicated transition to grassy tundra. In this case I used the Land Cover, circa 2000-Vector data available at Geogratis, and assigned progressively lighter colours to “coniferous dense” (which captures the main forest of spruce and fir), “coniferous open” and “broadleaf open” (which captures willow).

Seaton_sm

If timberline is something you just want to suggest, but don’t really need to accurately show, a method that can still look good is to style the digital elevation model in a series of fading greens. Have it become completely white at the elevation where timberline is typically encountered in the area. You get the effect of the “naked” mountains rising above forested slopes without introducing the complexity of avalanche tracks and the differences in where trees grow between north- and south-facing slopes.

RedRose_sm

Timberline in this area tends to be at 1500m

Sometimes the highest elevations on the map just touch timberline, in which case there can be a large zone where small meadows alternate with clumps of diminutive trees, an ecology often called parkland. I’ve tried representing this by scattering a host of little tree shapes, sort of trying to show how trees are still there, but not connected to the “mainland” of the forest.

QuickHills_sm

And this is just scratching the surface. The more you study timberline, the more you realize the folly of trying to show it accurately, yet the importance of indicating where we see it!

 

 

 

A new map of the Bulkley Valley

When I attended the ICA Mountain Cartography Workshop in Banff last spring, I was encouraged by a number of people there to consider printing and selling maps of my local area. Most of this last winter has been taken up by the production of this map, which was printed in April and is now available at Interior Stationary in Smithers, BC.Bulkley_Valley_poster_map_wordpressThere are a number of goals coming together here. From a cartography perspective, the idea was a map that would feature rich shaded relief generated in Blender. From a community perspective, the idea was to produce a map that represents the amazing topography of this area, that draws you in and makes you want to go exploring.

Every map is a story, and maps, being deliberately authored, contain bias. In this case, I wanted the story to be about the terrain — the mountains, rives and lakes — and my bias was to emphasize the romance of landscape at the expense of the technical ways we divide it up. So missing from this map are all those lines that chop up a landscape: town boundaries, ski area boundaries, recreation site boundaries, private land designations, First Nations reserves and Regional District subdivisions. Even provincial park boundaries got the axe, because while on the one hand a park indicates a piece of the landscape to be preserved from development, on the other hand is the implication that outside the boundary anything goes. Park boundaries on maps suggest that only within the line is the higher quality landscape. I wanted to avoid that suggestion.

As a concession to practicality, most roads are named and you can use it as a road map. But it’s not one of those maps where every road in the provincial database is present. I left out the vast majority of forestry roads, including only those mainline roads which are commonly driven. And here too bias plays a part, because although I’ve lived in the valley for 20 years, my notion about what roads are “commonly driven” will differ from someone else’s.

Bulkley_Valley_poster_map_detail2_wordpressAs another concession I put in some point markers for those small provincial parks, or recreation sites, where a person can camp or picnic. I did this because, like roads, these are important landmarks for people.

Trails are not shown on the map, partly because at this scale, 1:150,000, you can not use it as a trail map, but also because I didn’t want to suggest that trails are an important interpreter of the landscape. When you put trails on a map they stand out as a travel network, perceived in the mind not unlike a road network. We think that where trails go must be better than where trails don’t go. I didn’t want to bias the viewer in this way.

The colouring of the map is derived from land cover data, which is essentially a satellite or aerial photo that someone has looked at and classified into zones of different stuff. Coniferous and deciduous and mixed forest are all identified, as are ice and snow, bare rock, built-up areas and shrubland. I’m a big proponent of the idea that a map should faithfully give you the colours of a landscape before you ever arrive in the area, so the colours I chose were keyed to the summer landscape here — and consequently there’s a whole second possible project of a winter map of the valley. I think the colours work pretty well, although it’s something of an idealized summer landscape. For example, there are no lingering snowpatches such as one would actually see. It is quite clear how coniferous forests dominate upper slopes while the valley bottoms are full of cottonwood and aspen. You see how south facing slopes are different from north-facing. Given that it can’t account for how light changes throughout the day, it’s a pretty good representation of what we see.

Bulkley Valley detail 300dpiNaming is an interesting area in map-making. If you publish an article about how such and such peak should be named after your uncle Fred, people will agree with you or not, but you’ll be known as having proposed that name. However if a map-maker puts a name on a mountain or creek, no one really questions it. Maps are seen as authoritative.

But names can be seen as impositions on the landscape like other lines. Does it improve this river that we call it the Bulkley River? Would it be a better river if it was labelled the Wedzenkwe? Or is it best not labelled at all? I made a decision to go with the naming of features that reflects the British Columbia geographic names database — in other words, the names that appear on most maps and with which most map-readers are familiar — but I’m not entirely comfortable with this. For example, I’d like to see a version of this map with Wet’suwet’en and Gitxsan names on as many features as possible.

We have a lot of unnamed features in this area — or, I should say, features without official names. I hope that this map will act as a lightning rod for name information, that people will give me local names and I can add them on to future editions, and thus substantiate them.

Buffering (or haloing) text over complex backgrounds using the Screen blend mode

This is a common problem when you are working over a complicated, multi-colour background. Such backgrounds tend to swallow text!

Here’s an example, where text lies over a background with a lot of variability in value (darkness).text over complex background

Here’s the same image in greyscale, just to show the pixel values. Out of a total range of 0 to 255, the pixels here range from 6 (almost black) to 240 (almost white). Black text over light portions reads easily, but it’s hard to read the letters, M, N and T because of the way they overlap dark areas.value map of background

The classic, or perhaps I should say default solution, is to buffer the text in white. This works fine, in that you can read the text more easily, but it’s a crude solution and it’s a bit loud. It lacks subtlety. white 1mm buffer

A somewhat more elegant solution is to have the buffer match the background colour. In this example the text stretches across many background colours, so no matter what colour you select it’ll be a compromise. Still,  it looks better than white. The problem is that across the map as a whole each piece of text will probably get a different buffer colour, depending on what is dominant in its background.1mm colour buffer

What we really want is a way to pick up the background colour for each pixel, and colour the buffer that way. We can actually approximate this by using the blend mode of Screen.

Screen is the opposite of Multiply. Where Multiply transfers dark values onto another layer, Screen uses dark values to lighten another layer.

In this case, you set your buffer to something fairly dark, like 80% grey. Then set the blend mode of the buffer to Screen. You get this.screen 80 percent grey

Here the buffer colour echoes the underlying colour. The buffer around the T is orange-y, while the buffer around the O is grey. The buffer for the M uses both colours.

How would you do this in QGIS and Inkscape?

For text labels displayed in QGIS, this impressive piece of software can do this automatically. On the Buffering tab of the Labels dialogue, check Draw text buffer and specify an 80% grey for the buffer fill. (Note that 80% grey is RGB 51/51/51, HSV 0/0/20, or HTML #333333.)

Then set the blend mode to Screen. There are a couple things to note. One is that the sample shown will look ghastly; you can ignore that. The other is that the blend mode of the label is different from the blend mode of the buffer. The label itself probably has the default blend mode of Normal, and this can be checked on the Text tab.QGIS_buffer_screen

Here’s the result: Tiltusha example

In Inkscape, where you might be doing more complex text, text set on curved paths with adjusted kerning, you can also do this. Blend modes in Inkscape are applied to an entire layer, so you can’t give one piece of text a different blend mode than another. But by introducing a special layer that has a blend mode of Screen just beneath your text, you can put the buffers there.

The procedure goes something like this:

  1. Create a new layer below your text layer and set its blend mode to Screen.
  2. On the text layer, select your text and duplicate it with Ctrl-D. duplicate
  3.  Turn on Stroke for the duplicated text. (Text normally has fill but no stroke) and set the stroke colour to 80% grey. (You can do both of these steps simultaneously with a shift-click on the 80% grey swatch in the Palette at the bottom of the screen). Set its width as you see fit (1mm in this example). It should look pretty ghastly.buffer stroke on
  4. Press Shift-PgDown to send that duplicate piece of text to the layer below, the layer with blend mode of Screen. It should now look pretty nice. buffer moved to screen

I should point out one disadvantage of this two-layer approach: you can’t group the text with its buffer, a practice I normally do because it makes it easier to move labels around later.

Shaded relief with BlenderGIS, part 2

With the DEM read in, we’re now ready to adjust the many parameters in Blender which determine what we’ll get. Here again is the overall procedure (mentioned in the previous post) . We’re now at step 5.

  1. Select the Cycles Renderer
  2. Read in the DEM as DEM
  3. Create new material
  4. Adjust Z scaling
  5. Create and adjust the georef camera
  6. Correct final pixel dimensions to match the input DEM
  7. Adjust lamp location and properties
  8. Test render
  9. Adjust final parameters
  10. Render

Create and adjust the georef camera

With the plane selected, go to the the (new) GIS tab on the left, and, under BlenderGIS, click Georender.

geoRender.jpg

This creates a camera at an appropriate distance above the plane. It will appear in the tree of objects (upper right) as “Georef cam”. It has been placed high enough (Z value of 16791) that it is above the highest point on the plane surface. Note that it is important to have set your Z scaling before creating the Georef cam.

georef cam set up

If you look at the camera’s Object Data>Lens, you’ll see that BlenderGIS has automatically set it to Orthographic. It’s a good idea to note the camera’s clipping Start and End, because you may have to adjust these later.

geoRefCamDetails

As well, if you check out the Scene’s Scene panel, you will see that the current camera has been set to this new Georef cam.

Correct final pixel dimensions to match the input DEM

If you go to the Render tab, under Dimensions>Resolution you can see the pixel dimensions of the final image as suggested by BlenderGIS. In this case, as I noted before, one cell has been dropped in each dimension, so it is 5029 x 3276.

resolutionInitial

I correct these to the actual values (5030 x 3277, in this case) so that the dimensions of my final image match the DEM. That way I can pair the TFW world file with the image created by Blender, and read it into QGIS georeferenced.

dimensionsRepaired

While I’m at it, I’ll set the Resolution percentage to 25%. This way I can get a test image without doing a full render. While I’m here on the Render tab, I’ll go to the Light Paths panel and set it to Limited Global IlluminationlimitedGlobalNow we’ll check that all elevations in the DEM are visible to the camera. Pop into Camera Ortho view (Numpad 0 while the mouse is over the 3D View). You are looking for  black faces where the higher parts of the DEM have been cut off. It can look something like this:orthoView1BlenderGIS was updated in December 2014 to place the Georef Cam high enough so that this typically does not happen (Thanks!), but I’ll still discuss what to do in case you do encounter it. This could happen if you, say,

  • create the Georef cam before you increase Z scaling on the DEM.
  • create the Georef cam at a time when the View number in the Subsurf modifier is low (e.g., 6). In this case the displacement is quite generalized; the Georef cam will be placed high enough to clear it, but not the actual high points that pop out when you render at 11.

What’s happening here is that the camera only “sees” objects that fall between the Start and End clip distances we noted a few steps back. The Start and End are measured from the camera outward (or down, in this case).

Specifically, in the image above, I had a camera at 11,706m, with a Clip Start of 0.5 m, and a Clip End of 8011.5 m. If you do the math, you’ll see that this means that the maximum elevation at which it can see objects is 11,705.5 (that’s 11,706 – 0.5). The minimum elevation at which it can see objects is 3694.5 (that’s 11,706 – 8011.5).  Our DEM originally contained elevations from 739 to 2379, but now we’ve scaled it up 5 times, so its minimum is 3695 (= 739 x 5) and its maximum is 11,895 (= 2379 x 5).  The black patches are the elevations that are higher than 11,705.5 m.

If I do a test render with these settings I’ll get the image below; I’ve circled the black patches where the highest peaks were cut off.SRTM_2_3_cutoff_render

In this case we need to raise the camera to just above the max value in the DEM, and then adjust the Clip Start and End distances to encompass the whole DEM. Since the DEM ranges from 3695 to 11,895, let’s put the camera at 11, 896, leave Start at 0.5, and set End to 8201.5.

 Adjust lamp location and properties

This does not differ in any way from what you would do with the Subdivide method. Set the Lamp to Sun, and click Use Nodes. I like a sun size of 1, Cast Shadows checked, and a Strength of 4.

lampSettings

Then change to the Lamp’s Object tab, Transform section, and set its X, Y and Z Rotation. In this case I want an elevation of 40°and an azimuth of 337.5°, so I choose X=0, Y=40, and Z=112.5. Remember that Blender measures the sun azimuth counterclockwise from East. lampSettings2

 Test Render

Hit F12 to see a test render at 25%. It should look grainy but basically correct.

Final adjustments

At this point you can play with:

  • lamp rotation, strength, size
  • Z scaling on the model
  • the Render number under the Subsurf modifier on your plane

When you’re ready to go,

  • On the  Render tab>Dimensions, raise the Resolution from 25% to 100%
  • On the  Render tab>Sampling set the Samples Render number to 200. This high sampling number makes the rendering process take a long time, but it’s worth it. Check out this comparison of the same scene rendered with sampling 200 (left) vs sampling 90 (right):

    Note the grain

    Note the noise in the shadows when sampling only 90 times (right image) as opposed to 200 times (left image).

  • On the Render tab>Output, set the image type to TIFF (if you want to pair it with a pre-existing TFW after you render)
  • Save your session to a .blend file using File>Save As…

Render

Two options here. One is that you can hit F12 and go out for coffee.  When it’s done, hit F3 to save.  Or, you can close Blender and render that .blend file you just saved from the command line. This saves the GUI rendering, and so it takes less time — you might get 10-40% off the regular render time.

Let’s say your  .blend file is called tile_4_4.blend. The command line to render it would be

/path/to/blender -b tile_4_4.blend -o hillshade_4_4.tif -F TIFF -f 1

The switches are:

-b     run in background (no-GUI) mode
-o     the output filename
-f 1    render just one frame.

Pulling the result back into QGIS

Let’s say you have saved the rendered image as hillshade.tif. Rename a copy of the world file you created way back at the beginning to hillshade.tfw.

QGIS should be happy to read in hillshade.tif now as raster data, although it will probably ask for the CRS. (You know the CRS because you specified it when you re-projected the DEM, and when you read it into Blender.)

It works quite well in QGIS to set the Blend Mode of the hillshade to Multiply, and then overlay it with layers representing vegetation or elevation.

The shadows can be thunderingly deep in these hillshades…

QGIS1

…so it can be good to display them with increased brightness and decreased contrast.

QGIS2

 

Shaded relief with BlenderGIS, part 1

[The BlenderGIS add-on has now been improved to the point where generating shaded relief is downright easy. So, I have updated this post (April 2017) to reflect how one uses its latest version. Also the latest version of Blender (2.78c).]

You may have read Daniel Huffman’s great post about created shaded relief using the free, open-source animation software, Blender. As well, Daniel posted a video tutorial for his method, which is well worth watching.

In Daniel’s method, you subdivide a plane to create a grid of finely spaced points. You then apply a Displace modifier to the plane with a DEM, and get some very realistic results.

There are a few things about this method, which I’ll call the Subdivide method, that were difficult for me. For one thing my meager 4GB of RAM precluded subdividing the plane more than about 1500 times. So I went searching around, and found to my surprise that there’s a plug-in available for Blender called BlenderGIS that does much the same thing. But unlike the Subdivide method, BlenderGIS uses a Subsurf modifier. This saves a lot of memory.

Note that distinction: subdivide vs. subsurf. That’s important

As well, BlenderGIS offers a few interesting tools for the GIS user. One is the ability to create a GeoRef camera, correctly configured and positioned to shoot an image of your lit DEM. Another, which I haven’t really explored, is the ability to read in shapefiles to add features to the DEM surface.

In this series of posts I give a tutorial for how to use Blender GIS’s subsurf method to achieve results similar to those obtained with the subdivide method.

I’ll assume that you have watched Daniel’s tutorial, so finding your way around Blender for purposes of generating a shaded relief isn’t too hard for you; and this of course implies, as well, that you are the kind of GIS user who knows what a DEM is. My workflow mostly happens under Ubuntu 14.04 using QGIS and GDAL tools.

We’ll do this in two parts. Part 1 is about installing Blender GIS,  preparing your DEM, reading it in and adjusting the vertical exaggeration. Part 2 is about creating the georef camera, placing the lamp and adjusting some final parameters.

Installation of BlenderGIS

The master instructions are in the BlenderGIS wiki at https://github.com/domlysz/BlenderGIS/wiki/Install-and-usage  However, I will paraphrase them here.

Go to the BlenderGIS site on github and hit the Clone or Download button. You will receive a .zip file, which you can store pretty much anywhere.

Within Blender,  at File>User Preferences, the Addons tab, click Install From File, and point it at the .zip file you just downloaded. (You don’t have to unzip it.)

Once its installed, you still have to do two things:

  1. Check the box next to 3D view: BlenderGIS to enable it
  2. Click Save User Settings

That’s it. BlenderGIS is installed.

Preparing your DEM

I used to use divide my DEM up into square tiles, and feed these one by one through Blender. However BlenderGIS in no way requires this.

Projection and resolution

BlenderGIS claims that it can now read a DEM projected in “WGS84 latlon”, but I have not had very good results with this. I still prefer a DEM that is in a projection that uses metres. Typically for me this would be a UTM projection.

You may also want to re-sample your DEM, so the cell size is appropriate to the scale of the map you are making.  I use this guideline: the DEM resolution should be map’s [metric] scale divided by 10,000. For example, if the map scale will be 1:100,000, I want a DEM with 10 metre cells. For a 1:30,000 scale map, I want a 3 metre DEM. It’s just a general guideline.

Before you re-sample and re-project, however, be sure to float your DEM and create a world file.

Float DEM

Most DEMs come as integer values. Floating it means converting it to a DEM that can has floating point elevation values.  You can do that with a line like this with gdal_translate.

gdal_translate -ot Float32 myDEM.tif myFloatDEM.tif

where the -ot switch specifies output type.

World file

You’ll need the world file when you are done, to georeference your hillshade.

Almost all DEMs come as GeoTiff files, and the corresponding world file extension is .tfw. An easy way to create a TFW World file is with the GDAL command gdalwarp:

gdalwarp -co "TFW=YES" myGeoTiff.tif deleteme.tif

This just makes a second copy of myGeoTiff.tif called deleteme.tif, but in doing so it creates a world file called deleteme.tfw. (The switch -co “TFW=YES” means to use the creation option “TFW=YES.”) The world file is what you really want, so delete the deleteme.tif. The world file and the main image file need to share the same name, so rename the TFW to myGeoTiff.tfw.

Re-sample and re-project

I like to do these two operations simultaneously with gdalwarp. In the following example, I am re-projecting out of 4326 (the EPSG code for lat/long/WGS84) into 32609 (the EPSG code for UTM Zone 9 north/WGS84).

Where

-s_srs
the Spatial Reference System (projection & datum) of the original
-t_srs
the SRS of the result
-r
resampling method
-tr
resolution, in metres
-of
output format
gdalwarp  -s_srs EPSG:4326 -t_srs EPSG:32609 -r bilinear -tr 5 5 -of GTiff myFloatDEM.tif myFloatDEM_32609_5x5.tif

Note that the terms “projection/datum”, SRS (Spatial Reference System) and CRS (Coordinate Reference System) are interchangeable for our purposes.

Now we’re reading to open up Blender and get started.

Overall process

The process will essentially go like this (and you can use this as a checklist once you get familiar with it):

  1. Select the Cycles Renderer
  2. Read in the DEM as DEM
  3. Create new material
  4. Adjust Z scaling
  5. Create and adjust the georef camera
  6. Correct final pixel dimensions to match the input DEM
  7. Adjust lamp location and properties
  8. Test render
  9. Adjust final parameters
  10. Render

The green steps are probably already familiar to you from Daniel’s tutorial.

Select the Cycles Renderer

Open up Blender and delete the default cube that appears. Switch the selected renderer to the Cycles Renderer.

Cycles

Reading the DEM in

File>Import> Georeferenced rasteropenDEMAsPlane

Select your floated, re-projected, re-sampled DEM, and in the left margin, where it says Import georaster, set Mode to As DEM, Subdivision as Subsurf, and select the correct CRS for your DEM.

importGeoraster

If the CRS you want to use is not on the dropdown menu,you can add it by clicking the “+” button, checking the Search box, typing in the EPSG code for your CRS into the Query box, and hitting Enter. It should appear in the Results box, and then you just click OK.

addingANewCRS

With all these things set, you click Import Georaster, and the result is a somewhat flabby-looking, grey plane.

flabby-looking grey plane

Tap n (with the cursor over the 3D view) to bring up the 3D view numeric panel, and we’re ready to go for a tour of things that BlenderGIS has done automatically. You don’t ordinarily have to check these things, but they are worth knowing about.

  • On the 3D view numeric panel we can see that the dimensions of the plane are the real-world dimensions of the DEM. In this case I’m using a 5030 x 3277 DEM with 5 m cells. For some reason BlenderGIS drops one cell in each dimension, so the DEM is 5029 x 5 on the X dimension (25145), and 3276 x 5 on the Y dimension (16380). metres on a side. We also see the thickness (Z dimension) of the DEM, which is 1911.56 metres.

transfromPanel

  • Also on the 3D View Numeric Panel, under Display, we can see adjustments to the grid floor. It gives values for Lines and Scale. Ordinarily (e.g., with the default cube) these will be values like 16 lines and a scale of 1. This means lines extend out 16 from the origin and no more, and each square is 1. After BlenderGIS has read in a DEM,  these can be more like 25 lines with a scale of 1000. In this case we have 25 lines on each side of the origin and each represents 1000 units (corresponding to metres).

gridFloor

  • Again on the 3D View Numeric Panel, under View, the clip distances represent the nearest and farthest you can see the object in perspective view. Ordinarily (e.g., with the default cube) these will be values like 0.1 to 1000. With a DEM read in these can be more like 0.1 to 250,000.

clipDistances

  • No longer on the 3D View Numeric Panel, but under Scene>Custom Properties, BlenderGIS has added a crs X and crs Y which correspond to the UTM coordinates of the centre of the plane. This is how the plane can be centred at (0,0), but any additional georeferenced data placed over it will be in the correct place. There’s also an SRID, and even a latitude and longitude.

customProperties

To understand how the DEM was used to displace a plane, go to the Modifiers tab for the plane. You will see two modifiers have been created. The upper one, called DEM, is a Subdivision Surface, or Subsurf. This kind of modifier subdivides the plane on the fly, and you can specify separate subdivision factors for View mode (what we’re doing right now) and Render mode (when you render the image). The beauty of the Subsurf modifier is that you can leave the View mode number low (e.g., 6) but set the Render number high (e.g., 11).

The factor is used as an exponent of 2, so our View mode is 2^6 or 64 subdivisions along the side of the plane. This is a small number considering that the DEM itself is thousands of points on a side. Which is why the view looks lumpy. The second modifier, called DEM.001, is a Displace modifier, and it applies the DEM to the subdivided plane much as in the Subdivide method.

Try changing the View number on the Subsurf modifier, and watch the fineness of the plane increase. You’ll notice that it takes longer to adjust the 3D view each time you increase the View number by 1, because you are doubling the number of subdivisions. A view number of 6 corresponds to 64,  7 to 128, 8 to 256, 9 to 512, 10 to 1024 and 11 to 2048 subdivisions. increasingViewNumberBy the time you get to 10 or 11, you’ll notice the status bar at the top of Blender is reporting a lot more memory use, 1.2 GB in this case.

statusBarat11

Go ahead and set the View number back to 6, but increase the Render number to 10 or 11. Because this higher number of subdivisions will be performed on the fly during render, it takes no memory now.

Create new material

If I now go to the Material panel for the plane, I’ll see that at this point I have no material.

noMaterial

I click New to create a new one. I get the default Diffuse BSDF  material, with a roughness of 0. This is fine, and if I want I can change the material or some of its parameters later.

defaultMaterial

Adjust Z scaling

Last, we’ll adjust the Z scale factor in the 3D View Numeric Panel to reflect the vertical exaggeration we want. I’ll use 5 in this case. Note that the Z dimension of the plane jumps up to 5 times its previous value, to 9794.137 in this case.

ZScaled

Go on to Part 2.

Cutting large DEMs into square, partially overlapping tiles

I’ve been doing a lot of shaded relief using Blender, and one practice I’ve developed is feeding Blender square DEM tiles.

The DEMs I use in maps tend to be too big for my installation of Blender to handle. One I had recently was 6348 x 10899 pixels; another was 4155 x 3572.  So that’s one reason for tiling them: I can produce shaded relief of the tiles separately in Blender and then mosaic the results.

And as long as I’m tiling, I’d like to use square tiles. That’s because when you subdivide a plane in Blender you have to use the same factor on both X and Y axes. Call me superstitious, but I get uncomfortable knowing that I’ve mapped 1803 pixels of DEM onto 1500 cells along one axis, and 2403 pixels to the same 1500 cells on the other. It seems to me this should cause loss of precision.

By feeding Blender square tiles, I eliminate that. I can subdivide the plane 2048 times on both axes, and feed it a 2050 x 2050 DEM.

This is not the kind of tiling that you would use a tool like gdal2tiles.py to do. That tool produces tiles that do not overlap, but I want some overlap. This is for a couple of reasons.

One is that Blender sometimes produces a bit of an edge artifact on a tile’s down-light edge (that is, the edge of the tile away from the light source). I can trim this off, and if there’s a bit of overlap between adjacent tiles, I get a nice smooth mosaic.

Left: mosaic of two non-overlapping tiles -- note the "seam" down the centre. Right: the left tile trimmed, and the right tile (with an overlap of 200 pixels) laid under it. Light is from 337.5°

Left: mosaic of two non-overlapping tiles — note the “seam” down the centre. Right: the left tile trimmed, and the right tile (with an overlap of 200 pixels) laid under it. Light is from 337.5°

Another reason for overlap is that if a shadow-casting feature (a mountain, a cliff) is just off the up-light edge of a tile, Blender can’t know to cast the shadow from that feature onto the tile. Overlap cleans this up by putting a healthy margin of features on the up-light side.

How much you overlap seems to depend on the topography and the scale. With small scale mapping, you don’t need much overlap because the mountains can’t cast their shadows very far. At large scales they can cast them many, many pixels.

Comparison of shaded relief produced for 1:1,000,000 (left), 1:150,000 (centre), and 1:30,000 (right).

Comparison of shaded relief produced for 1:1,000,000 (left), 1:150,000 (centre), and 1:30,000 (right).

To make tiling-with-overlap easier, I’ve written a bash script, called maketiles.sh, to tile large images, which you can download. This script uses the gdal_translate tool, so you also need to have GDAL installed.

Usage looks like this (assuming linux here):

./maketiles.sh <image file> <desired tile size in pixels> <number of columns in image> <number of rows in image> <tile_filename_prefix> <overlap>

It’s the rare DEM that can be cut perfectly into tiles without the last row and column being a bit off-sized. This script resolves that by increasing overlap on the final row and column so that all tiles are the requested size.

Let’s run through an example.

I have a DEM called 093L03-4-5-6_dem_26909.tif which is 4155 x 3572 pixels.

(If I don’t initially know its dimensions, I can get them through gdalinfo:

gdalinfo 093L03-4-5-6_dem_26909.tif
Driver: GTiff/GeoTIFF
Files: 093L03-4-5-6_dem_26909.tif
093L03-4-5-6_dem_26909.tfw
Size is 4155, 3572)

Let’s say I want tiles that are 2000×2000 pixels with a 100 pixel overlap, and I want the resulting tiles (and their world files) to be named tile_1_1.tif, tile_1_2.tif, etc.

The command would be

./maketiles.sh 093L03-4-5-6_dem_26909.tif 2000  4155 3572 tile 100

The script first echos back your parameters, and then tells you how many tiles will be produced. In this case we’ll be getting 6 tiles in 2 rows and 3 columns. It also tells you how much extra overlap you’ll have on the last row and column.

image file is 093L03-4-5-6_dem_26909.tif
desired tile size is  2000 x 2000
image is  4155 columns by  3572  rows
prefix is tile
overlap is 100 pixels.

Estimating 3 columns of tiles (X), by 2 rows of tiles (Y).

Taking a closer look at columns.
When laying out 3 columns of tiles, it looks like the final column will fall at 5800 pixels.
That's bigger than the image, so we'll back off the final column of tiles to begin at 2155 and end exactly at 4155.
This final column of tiles will overlap the second-to-last column by 1745 pixels.

Taking a closer look at rows.
When laying out 2 row of tiles, it looks like the final row will fall at 3900 pixels.
That's bigger than the image, so we'll back off the final row of tiles to begin at 1572 and end exactly at 3572.
This final row of tiles will overlap the second-to-last row by 428 pixels.

Layout has 6 tiles in 2 rows and 3 columns.

Finally, it computes the maximum overlap you could have without increasing the number of tiles, and asks if you want to try that.

Would you like to go with an overlap of 428, which would spread the overlap out but keep the same number of tiles? [Y/n]

If you accept this increased overlap, it confirms the change, and then you get a final chance to approve the process or exit.

Overlap set to 428.
overlap is 428.

Hit Enter to proceed with cutting tiles, or Ctrl-C to exit.

Assuming you hit Enter/Return, it then calls gdal_translate many times as is necessary, giving you extent info for each tile.

Producing tile tile_1_1.tif (row 1, column 1).  Extents: (0, 0) to (2000, 2000).
gdal_translate  -co "TFW=YES" -srcwin 0 0 2000 2000 093L03-4-5-6_dem_26909.tif tile_1_1.tif

Input file size is 4155, 3572
0...10...20...30...40...50...60...70...80...90...100 - done.

Producing tile tile_1_2.tif (row 1, column 2).  Extents: (1572, 0) to (3572, 2000).
gdal_translate  -co "TFW=YES" -srcwin 1572 0 2000 2000 093L03-4-5-6_dem_26909.tif tile_1_2.tif

Input file size is 4155, 3572
0...10...20...30...40...50...60...70...80...90...100 - done.

The final tile in the row will be a bit too big.
To keep all tiles the same size, I'm backing off the upper left corner of this tile from (3144,0) to (2155,0).

Producing tile tile_1_3.tif (row 1, column 3).  Extents: (2155, 0) to (4155, 2000).
gdal_translate  -co "TFW=YES" -srcwin 2155 0 2000 2000 093L03-4-5-6_dem_26909.tif tile_1_3.tif

Input file size is 4155, 3572
0...10...20...30...40...50...60...70...80...90...100 - done.
ROW IS FINISHED.

Producing tile tile_2_1.tif (row 2, column 1).  Extents: (0, 1572) to (2000, 3572).
gdal_translate  -co "TFW=YES" -srcwin 0 1572 2000 2000 093L03-4-5-6_dem_26909.tif tile_2_1.tif

Input file size is 4155, 3572
0...10...20...30...40...50...60...70...80...90...100 - done.

Producing tile tile_2_2.tif (row 2, column 2).  Extents: (1572, 1572) to (3572, 3572).
gdal_translate  -co "TFW=YES" -srcwin 1572 1572 2000 2000 093L03-4-5-6_dem_26909.tif tile_2_2.tif

Input file size is 4155, 3572
0...10...20...30...40...50...60...70...80...90...100 - done.

The final tile in the row will be a bit too big.
To keep all tiles the same size, I'm backing off the upper left corner of this tile from (3144,1572) to (2155,1572).

Producing tile tile_2_3.tif (row 2, column 3).  Extents: (2155, 1572) to (4155, 3572).
gdal_translate  -co "TFW=YES" -srcwin 2155 1572 2000 2000 093L03-4-5-6_dem_26909.tif tile_2_3.tif

Input file size is 4155, 3572
0...10...20...30...40...50...60...70...80...90...100 - done.
ROW IS FINISHED.
IMAGE IS FINISHED.

Tiles are named with the prefix you gave, followed by _X_Y.tif, where X is the row position, and Y the column position, of the tile. So in this case we got tile_1_1.tif, tile_1_2.tif, etc., along with tile_1_1.tfw, tile_1_2.tfw, etc