# Finding a pole of inaccessibility relative to multiple points without using rasters

I want to find a 'point of inaccessibility' from a series of points, representing towns and villages in Scotland. Using FOSS GIS tools (QGIS, SAGA, GRASS, Postgres/PostGIS ... ).

This is the point furthest from any centre of population. This would be the centre of the largest circle which doesn't contain a town, but it will have a town on its edge.

There's a raster approach, which works - the GDAL Proximity filter.

Is there a VECTOR based way to do this? One that doesn't require rasterisation?

Here's my attempt with rasters...

• Start with a point layer in a suitable projection in meters (I'm using OSGB, 27700)
• Add a field called "one", and set to integer with constant value 1.
• Vector to raster , editing command to use -a one and using GEO rather than pixel. Make sure pixels are square!
• Create contours on the proximity raster (in this case, 1 km wide)

Here's what it looks like. I drew some candidate points by hand ;-)

I thought that dissolving repeated buffers (using the Multi ring buffer plugin) might work, but I suspect I don't have enough memory for this to work.

## 2 Answers

I guess it is the centre of triangle excircle with largest radius, that touches no more than 3 points. In the picture below first 11 such centeres shown. They are labelled by their ranking number.

It is enough to weed out ones that are outside triangles and define the champion, i.e. No3 in the picture.

UPDATE INSPIRED BY STEVEN FINDING:

Result above unlike first solution obtained in ArcGIS without scripting (but with advanced license) and it works! However it can produce questionable results on the edges of a study area, e.g. point No 1, i.e. second remotest point found...

• Wow, thanks sharing this! I tried to follow your workflow in QGIS by (1) Points to Delaunay (`Vector | Delaunay Triangulation`), (2) Select champion triangle (Calculated `\$area` on the attribute table of Delaunay polygons) and saved the largest one as a new layer (3) Find circumcenter (`Processing | GRASS | Vector | v.voronoi.skeleton` on the champion triangle). Do you think this method is correct? – Kazuhito Oct 1 '17 at 15:33
• First step is correct. I computed exradius centres for ALL triangles. To find centre I constructed inward perpendicular for midpoints of all 3 sides and find intersection ( combination 3,2). – FelixIP Oct 1 '17 at 19:21
• Circumscribed circle of triangle, not excircle. – FelixIP Oct 1 '17 at 19:39
• Thanks so much! (and thanks for nice picture, too. The more I lean the background about it, the more it becomes fascinating.) – Kazuhito Oct 1 '17 at 22:14
• Nice observation! I'll test it and update the answer. Unlike first solution, it will require no scripting at least in Arcgis. – FelixIP Oct 3 '17 at 20:50

Here's a follow-up. Thanks to @FelixIP for pointing me in the right direction!

Using the OSM data from Australia, I was able to find the "point of inaccessibility" on the Australian Mainland - I make it around 260km equidistant from Akarnenehe, Bedourie, and Mount Dare, at POINT(137.234888 -24.966466)

I found a fairly easy workflow in QGIS which uses a combination of the raster and vector approaches. I'm sure a similar approach would work in other GISes.

The 'a-ha' moment came from noticing this

The maxima all fall on tri-points on the voronoi mesh - points where adjacent triplets of voronoi polygons meet.

The approach was as follows :-

• work in a meter-based projection. I used 3857, not ideal but my OSM data was in that format :)
• get a layer representing points of interest place in ('city','village','town','hamlet')
• create a voronoi mesh from this layer
• use extract nodes to get the tripoints
• create the proximity raster (as shown in the question)
• clip the raster to the outline of the landscape, so sea pixels are set to 0.
• use Point Sampling tool on the extracted nodes against the raster

Then use Db Manager and Virtual Layers to find the node with the largest distance value on the raster.

``````select
rowid,
name,
proximity2,
st_astext(st_transform(geometry,4326)) as pt,
st_buffer(geometry, proximity2) as geometry
from
"samples"
order by
proximity2 desc
limit 1;
``````