Maths in code: distances on the Earth and rotations in 3D with Python

Reading about vectors, datums and precision is useful. Running a few lines of code that show them in action is better. In this tutorial you will use Python and numpy to measure distances on a curved Earth, rotate points in 3D, convert pixels to map coordinates and watch floating-point numbers misbehave. Each exercise is small, and each ends with a check that you can compare with the numbers in this article.

What you need

• Python 3.9 or newer.
• The numpy package.
• A text editor or notebook, and about thirty minutes.

pip install numpy

Exercise 1: distance between two cities

We will measure the distance between Kathmandu and Pokhara in two ways: the wrong way and the right way.

import math

EARTH_RADIUS_KM = 6371.0088

def haversine_km(lat1, lon1, lat2, lon2):
    phi1, phi2 = math.radians(lat1), math.radians(lat2)
    dphi = phi2 - phi1
    dlambda = math.radians(lon2 - lon1)
    a = math.sin(dphi / 2) ** 2 + math.cos(phi1) * math.cos(phi2) * math.sin(dlambda / 2) ** 2
    return 2 * EARTH_RADIUS_KM * math.asin(math.sqrt(a))

ktm = (27.7172, 85.3240)   # latitude, longitude
pkr = (28.2096, 83.9856)

good = haversine_km(*ktm, *pkr)
naive = 111.0 * math.hypot(pkr[0] - ktm[0], pkr[1] - ktm[1])   # treats degrees as flat

print(f"haversine: {good:.1f} km")
print(f"naive:     {naive:.1f} km  ({(naive - good) / good:+.1%})")

You should see about 142.4 kilometres for the haversine, and about 158.3 for the naive version, an error of roughly 11 percent. The naive version overestimates because a degree of longitude near 28 degrees north is about 98 kilometres, not 111.

Exercise 2: rotate a point and check its length

A rotation matrix turns a point around an axis without changing its distance from the origin. We build one for each axis and test it.

import numpy as np

def rot_z(deg):
    t = np.radians(deg)
    return np.array([[np.cos(t), -np.sin(t), 0],
                     [np.sin(t),  np.cos(t), 0],
                     [0,          0,         1]])

def rot_x(deg):
    t = np.radians(deg)
    return np.array([[1, 0,          0],
                     [0, np.cos(t), -np.sin(t)],
                     [0, np.sin(t),  np.cos(t)]])

point = np.array([1.0, 0.0, 0.0])
turned = rot_z(90) @ point

print(np.round(turned, 6))                                  # expect [0, 1, 0]
print(np.isclose(np.linalg.norm(turned), np.linalg.norm(point)))   # length preserved: True

A 90 degree turn about the vertical axis swings the point from the x axis to the y axis, and its length stays 1. Whenever you build a rotation, check both.

Exercise 3: the order of rotations matters

p = np.array([0.0, 1.0, 0.0])

a = rot_z(90) @ rot_x(90) @ p     # rotate about x first, then z
b = rot_x(90) @ rot_z(90) @ p     # rotate about z first, then x

print(np.round(a, 6))   # [0, 0, 1]
print(np.round(b, 6))   # [-1, 0, 0]

Same two rotations, different results. When you combine transformations in a game engine, a headset’s tracking code or a GIS, you are choosing an order, so choose deliberately. Matrix products apply from right to left: the one nearest the point acts first.

Exercise 4: from a pixel to a map coordinate

A georeferenced raster stores a few numbers that convert a pixel’s column and row to a position on the map. For an image with 10 metre pixels, whose top-left corner is at x 500,000 and y 3,100,000 metres, the conversion is as follows.

x0, y0, pixel = 500_000.0, 3_100_000.0, 10.0

def pixel_to_map(col, row):
    x = x0 + col * pixel + pixel / 2       # +pixel/2 gives the pixel centre
    y = y0 - row * pixel - pixel / 2       # minus: rows go DOWN, northings go UP
    return x, y

def map_to_pixel(x, y):
    col = int((x - x0) // pixel)
    row = int((y0 - y) // pixel)
    return col, row

x, y = pixel_to_map(250, 100)
print(x, y)                      # 502505.0 3098995.0
print(map_to_pixel(x, y))        # (250, 100): round trip works

The minus sign on y is the single most common source of upside-down maps. Always test a conversion by going there and back.

Exercise 5: when numbers run out of digits

Computers store real numbers with limited precision. Most graphics code uses 32-bit floats, which hold about seven significant digits. Watch what happens with large map coordinates.

big = np.float32(5_000_000.0)      # a UTM-style northing in metres
print(np.spacing(big))             # smallest step this number can take: 0.5
print(big + np.float32(0.01) == big)   # True: adding a centimetre does nothing

print(np.spacing(5_000_000.0))     # 64-bit float: about 9e-10

A 32-bit number at five million metres cannot represent anything finer than half a metre. In a 3D scene far from the origin, vertices jitter, and in VR that shaking is very noticeable. The cure is to subtract a local origin from your coordinates, so that numbers stay small, and keep the full coordinates in 64-bit precision.

Turn your checks into tests

These exercises are also a ready-made test suite. Collect them as automatic checks that you can run any time you change code.

def test_basics():
    assert abs(haversine_km(*ktm, *pkr) - 142.4) < 0.2
    assert np.allclose(rot_z(90) @ np.array([1.0, 0, 0]), [0, 1, 0], atol=1e-9)
    assert map_to_pixel(*pixel_to_map(250, 100)) == (250, 100)

test_basics()
print("All checks passed")

What to try next

• Compute the distance with a geodesic library, and compare it with the haversine. The difference is usually a few tenths of a percent.
• Build a function that converts latitude and longitude to a local east-north-up frame around a chosen point.
• Use degrees where you should have used radians, and see how wrong the output is. Knowing what the error looks like is half the cure.
• Ask an AI assistant to write the same functions, and run these tests against its code.

Frequently asked questions

Why use radians instead of degrees?

Because the mathematics of sine, cosine and their derivatives is simplest in radians, so most programming languages expect them. Convert at the edge of your code and stay consistent inside it.

Is the haversine accurate enough?

For most analysis, yes: it is usually within about half a percent. For surveying, use an ellipsoidal geodesic.

What is the difference between a rotation matrix and a quaternion?

Both represent orientation. Matrices are easy to understand and combine, while quaternions avoid gimbal lock and interpolate smoothly, so tracking systems prefer them.

Why do my 3D objects shake far from the origin?

Because single-precision numbers lose resolution at large values. Re-centre the scene near the area of interest.

None of these exercises needs more than a few lines. That is the point. The fundamentals are small, and once you can test them quickly, you can verify the larger, more complex code that AI tools generate for you.