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 6.

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

Create and adjust the georef camera

With the plane selected, go to the the (new) GIS tab on the left, and click Create/Update.

This creates a camera at an appropriate distance above the plane. It will appear in the tree of objects (upper right) as “Georef cam”. Notice that it has been placed high enough (Z value of 11,706) that it is above the highest point on the plane surface.

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.

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 they are 6150 x 6150.

However, I like the dimensions of my final image to match the DEM. That way I can pair the TFW world file from the DEM with the image created by Blender, and read it into QGIS georeferenced. So I’ll just change these to 2050. 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, here we have 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. lampSettingsThen 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. In this case my 2050 x 2050 tile takes just over an hour on my machine. 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.

Better yet, you can render many tiles at once using a script. Having saved all the .blend files, just make a script that includes the line above many times with the parameters changed on each one as necessary. E.g.,

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

Which brings us to our last section…

Saving the .blend file for another tile with the same parameters

If you have tiled your DEM, you may want to save the .blend file for this render and then pull in the next tile. You can save its .blend file with the same parameters — not having again to spend the time for setting up the lamp and camera. This allows faster production of .blend files, which can then be batch-rendered through a command line script.

This assumes all your tiles are the same size and share the same projection, but, of course, are in different locations. It should take less than two minutes per tile.

  1. Open a previous blender file
  2. Re-save with a new name for this tile.
  3. Delete the plane you have been using
  4. Read in the new DEM (File>Import>georeferenced raster.) Select On Plane and check Consider georeferencing and Adjust 3D view.
  5. Read in the new DEM again As DEM.
  6. On the new plane
    • Material>delete the default material and make a new one,
    • Texture>Image Mapping>Extend
    • Modifier>set the Render subsurf to the value you used with the previous tile
    • Modifier>add a Smoothing modifier if you had one in the last tile
    • set the Z scaling to the value you used with the previous tile
  7. Set Origin>Geometry to Origin
  8. Render>Light Paths>Limited Global Illumination
  9. Test camera Ortho view — see if there are any black patches — if so, adjust camera clipping distances
  10. Save again

Shaded relief with BlenderGIS, part 1

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 propose to give a tutorial for how to use Blender GIS’s subsurf method to achieve results similar to those obtained with the subdivide method. The Subsurf method isn’t always better, but it’s a good tool to have in your pocket. The method I use depends on the resolution of the DEM I’m using, how noisy the data is and, well, how much time I have.

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. Part1 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 download:

These four items go in Blender’s install tree, in the scripts/addons subdirectory. For me this is ~/blender-2.70a-linux-glibc211-x86_64/2.70/scripts/addons

Then, within Blender,  at File>User Preferences, the Addons tab, enable the four new features by putting a checkmark in the checkbox next to each:

  • 3D view: Set georef cam
  • Import-Export: Import raster georeferenced with world file
  • Import-Export: Import/export from ESRI shapefile file format (.shp)
  • Mesh: Delaunay Voronoi

That’s it. BlenderGIS is installed.

Preparing your DEM

In theory, whether you’re using the Subdivide or the Subsurf method, you can work with DEMs of any size — until you exceed your RAM. However, as I posted previously, I like to use  square DEM tiles. If I have a DEM whose pixel dimensions are 3805 x 2456, I’ll cut it into four  2050 x 2050 tiles.  For my workflow in doing this, as well as the bash script is use to cut large images into tiles, see that previous post.

Later, I’ll explain why I choose a size of 2050 x 2050 for these tiles.

Projection and resolution

Blender needs a DEM that is projected in a projection that uses metres (or feet, I suppose), not one that uses degrees. Typically for me this would be a UTM projection. Describing how to reproject is beyond the scope of this post, so I’ll assume you know how to do that.

You will probably want, at this point, to check the resolution of your DEM, because the scale of the final map will require a certain resolution. I tend to 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 10 metre DEM. For a 1:30,000 scale map, I want a 3 metre DEM. It’s just a general guideline.

Canadian DEMs in the province of British Columbia come in a lat/long projection at a resolution of 1 second. At my latitude (about 55° north), this means each cell is roughly 12m x 24m. If I reproject to UTM without specifying a cell size, the default resampling gives me a compromise dimension, maybe something like 19m x 19m. I often then wind up resampling to a different resolution, such as 15m x 15m.

World file

