Shaded relief with BlenderGIS (2019), part 3

[Back to Part 2]

Why do this?

The shaded relief images (also often called a hillshade) that GIS software produces are pretty good. Here’s one made with QGIS.

QGIS 2-up

The hillshade by itself (left), and blended with hypsometric tinting (right)

Hillshading algorithms typically take three parameters: the elevation of the sun, the azimuth (direction) of the sun, and the vertical exaggeration of the landscape. (60°, 337° and 1X, in this case.) They do some ray-tracing to figure out what parts of the landscape are being hit by the sun and what parts are in shadow, and then produce a greyscale image.

You can tweak the brightness and contrast in your GIS software, or stretch the histogram to your liking. This can do a lot to lighten up the darkest shadows, or to turn the hillshade into a ghostly wash that merely suggests relief without dominating everything else on the map.

Much of the time what we can get out of GIS software meets our needs for a hillshade, which, after all, is not the final map but a merely a layer of the map. (Although occasionally it is the pièce de résistance…)

But perhaps you don’t like the glossy, shiny quality of that hillshade above. You think the mountains look like they were extruded in plastic. They remind you a little too much of Google maps Terrain layer. Maybe you’d like something that looks more like it was chisled from stone, like this…

Blender 2-up

Or perhaps you are interested in…

Blender 2-up 2 suns

A second sun, shinning straight down to eliminate the darkest shadows

Blender 2-up warm and cool light

A warm (yellow) sun in the northwest, cool (blue) shadows

Blender 2-up denoised

Denoising performed on the image after rendering

This might be why you are investigating Blender.

The big settings that make a difference (besides the conventional settings of azimuth, elevation and vertical exaggeration) are…

  • material
  • multiple lights and their colours
  • amount of light bounce
  • denoising the render

Material

In animation modelling, Material is what reflects, absorbs and scatters light. Blender spends much of its time tracing light rays and deciding how the surfaces they encounter affect them.

With your plane selected, go to the Material tab icon Material tab and hit New. The default material that comes up has these Surface properties. (And this is what Blender used for your basic render.)

Principled BSDF initial properties

Principled BSDF is a sort of super-versatile material that allows you to have all of these properties (subsurface scattering, metallic look, sheen, transmission of light) that in earlier versions of Blender were assigned to specific surface types, like “Diffuse BSDF,” “Glossy BSDF” or “Subsurface scattering.”

(BSDF stands for “bidirectional scattering distribution function.”)

These various surfaces are actually shaders, which are pieces of software that render the appearance of things (and may actually run on your GPU, not your CPU). If you click on “Principled BSDF” next to the word “Surface” a list of all of the possible shaders comes up.

shader list

You can learn a lot more about shaders in the Blender manual, but Blender essentially wants to give you the tools to be able to simulate any material, from water to hair, and some of the effects you can get applying this to shaded relief are pretty weird.

8 kinds of material

The same piece of terrain rendered with different shaders. Top, left to right: Diffuse BSDF with a wave texture, Glass BSDF, Hair BSDF, Principled BSDF. Bottom, left to right: Toon BSDF, Translucent BSDF, Translucent BSDF with Principled volume emission, Velvet BSDF.

The main shader you probably want to play with is the Mix Shader, which allows you to mix the effects of two different shaders. The Mix Shader’s  factor (from 0 to 1) determines how much the results are influenced by the second shader.

Mix Diffuse Glossy

On left, the original render; on right, a Mix shader (Fac=0.2) of Diffuse BSDF (defaults) and Glossy BSDF (IOR=2). This adds just a bit of glossiness to the surface.

 

Mix Diffuse Toon

On left, the original render; on right, a Mix shader (Fac=0.3) of Diffuse BSDF (defaults) and Toon BSDF (defaults). This brings in bright highlights that tend to wash out flat surfaces.

The other aspect of material that is worth experimenting with is value. Blender’s default material, the Principled BSDF, has a default colour of near-white, and its value is 0.906 (on a scale of 0 to 1).

principled BSDF default base colour

By darkening this you create a more light-absorbing material.

Principled BSDF base colour 0.5 lamp 45 337 str 3 v2

On left, the original render; on right, Principled BSDF (base color value turned down to 0.5).

You might think that the effect of a darker material is just to shift the distribution of pixel values down in the render, and that you could get the same effect by turning down the brightness of the original. But if you look at the histogram of each image you see the pixel values are distributed differently, and even when both are displayed with “histogram stretch,” they are distinct from each other.

base color 0.906 versus 0.5 both histogram stretched

Left: the value of the base color of material (Principled BSDF) is 0.906 (the default). Right: the value of the base color is 0.5. Both are displayed with histogram stretch.

And of course if it works for your map you can give the material real colour. (Don’t forget to save the render as RGB rather than BW). This, however, is very similar to colorizing your shaded relief in GIS software.

sandstone yellow material

Multiple lights and their colours

To add another light to your scene, go Add>Light>Sun. You can also experiment with adding other types of lights: point lights (which shine in all directions from a specific place), spot lights (which send out a specific cone of light) and area lights.

A second light can bring out features in a very nice way.

second sun from the east

On left: vertical exaggeration 2x, one sun in the NNW, elevation 45°, angle= 10°, strength = 3. On right: the same scene plus a second sun in the east, elevation 45°, angle = 1°, strength 0.5.

By giving colour to lights, you can create differential lighting and colour.

 

Diffuse 2x BSDF default 45 337 str 3 color orange LAMP2 str 0.5 az east color yellow

Same scene as above, but now the sun in the NNW is orange, and the sun in the east is yellow.

There is one more light in your scene that often goes unnoticed, and this is what Bender called the world background colour. You can think of this as a very far away surface that nonetheless does contribute a bit of light to your scene—from all directions. You can see the world colour in action if you render your scene with the sun strength turned down to 0.

world background only

No lights in the scene: mysteriously there’s still something there. This is the effect of the world background colour.

If there are places in your scene that sunlight does not reach, the world background still contributes to their illumination—and colour.

By default the word background colour is 25% grey (value = 0.25), but you can change this to a dark blue if you would like dark blue light to collect in your shadows.

world background blue

World background colour set to blue (Hex 414465), strength 1. Secondary pale yellow (Hex FFE489) light in east, strength 0.5

 

Amount of light bounce

Part of the charm of Blender is that it calculates how much light reflects off surfaces and then hits other surfaces, which is why even shadowed areas have some light. You can control, however, the number of bounces a light ray has before it expires. This is on the Render tab icon Render tab, in the Light Paths section.

light paths default

By default, the light paths settings are as above. With the material I am using, the Max Bounces for Diffuse materials is the important part, and by default it is 4.

The presets exist icon icon to the right of “Light Paths” indicates that there are presets available. Clicking here reveals three key presets:

  • Direct Light: light gets few bounces, and none off diffuse material
  • Limited Global Illumination: light gets one bounce off diffuse material
  • Full Global Illumination: light gets 128 bounces off diffuse material

Changing among these does not make a huge difference, but in scenes with deep shadows it is visible.

light paths compared

Clockwise from top left, the presets are: Direct Light, Limited Global Illumination, Full Global Illumination, and the default settings. (With vertical exaggeration at 2x)

Denoising the render

Blender offers the post-processing feature of denoising the render. I liken this to a blanket of snow on your landscape: it erases tiny details.

denoise comparison

Left: 2x vertical exaggeration, sun in NNW, elevation 45° strength 3. Right: the same, plus denoising (Strength = 1, Feature Strength = 1)

To activate denoising, go to the View Layer tab icon View Layer tab, scroll to the very last section, Denoising, and check the Denoising box.

denoising panel

I do not find that denoising does much if you keep the default settings. I tend to turn both Strength and Feature Strength up to 1.0.

In conclusion

This is only been the tip of the iceberg. Some things I have not talked about are

  • Applying a subdivision surface modifier to your mesh so that Blender is interpolating and smoothing it on the fly.
  • BlenderGIS’s ability to read your DEM As DEM Texture, which creates a plane with subdivision and displace modifiers instead of a plane whose every vertex corresponds to a DEM cell.
  • How to focus the orthographic camera down on one small part of your DEM where you can do quick test renders while you work out lighting, material, etc.
  • The node editor for complex materials
  • Using a perspective camera to shoot a scene that is not straight down
  • Cutting DEMs into tiles for Blender to work with, when the whole DEM is simply too much for the RAM on your computer.

