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.
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.
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.
The process will essentially go like this (and you can use this as a checklist once you get familiar with it):
- Select the Cycles Renderer
- Read in the DEM as plane
- Delete default material and replace with new material
- Read in the DEM as DEM
- Adjust Z scaling
- Create and adjust the georef camera
- Correct final pixel dimensions to match the input DEM
- Adjust lamp location and properties
- Test render
- Adjust final parameters
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.
Reading the DEM in
File>Import> Georeferenced raster
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.
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. By 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.
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.