import math

Distance Between Points

When you’re dealing with long distances—like planning infrastructure between cities or analyzing migration patterns—great-circle distance is the gold standard. It accounts for Earth’s curvature, giving you the shortest path over the surface between two lat/lon points.

def calculate_distance(lat1, lon1, lat2, lon2):
    R = 6371
    lat1_rad = math.radians(lat1)
    lon1_rad = math.radians(lon1)
    lat2_rad = math.radians(lat2)
    lon2_rad = math.radians(lon2)
    
    a = math.sin((lat2_rad - lat1_rad) / 2) * math.sin((lat2_rad - lat1_rad) / 2) + math.cos(lat1_rad) * math.cos(lat2_rad) * math.sin((lon2_rad - lon1_rad) / 2) * math.sin((lat2_rad - lat1_rad) / 2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
    
    return R * c

Earths radius can be in kilometers (6371), statute miles (3959), or nautical miles (3444) to give a distance in that unit. As a best practice, use the same unit of measurement in all of your formulas, regardless of which unit you choose.

Result

Intermediate Points

Sometimes, you don’t just care about where something starts and ends—you need to know what’s in between. Being able to calculate intermediate points it’s a practical skill when needing to interpolate travel paths, generating regular sampling points along a flight path, or even animating movement across a globe. GIS users can use this to create evenly spaced waypoints or to understand where an object would be at a specific distance along a path.

def calc_midpoint(lat1, lon1, lat2, lon2, interval):
    R = 6371
    d = calculate_distance(lat1, lon1, lat2, lon2)
    lat1_rad = math.radians(lat1)
    lon1_rad = math.radians(lon1)
    lat2_rad = math.radians(lat2)
    lon2_rad = math.radians(lon2)
    
    f = interval / d
    a = math.sin((1 - f) * (d / R)) / math.sin(d / R)
    b = math.sin(f * (d / R)) / math.sin(d / R)
    x = a * math.cos(lat1_rad) * math.cos(lon1_rad) + b * math.cos(lat2_rad) * math.cos(lon2_rad)
    y = a * math.cos(lat1_rad) * math.sin(lon1_rad) + b * math.cos(lat2_rad) * math.sin(lon2_rad)
    z = a * math.sin(lat1_rad) + b * math.sin(lat2_rad)
    
    return to_latlon(x, y, z)

While this calculation produces a point in three-dimensional Cartesian space (X, Y, Z), the return value is that of a standard latitude and longitude. This conversion is done via a helper function that can be reused as needed.

def to_latlon(x, y, z):
    new_lat = math.degrees(math.atan2(z, math.sqrt(x * x + y * y)))
    new_lon = math.degrees(math.atan2(y, x))
    
    return {"lat": new_lat, "lon": new_lon}

Result

Great Circle Intersections

For global-scale spatial modeling, finding the intersection of two great circle arcs is a powerful technique. It’s particularly relevant in airspace deconfliction, transcontinental logistics, or anywhere else you need to identify the precise point where two curved paths cross.

We'll start by making another helper function, this time to convert from latitude and longitude to cartesian coordinates.

def to_cartesian(lat, lon):
    lat_rad = math.radians(lat)
    lon_rad = math.radians(lon)
    
    x = math.cos(lat_rad) * math.cos(lon_rad)
    y = math.cos(lat_rad) * math.sin(lon_rad)
    z = math.sin(lat_rad)

    return {"x": x, "y": y, "z": z}

Now each point can easily be converted into Cartesian coordinates. From there, we can calculate the normal vectors of the planes defined by each arc using cross products. These vectors represent the orientation of each great circle. The function then finds the line where these two planes intersect by computing the cross product of the two normals. That gives us a 3D point, which is then converted back into latitude and longitude coordinates to find the actual location of the intersection on the globe.

def calculate_intersection(lat1, lon1, lat2, lon2, lat3, lon3, lat4, lon4):
    point1 = to_cartesian(lat1, lon1)
    point2 = to_cartesian(lat2, lon2)
    point3 = to_cartesian(lat3, lon3)
    point4 = to_cartesian(lat4, lon4)

    N1x = point1["y"] * point2["z"] - point1["z"] * point2["y"]
    N1y = point1["z"] * point2["x"] - point1["x"] * point2["z"]
    N1z = point1["x"] * point2["y"] - point1["y"] * point2["x"]
    N2x = point3["y"] * point4["z"] - point3["z"] * point4["y"]
    N2y = point3["z"] * point4["x"] - point3["x"] * point4["z"]
    N2z = point3["x"] * point4["y"] - point3["y"] * point4["x"]

    N1 = math.sqrt(N1x ** 2 + N1y ** 2 + N1z ** 2)
    N2 = math.sqrt(N2x ** 2 + N2y ** 2 + N2z ** 2)

    S1x, S1y, S1z = N1x / N1, N1y / N1, N1z / N1
    S2x, S2y, S2z = N2x / N2, N2y / N2, N2z / N2

    P1x = S1y * S2z - S1z * S2y
    P1y = S1z * S2x - S1x * S2z
    P1z = S1x * S2y - S1y * S2x

    return to_latlon(P1x, P1y, P1z)

However, this point is one of two possible points. It might be on the great circle but beyond the arc segment. The only way to verify is to test the point. Replace the return line with the code below.

    possible = to_latlon(P1x, P1y, P1z)

    arc1 = on_arc(point1, point2, possible["lat"], possible["lon"])
    arc2 = on_arc(point3, point4, possible["lat"], possible["lon"])

    if arc1 and arc2:
        return {"lat": round(possible["lat"], 4), "lon": round(possible["lon"], 4)}

But what if the point isn't on those arcs? In that case, we'll need to test the other point, the antipode (the directly opposite point on the globe.) Add the following code to create and test this other point.

    anti_lat = -possible["lat"]
    anti_lon = possible["lon"] + 180 if possible["lon"] < 0 else possible["lon"] - 180
    anti_point = {"lat": anti_lat, "lon": anti_lon}

    arc1 = on_arc(point1, point2, anti_point["lat"], anti_point["lon"])
    arc2 = on_arc(point3, point4, anti_point["lat"], anti_point["lon"])

    if arc1 and arc2:
        return {"lat": round(anti_point["lat"], 4), "lon": round(anti_point["lon"], 4)}

These if-statements are dependent on new function, on_arc, which tests whether the point is on the arc. In this function, a series of dot products are computed then turned into angles that let compare the angles between the test point, and the two arc points. If the angles match (within a tiny tolerance for floating point rounding errors), then the point is on the arc and the function returns a True value. If they don’t match, then the point is somewhere else and a False value is returned.

def on_arc(arc_start, arc_end, test_lat, test_lon):
    test_point = to_cartesian(test_lat, test_lon)
    
    ab = arc_start['x'] * arc_end['x'] + arc_start['y'] * arc_end['y'] + arc_start['z'] * arc_end['z']
    at = arc_start['x'] * test_point['x'] + arc_start['y'] * test_point['y'] + arc_start['z'] * test_point['z']
    tb = test_point['x'] * arc_end['x'] + test_point['y'] * arc_end['y'] + test_point['z'] * arc_end['z']
    
    ang_ab = math.acos(min(1, max(-1, ab)))
    ang_at = math.acos(min(1, max(-1, at)))
    ang_tb = math.acos(min(1, max(-1, tb)))
    
    if abs((ang_at + ang_tb) - ang_ab) < 0.000001:
        return True
    else:
        return False

Now we can determine where the intersection of these two arcs are. However, this can still fail if the two arcs don't intersect, if they're the same arc, or in some other abstract edge case, and there is no return case in those scenarios. To be truly safe, you'll need to come up with some sort of way to handle those errors, the specifics dependent on your specific application.

Result

Azimuth and Elevation

Whether you’re setting up satellite TV, configuring a Broadband Global Area Network terminal, or determining the availability of satellite communications for aircraft while in flight, azimuth and elevation tell you exactly how to aim your antennas and ensure your equipment is aligned correctly with geostationary satellites for the strongest possible signal. In GIS, understanding these angles are essential to modeling satellite visibility.

Websites such as https://www.n2yo.com/ can help identify a satellite's longitude

Calculating Elevation

Elevation tells you how high above the horizon a satellite appears from a specific location on Earth. A higher elevation angle typically means better visibility and less atmospheric interference. If the elevation is too low—or negative—it’s a quick way to determine that the satellite is below the horizon and out of view

def calculate_elevation(lat, lon, sat_lon):
    my_lat_rad = math.radians(lat)
    my_lon_rad = math.radians(lon)
    sat_lon_rad = math.radians(sat_lon)
    
    g = sat_lon_rad - my_lon_rad
    
    return to_degrees(math.atan((math.cos(g) * math.cos(my_lat_rad) - 0.1512) / math.sqrt(1 - (math.cos(g) * math.cos(g)) * (math.cos(my_lat_rad) * math.cos(my_lat_rad)))))

At first glance, it might seem that there's nothing in this formula that corresponds to the satellite's altitude. However, the 0.1512 is a constant representing the radius of the Earth (6371 km), divided by the radius of the satellite's orbit (42,164 km from the Earth's center.) If the satellite was at a different altitude, you would calculate a new value using 6371 / <Orbit radius>

Calculating Azimuth

While elevation tells you how high, azimuth tells you which direction. It’s the compass bearing—measured in degrees clockwise from true north—that you would point an antenna to reach the satellite.

def calculate_azimuth(lat, lon, sat_lon):
    my_lat_rad = math.radians(lat)
    my_lon_rad = math.radians(lon)
    sat_lon_rad = math.radians(sat_lon)
    
    return math.degrees(math.pi - math.atan(math.tan(sat_lon_rad - my_lon_rad) / math.sin(my_lat_rad)))

Result

NOTE: Websites such as https://www.n2yo.com/ can help identify a satellite's longitude