Using square tiles is optional, but BlenderGIS definitely needs an DEM file paired with a world file. The DEM can be TIF with TFW, JPG with JGW, PNG with PGW or BMP with BPW. I tend to work with Geotiffs, and 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.

Image Type

BlenderGIS wants DEMs that are either Unsigned 16-bit Integer (UInt16) or Signed 16-bit Integer (Int16). If your DEM is Float32, you’ll need to convert it to Int16 or UInt16. An easy GDAL line for this is

gdal_translate -ot Int16 myFloatDEM.tif myIntDEM.tif

where the -ot switch specifies output type.

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 plane
  3. Delete default material and replace with new material
  4. Read in the DEM as DEM
  5. Adjust Z scaling
  6. Create and adjust the georef camera
  7. Correct final pixel dimensions to match the input DEM
  8. Adjust lamp location and properties
  9. Test render
  10. Adjust final parameters
  11. 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 Cycles Renderer.Cycles

Reading the DEM in

File>Import> Georeferenced rasteropenDEMAsPlane

Select your DEM, and in the left margin, where it says Import georaster, set Mode to As Plane and make sure Adjust 3D view is checked.

The result is a large black plane.

You’re zoomed in a little too far to comfortably see it, so, with the cursor hovering over the 3D view, hit period (.) on the NumPad to zoom to view the whole selected object.fullPlaneBlack
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 2050 x 2050 DEM with 15 m cells, so the DEM is 30,750 metres on a side. This is quite different from the Subdivide method, where our plane might be 2 or 3 units wide.
  • 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 31 lines on each side of the origin and each represents 1000 units (corresponding to metres).
  • 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.
  • No longer on the 3D View Numeric Panel, but under Scene>Custom Properties, BlenderGIS has added a Georef X and Georef 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.

Delete default material and replace with new material

If you now go to the Material panel for the plane, you’ll see BlenderGIS has created a material called “RastMat.” I would rather choose my own material, so I will delete this by hitting the “-“ to the right of it. At this point I have no material.

And then I click New to create a new one. I’ll 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.

Read in the DEM as DEM

Now it’s time to displace the plane with the DEM. Again, I go File>Import>georeferenced raster, but this time, when I choose the (same) DEM file, I make sure that in the left margin, where it says Import georaster, I set Mode to As DEM and adjust Image depth to 16 bits signed. (If you use 16-bit Unsigned Integer DEMs, you won’t need to change this.)

What you will get is a plane that has been displaced, but looks somewhat lumpy and generalized.

To understand how the displacement worked, 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 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 Z dimension of the plane has tightened up to match the actual dimensions of the DEM. In this case, the minimum value of the DEM is 739, and the maximum is 2379. The distance between the two is 1640, and this is what the Z dimension now shows.1640

At 11, the status bar at the top of Blender is reporting a lot more memory use, 1.7 GB in this case.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.

With a tile that is 2050 x 2050 I will set Render number to 11. Render 11 corresponds to a plane subdivided 2048 times on a side. If you picture a plane with 1 division, it has 3 vertices per side: one at each corner and one at the midpoint of each side. If subdivided twice, it has four vertices on each side, and so on. At 2048 subdivisions, it has 2050 vertices per side. This is how many vertices I want to displace, hence the DEM tiles of 2050 x 2050.

I have found that using tile sizes that do not correspond to 2^n + 2, I get a faint grid of lines across my final rendered image. The trustworthy, 2^n + 2, sizes are 258 (256 + 2), 514 (512 + 2), 1026 (1024 + 2) and 2050 (2048 + 2). I suspect these lines are an aliasing effect that has to do with Blender resampling to image to apply it to the number of vertices in the subdivided plane. Since the subsurf method restricts me to subdividing the plane by powers of 2, I match it by using tiles of 1026 (with Render 10) or 2050 (with Render 11).

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 if you still have render set to 10 or 11, the Z dimension of the plane jumps up to 5 times its previous value, to 8200.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

Moving Mountains

I’ve been reading Gottfried Merzbacher’s Central Tian Shan which was published in 1905. (Archive.org has a free copy here.) This German geologist spent 1902-3 in the Tian Shan exploring valleys, noting strike and dip, and, among other things, trying to figure out the actual location of the peak Khan Tengri. This 7000m peak was not marked correctly on the 40-verst Russian map he was using. (40-verst maps were a series at 40 versts to the inch, or 1:1,680,000.)

While reading Merzbacher’s account, I looked at my map of Charles Howard-Bury’s journey in the Tian Shan (posted previously) and discovered I too had Khan Tengri in the wrong place!

So I did a little homework to make sure I had it moved it to the right place, because not being sure where Khan Tengri is located should be the kind of problem they only had 100 years ago, right? Well, maybe not.  Let’s compare Google maps, on the left, and Bing maps on the right.

Google-Bing Khan Tengri comparisonWow, that’s some pretty serious difference. Not only is Khan Tengri in a different place (relative to the glaciers, for example), but the Chinese-Kyrgyz border is also in a different place. In Bing’s world the border runs atop the highest peaks, while in Google’s the upper portion of the South Inyl’chek glacier is actually in China.

Bing, by the way, is using Russian peak names, like Pobedy Peak, while Google has its Kyrgyz name, Jengish Chokisu.

Turns out that OpenStreetMap has the same border Google does; but Yahoo and Here have the border that Bing does!

A quick web search reveals that China and Kyrgyzstan settled their border disputes in 2009. It’s likely that one of the sets of lines we see above dates from prior to that settlement, and one from after. But which is which?

Administrative boundaries downloaded from GADM and DIVA indicate that Bing/Yahoo/Here are actually more up to date.

But as far as the location of Khan Tengri is concerned, I’m afraid Google/OSM are correct.

Looking for Lepzinsk

In chapter 3 of Charles Howard-Bury’s Mountains of Heaven, written in 1913 but published in 1990, the author details a ten-day journey on Russian post roads. (I discussed a map I made of the central part of his journey in this previous post.)

He arrives in the town of Semipalatinsk by steamer (coming up the Irtish river from Omsk) on June  7th,1913, and then he and his servant John Periera prepare for this overland journey by buying a tarantass, a type of 4-wheeled cart. Although in the post road system one would typically receive a new horse and tarantass at each post station (about 20 to 30 versts apart — a verst being 1067 metres or slightly longer than a kilometre), Howard-Bury has strategized that by providing his own tarantass he’ll avoid repacking the heaviest boxes at each change of horses.

The list of luggage is revealing. A valise and bedding, camp bed, bath, table and chair, baking powder and jam in tins, photographic films, a few varieties of soups, some condensed milk and “one or two delicacies such as sardines and potted meats;” two rifles, a shotgun, and ammunition; tea, bought in Moscow; and “a plentiful supply of literature.”

They depart Semipalatinsk on June 9th in the late afternoon. The first three days of post-road travel takes them to the town of Sergiopol, covering on those days, by his reckoning, 47, 101 and 124 versts respectively (53, 108 and 132 km). Mountains of Heaven does not give us a detailed map, but it’s not hard to determine the route to Sergiopol, which is an old name for the town of Ayagoz, Kazakhstan. By following the most direct road between the two towns on Soviet maps made in the 1980s, I can get just these distances that Howard-Bury reports.

After Sergiopol, things get a bit trickier. Howard-Bury and Pereira stop sleeping at post stations and the author takes to giving distances from one midnight to the following midnight. At the end of day 6, they are one stage beyond the town of “Kapal,” which it’s easy to be confident is the still-extant town of Kapal, Kazakhstan. But by which route did they get from Sergiopol to Kapal?

Howard-Bury mentions all sorts of landmarks over these three days, which we can try to match up with the terrain. He gives the travel distances as 140, 153, and 105 versts, respectively (150, 163, 112 km). The most crucial landmark is a town named “Lepzinsk,” which they pass through towards the end of Day 5. If we can locate this Lepzinsk, we’ll be close to drawing the right post-road connection from Sergiopol to Kapal.

Here are the clues which Howard-Bury gives us.

Day 4 (140 versts, starting at Sergiopol):

  • For some distance we followed the valley of the Aksu, crossing and recrossing the river…
  • Towards evening we got into flatter country again, where the plains were covered with wormwood as far as the eye could see.

Day 5 (153 versts):

  • We continued travelling all night but at eight o’clock in the morning had only succeeded in covering 32 versts since midnight, owing to the sandy nature of the soil…
  • We had to cross several arms of the desert which stretched across the road from the west in broad rivers of sand several miles wide…
  • We now had to cross a range of rocky hills…
  • We are now not far from the shores of Lake Balkhash… there is good fishing to be got at this end of the lake where the Sergiopol River flows into it.
  • Low rocky ridges alternate with broad grey and yellow plains…
  • Further on the plain became at times white as snow from the alkali in it that had come to the surface…
  • Towards evening we got our first glimpse of lofty snowy mountains to the south…
  • The going was very sandy and heavy. We kept crawling slowly up and down over big sand dunes, when all of a sudden in the midst of the sand, we came across a fair sized river on the far side of which, across a wooden bridge, stood the town of Lepzinsk, set in a bower of green. What a delightful contrast!… There were pleasant grassy meadows all around, well irrigated and with many a wild flower growing on them. The contrast to the barren country through which we had been passing made Lepzinsk look doubly beautiful. …From now on to the mountains, the rainfall is greater and the climate milder, so that vegetation becomes everywhere more luxuriant.

Day 6 (105 versts):

  • It took us ten hours to do 47 versts, as we had to cross a range of rocky hills that formed a kind of buttress to the snowy Ala-tau mountains. The road… followed the bottom of a narrow and stony gorge…
  • On arrival at the top of the pass, which was only about 4,000 feet in height, there was a glorious view of the snowy mountains, only separated from us by a narrow grassy valley.
  • From now onward, the road kept at the height of over 4,000 feet and we drove across beautiful grassy meadows…
  • Soon after mid-day we arrived at Kapal… Except for a few trees at Sergiopol and Lepzinsk, this is the first place in the 650 versts that we have traversed since leaving Semipalatinsk, where they are at all plentiful.

So there are some wonderful clues here: the arms of a sandy desert coming in from the west, the proximity of Lake Balkhash, verdant Lepzinsk in the midst of sandy dunes and first seen on the far side of a “fair sized river,” a pass at 4,000 feet in a range of rocky hills “that formed a kind of buttress to the snowy Ala-tau mountains,” and a plateau over 4,000 feet that leads on to Kapal.

To make initial sense of all this, let’s take a look at this map. This is a portion of a Soviet topographic map at 1:1 million produced in 1973.

Locales that are easy to identify from Howard-Bury's narrative

Locales that are easy to identify from Howard-Bury’s narrative

As you can see, Kapal lies near the large northeast-trending mountain range at the bottom of the map, the Dzhungarskiy Ala-tau, which Howard-Bury calls the Ala-tau mountains.

Starting from the end and working backwards, Let’s look at the area around Kapal for a traveller arriving from the north.

L-44 version Kapal

This topography seems to correspond well with what Howard-Bury describes on Day 6. There’s a line of rocky hills forming “a kind of buttress to the snowy Ala-atu mountains,” and there a pass through them marked, on this Soviet map, with an elevation of 1052 metres, or 3500 feet. From the pass to Kapal the route continues to gently ascend.

The town at the foot of that 1052 m pass was marked as Джансугуров/Dzhansugurov in Russian. Now it is Жансугиров/Zhansugirov, Kazakhstan. We have reduced our task to determining the route from Sergiopol to Zhansugirov.

Let’s back up for a moment and look at rivers in the area.

Sergiopol to Dzhangusurov rivers2

There are three main rivers that feed Lake Balkash in this area. From north to south today they are called the Ayaguz (Аягуз), Lepsi (Лепсы) and Aksu (Аксу). These are three of the famous Seven Rivers that give the region its name in Kazakh (Jetyysu) and Russian (Semirechinsk).

Howard-Bury mentions following the Aksu as he leaves Sergiopol, but confusingly the Aksu is the southernmost of the three rivers we see here. It’s possible that after staying up all night for so many days became confused about which river he saw on which day. But it’s also possible that someone told him the Ayaguz was called “Aksu,” because Aksu is a very common name for rivers in Central Asia (meaning, essentially, white water or clear water).

He also refers to a Sergiopol river, which, from his reference to it entering Lake Balkhash we can assume is the Ayaguz.

We might expect a town called Lepzinsk to be on the Lepsi river, since “Lepzinsk” is essentially a Russian construction meaning “of or pertaining to Lepsi.” Surprisingly, there are three candidates for this town, and, correspondingly, three main route choices.

Sergiopol to Dzhangusurov 2

Western Route

