OpenStreetMap

Code for reading a survey description file

Posted by BaddDadd on 11 September 2020 in English (English)

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] = readsurvey(boundaryfile);
%  [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 = [0];
y = [0];
npoints = 1;
ipb = [];

nlines = 0;
line = fgetl(fid);

arcrad = 0;
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
    arcrad = data(1);
%
%   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 ', ...
          'radius ', num2str(arcrad), cstr])
    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;
    dev = arcrad - sqrt(arcrad^2 - l^2);
    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);
      e = arcrad - sqrt(arcrad^2 - 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)
    arcrad = 0;
  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

Login to leave a comment