0

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

3
  • Not criticising, just wondering why you want to roll your own reprojection code instead of using something like pyproj?
    – user2856
    Mar 1 at 2:11
  • I would be happy to use pyproj or any other library that can do this out of the box. I will research pyproj but am currently unaware of how that would be done. I would certainly see any library-based implementation as a valid solution to this question. Mar 1 at 19:28
  • I realize that my question is missing some key information and that I have mixed up some terminology. My points are Euclidean, not Cartesian. That is to say, I don't have a projection system for the points, only their relative location to the anker point. I will edit the question to reflect this. Mar 1 at 19:35

1 Answer 1

0

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])

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.