Converting DEM rasters into OBJ meshes for use in Blender

This is the technique I use to get Digital Elevation Models (DEMs) into Blender from QGIS. Here's an example using this technique:-

Indian Ocean. Sea was added as an additional, transparent mesh at sea level (elevation 0)


For this, you’ll need

- Blender (I use 2.79)
- QGIS
- Python
- dem2obj (here on Github)
- git

First of all, fork the project in Github. That'll give you a fresh copy to play with (and enhancements/fixes are welcome)

cd /path/to/project
git clone https://www.github.com/your_user_name/dem2obj
pyvenv venv
source ./venv/bin/activate

You should now see the command line prompt change to (venv).

(venv) pip install -r requirements.txt

That might take a while (a few minutes, potentially)

Running dem2obj

When that's done, you should be able to run this. If you're not seeing the (venv) prompt, you'll need to reactivate the virtual environment, as follows:-

cd /path/to/project
source ./venv/bin/activate

You can now run the tool. You need to run AT LEAST --input FILENAME and --output FILENAME.

If your DEM is in degrees (e.g. WGS84) and elevation in meters, pass in the -w (or --wgs84) flags

(venv) python rastertomesh.py --input /path/to/input.tiff --output /tmp/foo.obj

Importing into Blender

Blender supports OBJ file import as a standard, no need for plugins.

In Blender, File > Import > Wavefront (.obj)

You need to make sure your settings look like this. I'm using Y forward, Z up. You can save these settings and recall them later using the [+] button. You'll need to reselect this each time you import a new mesh!

Screen Shot 2017-11-18 at 20.53.54.png

Styling Relief Maps in Blender

There are lots of ways of styling in Blender. I find the traditional way (lighting from the top-left, and viewing straight down with an Orthogonal camera) works really well.

This is the Nodes recipe I use (with Cycles renderer). This applies a relief-map colour ramp to the z-axis to give cartographic-looking results...

an example configuration

Troubleshooting

Out-of-memory (Blender crashes). I find a 4Gb machine can handle 1000x1000 fine; 8Gb RAM is good for 2000x2000. If your raster is too large Blender will crash during loading. You can use QGIS or gdalwarp to resize your raster to something your machine can cope with.

If the tool runs, and Blender imports it but you can't see anything, it may be that the scale means the mesh is too large (clipping means that very large models can disappear). You can pass this in using the --scale VALUE option. The default is 1.0, which means 1 unit of size in the DEM equals 1 Blender unit. You can scale it up/down in Blender, or run it again with this option.

If the model imports but is rotated on its side, check you imported using the suggested settings. Remember, y forward, z up.

Finding a point of Inaccessibility using QGIS and FOSS tools

Finding points of inaccessibility

Question : What is the point on the Scottish mainland which is furthest from a building (as the crow flies)? This is something called a point of inaccessibility

This is something that is possible to do with rasters, in particular the gdal_proximity tool.

This was done using QGIS 2.18.10

Data sources

I used Ordnance Survey Open Data for this. In particular, a shapefile courtesy of Alasdair Rae. This has all building outlines clipped to the border. I recommend getting this into sqlite or postgres before attempting further processing - it'll save hours of processing time ;-)

Also used a sqlite database with an outline of Scotland for clipping. I created this a long time ago, using data from openstreetmapdata.com, namely the split land polygons, and dissolved them. This is optional, but recommended, as most of the furthest points from buildings are in the North Sea. You could also use Natural Earth data 1:10 admin level 1 for this... the coastline is much coarser, though.

Process

Throughout, I've used EPSG:27700, the standard UTM projection for the UK, in meters.

1. shp2pgsql to load buildings into Postgres table

2. Get centroid layer using DB Manager. This took 4 secs; I gave up on the Shapefile after 5 minutes when it hadn't even got around to 1% done :-) The "1 as one" field will be used later, it is needed!

select st_centroid(geom) as geometry, 1 as one from buildings;

3. Save layer as a permanent OGR source (I used .sqlite). This is needed as the following step doesn't (currently) work with virtual layers.

4. Rasterise to 100m resolution grid. By using the "one" field, all 100x100m cells with a building take the value 1, and all other cells are 0. (I could also use the -burn option instead)

Dundee and environs, as a binary raster:-

Dundee, with 100x100 cells. Cells with buildings are white.

Dundee, with 100x100 cells. Cells with buildings are white.

5. Proximity. We want to find the distance to the nearest "one" pixel.

gdal_proximity.py /tmp/raster100m.tiff /tmp/proximity100m.tiff -values 1 -distunits GEO -of GTiff

if all goes well, we should see something like this:-

proximity raster

proximity raster

6. Clip proximity to a cutline so that sea pixels are set to 0

gdalwarp -q -cutline /tmp/scotland_dissolved.sqlite -tr 100.0 100.0 -of GTiff /tmp/proximity100m.tiff /tmp/proximity100m_clipped.tiff

7. Contourise (Extraction > Contours) using 200 meter spacing.

8. Use rule-based rendering to hide contours less than 6km from a building. This instantly shows candidates - only a handful of areas of the mainland (and several remote uninhabited islands) show up.

This image shows contour lines > 6.6 km from the nearest building:-

areas more than 6600 meters from a building

areas more than 6600 meters from a building

And the winner is...

I make the point of inaccessibility somewhere around OS grid square NN8682, which is ~ 7150m away from the nearest 3 buildings. The limited resolution of the grid (100m) rules out anything more accurate than that :) According to the Lat Lon tools plugin the coords are approximately latitude 56.9211, longitude -3.8661 (rounded down to reflect lack of accuracy)

This is close just to the NW of Beinn Bhreac

it's in the Cairngorms, there's a surprise :-)

it's in the Cairngorms, there's a surprise :-)

To get the circle, I used a single point layer in 27700, and used a Geometry Generator symbol using

buffer($geometry,7150)

Other reading

This came out of asking a question on GIS Stackexchange.

One feature of these points is that the lie on the intersection points of the voronoi mesh, which provides another technique for finding maxima. In this case the number of points was quite high (300k+) so I decide to go for a more pragmatic approach :)

Detecting adjacent polygons using QGIS Virtual Layers

Challenge:-

- to take a series of polygons in a shapefile
- find which pairs of polygons share a border.
- draw a line between the centroids of adjacent polygons.
- Do this without manually loading into a database.

This is possible in QGIS using DB Manager with Virtual Layers. This uses Spatialite under the hood - I'm more used to PostGIS, but many of the functions have the same name.

finding polygons which share a border, excluding those which touch at a point

finding polygons which share a border, excluding those which touch at a point

Using the following query in the DB Manager SQL window, and adding the view as a new layer...

select
    a.gid,
    b.gid,
    make_line(st_pointonsurface(a.geometry), st_pointonsurface(b.geometry)) as edge
from
    areas20170628201438984 as a,
    areas20170628201438984 as b
where
    a.gid != b.gid and
    a.gid < b.gid and
    st_relate(a.geometry,b.geometry,'FF2F11212');

The st_relate value of FF2F11212 will find polygons whose intersection is a line. This means that polygons which intersect at a point are ignored.

The a.gid < b.gid makes sure there's only one line between each polygon pair. (Similar to filling in one diagonal half of a matrix)

st_pointonsurface is similar to st_centroid, but guarantees that the points are inside the polygon... given the irregular shapes of UK Census Output Areas, this is more useful than st_centroid.

I find virtual layer queries like this are painfully slow to render when the source layer is a shapefile; it's probably better to do this on a layer from a database ;-)

TIP: I also found that I had to run Execute (F5) twice to allow the 'edge' geometry field to be selected when adding the layer. (QGIS 2.18.7).