Have fun exploring Blender’s many features!

Shaded Relief with BlenderGIS (2019), part 2

[Back to Part 1]

In the overall process we are now at step 3, but things will go faster now.

  1. Prepare your DEM
  2. Read the DEM into Blender as Raw DEM
  3. Render tab: select the Cycles rendering engine
  4. Adjust Z scaling (vertical exaggeration)
  5. Create and adjust georef camera
  6. Output tab: correct the final pixel dimensions of the output image to match the DEM. Set final image type to be TIFF
  7. Turn the light into a Sun, and adjust its properties
  8. Do a test render (F12)
  9. Render

At the end I’ll cover some of the variations on this process and extra tweaks you can do.

Render tab: select the Cycles rendering engine

Go to the Render tab icon Render tab and observe the Render Engine setting:

Render tab engine and sampling

Blender 2.8 comes with the Eevee render engine as the default, but Eevee is just a basic  engine for previews. Change it to the Cycles render engine.

Render tab cycles selected

If you have a powerful GPU (i.e., video card), you can set Device to GPU Compute.

While you are here, adjust Sampling. To the right of the word “Sampling,” note the menu icon presets exist icon. This indicates that presets exist for this setting. Clicking on it, select Final. This increases the sampling numbers for the render, which will improve image quality.

Adjust Z scaling (vertical exaggeration)

The first thing to adjust on your plane is what we usually call the vertical exaggeration, but Blender thinks of as the Z scaling.

With the plane selected, go to the Object tab icon Object tab. Here you should see a number of sections, the first of which is the Transform section.

Object Transform section

Under Scale increase the Z value if you want to create vertical exaggeration. You should see your plane change shape.

[It’s important to increase Z scaling before setting up the georef camera, because if you increase Z scaling afterwards, you may accidentally have set the camera below the tops of the highest features.]

In my case the terrain has plenty of relief, so I’ll leave Scale Z at 1.000.

Create and adjust georef camera

Make sure the plane is still selected (in the Outliner) and go GIS>Camera>Georef.

A new camera (“Georef cam”) will appear in the Outliner. You do not need to delete the old camera.

When the georef cam is selected, then in the 3D Viewport you’ll see something like this (you may have to pull back).

georef cam createdThis is an orthographic camera (meaning, in effect, that all points of its lens look directly down: there is no perspective distortion) placed over the plane.

To check that your camera is properly pointed at your landscape, and sees all of it, hover the mouse over the 3D Viewport, and hit 0 on the numeric keypad (or go View>Viewpoint>Camera). You should see what the camera sees (“camera orthographic view”)

camera orthographic view

Hit NumPad-0 again (or go View>Viewpoint>Camera again) to go back to regular (“User perspective”) view.

Output tab: Correct the final pixel dimensions of the output image to match the DEM, and set final image type to be TIFF

Go to the Output tab icon Output tab, and under Dimensions check the Resolution X and Y that were set up when you created the georeferenced camera. Also note the percentage (%), which is 100% at this point.

output resolution

Change Resolution X and Resolution Y to match the pixel dimensions of your DEM. In my case, the initial DEM was 839 x 702 pixels, so I enter these two numbers for X and Y. With the dimensions of the hillshade matching those of the DEM I can, at the end, apply the DEM’s world file to the Blender output, and georeference it.

It’s handy at this stage in the process to set percentage to something small, like 20%. This way you can do some quick test renders: Blender will produce a rendering whose pixel dimensions are only 20% of the overall Resolution numbers you set. (Before the final render you’ll be back here and set this to 100% again.)

Lower down on the same tab, under Output, set the parameters for the type of image you want.

Output output parameters

I typically set these to File Format = TIFF, Color=BW and Color Depth = 8, but you can get JPG, PNG, etc. If you choose later to assign colour to the world background and the sun,  other than grey or white, you will probably want to set Color to RGB or RGBA.

Turn the light into a Sun, and adjust its properties

Your plane is built, your camera is set: now all that is left to arrange is the Sunlight.

Select the Light in the outliner, and go to the Light object data tab icon Object Data tab.

Initially you should see something like this under Light:

Light object data Light initial

Click on Sun, and then change Strength to 1.

Click on the Use Nodes button below, and set its Strength to 2. It should now look like this.

Light object data Light final

You will want to play around with the Strength setting (in the Nodes section) when you do test renders. The strength of the light affects how bright your hillshade is.

The Angle setting for the sun is important. The Angle is how many degrees across the sun’s face looks from the ground. A sun with a smaller angle produces shadows with sharper edges.

sun angle comparison

Sun angle of 1° (left) and 10° (right)

A 1° sun is essentially a point source, and it casts sharp shadow. At 10° sun is bigger and the edges of shadows are diffuse. Take your pick.

Note that if you decide to use a coloured sun later (see “Part 3“) this is where you will change its colour: Light>Object Data>Nodes>Color.

Now go to the Object tab icon Object tab, and consider the location, rotation and scale of the light.

Light object transform initial

The location of the light does not matter: a “Sun” type light is considered to shine from an infinite distance regardless of its location. Scale should be left as all 1’s.

The rotation of the light however is all-important to us. It is here that we set the sun elevation and azimuth. These are both counter-intuitive in Blender, so here’s how it works.

Rotation X is always 0

Rotation Y is the complement of the familiar elevation angle (that is, 90 minus elevation). E.g., if you want a sun 60° above the horizon, set Rotation Y to 30°. Think of a light in a theatre clamped to a lighting pipe suspended over the stage. The pipe is the Y axis. The light is initially pointing straight down (0°) and you are swinging it up by 30°.

Rotation Z is the azimuth of the light, but instead of beginning at 0° for North and increasing clockwise (as we are familiar with from compass bearings), it begins at 0° for East and increases counterclockwise. (Think of angles measured on a coordinate plane.) It’s generally quickest to figure out by sketching on the back of an envelope, but the actual formula is that the Rotation Z = 450 – Azimuth, and subtract 360 from your answer if it is greater than 360.

A few common angles: common rotations ZIf I want a sun 50° above the horizon, coming from azimuth 337° (North-northwest), my light’s Rotation numbers will look like this:

light rotation for 50 at 337

Do a test render (F12)

Hit F12 to get a test render at 20% of the size of the final render. A little Render window opens and after a pause you (hopefully!) get something like this.

test render basic

The test render often catches common mistakes, like forgetting to set up your light or camera. It also reveals how the overall image will look, and so you can adjust your light strength, azimuth and elevation to increase or decrease shadows and the overall brightness of the image.

Render

Go back to the Output tab icon Output tab, and set percentage (%) to 100%.

Hit F12 to render.

Note that you can interrupt a render by hitting ESC, or by closing the render window.

I get this result:

basic render

Pressing Shift-S in the Render window (or Image>Save As) allows you to save the image. If I saved this as “Blender shaded relief 1x 50 337.tif” I would then make a copy of the TFW world file I created back at the beginning and name it “Blender shaded relief 1x 50 337.tfw” This makes something georeferenced that I can pull into my GIS.

That’s it! You made a hillshade in Blender!

Now, to consider the many other things we can fiddle with in Blender—material, number of lights, denoising, etc.—go to Part 3.

Shaded relief with BlenderGIS (2019), part 1

This tutorial replaces my tutorial from five years ago,  “Shaded Relief with BlenderGIS” tutorial. At this point (November 2019) I am using Blender 2.8 and the most recent BlenderGIS addon.  Daniel Huffman’s tutorial, which uses a different technique, has also been updated.

Blender 2-up denoised

Blender is tremendously complex. To keep things from getting out of hand, we’re going to take a straight path right through the centre of the software, with a focus on getting out the other side with a completed render of shaded relief. But at the end, we will return and explore some of the interesting side trails that lead to the features that make the Blender hillshades so interesting.

Here’s the basic structure we will follow to create a hillshade, once you have Blender with the BlenderGIS addon installed. This can function as a checklist when you are familiar with the process:

  1. Prepare your DEM
  2. Read the DEM into Blender as Raw DEM
  3. Render tab: select the Cycles rendering engine
  4. Adjust Z scaling (vertical exaggeration)
  5. Create and adjust georef camera
  6. Output tab: correct the final pixel dimensions of the output image to match the DEM Set final image type to be TIFF
  7. Turn the light into a Sun, and adjust its properties
  8. Do a test render (F12)
  9. Render