This town of Lepsy is absent on historical maps, but it exists today and can be found, for example, on Google maps. We’ll call this one “Low Lepzinsk” because it is lowest on the Lepsi river. Its satellite image shows that one would indeed come across this river “all of a sudden in the midst of the sand.”

Lepsy satellite

This route is initially appealing because it mimics the railway corridor shown on the 1973 map. It passes close to Lake Balkhash. But, there are a few problems. The town of Lepsy is on the north side of the river, and Howard-Bury says that he sees the river “on the far side of which, across a wooden bridge, stood the town of Lepzinsk” As well, prior to getting close to Lake Balkhash, arms of desert crossing this route would come from the east.

Eastern Route

This Lepzinsk still exists, and we can see it shown as Lepsinsk/Лепсинск on the 1973 Soviet topo. It’s much higher on the Lepsi river, near its headwaters, so we’ll call it “High Lepzinsk.”

Lepsinsk/Лепсинск is right in the centre of this detail

Lepsinsk/Лепсинск is right in the centre of this detail

There are a number of problems with this route as well. Overall the eastern route is too long to match the distances Howard-Bury gives. It goes nowhere near Lake Balkhash. A satellite picture indicates that rather than travelling “up and down over big sand dunes” before coming upon the Lepsi river in this area, one is in a well-vegetated region.

approach to Lepzinsk3

Central Route

This town of Lepzinsk, which we’ll call Mid Lepzinsk because it is mid-way along the Lepsi river,  is shown on a 1912 map of Central Asia, by Alexander Keith Johnson:

Detail of Central Asia by the London mapmaker Alexander Keith Johnson, 1912

Detail of Central Asia by the London mapmaker Alexander Keith Johnson, 1912

Notice he actually has two Lepzinsks. There is Mid Lepzinsk, labelled ‘Lepsinski,” and also High Lepzinsk, here labelled ‘Lepsinskaia (Lepsa).” Also, importantly, he shows the post road from Sergiopol to Kopal running through Mid Lepzinsk.

Other period maps corroborate. Here’s Edward Stanford, 1901:

Detail of Central Asia by London mapmaker Edward Stanford, 1901

Detail of Central Asia by London mapmaker Edward Stanford, 1901

Again we have the two Lepsinsks, and the same post road passing through Mid Lepzinsk.

Here’s Vivien de Saint-Martin & Schrader, Atlas Universel de Geographie, published in 1935:

Detail from Turkestan by Louis Vivien de St Martin, Paris, 1935

Detail from Turkestan by Louis Vivien de St Martin, Paris, 1935

They too show both Lepsinks, although High Lepzinsk is bigger, and Mid Lepzinsk is ‘Lepsinskoïe.” At this point the railway has been built to the west and the older post road is not shown

It’s also probably worth looking at the map that comes with Mountains of Heaven.

Detail from the map in Mountains of Heaven. (Original is greyscale.)

Detail from the map in Mountains of Heaven. (Original is greyscale.)

Lepzinsk is there, but it is hard to say from this map whether Mid or High Lepzinsk is intended. Looking at other maps we can see that a bearing from Kapal to Mid Lepzinsk is more north than east, whereas one to High Lepzinsk is more east than north. Without longitude lines on this map we’re making a guess, but the Lepzinsk shown seems to be north-northeast of Kopal, and therefore would be Mid Lepzinsk.

What’s confusing is that the only river shown is the Aksu, which we know to be south of the Lepsi, and in fact the author’s route is shown as going to Lake Balkhash and then up the Aksu to Lepsinsk. This seems an unlikely route choice, as if one actually went up the Aksu he would have to turn north to get to a town on the Lepsi.

An argument against Mid Lepzinsk is that it no longer exists by that name. At this site today there is a town and a crossing of the Lepsi (and it is suitably situated with sandy areas north of the river, as the satellite image below shows) but it is called Kokterek (or “Blue Tree”), Kazakhstan. I can’t find any evidence it was ever named Lepzinsk.

Modern day Kokterek, Kazakhstan -- formerly known as Lepzinsk?

Modern day Kokterek, Kazakhstan — formerly known as Lepzinsk?

In favour of Mid Lepzinsk and the Central Route option is that on this route it is easy to make distances agree with the daily distances Howard-Bury gives.

So I vote that Howard-Bury’s Lepzinsk was the town today known as Kokterek, Kazakhstan, and that he took what I have shown above as the central route from Sergiopol to Kapal.

