Create a longer line segment from a shorter one?

Question

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?

Answer 1

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 then turf.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):

var a = turf.point([-46.84417083099993, 37.821055535000085]);
var b = turf.point([-46.64967213299991, 37.89112389400009]);

var distance = turf.distance(a, b);
var bearing = turf.bearing(a, b);

var a2 = turf.destination(a, (distance / 4), (bearing + 180));
var b2 = turf.destination(b, (distance / 4), bearing);

var ap1 = turf.destination(a, (distance / 4), (bearing + 90));
var ap2 = turf.destination(a, (distance / 4), (bearing - 90));

var bp1 = turf.destination(b, (distance / 4), (bearing + 90));
var bp2 = turf.destination(b, (distance / 4), (bearing - 90));

That's how this would look like on the map:

Answer 2

You are using planar geometry equations (Euclidean distance measurement, etc.) on spherical coordinates, which will introduce errors.

If your data only cover a limited geographical extent, try converting the coordinates to a planar coordinate system (such as WGS84/UTM37S - check the linked page if the extent is sufficient), doing the calculations in those coordinates and then converting back. Geopandas make the conversions easy to do in a single function call.