Import pyproj, numPy, math module

In [1]:
import pyproj
import numpy as np
import math

import pyproj as ppj

print("# numpy version = ",np.__version__)
print("# pyproj version = ",ppj.__version__)

# numpy version =  1.19.1
# pyproj version =  2.6.1.post1

Coordinate center

In [2]:
# M3+ events
lond = -115.6000000
latd = 33.1750000

# where latd and lond at xy coordinate
x0 = 0
y0 = 0

xy locations

This example assumes: x -> longitude, y-> latitude. if xy coordinate is rotated with some degrees, we need to use "azimuth" or "invAngle" when we have latd2 and lond2.

In [3]:
x1s = [-20, 10] # distance (km) from x0. will change longitude
y1s = [0,-10]   # distance (km) from y0. will change latitude

Convert from x,y coordinates to longitude, latitude

This example uses WGS84 system. latd2 and lond2 are projected values in lat & lon system.

In [4]:
geod = pyproj.Geod(ellps='WGS84')

latd2out = np.zeros(len(x1s))
lond2out = np.zeros(len(x1s))
for i, x1 in enumerate(x1s):
    #print (x1)
    # do get x1 and y1 in meter
    x1 = x1 * 1000
    y1 = y1s[i] * 1000
    #print ("# x1 = ",x1," y1 = ",y1)
    distance = math.sqrt((x1-x0)**2+(y1-y0)**2)
    azimuth = math.degrees(math.atan2((x1-x0),(y1-y0))) # from x0,y0
    lond2,latd2,invAngle = geod.fwd(lond, latd, azimuth, distance)
    #print (x1, y1,latd2, lond2, invAngle)
    print("# latd2 = ", latd2, " y1 = ", y1)
    print("# lond2 = ", lond2, " x1 = ", x1)
    print("# distance = ", distance, " azimuth = ", azimuth)

# latd2 =  33.17481534720177  y1 =  0
# lond2 =  -115.81443485643108  x1 =  -20000
# distance =  20000.0  azimuth =  -90.0
# latd2 =  33.084788099818695  y1 =  -10000
# lond2 =  -115.49289204670026  x1 =  10000
# distance =  14142.13562373095  azimuth =  135.0
In [ ]: