I do not believe my math is incorrect. This is seems like it could be a projection issue.
I start with the following two coordinates ( lon, lat in degrees )
a: [37.821055535000085, -46.84417083099993]
b: [37.89112389400009, -46.64967213299991]
When visualized with geopandas, the lineSegment
defined by the points looks like:
path = 'data.geo.json'
data = json.load( open( path, "r" ) )
gdf = geopandas.GeoDataFrame.from_features( data, crs = 4326 )
What is drawn here is correct.
I want to extend this line segment in both directions and draw the perpendicular line segment.
To do this, I am using turf and some math:
const direction = [b[0] - a[0], b[1] - a[1]]
const length = Math.sqrt(Math.pow(direction[0], 2) + Math.pow(direction[1], 2))
direction[0] = direction[0] / length
direction[1] = direction[1] / length
const b2 = [a[0] + direction[0] * length, a[1] + direction[1] * length]
const lineSegment = turf.lineString([a, b2])
const perpendicularDirection = [-direction[1], direction[0]]
const boxDiagonalLength = length * 2
let boxPoint = a
const iFactor = direction[0] * boxDiagonalLength
const jFactor = direction[1] * boxDiagonalLength
const normalLineEndA = [boxPoint[0] + iFactor, boxPoint[1] + jFactor]
const normalLineEndB = [boxPoint[0] - iFactor, boxPoint[1] - jFactor]
const normalLine = turf.lineString([normalLineEndA, normalLineEndB])
const piFactor = perpendicularDirection[0] * boxDiagonalLength
const pjFactor = perpendicularDirection[1] * boxDiagonalLength
const pNormalLineEndA = [boxPoint[0] + piFactor, boxPoint[1] + pjFactor]
const pNormalLineEndB = [boxPoint[0] - piFactor, boxPoint[1] - pjFactor]
const pNormalLine = turf.lineString([pNormalLineEndA, pNormalLineEndB])
When I visualize the extended line ( normalLine
) and the perpendicular line to the extended line ( pNormalLine
) with geopandas, it produces:
The extended line is shifted to the right. The line that should be perpendicular isn't.
I expected the extended line to overlap the line segment and for the perpendicular line to be perpendicular. Additionally, I would expect the point of intersection to be at a
.
I am not sure what is going wrong. I am guessing there is a fundamental concept I am missing.
What am I missing?
Best Answer
Proposal from @Jan Šimbera to convert coordinates to projected CRS EPSG:32737 (UTM zone 37S) and then use cartesian math would in this case work quite OK, since line length is only about 20km and central meridian for EPSG:32737 is 39, very close to the line coordinates.
But much simpler and generally valid solution would be to use
turf.bearing
method to get line bearing and thenturf.destination
method to get coordinates of desired points.Code could then look something like this (point offsets are one quarter of line length, see image below for meaning of points; code is JS, map by Leaflet):
That's how this would look like on the map: