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

Well, really just metes…

Anyway, I had a description of the property, which can be hard to follow for most of us. Here’s an example description of part of Buttonbush Nature Area:

Commencing at the South ¼ corner of Section 10, T2S, R6E, City of Ann
Arbor, Washtenaw County, Michigan; thence N01°57’53”W 60.00 feet
along the North-South ¼ line of said Section 10 and along the East line of
Foxfire Condominium, Washtenaw County Condominium Subdivision Plan
No. 136 for a PLACE OF BEGINNING; thence continuing N01°57’53”W
2208.89 feet along said North-South ¼ line and along the East line of
Foxfire Condominium, Washtenaw County Condominium Subdivision Plan
No. 136 and along the East line of Fox Ridge Commons Condominium,
Washtenaw County Condominium Subdivision Plan No. 176; thence 96.56
feet along the arc of a 7442.41 foot radius circular curve to the right, chord
bearing N81°07’45”E 96.56 feet along the South right-of-way line of US-23
(variable width); thence S07°00’18”E 195.59 feet; thence S48°50’07”E
5.51 feet; thence S01°57’53”E 2022.40 feet; thence S88°16’54”W 117.07
feet along the proposed north right of way line of Dhu Varren Road (60
feet proposed half width) to the Place of Beginning, being a part of the
Southeast ¼ of said Section 10 and containing 5.90 acres of land, more or
less.

My eyes glazed over the first time I read that. But after looking at it for a while, I realized that i only needed a little bit of that to figure out the boundary. Most of it was just prose.

To get the data I needed, the first step was to eliminate all the prose, and put it in a standard form that could be more easily read in. I copied the text and just started deleting the extra words. This became:

# Commencing at the South 1/4 corner of Section 10, T2S, R6E,
# City of Ann Arbor, Washtenaw County, Michigan
# Containing 5.90 acres of land
N01,57,53,W 60.00 feet
PLACE OF BEGINNING
N01,57,53,W 2208.89 feet
Arc 7442.41 foot radius
N81,07,45,E 96.56 feet
S07,00,18,E 195.59 feet
S48,50,07,E 5.51 feet
S01,57,53,E 2022.40 feet
S88,16,54,W 117.07 feet
Place of Beginning

This is almost entirely just deleting unnecessary parts of the description. It’s easier to do that manually than to write a program that does it. What’s left is the description of the starting point in comments in the first couple lines, a pair of lines specifying the starting/ending corner of the boundary, and one additional line where the boundary is a circular arc rather than a straight line. I did change the degree, minute, and second symbols to commas [Edit: No longer necessary!], just so I didn’t have to search for each symbol separately. (Using spaces would also work.)

The comment lines at the top aren’t used. I kept the description of the starting location, and also moved the “Containing 5.90 acres” description there, because I like keeping that information.

When I read in the data, I look at the first character of each line to determine what the line is. N, S, A, P, and # are recognized, and a line beginning with any other character is ignored. N or S means that line is one part of the boundary, usually a straight line. A for Arc means the following line is a circular arc, and gives the radius of the circle. I report the maximum deviation from a straight line. P for Place of Beginning denotes the beginning/ending corner of the property boundary. Any lines prior to the first Place are just to get from a known reference location to a corner. The segments between the two Places should form a closed loop, with only a tiny error. If they don’t, something is wrong in the description, either in the source document, or in the process of converting it to the standard form.

I should add the caveat that I don’t know if the way the prose is written is fairly standard, or if it varies between jurisdictions. For example, Place of Beginning might be a different term elsewhere.

I’ll post the code in the following entry.