I'm looking through a variety of land surveys from construction permits in Florida, and I'm trying to generate a set of lat/lon points from them, in order to map them out. The surveys list a PLSS reference as a starting point, and give bearing and distance all around the shape.
From what I've read, I think I should be able to use the Vincenty formula to generate points based on that information. I've got PLSS section corner points from Geocommunicator to use as starting points, and I'm using the WGS-84 constants for the formula.
I've tried Geopy's implementation as well as a Python port of a JavaScript implementation, and they both suffer from the same problem: the overall shape looks right, but it's not aligned correctly. Here's what it looks like in JOSM, overlaying the resulting shape (including the path from the PLSS point to the survey's Point of Beginning) on aerial photography for reference.
Why is this happening? The shape looks right and the size fits, but it's mis-aligned. My first thought is that I had wrong PLSS data for my starting point, but the alignment problem isn't uniform. It's usually off by within 200-300 feet, and mostly east of where it should be, but sometimes a little west instead, and the north/south offset varies as well, but usually it's not as dramatic as the east/west misalignment.
Is it possible that I just need better PLSS points to start from? Or should I be using some other approach to calculate the points, besides the Vincenty formula? Do I need to reproject something? I'm still learning all the tools here, and I'm sure I'm just missing some step of the process.
Update:
In response to some requests for more information, the PLSS shapefile I have is in NAD83 coordinates (prj2epsg identifies it as EPSG:4269), which identifies the corners of each section within the county. I then used simple arithmetic to identify the midpoints where appropriate, as some of the surveys use those instead of the outer corners. I've tried reprojecting into EPSG:4326 using ogr2ogr, but it doesn't seem to change the coordinates at all.
Also, the survey data consists of bearings and distances, which I think is different from metes and bounds. For example, the permit identified by the above image is as follows:
Point of reference: South quarter corner of Section 11 of Township 24 South, Range 27 East off the Talahassee Meridian. According to my PLSS shapefile and midpoint math, that's 28.4056197, -81.58239195.
From there, here are the measurements:
N 89° 41' 52" E 237.44 ft
N 00° 00' 00" E 31.87 ft
Point of Beginning
N 74° 54' 56" W 160.57 ft
N 15° 05' 04" E 14.09 ft
N 74° 54' 56" W 239.99 ft
N 15° 05' 04" E 92.44 ft
S 74° 54' 56" E 213.61 ft
N 15° 05' 04" E 13.29 ft
S 74° 54' 56" E 36.43 ft
S 15° 05' 04" W 27.38 ft
S 74° 54' 56" E 150.53 ft
S 15° 05' 04" W 92.44 ft
I've also tried making sure I'm using survey feet, rather than standard feet, but the difference is minuscule at these distances.
Also, the image comes from Bing aerial photography, so I assume it's EPSG:3857. JOSM uses that as a fairly standard layer for tracing, so I also assume it's done whatever work is necessary to line it up with lat/lon points properly.