# OpenStreetMap  ## Code for reading a survey description file

This is the code referred to in my previous entry. This is written for Matlab. It should work in Octave as well, but I haven’t tested it.

Edit: Updated the code to return arc segments split into two or four pieces, depending on deviation from a straight line. Also, code now supports many bearing delimiters for degrees, minutes, and seconds.

``````function [x, y, err] = readsurvey(file)
%   Read a survey description file.
%
%  [x, y, err] = readsurvey(boundaryfile);
%
%  boundaryfile is the survey description file.
%
%  x, y are the east and north coordinates
%  err is distance between first and last points of closed loop, useful as an
%    accuracy check.
%
%   # means comment line, e.g. starting reference location
%   N means line direction referenced to North
%   S means line direction referenced to South
%     W,E mean angle is to West or East of North or South.  Numbers are
%     degrees,minutes,seconds, then length.  Many delimiters are recognized.
%   P means Place of Beginning, a corner of property.  Should be listed twice.
%     Lines before first instance are from reference location not on boundary.
%   Arc means line gives the radius for a circular boundary, and next line is
%     the chord for that boundary.  Up to four segments will be returned for
%     circular arcs, depending on deviation from a straight line.

%
%   Delimiters that may be used to specify or separate degrees, minutes, and
%  seconds of direction angle.
%
delimiters = {'''', '"', '`', '°', '’', '”', ','};

fid = fopen(file, 'r');
if (fid == -1)
disp(['Error opening ', file])
disp(['Current directory is ', pwd])
return
end
%
%   Set cutoffs for making additional points along arcs.
%
devmin1 = 2;
devmin2 = 8;
devrel = 0.05;

x = ;
y = ;
npoints = 1;
ipb = [];

nlines = 0;
line = fgetl(fid);

while((ischar(line)) & (nlines < 100000))
nlines = nlines + 1;

if ((strcmp(upper(line(1)), 'N')) | (strcmp(upper(line(1)), 'S')))
%
%   Determine angle direction.
%
ii = find(upper(line) == 'W');
jj = find(upper(line) == 'E');
if (isempty(ii))
if (isempty(jj))
error('Could not determine angle direction')
else
angsign = -1;
line(jj(1)) = ' ';
end
else
if (isempty(jj))
angsign = 1;
line(ii(1)) = ' ';
elseif (ii(1) < jj(1))
angsign = 1;
line(ii(1)) = ' ';
else
angsign = -1;
line(jj(1)) = ' ';
end
end
end

if (strcmp(upper(line(1)), 'N'))
%
%   Can recognize minutes and seconds, but not degree symbols, or other weird
%  symbols.  Make any found just be spaces instead.
%
line = regexprep(line, delimiters, ' ');
line(1) = ' ';

data = sscanf(line, '%f');
if (length(data) ~= 4)
warning('Should be exactly four data values...')
end

ang = data(1) + data(2)/60 + data(3) / 3600;
ang = 90 + angsign * ang;
len = data(4);

dataline = 1;
elseif (strcmp(upper(line(1)), 'S'))
%
%   East and west have opposite sign meanings for southerly boundary lines.
%
angsign = -angsign;

line = regexprep(line, delimiters, ' ');
line(1) = ' ';

data = sscanf(line, '%f');
if (length(data) ~= 4)
warning('Should be exactly four data values...')
end

ang = data(1) + data(2)/60 + data(3) / 3600;
ang = 270 + angsign * ang;
len = data(4);

dataline = 1;
elseif (strcmp(upper(line(1)), 'P'))
%
%  Place of beginning
%
ipb = [ipb,npoints];
dataline = 0;
elseif (strcmp(upper(line(1)), 'A'))
%
%   Next line is the chord of an Arc.  Read the radius now.
%
line(1:3) = ' ';
data = sscanf(line, '%f');
if (length(data) < 1)
warning('Should be at least one data value...')
end
%
%   Find direction arc is curving.
%
ii = strfind(lower(line), 'right');
jj = strfind(lower(line), 'left');
if (~isempty(ii) & isempty(jj))
curve = 1;
cstr = ', curving to the right.';
elseif (~isempty(jj) & isempty(ii))
curve = 2;
cstr = ', curving to the left.';
else
curve = 0;
cstr = ', curving in unknown direction.';
end

disp(['Segment ', num2str(npoints), ' is a circular boundary line of ', ...
dataline = 0;
else
%
%   Comment or unrecognized line.
%
dataline = 0;
end

if ((arcrad > 0) & dataline)
%
%   This is the chord of a circular arc.  Report the maximum deviation from
%  straight.
%
l = 0.5 * len;
disp(['Maximum deviation from straight line = ', num2str(dev)])

if ((dev > devmin2) | (dev > len * devrel))
%
%   Get two additional points along arc.
%
h2 = 0.25 * (l^2 + dev^2);
h = sqrt(h2);
er = l / (2 * h);
el = dev / (2 * h);
else
e = 0;
end
else
dev = 0;
end

if (dataline & (dev > 0) & (curve > 0))
dx = len * cosd(ang);
dy = len * sind(ang);

if ((dev > devmin2) | (dev > len * devrel))
%
%   Add three points on circular boundary if deviation is large.
%
ex = dev * cosd(ang+90);
ey = dev * sind(ang+90);
fx = e * cosd(ang+90);
fy = e * sind(ang+90);

if (curve == 1)
%
%   Points along the curve.
%
xmid = x(npoints) + 0.5 * dx + ex;
ymid = y(npoints) + 0.5 * dy + ey;

xend = x(npoints) + dx;
yend = y(npoints) + dy;

x1q0 = 0.5 * (x(npoints) + xmid);
y1q0 = 0.5 * (y(npoints) + ymid);

x3q0 = 0.5 * (xend + xmid);
y3q0 = 0.5 * (yend + ymid);

x(npoints+1) = x1q0 + (fx * er - fy * el);
y(npoints+1) = y1q0 + (fy * er + fx * el);
npoints = npoints + 1;

x(npoints+1) = xmid;
y(npoints+1) = ymid;
npoints = npoints + 1;

x(npoints+1) = x3q0 + (fx * er + fy * el);
y(npoints+1) = y3q0 + (fy * er - fx * el);
npoints = npoints + 1;

x(npoints+1) = xend;
y(npoints+1) = yend;
npoints = npoints + 1;
else
xmid = x(npoints) + 0.5 * dx - ex;
ymid = y(npoints) + 0.5 * dy - ey;

xend = x(npoints) + dx;
yend = y(npoints) + dy;

x1q0 = 0.5 * (x(npoints) + xmid);
y1q0 = 0.5 * (y(npoints) + ymid);

x3q0 = 0.5 * (xend + xmid);
y3q0 = 0.5 * (yend + ymid);

x(npoints+1) = x1q0 - (fx * er + fy * el);
y(npoints+1) = y1q0 - (fy * er - fx * el);
npoints = npoints + 1;

x(npoints+1) = xmid;
y(npoints+1) = ymid;
npoints = npoints + 1;

x(npoints+1) = x3q0 - (fx * er - fy * el);
y(npoints+1) = y3q0 - (fy * er + fx * el);
npoints = npoints + 1;

x(npoints+1) = xend;
y(npoints+1) = yend;
npoints = npoints + 1;
end
elseif (dev > devmin1)
%
%   Add midpoint of circular boundary if deviation is small but greater than
%  devmin1.
%
ex = dev * cosd(ang+90);
ey = dev * sind(ang+90);

if (curve == 1)
x(npoints+1) = x(npoints) + 0.5 * dx + ex;
y(npoints+1) = y(npoints) + 0.5 * dx + ey;
npoints = npoints + 1;
else
x(npoints+1) = x(npoints) + 0.5 * dx - ex;
y(npoints+1) = y(npoints) + 0.5 * dx - ey;
npoints = npoints + 1;
end
else
x(npoints+1) = x(npoints) + dx;
y(npoints+1) = y(npoints) + dy;
npoints = npoints + 1;
end

curve = 0;
elseif (dataline)
dx = len * cosd(ang);
dy = len * sind(ang);

x(npoints+1) = x(npoints) + dx;
y(npoints+1) = y(npoints) + dy;
npoints = npoints + 1;
end

if ((arcrad > 0) & dataline)
end

line = fgetl(fid);
end

if (nargout >= 3)
if (length(ipb) ~= 2)
disp('Assuming end point of data is also the Place of Beginning.')
elseif(ipb(2) ~= npoints)
disp('Warning, second Place of Beginning is not end of data.')
end
ipb1 = ipb(1);
%
%   Third argument is error at end point.
%
err = sqrt((x(end) - x(ipb1))^2 + (y(end) - y(ipb1))^2);
end
``````

Following are the results returned when run using the data in my previous post, saved as file ‘buttonbush.dat’. Positive X is east of the starting reference point, and negative X is to the west. Similarly, positive Y is north of the starting point. The Place of Beginning is the second point, and it’s replicated as the last point, with only a small error, around 1/15 of an inch. Also, the output returned for a circular arc is shown. For this short piece, the deviation from a straight line is only a couple inches.

``````>> [x,y,err] = readsurvey('buttonbush.dat')
Segment 3 is a circular boundary line of radius 7442.41, curving to the right.
Maximum deviation from straight line = 0.1566
x =
0       -2.057      -77.787       17.618       41.471       45.619       114.96      -2.0617
y =
0       59.965       2267.6       2282.4       2088.3       2084.7       63.479       59.968
err =
0.0056768
``````