solve two simultaneous equations: one contains a Python function

Question

On the below map, I have two known points (A and B) with their coordinates (longitude, latitude). I need to derive the coordinates of a point C which is on the line, and is 100 kilometres away from A.

First I created a function to calculate the distances between two points in kilometres:

# pip install haversine
from haversine import haversine

def get_distance(lat_from,long_from,lat_to,long_to):
    distance_in_km = haversine((lat_from,long_from),
                               (lat_to, long_to),
                                unit='km')

    return distance_in_km

Then using the slope and the distance, the coordinates of point C should be the solution to the below equations:

# line segment AB and AC share the same slope, so
# (15.6-27.3)/(41.6-34.7) = (y-27.3)/(x-34.7)

# the distance between A and C is 100 km, so
# get_distance(y,x,27.3,34.7) = 100

Then I try to solve these two equations in Python:

from sympy import symbols, Eq, solve

slope = (15.6-27.3)/(41.6-34.7)
x, y = symbols('x y')
eq1 = Eq(y-slope*(x-34.7)-27.3)
eq2 = Eq(get_distance(y,x,34.7,27.3)-100)
solve((eq1,eq2), (x, y))

The error is TypeError: can't convert expression to float. I may understand the error, because the get_distance function is expecting inputs as floats, while my x and y in eq2 are sympy.core.symbol.Symbol.

I tried to add np.float(x), but the same error remains.

Is there a way to solve equations like these? Or do you have better ways to achieve what is needed?

Many thanks!

# there is a simple example of solving equations:
from sympy import symbols, Eq, solve
x, y = symbols('x y')
eq1 = Eq(2*x-y)
eq2 = Eq(x+2-y)
solve((eq1,eq2), (x, y))

# output: {x: 2, y: 4}

I'm not surprised. I don't know anything about `havesine`, but assuming it uses some sort trig calculations, it either uses `math.sin` or `np.sin`. `math.sin(x)` produces this error. It would be interesting to see the full traceback, but I don't think it will help. `scipy` has some numerical solvers that might work better. — hpaulj, Dec 14 '21 at 00:24
does the data type `sympy.core.symbol.Symbol` contain some sort of field/attribute from which you can fetch float values? — Ryan Millares, Dec 14 '21 at 00:42
If you want a symbolic solution, you need to rewrite your function to only use authentic sympy functions. See e.g. [this post](https://stackoverflow.com/questions/4913349/haversine-formula-in-python-bearing-and-distance-between-two-gps-points) for a formula, where you should change all `sin`, `cos`, `asin`, `sqrt` and conversion to radians to their sympy equivallent. — JohanC, Dec 14 '21 at 00:43
This can be calculated and doesn't need to be solved with solver. That is, if there is a solution. If both points are less than 100KM away from each other, there is no solution — Willem Hendriks, Dec 14 '21 at 08:09
https://stackoverflow.com/questions/38767074/intermediate-points-between-2-geographic-coordinates?rq=1 is related. Let me know if you want a python example — Willem Hendriks, Dec 14 '21 at 08:11

score 3 · Accepted Answer · answered Dec 14 '21 at 09:53

You can directly calculate that point. We can implement a python version of the intermediate calculation for lat long.

Be aware this calculations assume the earth is a sphere, and takes the curve into account, i.e. this is not a Euclidean approximation like your original post.

Say we have two (lat,long) points A and B;

import numpy as np

A = (52.234869, 4.961132)
B = (46.491267, 26.994655)

EARTH_RADIUS = 6371.009

We can than calculate the intermediate point fraction f by taking 100/distance-between-a-b-in-km

from sklearn.neighbors import DistanceMetric
dist = DistanceMetric.get_metric('haversine')

point_1 = np.array([A])
point_2 = np.array([B])

delta = dist.pairwise(np.radians(point_1), np.radians(point_2) )[0][0] 

f = 100 / (delta * EARTH_RADIUS)

phi_1, lambda_1 = np.radians(point_1)[0]
phi_2, lambda_2 = np.radians(point_2)[0]

a = np.sin((1-f) * delta) / np.sin(delta)
b = np.sin( f * delta) / np.sin(delta)

x = a * np.cos(phi_1) * np.cos(lambda_1) + b * np.cos(phi_2) * np.cos(lambda_2)
y = a * np.cos(phi_1) * np.sin(lambda_1) + b * np.cos(phi_2) * np.sin(lambda_2)

z = a * np.sin(phi_1) + b * np.sin(phi_2)
phi_n = np.arctan2(z, np.sqrt(x**2 + y**2) )
lambda_n = np.arctan2(y,x)

The point C, going from A to B, with 100 km distance from A, is than

C = np.degrees( phi_n ), np.degrees(lambda_n)

In this case

(52.02172458025681, 6.384361456573444)

Thanks for your answer, I am still thinking about your answer, and trying your codes with my real research data, I will approve it as soon as I can. — Jeremy, Dec 14 '21 at 13:32
Initially I had concerns with precisions due to assumptions made in this solution or any other ones, like the fixed value of earth radius. But then I think these are inevitable. And there are already large uncertainties in my research data (i.e. the satellite retrievals). So for now, I decide to carry on. Thanks for your solution! — Jeremy, Dec 16 '21 at 17:54
thx and happy it (somewhat) works. There are indeed precision errors, good to read https://www.movable-type.co.uk/scripts/latlong.html#ellipsoid on this. This also mentioned the Vincenty alternative. It is more involved but more precise. I am not aware of a (elegant) form to get intermediate point for vincenty — Willem Hendriks, Dec 16 '21 at 21:49

solve two simultaneous equations: one contains a Python function

1 Answers1