I'm trying to replicate a nasty corner case in a LRS implementation that doesn't agree with ArcMAP measure tool. I've tried measuring distance in its native shp coordinate system NAD83 (27916) to no avail. (basically same results)
My question is how are the distances calculated along a line for the LRS system. My current distance function (for 4326) is the following.
EARTHRADIUS = 3963.190592 # miles ''' used in distance ''' def deg2rad(d): return d * math.pi / 180.0 # in miles def distance(pt1, pt2): dLat = deg2rad(pt2 - pt1) dLng = abs(deg2rad(pt2 - pt1)) dLat2Sin = math.sin(dLat / 2) dLng2Sin = math.sin(dLng / 2) a = dLat2Sin*dLat2Sin + math.cos(deg2rad(pt1))*math.cos(deg2rad(pt2))*dLng2Sin*dLng2Sin return 2.0 * EARTHRADIUS * math.atan2(math.sqrt(a), math.sqrt(1-a))
For most geometries the distance between the ESRI LRS and mine are nominal but this one is different. I've dealt with larger linestrings before (30+ miles) with no issues. Previously these other larger linestrings have worked fine throughout the alignment, however with this edge case I'm getting ~.5 miles difference (my implementation's distances are always larger) between the two implementations by the end.
Is distance calculated a different way via the LRS implementation, is my geometry duplicating previous points possibly, or is it possible coordinate system error is actually sneaking in?
For the ground truth measurement, I'm just hitting a backend ESRI LRS server I think? After playing with the code from this webmap. Basically all I'm doing is iterating through points in the subject linestring, adding my distance implementation to the distance along the line for each point then comparing the two. I did modify the API parameters for the lowest tolerance as I'm using the actually constituent point geometry and the shp's original coordinate system. I have no idea about how the LRS was implemented to be honest, other than the units are miles.
I tried to debug this for quite some time yesterday, I've come to conclusion I don't think its within the distance formula (although possible) more so some abstract set of calibration points that ESRI's LRS is defined by. The road's geometry represents a parallel set of vector for a good portion of it, and apparently some LRS systems use that original centerline to derive the LRS???
On top of that I noticed points are very dense multiple on top of one another (never shooting back though) I think underlying simplification may contribute as well. Also this geometry is about to be changed massively due to construction of a new roadway, so I think a small change to the real feature maybe potentially didn't bubble up to the geometry that the LRS is based on for this roadway. However I am sure I have the most up to date data set as it was dumped just to debug the same issues in ArcDesktop.
The real confirmation of my suspicion is when I dropped the raw shp into Earth Pro with the original CRS and opened up the elevation profile, and it had the same milepoint I did for my debugging landmark. So I feel there dealing with either a simplified or entirely different geometry, which is fun.
Like I said I could be entirely wrong but I tried every distance formula (geodesic, viscinity?, normal, euclidean) I know at the OG projection, and raw 4326 nothing worked. I hope that this will be of use to someone else and if I can I'll hunt up the ESRI documentation talking about abstract LRS systems from a centerline, and post a picture with the earth versus webmap output.