latlon2cartesian-new.py

latlon2cartesian-new.py

EARTH_RADIUS = 6371
def haversine_distance(point1, point2, radius=EARTH_RADIUS):
    assert point1.shape[1] == 2
    assert point2.shape[1] == 2
    
    latitude1, longitude1 = np.radians(point1).T
    latitude2, longitude2 = np.radians(point2).T
    
    dlongitude = longitude2 - longitude1 
    dlatitude = latitude2 - latitude1 
    a = np.sin(dlatitude/2)**2 + np.cos(latitude1) * np.cos(latitude2) * np.sin(dlongitude/2)**2
    c = 2 * np.arcsin(np.sqrt(a)) 
    km = radius * c
    return km

def cartesian2latlon(xyz):
    assert xyz.shape[1] == 3
    
    xyz /= np.sqrt((xyz ** 2).sum(axis=-1))
    x, y, z = xyz[:,0], xyz[:,1], xyz[:,2]
    latlon = np.column_stack([
        np.pi / 2 - np.arccos(z), # lat
        - np.arctan2(y, x),       # lon
    ])
    return np.degrees(latlon)


def latlon2cartesian(latlon, d=1):
    assert latlon.shape[1] == 2
    
    latlon = np.radians(latlon)
    return np.column_stack([
        d * np.cos(latlon[:,0]) * np.cos(-latlon[:,1]), # x
        d * np.cos(latlon[:,0]) * np.sin(-latlon[:,1]), # y
        d * np.sin(latlon[:,0]),                        # z
    ])