How to find the point in a country furthest from a boundary?

i'm trying to find the point in south africa furthest from the sea or a landward boundary.

my method thus far:

1. i obtained a first guess at the centre of the inscribed circle and the three segments of boundary this circle might touch using a compass and 1:2 500 000 map.

2. i then clipped these three segments of boundary from a data set of boundary coordinates supplied by national geo-spatial information (WGS84; decimal degrees; 9 decimal places; roughly one point for every 40 metres of boundary).

3. i converted the boundary coordinates to radians and used cartesian equations (distance between two points and centre of a circle through three points) to iteratively improve the centre of my circle and the three boundary points it touches. (obviously, i'm ignoring curvature; i'm effectively finding the centre of the base of the spherical cap defined by the inscribed circle -- i.e. a point approximately 16km underground.)

my questions:

1. is it reasonable to plug lat-long coordinates expressed in radians into cartesian equations to find the centre of the base of the spherical cap defined by the inscribed circle (given that the surface distance represented by one radian varies with latitude, and that the base of the spherical cap in question is roughly 800km in diameter and spans latitude 30S)? what order of error might be involved here?

2. is it reasonable to assume that the three boundary points tangential to the base of the spherical cap will likewise be the three boundary points closest to the actual point in the country furthest from a boundary?

3. i have the feeling i'm not even asking the right questions here. i guess there must be a better way to find the point in a country furthest from a boundary ...

cheers, william

• Maybe you can use a negative buffer. Decrease the buffer until there is no area left and use the centroid of the last positive area. In POstGIS it should be possible to do that with a recursive query. You might get problems that you get many small areas instead of just one. – Nicklas Avén Jul 11 '12 at 15:16
• This question is asked at gis.stackexchange.com/questions/13871/… and answered there with a raster solution. – whuber Jul 11 '12 at 15:32
• thank you for the input. i'm going to have to learn more about gis software to fully understand some of these answers (which i will certainly do). several avenues have opened up ... – user8722 Jul 11 '12 at 19:50

I think Voronoi polygons are the best approach. I enclose a Matlab function that does this quite efficiently. As input I use a set of border coordinates, x and y, including enclaves (i.e. Lesotho). In addition I include the external border (boundary).

``````function [xfar,yfar,distance,nearest,polygons] = farthest(x,y,boundary)
%FARTHEST Find the position that is farthest from any points in (x,y)
%  and also inside or at the boundary
%  Input:
%   x,y: arrays of coordinates (must be within the boundary)
%   boundary: [xb,yb] where (xb(i),yb(i)) is corner no i in the boundary
%     polygon.
%     default: [0 1 1 0;0 0 1 1]'
%  Output:
%   xfar,yfar : coordina4tes of most distant position
%   distance : distance from nearest point
%   nearest  : nearest point is (x(nearest),y(nearest))
%   polygons : cell array of Voronoi polygons

% Author Are Mjaavatten

if nargin < 3
boundary = [0 1 1 0;0 0 1 1]';
end
if boundary(1,:) == boundary(end,:)
boundary(end,:) = []; % Convert to open polygon if needed
end

% Add new points to make sure the convex hull is outside
% the boundary.  The new points must be distant enough that the
% midpoints between these and original points are outside the boundary
factor = 5;
center = [mean(x),mean(y)];
external = boundary;
for i = 1:length(boundary)
external(i,:) = center + (boundary(i,:)-center)*factor;
end
pointlist = [x(:),y(:)];
X = [pointlist;external];
DT = DelaunayTri(X);
[V,C] = voronoiDiagram(DT);
npoints = length(C);
distance = 0;
polygons = cell(npoints,1);
for i = 1:npoints
poly = V(C{i},:);   % Voronoi polygon around point i
for j = 1:size(poly,1)
if inpolygon(poly(j,1),poly(j,2),boundary(:,1),boundary(:,2))
d = norm(pointlist(i,:)-poly(j,:));
if d > distance
distance = d;
xfar = poly(j,1);
yfar = poly(j,2);
nearest = i;
end
end
end
polygons{i} = poly;
end
end
``````

Testing, simplistically using longitude and latitude coordinates, I get, for about 2000 points:

``````xfar =
23.3188
yfar =
-30.2397
distance =
3.7503
``````

The handling of enclaves is not entirely robust, but will work in most practical cases. A robust, but slower, version named remotest has been submitted to the Matlab File Exchange.

To be honest, it sounds like you know what you're doing more than me, but my approach would be:

• Rasterize geopoltical/sea boundaries when the polylines are 1 and everything else is 0.

• Use the Euclidean Distance tool in Arc to produce raster of 'distance to the nearest 1' (i.e. nearest boundary/sea pixel)

• Follow the advice from this Stackexchange explaining how to get the location of the pixel with the highest value in a raster.

Here's another approach, if you're going to be doing this manually rather than algorithmically, which keeps things in a vector form:

I'm working out of QGIS.

• Extract the nodes of the country polygon.
• Create voronoi polygons from these nodes.
• Follow the acute angles until you find an intersection that has bordering polygons extending outwards to split the country roughly in three.

I'm not 100% sure how to perform this last step without human intervention (perhaps something with measuring the principal axis?). I can say that it will be something other than the centroid of the polygon in most cases.

Disclaimer: Determine very deliberately how far upstream into a bay/river you should go before you consider yourself "in the country". Also won't work for countries where a node-less line segment is one of the three nearest boundaries to the Farthest Point.

EDIT:

Instead of Voronoi polygons, use the related Delaunay Triangulation tool, and select the triangle with the largest area. Take the centroid of that triangle, and you will have your point.

EDIT2:

I'm not 100% certain if the 'largest area' criteria will always align with the farthest point, mathematically. Interesting problem... In the end, if you want a strong scriptable solution you may need to use iterative buffers, or resort to the aforementioned raster.