Initially, it was tempting to write this tutorial for an audience that uses GIS software but has never made its own hillshades. But this would involve explaining digital elevation models (DEMs), cell sizes, nodata values and projections, plus the art of combining shaded relief with other layers, stretching histograms and adjusting brightness and contrast—all too much to include here. So instead I am going to assume that you already use shaded relief, and have even tried producing your own in your favourite GIS software. You already know that, at its most basic level, creating shaded relief involves specifying a sun elevation and azimuth, and selecting a vertical exaggeration for your terrain.

I use the free GDAL command-line tools gdal_translate and gdalwarp to do re-projection and re-sampling, as well as produce world files. If the command line makes you queasy, your GIS software should offer these abilities as well.

Your first step will be to go to https://www.blender.org/ and obtained the free animation software Blender 2.8.

Installation of BlenderGIS

You only have to do this once, and then Blender is prepared to accept geographic data.

The BlenderGIS addon has installation instructions in its wiki at https://github.com/domlysz/BlenderGIS/wiki/Install-and-usage. Basically, they go as follows.

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

Within Blender,  go Edit>Preferences, and select the Addons tab. Click Install…, select the BlenderGIS-master.zip file, and click Install Add-on from File.

Once it is installed, be sure to check the box next to 3D view: BlenderGIS to enable it.

BlenderGIS installed

Note that a new menu appears in Blender: the GIS menu.

BlenderGIS menu

Preparing your DEM

BlenderGIS can read DEMs whose data type is Int16, UInt16 or Float32. It works best with DEMs that are in a projection measured in metres (as opposed to degrees).

You’ll notice later that BlenderGIS claims that it can read a DEM projected in “WGS84 latlon” —in other words, EPSG 4326. I have found this works only when you are reading that DEM into a Blender scene already georeferenced in a metric projection.  You have to do a kind of complicated head-stand to make this occur, so we’ll take the simpler route and feed to BlenderGIS a DEM that is projected in metres. Typically for me (in northern British Columbia) this would be UTM Zone 9 North/WGS84 (EPSG 32609) or the BC Albers Equal-Area projection/WGS84 (EPSG 3005).

I’m not going to suggest converting your DEM to Float32 data type. This makes a difference only with large scale mapping (e.g., 1:25.000 or larger), where you might actually see “steps” in your hillshade where the elevations change by whole metres. If you are doing such mapping, consider floating your DEM first.

You will also want a world file for the DEM, so you can georeference the hillshade Blender produces. World files used to be pretty common but now raster data so often comes with its georeferencing built in, I will explain what they are.

The world file was a clever invention: a six-line text file that contains the cell dimensions and the coordinates of the image origin (upper left corner). It is paired with the image file by having the exact same filename, but with a TFW extension (for TIFF images). (E.g., you could have myDEM.tif and myDEM.tfw.) It allows the georeferencing information to be held externally to the image file.

world file

When an image file has no georeferencing (i.e., it is an ordinary TIFF image, not a GeoTiff), the world file is all your GIS software needs to correctly place and scale the image. The only additional piece that it will have to ask you for is the projection of the image, in order to make sense of the world file.

For JPG and PNG images, the world file extensions are JGW and PNW, respectively.

We can use gdalwarp to re-project, re-sample and produce a world file. In the following example, I am re-projecting out of  lat/long/WGS84 (EPSG code 4326) into UTM Zone 9 north/WGS84 (EPSG code 32609), and re-sampling to 16 metre square cells. I will assume your DEM is a geotiff file.

gdalwarp  -s_srs EPSG:4326 -t_srs EPSG:32609 -r bilinear -tr 16 16 -of GTiff -co "TFW=YES" myDEM.tif myDEM_32609_16x16.tif

Where

-s_srs
the Spatial Reference System (projection & datum) of the original (the “source”)
-t_srs
the SRS of the result (the “target”)
-r
resampling method (bilinear, cubic, etc.)
-tr
resolution, in metres (give two values, one for X and one for Y)
-of
output format
-co
additional creation options (TFW=YES means create a world file)

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