Mapping the travels of Charles Howard-Bury

[Update: an expanded version of this post has been published in Ripcord Adventure Journal, vol. 1, no. 1]

I’ve just finished a map for a book that had no map: Mountains of Heaven: Travels in the Tian Shan, 1913. This version of the Charles Howard-Bury’s journals from his travels in what are today eastern Kazakhstan and the adjacent region of China, was published in 1990 and although it’s out of print, it’s generally available in the used book market.

Here’s the map.

Charles Howard-Bury in the Tian Shan 1913 - map only - smaller

What I set out to do here was to resolve all the geographic references in the text. That is, you’re reading along, and he mentions the Yulduz plains or the village of Karachekinskaia, and you don’t know where that is. The various maps, atlases and website you can check don’t have it. So the book needs a supporting map. Pretty simple.

But it evolved into something else. Initially I planned a 24″x24″ sheet of paper, a size typical of folded maps that used to be glued inside the back covers of books. It would contain three maps. There would be a 1:1 million map detailing the central portion of his journey, the area where he spends the most time and mentions the most local details. There would be a 1:3 million map showing how he got from Semipalatinsk to Jarkent, the ten-day journey along Russian post roads for whose specific route he gives intriguingly few clues. Last, there would be a map at 1:7 million showing his exit route from Central Asia, from Jarkent to Tashkent, and then on by rail to the Caspian Sea.

Then things happened which illustrate some general hazards about mapping for old books.

The first map became the main project, and eventually pushed both other maps off the sheet. It’s what you see above. There were simply too many details to squeeze it in at a smaller scale onto part of the sheet. I let go of the idea that all the mapping one would need could fit into 24″x24″. I began building a hillshade of the area in Blender, and figured out an hypsographic colouring scheme inspired by an old Soviet topographic map of the area (К-44-Б): low elevations are white, and then rise through oranges and brown to white summits.

At the same time, as I tried to figure out just where Howard-Bury had gone, I found myself buried in old maps, wikipedia articles, and all sorts of other documents such as you get when you search on obscure terms like “Kuitun River.” I was reading maps in Russian and putting German websites through Google Translate. With each Aha! answer to a question, I felt more convinced that these puzzles — for the map had come to represent a series of puzzles — would make good stories.

At this point the map became a poster, the bottom half of which contained explanations of stuff. Here’s the poster:

Charles Howard-Bury in the Tian Shan, 1913_smaller

The sections in the bottom half cover the reason for the map, the plethora of alternate place names in Central Asia, a glossary of obscure terms Howard-Bury uses, an explanation of the Soviet topographic mapping system, some Turkic geographical equivalents, and a list of references. It is, in short, a brief tour of all the places I had to go to produce the map.

I’ll do another post here to present an example of one of the interesting puzzles: where was Lepzinsk?

Moonlight Mountain

Unnamed ridge

Unnamed ridge of Kispiox Mtn.

Yesterday I went up to Kispiox and hiked the Moonlight Mountain trail. I’m not sure the folks who showed me this trail would like its location publicized yet, so I won’t disclose precisely where it is, but it approaches Kispiox Mountain from the north.

The name “Moonlight” is a bit of a mystery, since there is no gazetted Moonlight Mountain, nor is there a Moonlight Creek, or any other “Moonlight” feature in the area. My theory is that it is a corruption of Moonlit Creek, which runs west from here into the Kitwanga River near Kitwancool Lake.

As the photo below shows, the rock here is composed of tilted – seriously tilted – strata. They vary from grey to deep black. So black, in fact, that I thought of oil, and then it occurred to me that these might be Bowser Basin sediments.

North Ridge of Moonlight

North Ridge of Moonlight — zoom in to see the tilted strata.

Bowser Basin sediments, for those of you who don’t live in northern BC, are the infamous rock unit that hosts quite a bit of coal, gas and oil — various kinds of hydrocarbons — and are responsible for much of the unwelcome interest that oil and coal companies have in northwest BC.

And sure enough, these rocks are Bowser Basin sediments. A map of that rock group (try here) shows that once you go north of Hazelton you are in Bowser territory. However, this area does not appear, thankfully, to host coal. (It probably does host recoverable oil or gas were one to use something like fracking.)

These beds are so beautifully up-ended and even perhaps overturned: you can see BC’s tectonic assemblage process in action.