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 :)

Toponymy - mapping place name patterns in Great Britain with QGIS

This project was to study the distribution of certain patterns in place names of Great Britain. To do this, I used the Ordnance Survey Open Gazetteer (crown copyright and database rights). These maps were created in QGIS.

Place names in Britain come from a number of sources; Old English, Brythonic (Old British), Old Norse (via the Vikings), Pictish and Gaelic. Wikipedia has a write up here; I used this as inspiration for the queries.

The maps use a hexagonal grid of approx 20x20km. The colour represents the number of places with that pattern within the cell.

Danelaw and the Viking Influence

This map shows the distribution of places ending in -by, -thorpe, -thwaite, -toft, -kirk, -ness. In England, it shows the diagonal band from Cumbria to Lincolnshire which was under DaneLaw. In Scotland, it shows most influence in Orkney and Shetland. Most of the remaining ones in Scotland are 'kirk'. I did a map like this a few years back, but the cartography was pretty basic ;-)

places with Viking imports

Legend

These maps all use the following legend. I originally used Jenks Breaks to categorise the regions; they almost always came out with a distribution of 1-2, 2-4, 4-8 etc. So that's what I used.

legend.png


Pictish place names

One of the few Pictish place names I could find was -carden- (or -cardine). For example, Kincardine.

-carden or -cardine. The ones on the Welsh border are probably not of pictish origin!

Places beginning with Llan-

If there's an obvious Welsh place name pattern, it's the prefix Llan- (meaning Church). This is borne out by the map, which shows almost all such places in Wales (or near the English/Welsh border)

places beggining with Llan-

Influence of Brythonic - combe, coombe, cwm

Brythonic is one of the old British languages. Its word for valley is spelt differently; -coome or -coombe in England, cwm in Welsh. The sphere of influcence is felt most strongly in Wales, Devon and Somerset. Brythonic names are certainly found as far north as Edinburgh (for example, Penicuik is derived from Pen Y Cog, a name which wouldn't look out of place in Wales).

-combe, -coombe, -cwm, Brythonic

-combe, -coombe, -cwm, Brythonic

Making of the Map

In QGIS, I used MMQGIS to make the hexagonal grid. I removed cells which didn't overlap land, and manually removed some more hexes for aesthetic reasons (e.g. to separate the Western Isles from the mainland).

Each map was done by querying the gazetteer, saving as a new points layer, and doing a points-in-polygon analysis of this layer with the hex grid.

Finally, I buffer/dissolved the hex cells into a new layer to get the outline. The dashed outline was from applying Multi Ring Buffer plugin to this layer.

 

Making a straightened railway strip map

This map shows the first few Km of the recently re-opened Borders Railway which runs between Edinburgh and Galashiels. This was created using data from OpenStreetMap (OSM), data copyright OpenStreetMap contributors.

Straightening the map

To save space, the main railway line has been straightened, and the space around it has been distorted.. Strip Maps like this are nothing new; centuries ago, they were used to illustrate stagecoach routes and roads between towns.

"Up" in this map changes direction as you read from left to right; it goes from N in Edinburgh, to NE in the middle, and E at the right.

Preparing the data

In QGIS, I fetched the OSM relation for the Borders railway using QuickOSM. This was buffered by 1km, and the buffer was used to clip the OSM road network. Nodes were extracted (using Extract Nodes) and all of these were loaded into Postgres using shp2pgsql.

Warping the map

The points were warped/reprojected using PostGIS. The basic methodology is shown in the diagram below.

path7858.png
  • The shortest line between each point and the railway is drawn (L)
     
  • Where that line meets the railway, we measure the cumulative distance from the start of the railway (A). This is the x coordinate.
     
  • The y coordinate is the length of L.
     
  • The azimuth (bearing) of L is computed. For certain ranges of angles, we negate the Y coordinate. This has to be done, otherwise all points end up on one side of the railway!

These are fairly easy to compute in PostGIS, using ST_ShortestLine() and ST_LineLocatePoint().

Problems Encountered - OSM and Geometry

One of the things that took longest was trying to get it to work with the existing OSM route.

Eventually I took a pragmatic approach, and traced over the OSM route manually using QGIS.

The reason? OSM and geometry.

This is not a criticism of OSM mappers, or of OSM. OSM data often needs processing to make it amenable to coded approaches like this one.

I mention this as it's tripped me up a few times before. But there are always ways around it. Here were the hurdles I had to overcome...

  • There were gaps (mostly caused by tunnels), which meant the route was broken into segments.
     
  • These segments were not always aligned head-to-tail, and varied in direction. This made it difficult to create a single, unified linestring.
     
  • This railway has been under construction for a few years now, so it will have been mapped incrementally. OSM tends to reflect "the truth on the ground".
     
  • There are no rules in OSM about geometry. Mappers don't need to worry about whether they draw a line from North to South, or from South to North. They're free to draw a polygon clockwise or anti-clockwise. That's not their concern, nor should it be.

In the end, I used the pragmatic approach, and traced a single line over the OSM route using QGIS. This made a single line segment with consistent direction. There are QGIS Plugins ("Join Lines", "QChainage") which I've used in the past to fix such problems.

Finally, thank you to the OSM mappers who made this map possible, and in such a timely manner - the route was mapped by the time the line was opened.

Credits

Uses data copyright OpenStreetMap contributors.

Font used is "Transport Heavy" by Nathaniel Porter, based on the UK Road Sign font (free for non-commercial use)

See a larger version

You can see/download a larger version on my flickr stream.