For example, I have a list of coordinates (lat, lon):
Coords_of_points = [(57.769999999999996, 108.06),
(60.271336095306005, 119.70058684113518),
(55.44781555321857, 121.56598882454766),
(42.39597453807919, 141.90642068006264),
(47.416832694979675, 126.65598006176236)]
Center_of_coordinates = (48.528055555555554, 135.18805555555556)
For calculating X and Y of points from Center_of_coordinates
, I use this function:
def transform(Coords_of_points, Center_of_coordinates):
for i in range(len(Coords_of_points)):
line = Coords_of_points[i]
rad = 6372795
lat_point = math.radians(line[0])
lon_point = math.radians(line[1])
lat_0 = math.radians(Center_of_coordinates[0])
lon_0 = math.radians(Center_of_coordinates[1])
cl1 = math.cos(lat_0)
cl2 = math.cos(lat_point)
sl1 = math.sin(lat_0)
sl2 = math.sin(lat_point)
delta = lon_point - lon_0
cdelta = math.cos(delta)
sdelta = math.sin(delta)
y = math.sqrt(math.pow(cl2 * sdelta, 2) + math.pow(cl1 * sl2 - sl1 * cl2 * cdelta, 2))
x = sl1 * sl2 + cl1 * cl2 * cdelta
ad = math.atan2(y, x)
dist = ad * rad
x = (cl1 * sl2) - (sl1 * cl2 * cdelta)
y = sdelta * cl2
z = math.degrees(math.atan(-y / x))
if (x < 0):
z = z + 180
z2 = (z + 180.) % 360. - 180
z2 = - math.radians(z2)
anglerad2 = z2 - ((2 * math.pi) * math.floor((z2 / (2 * math.pi))))
x = round(((math.sin(anglerad2) * dist) / 1000), 3)
y = round(((math.cos(anglerad2) * dist) / 1000), 3)
yield x,y
This function works good when points located near with Center_of_coordinates
, but if they located to far from Center_of_coordinates
, X and Y has been calculating with mistake, I guess this is happening due to Earth Curvature.
So, if i have Point (57.769999999999996, 108.06)
, its X and Y from Center_of_coordinates
will be calculated as (-1577, 1326.651)
. The coorrect answer should be: (-1590, 1338.34)
.
Anybody know how to modify or maybe change the whole algoritm to get correct X and Y of points using WGS84?