Euclidean Coordinates to Lat Long Based on Reference Location in Python

Question

Goal

Convert Euclidean Coordinates to Latitude and Longitude given the Latitude and Longitude of euclidean origin using Python. To give a sense of scale and accuracy for my application, the coordinates will always be within 2,000 meters of the origin.

Data on hand.

#A 'Anker Point' or the earth location for the cartesian origin (0,0)
euclidean_anker_point = [37.326695, -121.902583]

#The list of points to convert to lat lon
euclidean_points = [[1450.3, 400.8],[1460.2, 460.6],[1470.2, 470.2]]

The euclidean points currently have no projection system, have the units of Meters, and are only understood in euclidian space relative to the Anker point. I have no projection system for these coordinates.

First Pass

The below code does not provide accurate results but does seem to be in the ballpark. I believe how the bearing is calculated is not correct given my context.

Does anyone have a suggestion on how to improve this approach?

import numpy
import math
from pygeodesy.ellipsoidalVincenty import LatLon 

#create pygeodesy anker point
origin_point = LatLon(euclidean_anker_point[0], euclidean_anker_point[1])

#create ouput array
geographic_points = []

#for each point to convert to lat lon
for index, point in enumerate(euclidean_points):

    #get distance from anker point using Pythagorean theorem
    dist = math.sqrt((point[0]*point[0])+(point[1]*point[1])) 

    #get bearing using function (see below) 
    bearing = get_bearing(0,0,point[0],point[1])

    #use pygeodesy to calculate desination lat lon
    dest = origin_point.destination(dist, bearing)

    #add lat lon version of point to output array
    geographic_points.append([dest.lon,dest.lat])


def get_bearing(lat1, long1, lat2, long2):
    dLon = (long2 - long1)
    x = math.cos(math.radians(lat2)) * math.sin(math.radians(dLon))
    y = math.cos(math.radians(lat1)) * math.sin(math.radians(lat2)) - 
    math.sin(math.radians(lat1)) * math.cos(math.radians(lat2)) * 
    math.cos(math.radians(dLon))
    brng = numpy.arctan2(x,y)

    return brng

Credits

I am trying the pygeodesy approach listed by Antonia Falciano and using Get_Bearing from Aliff Daniel

Answer 1

Here is the solution I have gone with.

Brief

The fastest way to deal with projections is using one of the many robust projection libraries out there. Most of these expect that your original coordinates have a projection system, so the best thing to do is to focus on converting the Euclidean points to Cartesian points, then re-project those into Latitude and Longitude.

Step 1

Get the anker_point into a projection system with units that match your coordiante units.

#A 'Anker Point' or the earth location for the cartesian origin (0,0)
anker_point = [37.326695, -121.902583]

#Cartesian Projection System
inProj = Proj(init='epsg:2227')

#Anker Projection System
outProj = Proj(init='epsg:4326')

#using the pyproj library (thanks user2856!)
from pyproj import Proj, transform

#produce a anker point with same units as euclidean points
cartesian_anker_point= transform(outProj,inProj,anker_point[1],anker_point [0])

Step 2

Now we can use the cartesian anker point to locate the Euclidean coordinates in our cartesian space.

#euclidian point to locate
euclidan_point = [1450.3, 400.8]

#located point in  inProj system
cartisan_point = [euclidan_point[0]+cartesian_anker_point[0],euclidan_point[1]+cartesian_anker_point[1]]

Step 3

With a located Cartesian point we can now project to lat lon

#located point in  outProj system
latLon_Point = transform(inProj,outProj,cartisan_point[0],cartisan_point [1])