So at the end of this preparation step I have three things:

  • a DEM in a metric projection (e.g., myDEM_32609_16x16.tif)
  • a corresponding world file (e.g., myDEM_32609_16x16.tfw)
  • I know the EPSG projection code for the projection the DEM is in (e.g., 32609)

Read the DEM into Blender as Raw DEM

Open Blender, and note the default cube that is always there in a new workspace. To delete it, hover the mouse over the cube and hit the Delete key.

On the GIS menu, go Import>Georeferenced raster, and navigate to your DEM file.

In the left margin, where it says Import georaster, set Mode to Raw DEM and select the correct CRS for your DEM. (If you do not check the Build faces box, you will get a point cloud.)

importGeoraster options

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 check Save to addon preferences and click OK.

addingANewCRS

After a pause, during which Blender builds a plane mesh out of your DEM, you hopefully see a grey, 3D rendering of your terrain, floating in a coordinate space.

initial screen

A quick tour of Blender

Let’s take a brief detour here and learn some of the the language of Blender, as well as some of its controls and unexpected behaviours.

Blender, being 3D animation software, is a world of vertices (points in space), edges (which connect vertices) and faces (which fill in between edges). In this case, the cells of the DEM have been translated into a square array of vertices, each with the appropriate height for the DEM cell it represents. These vertices have been connected with edges that form a square lattice like the cells of the DEM. Each square of four edges has a face.

The whole set of vertices, etc. is a mesh, specifically a plane mesh.

In the lower right  corner, you will notice a status bar:

initial status bar

This tells us that we have 588,978 vertices, 587,438 faces. My DEM happens to be 839 columns by 702 rows, and 839 x 702 = 588,978, so this is the right number of vertices. The number of faces is slightly less because no faces are generated by the final row and final column (839 + 702 = 1541, minus 1 face which would be generated by both the final row and final column, hence 1540 fewer faces than vertices.)

The status bar also tells us we are presently using 291MB of RAM. This will go up during the final render. Near as I can tell, available RAM is the only limit to the size of DEMs that Blender can handle.

You can see the vertices, edges and faces of your plane mesh if you zoom in and press TAB to toggle into Edit mode.

initial mesh in edit mode

Don’t worry if you can’t figure out how to do this just yet. If you do though, be sure to press TAB again to leave Edit mode before continuing on.

The upper right corner of Blender has a panel called an outliner that shows a tree of objects in the Blender world. At present it looks like this.

outliner

Notice that the plane (called “Thomlinson 32609 16m clip” in my case) is selected. In addition there is a camera and a light.

Below the outliner is a Properties panel. On its left margin you will see a series of tabs identified by these icons. Hovering over them reveals their names.

Blender tab selector

We won’t have anything to do with some of these (Particles, Physics, etc.) The important ones for our purposes are Render, Output, Object and Object Data.

Note: sometimes the collection of tabs changes, depending on what object in the Blender world is selected. At this moment the plane is selected, but if the light were selected some of these tabs would disappear, and the Object Data tab would display a Light object data tab icon green light bulb icon.

The rest of the screen is taken up by the large 3D Viewport panel. Navigating in this is similar to using Google Earth, but a little different.

  • Rolling the mouse wheel zooms you in and out.
  • Holding down the mouse wheel button while moving the mouse in this panel allows you to rotate around objects.
  • To pan, hold down Shift and the mouse wheel button while moving the mouse to either side.

The other thing to know about Blender is that some keys work only when the cursor is over the 3D Viewport.

  • Numpad -. (period) centres the view on the selected object.
  • NumPad-0 toggles between this view and camera view.
  • TAB toggles the selected mesh in and out of Edit mode (and we won’t be using this).

Remember that Blender is a studio for animators. It is designed to allow you to sculpt new objects, and the generate many frames of animation in which various objects are lit by lights and shot from a specific camera. We are not going to create any new objects, and are going to shoot only a single frame; so there are many controls we will not use.

However, setting up our light and camera are crucial, and that’s what we do in Part 2.

 

 

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.

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

Note: This tutorial has now been superseded by a new one which uses Blender 2.8.

[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 tutorial about created shaded relief using the free, open-source animation software, Blender. This is an update of Daniel’s video tutorial for his method.

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