solve two simultaneous equations: one contains a Python function

user_9851307 0
Viewed 308

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. enter image description here

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}
sympy python
1 Answers

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)
user_31849500
Related