GIS & Geospatial Data (QGIS, GeoPandas, PostGIS) Module 2 – Coordinate reference systems

Reprojecting data without silent errors

Transform layers between CRSs and catch the silent errors, from set_crs against to_crs and swapped axes to a missing operation and a doubled step.

Intermediate
20 minutes
Examples run with Python 3.14.8 and Pyodide 314.0.7
By MySmartCoPilot

Free: the whole lesson. With a pass: the model answers to its 2 interview questions.

What you will learn

Distinguish labelling a layer's CRS from transforming it, and repair a layer labelled with the wrong CRS
Pass coordinates in the axis order that pyproj, DuckDB, GeoPandas and PostGIS each expect
Choose among the operations between two CRSs by area of use and accuracy, and refuse when none applies
Explain on-the-fly reprojection in QGIS and why analysis works on data saved in one projected CRS
Show that a round trip is safe and that a transformation applied twice is not

Before you start

EPSG codes, WKT2 and PROJJSON in practice

What you’ll learn · Before you start

What you will learn

Distinguish labelling a layer's CRS from transforming it, and repair a layer labelled with the wrong CRS
Pass coordinates in the axis order that pyproj, DuckDB, GeoPandas and PostGIS each expect
Choose among the operations between two CRSs by area of use and accuracy, and refuse when none applies
Explain on-the-fly reprojection in QGIS and why analysis works on data saved in one projected CRS
Show that a round trip is safe and that a transformation applied twice is not

Before you start

EPSG codes, WKT2 and PROJJSON in practice

On this page

Label a layer of Pune’s clinics, stored in degrees, as UTM zone 43N instead of transforming it, and every clinic moves to the Indian Ocean near the equator, 2,081 km away. Pass the same coordinates longitude first to a transformation that expects latitude first, and they land 7,000 km away, off Norway. Neither mistake raises an error.

This lesson separates the steps where such errors hide: labelling a CRS and transforming to one, the axis order each engine expects, the choice among several operations or none, on-the-fly reprojection in QGIS, and a transformation applied twice. Each is reproduced by a program or stepped through in a figure, with a check that catches it.

Labelling a CRS is not transforming it

Every library has two operations that sound alike. GeoPandas puts it plainly: “Setting a projection is how one tells GeoPandas how to interpret coordinates”, while “re-projecting is the process of changing the representation of locations from one coordinate system to another” (GeoPandas). set_crs labels; to_crs transforms. PostGIS warns that ST_Transform “is often confused with ST_SetSRID”: the first changes the coordinates, the second “simply changes the SRID identifier” (PostGIS). The program does both to Pune:

Pune labelled as UTM, and Pune transformed to UTM Python · reproject/label_or_transform.py

"""Labelling coordinates with a CRS (set_crs) against transforming them (to_crs), for a point in Pune."""
from crs import distance, from_utm, to_utm

lon, lat = 73.8567, 18.5204  # Pune, in degrees on WGS 84 (EPSG:4326)
x, y = to_utm(lat, lon, 43)
print(f"Transformed to UTM 43N: easting {x:,.2f} m, northing {y:,.2f} m")

# Labelled UTM 43N without transforming: the degree values are now read as metres.
wrong_lat, wrong_lon = from_utm(lon, lat, 43)
km = distance(lat, lon, wrong_lat, wrong_lon) / 1000
print(f"Labelled UTM 43N as they are: easting {lon} m, northing {lat} m")
print(f"  which is the point {wrong_lat:.6f}° N, {wrong_lon:.4f}° E, {km:,.0f} km from Pune")


def label_problem(e, n):
    """Why an easting and northing cannot be UTM north metres, or None if they could be."""
    if not 100_000 <= e <= 900_000:
        return f"an easting of {e:,} m is outside 100,000 to 900,000 m"
    if not 0 <= n <= 9_400_000:
        return f"a northing of {n:,} m is outside 0 to 9,400,000 m"
    return None


print("Check of the UTM label:", label_problem(lon, lat) or "plausible")
print("Check of the transformed point:", label_problem(round(x, 2), round(y, 2)) or "plausible")

Output

Transformed to UTM 43N: easting 379,320.82 m, northing 2,048,144.62 m
Labelled UTM 43N as they are: easting 73.8567 m, northing 18.5204 m
  which is the point 0.000167° N, 70.5119° E, 2,081 km from Pune
Check of the UTM label: an easting of 73.8567 m is outside 100,000 to 900,000 m
Check of the transformed point: plausible

Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 label_or_transform.py

It imports its operations from a module that every program of this lesson shares: UTM zone 43N both ways, with the transverse Mercator series of the previous lessons, a datum shift by geocentric translation, and Vincenty’s distance on WGS 84:

crs.py: the coordinate operations of this lesson's programs Python · reproject/crs.py

"""Coordinate operations with the standard library, for the other programs of this lesson: UTM both ways, datum
shifts by geocentric translation (EPSG method 9603), and distances on WGS 84."""
import math

ELLIPSOIDS = {"WGS 84": (6_378_137.0, 298.257223563), "Everest 1830 (1937 Adjustment)": (6_377_276.345, 300.8017),
              "Everest 1830 (1975 Definition)": (6_377_299.151, 300.8017255)}
A, INV_F = ELLIPSOIDS["WGS 84"]
F = 1 / INV_F
N = F / (2 - F)
E = math.sqrt(F * (2 - F))
B = A / (1 + N) * (1 + N ** 2 / 4 + N ** 4 / 64)
FWD = (N / 2 - 2 * N ** 2 / 3 + 5 * N ** 3 / 16 + 41 * N ** 4 / 180,
       13 * N ** 2 / 48 - 3 * N ** 3 / 5 + 557 * N ** 4 / 1440,
       61 * N ** 3 / 240 - 103 * N ** 4 / 140,
       49561 * N ** 4 / 161280)  # the 4th-order series of IOGP report 373-07-2, forward
INV = (N / 2 - 2 * N ** 2 / 3 + 37 * N ** 3 / 96 - N ** 4 / 360,
       N ** 2 / 48 + N ** 3 / 15 - 437 * N ** 4 / 1440,
       17 * N ** 3 / 480 - 37 * N ** 4 / 840,
       4397 * N ** 4 / 161280)  # and back


def to_utm(lat, lon, zone):
    """Latitude and longitude on WGS 84 to easting and northing in a northern UTM zone (EPSG method 9807)."""
    p = math.radians(lat)
    beta = math.atan(math.sinh(math.asinh(math.tan(p)) - E * math.atanh(E * math.sin(p))))
    eta0 = math.atanh(math.cos(beta) * math.sin(math.radians(lon - (6 * zone - 183))))
    xi0 = math.asin(math.sin(beta) * math.cosh(eta0))
    xi = xi0 + sum(h * math.sin(2 * j * xi0) * math.cosh(2 * j * eta0) for j, h in enumerate(FWD, 1))
    eta = eta0 + sum(h * math.cos(2 * j * xi0) * math.sinh(2 * j * eta0) for j, h in enumerate(FWD, 1))
    return 500_000 + 0.9996 * B * eta, 0.9996 * B * xi


def from_utm(x, y, zone):
    """Easting and northing in a northern UTM zone back to latitude and longitude on WGS 84."""
    eta_, xi_ = (x - 500_000) / (0.9996 * B), y / (0.9996 * B)
    xi0 = xi_ - sum(h * math.sin(2 * j * xi_) * math.cosh(2 * j * eta_) for j, h in enumerate(INV, 1))
    eta0 = eta_ - sum(h * math.cos(2 * j * xi_) * math.sinh(2 * j * eta_) for j, h in enumerate(INV, 1))
    beta = math.asin(math.sin(xi0) / math.cosh(eta0))
    q = q0 = math.asinh(math.tan(beta))
    for _ in range(15):
        q = q0 + E * math.atanh(E * math.tanh(q))
    lon = 6 * zone - 183 + math.degrees(math.asin(math.tanh(eta0) / math.cos(beta)))
    return math.degrees(math.atan(math.sinh(q))), lon


def shift(lat, lon, source, shift_m, target="WGS 84"):
    """A datum shift by geocentric translations: X, Y, Z on the source ellipsoid, plus shift_m, back on the target."""
    a, inv_f = ELLIPSOIDS[source]
    e2 = (2 - 1 / inv_f) / inv_f
    p, l = math.radians(lat), math.radians(lon)
    nu = a / math.sqrt(1 - e2 * math.sin(p) ** 2)
    x = nu * math.cos(p) * math.cos(l) + shift_m[0]
    y = nu * math.cos(p) * math.sin(l) + shift_m[1]
    z = nu * (1 - e2) * math.sin(p) + shift_m[2]
    a, inv_f = ELLIPSOIDS[target]
    e2 = (2 - 1 / inv_f) / inv_f
    r = math.hypot(x, y)
    lat2 = math.atan2(z, r * (1 - e2))
    for _ in range(10):
        lat2 = math.atan2(z + e2 * a / math.sqrt(1 - e2 * math.sin(lat2) ** 2) * math.sin(lat2), r)
    return math.degrees(lat2), math.degrees(math.atan2(y, x))


def distance(lat1, lon1, lat2, lon2):
    """Metres between two points on WGS 84 (Vincenty's inverse formula)."""
    b = A * (1 - F)
    u1, u2 = math.atan((1 - F) * math.tan(math.radians(lat1))), math.atan((1 - F) * math.tan(math.radians(lat2)))
    lam = big_l = math.radians(lon2 - lon1)
    for _ in range(200):
        sin_s = math.hypot(math.cos(u2) * math.sin(lam),
                           math.cos(u1) * math.sin(u2) - math.sin(u1) * math.cos(u2) * math.cos(lam))
        if sin_s == 0:
            return 0.0
        cos_s = math.sin(u1) * math.sin(u2) + math.cos(u1) * math.cos(u2) * math.cos(lam)
        sigma = math.atan2(sin_s, cos_s)
        sin_al = math.cos(u1) * math.cos(u2) * math.sin(lam) / sin_s
        cos2 = 1 - sin_al ** 2
        c2m = cos_s - 2 * math.sin(u1) * math.sin(u2) / cos2 if cos2 else 0.0
        c = F / 16 * cos2 * (4 + F * (4 - 3 * cos2))
        last = lam
        lam = big_l + (1 - c) * F * sin_al * (sigma + c * sin_s * (c2m + c * cos_s * (-1 + 2 * c2m ** 2)))
        if abs(lam - last) < 1e-12:
            break
    u_sq = cos2 * (A * A - b * b) / (b * b)
    big_a = 1 + u_sq / 16384 * (4096 + u_sq * (-768 + u_sq * (320 - 175 * u_sq)))
    big_b = u_sq / 1024 * (256 + u_sq * (-128 + u_sq * (74 - 47 * u_sq)))
    d_s = big_b * sin_s * (c2m + big_b / 4 * (cos_s * (-1 + 2 * c2m ** 2)
                                            - big_b / 6 * c2m * (-3 + 4 * sin_s ** 2) * (-3 + 4 * c2m ** 2)))
    return b * big_a * (sigma - d_s)

Transformed, Pune is at easting 379,320.82 m and northing 2,048,144.62 m, the values GeographicLib’s converter gives (GeoConvert). Labelled as UTM 43N without transforming, the degrees are read as metres, an easting of 73.86 m and a northing of 18.52 m, which is a point 0.000167° north of the equator at 70.51° E, 2,081 km from Pune. A range check exposes the label: UTM eastings in a zone lie between about 100,000 and 900,000 m, and 73.86 m does not. To repair such a layer, put back the label that matches the numbers, here EPSG:4326, and then transform it.

Fixing a projection with set_crs

A layer drawn in the wrong place does not need a new label unless its label is wrong. If the numbers are degrees on WGS 84 and the label says so, set_crs(32643, allow_override=True) only makes it worse; to_crs(32643) moves the numbers into the UTM grid. Change a label with allow_override=True only when the old label was wrong.

Every engine has an axis-order rule

EPSG defines its geographic CRSs with latitude first, while most files and code write longitude first. Each engine resolves this in its own way, and each documents it:

Engine	Order for EPSG:4326 coordinates	Its documentation
pyproj `Transformer`	Latitude first, unless created with `always_xy=True`	“the traditional GIS order” with always_xy (pyproj)
DuckDB `ST_Transform`	The CRS definition’s order, unless `always_xy` is true	an “[easting, northing] coordinate axis order” with always_xy (DuckDB)
GeoPandas	Always x, y: longitude, latitude	“always stored as (x, y), and thus as (lon, lat) order” (GeoPandas)
PostGIS	x, y: longitude, latitude	“For geodetic coordinates, X is longitude and Y is latitude” (PostGIS)

A transformation that follows the authority order and receives longitude first reads 73.8567 as a latitude. Step through one such call, then run a stand-in that follows the same rule as pyproj’s Transformer:

Interactive · message by message Pune's longitude and latitude through a transformer, twice Step through the same call without and with the option that makes the transformer read x before y.

The drawing shows the whole scene; its 5 steps are listed below.

Every step as text

The pair. The code holds Pune as 73.8567, 18.5204, longitude first, the order of GeoJSON and of most GIS code.
Authority order. The transformer is created without options, so it follows the axis order EPSG defines for EPSG:4326, latitude first.
A point off Norway. It transforms the point 73.86° N, 18.52° E and returns easting −1,009,879.90 m, northing 8,984,276.55 m, 7,000 km from Pune. Nothing reports an error.
Ask for x, y. Created with always_xy=True, the same transformer reads the first number as longitude and the second as latitude.
Pune. It returns easting 379,320.82 m, northing 2,048,144.62 m, the same UTM coordinates as GeographicLib's converter gives for Pune.

One pair of numbers, two axis orders Python · reproject/axis_order.py

"""The same pair of numbers through a transformation from EPSG:4326 that follows the authority's axis order
(latitude first) unless told to use x, y order, as pyproj's Transformer and DuckDB's ST_Transform do."""
from crs import distance, to_utm


def transform(first, second, always_xy=False):
    """EPSG:4326 to UTM 43N: (latitude, longitude) by default; (longitude, latitude) with always_xy=True."""
    lat, lon = (second, first) if always_xy else (first, second)
    return to_utm(lat, lon, 43)


pune = (73.8567, 18.5204)  # longitude, latitude, the order of GeoJSON and of most GIS code
right = transform(*pune, always_xy=True)
wrong = transform(*pune)
print(f"always_xy=True:  easting {right[0]:,.2f} m, northing {right[1]:,.2f} m")
print(f"authority order: easting {wrong[0]:,.2f} m, northing {wrong[1]:,.2f} m")
print(f"The wrong pair is the point {pune[0]}° N, {pune[1]}° E,"
      f" {distance(18.5204, 73.8567, 73.8567, 18.5204) / 1000:,.0f} km from Pune")
print(f"Latitude first, as the authority defines EPSG:4326: {transform(18.5204, 73.8567)[0]:,.2f} m easting")

Output

always_xy=True:  easting 379,320.82 m, northing 2,048,144.62 m
authority order: easting -1,009,879.90 m, northing 8,984,276.55 m
The wrong pair is the point 73.8567° N, 18.5204° E, 7,000 km from Pune
Latitude first, as the authority defines EPSG:4326: 379,320.82 m easting

Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 axis_order.py

The wrong call returns easting −1,009,879.90 m and northing 8,984,276.55 m, the UTM coordinates of a point off Norway, 7,000 km from Pune. Both readings are valid numbers, so nothing fails. Set the order once where the data enters the code, always_xy=True in pyproj, always_xy := true in DuckDB, and keep longitude first everywhere after.

Between two datums there may be several operations, or none

A change of datum needs a coordinate operation: A method with parameters that changes coordinates from one CRS to another: a conversion such as a map projection, or a transformation between datums, which has an area of use and a stated accuracy. with parameters, and the EPSG dataset can hold several for one pair of datums, each with its own area of use and accuracy, or none at all. pyproj’s TransformerGroup lists the candidates, “sorted with the most relevant ones first”, each with an accuracy that is “-1 if unknown” (pyproj); QGIS can ask which to use, and shows each one’s accuracy and area of use (QGIS).

Kalianpur 1937 shows the problem. EPSG has two operations to WGS 84: EPSG:1155 for onshore Bangladesh, with a stated accuracy of 18 m, and EPSG:1533 for the Moattama area of Myanmar, with 5 m (EPSG:1155, EPSG:1533). The figure draws their bounding boxes; the program applies both to a point in Kolkata and a point in Pune:

The bounding boxes of the two Kalianpur 1937 operations, with Kolkata and Pune

Text description of the diagram

A plot of longitude 70° E to 100° E across and latitude 8° N to 28° N up, with grid lines every 5 degrees and no map outline. A solid blue rectangle from 20.52° to 26.64° N and 88.01° to 92.67° E is the bounding box of EPSG:1155, the operation for Bangladesh onshore, with a stated accuracy of 18 m. A dashed black rectangle from 9.48° to 17.87° N and 93.94° to 99.66° E is the bounding box of EPSG:1533, the operation for the Moattama area of Myanmar, with a stated accuracy of 5 m.

A dot for Kolkata, at 22.57° N, 88.36° E, lies just inside the western edge of EPSG:1155's rectangle, though Kolkata is not in Bangladesh. A dot for Pune, at 18.52° N, 73.86° E, lies outside both rectangles. A note says that the rectangles are bounding boxes, not the areas of use themselves.

Kalianpur 1937 to WGS 84: two operations and none Python · reproject/operations.py

"""Kalianpur 1937 to WGS 84: the operations the EPSG dataset offers, where each applies, and what each does."""
from crs import distance, shift

OPERATIONS = [  # code, area of use, stated accuracy (m), bounding box (south, west, north, east), shift (m)
    ("EPSG:1155", "Bangladesh - onshore", 18, (20.52, 88.01, 26.64, 92.67), (282.0, 726.0, 254.0)),
    ("EPSG:1533", "Myanmar - Moattama area", 5, (9.48, 93.94, 17.87, 99.66), (214.0, 804.0, 268.0)),
]
PLACES = {"Kolkata": (22.5726, 88.3639), "Pune": (18.5204, 73.8567)}

for place, (lat, lon) in PLACES.items():
    print(f"{place}, {lat}, {lon} on Kalianpur 1937:")
    results = {}
    for code, area, accuracy, (s, w, n, e), shift_m in OPERATIONS:
        inside = s <= lat <= n and w <= lon <= e
        results[code] = shift(lat, lon, "Everest 1830 (1937 Adjustment)", shift_m)
        moved = distance(lat, lon, *results[code])
        print(f"  {code}, {area}, accuracy {accuracy} m")
        print(f"    inside its bounding box: {'yes' if inside else 'no'}; moves the point {moved:,.1f} m")
    apart = distance(*results["EPSG:1155"], *results["EPSG:1533"])
    print(f"  the two put the point {apart:,.1f} m apart; a ballpark operation does not move it")

Output

Kolkata, 22.5726, 88.3639 on Kalianpur 1937:
  EPSG:1155, Bangladesh - onshore, accuracy 18 m
    inside its bounding box: yes; moves the point 274.0 m
  EPSG:1533, Myanmar - Moattama area, accuracy 5 m
    inside its bounding box: no; moves the point 202.2 m
  the two put the point 72.1 m apart; a ballpark operation does not move it
Pune, 18.5204, 73.8567 on Kalianpur 1937:
  EPSG:1155, Bangladesh - onshore, accuracy 18 m
    inside its bounding box: no; moves the point 125.5 m
  EPSG:1533, Myanmar - Moattama area, accuracy 5 m
    inside its bounding box: no; moves the point 101.9 m
  the two put the point 87.1 m apart; a ballpark operation does not move it

Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 operations.py

At Kolkata the two operations move the point by 274.0 and 202.2 m and disagree by 72.1 m. Kolkata even lies inside EPSG:1155’s bounding box, although Bangladesh is its area of use: a rectangle is a quick test, not the area itself. At Pune neither operation applies. A tool that finds none may use a ballpark transformation: PROJ's fallback between two geographic CRSs when it knows no operation: it handles only axis order, units and the prime meridian and attempts no datum shift, so its error can reach a few hundred metres., which, PROJ says, “does not attempt any datum shift” and can be wrong by a few hundred metres (PROJ glossary): here it leaves Pune where the numbers say, 101.9 and 125.5 m from where the two operations would put it. pyproj’s allow_ballpark=False makes it refuse instead. Choose the operation whose area of use contains the data, prefer the better accuracy, record its code, and stop when none fits.

Checkpoint

Check: Choosing a coordinate operation

Question 1 of 2 Which rule picks a datum transformation for a layer?
Choose one answer.

Whichever the tool picks, since tools know best An operation whose area of use contains the data, then the best stated accuracy; record its code The first operation in the list The most accurate operation, wherever its area of use is

Show the answer: Question 1 of 2

Answer: An operation whose area of use contains the data, then the best stated accuracy; record its code

An operation applies inside its area of use, with its stated accuracy. Name the chosen one so the result can be reproduced.
Question 2 of 2 A layer in Pune is in Kalianpur 1937, for which EPSG has no operation to WGS 84 in India. A tool transforms it without complaint. What has most likely happened?
Choose one answer.

It used a ballpark operation, with no datum shift, leaving errors of about 100 m or more It used EPSG:1156, the Kalianpur 1975 operation It used EPSG:1155, the operation for Bangladesh Nothing is wrong; the transformation is exact

Show the answer: Question 2 of 2

Answer: It used a ballpark operation, with no datum shift, leaving errors of about 100 m or more

When no operation applies, a tool may still return numbers. Refuse instead, for example with allow_ballpark=False in pyproj, and ask the provider for control points or a datum the data can be moved from.

QGIS reprojects for display; analysis needs one CRS

QGIS “transparently reprojects all layers contained within your project into the project’s CRS”, so layers in different CRSs line up on the canvas (QGIS). That on-the-fly reprojection: Reprojecting layers in other CRSs into a project's CRS only for drawing them, as QGIS does, without changing the layers' data or the CRS they are measured in. changes only the drawing. Measurements and processing still use each layer’s own CRS: the expression area($geometry) is “always planimetric in the Spatial Reference System (SRS) of this geometry”, so on a layer in EPSG:4326 it returns square degrees, while $area follows the project’s ellipsoid setting (QGIS). Step through two layers in one project:

Interactive · stage by stage Two layers in one QGIS project, drawn and measured Step through what QGIS reprojects for you, what it does not, and the step that makes measurements consistent.

The drawing shows the whole scene; its 6 steps are listed below.

Every step as text

Two CRSs. A layer of wards arrives in EPSG:4326, in degrees, and a layer of roads in EPSG:32643, UTM 43N, in metres.
One project CRS. The project uses EPSG:3857 to match a web basemap. QGIS transparently reprojects every layer into the project's CRS for display.
They line up. On the canvas the wards and the roads sit where they should. Nothing about the drawing shows that the two layers use different CRSs.
Measuring. An expression such as area($geometry) is computed in the layer's own CRS, so on the wards it returns square degrees, while $area follows the project's ellipsoid setting.
One CRS for analysis. Reproject the wards with the Reproject layer algorithm, native:reprojectlayer, into EPSG:32643, choosing the coordinate operation if several exist, and save the result.
Consistent results. Now both layers share one projected CRS, so areas, lengths, buffers and joins all work in metres in the same grid.

Before analysis, reproject the layers into one projected CRS and save them. QGIS’s Reproject layer algorithm, native:reprojectlayer, takes an optional coordinate operation for “control over the exact transformation pipeline” (QGIS), so the operation chosen in the previous section can be named there too.

Round trips are safe; doubled steps are not

Transformations themselves are precise. The program sends 1,000 random points in zone 43N’s range to UTM and back, then makes the mistake that pipelines make: it applies the Kalianpur 1975 shift, EPSG:1156, to data that has already had it:

A round trip, and a transformation applied twice Python · reproject/round_trip.py

"""What a round trip costs, and what a transformation applied twice costs."""
import random

from crs import distance, from_utm, shift, to_utm

rng = random.Random(43)
points = [(rng.uniform(8.0, 35.0), rng.uniform(72.0, 78.0)) for _ in range(1000)]  # inside UTM zone 43N
worst = max(distance(lat, lon, *from_utm(*to_utm(lat, lon, 43), 43)) for lat, lon in points)
print(f"1,000 points to UTM 43N and back: largest error {worst * 1e9:.1f} nanometres")

KALIANPUR_1975 = "Everest 1830 (1975 Definition)"
SHIFT = (295.0, 736.0, 257.0)  # EPSG:1156, Kalianpur 1975 to WGS 84
once = shift(18.5204, 73.8567, KALIANPUR_1975, SHIFT)
twice = shift(*once, KALIANPUR_1975, SHIFT)  # the output treated as Kalianpur 1975 again
print(f"Pune moved once: {distance(18.5204, 73.8567, *once):.1f} m")
print(f"moved twice: a further {distance(*once, *twice):.1f} m, with nothing to say it happened")

Output

1,000 points to UTM 43N and back: largest error 37.8 nanometres
Pune moved once: 130.0 m
moved twice: a further 130.0 m, with nothing to say it happened

Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 round_trip.py

The round trip loses at most 37.8 nanometres. Applying the shift a second time moves Pune another 130.0 m, and nothing reports it, because the output of the first step is a valid input for the second. Transform once, from the CRS the data is documented in, write the operation’s code into the layer’s metadata, and let every later step read that label rather than guess.

Key idea

Errors in reprojection come from the steps around the formulas: a wrong label, a swapped pair, a missing or wrong operation, a step done twice. Check the label, fix the order at the boundary, name the operation, and transform once.

Checkpoint

Check: On-the-fly reprojection in QGIS · Round trips and doubled steps

Question 1 of 2 In a QGIS project in EPSG:3857, a layer of wards in EPSG:4326 lines up with the basemap. The expression area($geometry) on a ward gives 0.0003. Why?
Choose one answer.

Web Mercator distorts the area The ward really is 0.0003 m² The project's on-the-fly reprojection shrank the ward area() works in the layer's own CRS, so it returns square degrees

Show the answer: Question 1 of 2

Answer: area() works in the layer's own CRS, so it returns square degrees

Layers line up on screen because QGIS reprojects them for display. Measure after reprojecting and saving the data in a projected CRS, or use $area with an ellipsoid set.
Question 2 of 2 According to the program, how many metres further does Pune move when EPSG:1156 is applied a second time? Give one decimal.
Type a number.
Your answer m

Show the answer: Question 2 of 2

Answer: 130 m (anything from 129.95 to 130.05 counts)

The second application moves the point another 130.0 m. The output of the first step is a valid input for the second, so nothing reports the mistake.

Practise: check a layer before you reproject it

Build the two checks this lesson relies on: one that catches a layer labelled with the wrong CRS, and one that picks the operation for a point, or refuses when none applies:

Code lab · 2 steps

Check a layer before you reproject it

Catch a layer whose numbers cannot be in the CRS it claims, then pick a datum transformation by its bounding box and accuracy, or refuse.

Catch a layer labelled with the wrong CRS To do
Pick a datum transformation, or refuse To do

Exercise · Easy · Python

Catch a layer labelled with the wrong CRS

Write label_problem(points, label) in labels.py. points is a list of (x, y) pairs as a layer stores them, and label is the CRS the layer claims: "EPSG:4326", or a northern UTM zone such as "EPSG:32643". Return:

None when every point could be in that CRS;
"not degrees" when the label is EPSG:4326 and any x is outside −180 to 180 or any y outside −90 to 90;
"not UTM metres" when the label is a UTM zone, "EPSG:326" and two digits, and any easting is outside 100,000 to 900,000 m or any northing outside 0 to 9,400,000 m.

Values on the limits are inside. Raise a ValueError for any other label. A range test cannot prove a label right, but it catches the two mistakes of the lesson: degrees labelled as metres, and metres labelled as degrees.

Starter code · labels.py

def label_problem(points, label):
    """None if every (x, y) could be in the CRS the label names, else "not degrees" or "not UTM metres"."""
    # Replace this line: test every point against the ranges of EPSG:4326 or of a UTM north zone, EPSG:326zz.
    return None

The sample tests · test_labels.py

from labels import label_problem

PUNE_DEGREES = [(73.8567, 18.5204), (73.8601, 18.5310)]
PUNE_UTM = [(379320.82, 2048144.62), (379684.15, 2049318.40)]


def test_right_labels_pass():
    """accepts degrees labelled EPSG:4326 and UTM metres labelled with their zone"""
    assert label_problem(PUNE_DEGREES, "EPSG:4326") is None
    assert label_problem(PUNE_UTM, "EPSG:32643") is None


def test_degrees_labelled_utm():
    """catches degrees labelled as a UTM zone, the set_crs mistake of the lesson"""
    assert label_problem(PUNE_DEGREES, "EPSG:32643") == "not UTM metres"


def test_metres_labelled_degrees():
    """catches metres labelled EPSG:4326"""
    assert label_problem(PUNE_UTM, "EPSG:4326") == "not degrees"


def test_one_bad_point_and_the_limits():
    """flags a layer when any one point is out of range, and accepts values on the limits"""
    assert label_problem(PUNE_UTM + [(73.8567, 18.5204)], "EPSG:32643") == "not UTM metres"
    assert label_problem([(180.0, 90.0), (-180.0, -90.0)], "EPSG:4326") is None
    assert label_problem([(100000.0, 0.0), (900000.0, 9400000.0)], "EPSG:32644") is None
    assert label_problem([(99999.0, 10.0)], "EPSG:32644") == "not UTM metres"


def test_unknown_label():
    """refuses a label it cannot check, including a UTM code with a digit missing"""
    for label in ("EPSG:7755", "EPSG:3264"):
        try:
            label_problem(PUNE_UTM, label)
        except ValueError:
            continue
        assert False, f"{label} should raise ValueError"


def test_one_axis_out_of_range():
    """flags a point when only one of its two numbers is out of range"""
    assert label_problem([(73.8567, 2048144.62)], "EPSG:4326") == "not degrees"
    assert label_problem([(379320.82, -10.0)], "EPSG:32643") == "not UTM metres"

A hint

all(-180 <= x <= 180 and -90 <= y <= 90 for x, y in points) tests every point at once. For the UTM label, check label.startswith("EPSG:326") and its length, then test the easting and northing ranges the same way.

The sample tests run on this device, in your browser (Pyodide): nothing is sent to mysmartcopilot.com. The first run downloads Python (about 13.5 MB), which is kept for the next runs. A check in your browser is feedback for you, not proof that the code is right for every input.

Exercise · Medium · Python

Pick a datum transformation, or refuse

Write pick_operation(lat, lon, operations) in operations.py. Each operation is a dict with a "code", an "accuracy" in metres (or None when it is unknown) and a "box", its bounding box as (south, west, north, east) in degrees. Return the code of the operation to try first for the point:

keep only the operations whose box holds the point (a point on an edge is inside);
among those, return the one with the smallest known accuracy;
if the boxes that hold the point all have an unknown accuracy, return the first of them;
if no box holds the point, return None: refuse, rather than fall back to an operation without a datum shift.

This is the quick test of the lesson: a box that does not hold the point rules an operation out, while a box that holds it still leaves the area of use to check, as the lesson found for Kolkata, inside the box of the operation for Bangladesh. The sample tests use the two Kalianpur 1937 operations of the lesson, at Dhaka and at Pune.

Starter code · operations.py

def pick_operation(lat, lon, operations):
    """The code of the most accurate operation whose box (south, west, north, east) holds the point, or None."""
    # Replace this line: keep the operations whose box holds the point, then pick the smallest known accuracy.
    return operations[0]["code"] if operations else None

The sample tests · test_operations.py

from operations import pick_operation

KALIANPUR_1937 = [
    {"code": "EPSG:1155", "accuracy": 18, "box": (20.52, 88.01, 26.64, 92.67)},
    {"code": "EPSG:1533", "accuracy": 5, "box": (9.48, 93.94, 17.87, 99.66)},
]


def test_one_operation_fits():
    """picks the only operation whose box holds Dhaka, in Bangladesh"""
    assert pick_operation(23.8103, 90.4125, KALIANPUR_1937) == "EPSG:1155"


def test_none_fits():
    """refuses, with None, when no box holds the point, as for Pune"""
    assert pick_operation(18.5204, 73.8567, KALIANPUR_1937) is None
    assert pick_operation(18.5204, 73.8567, []) is None


def test_most_accurate_wins():
    """prefers the smaller stated accuracy when several boxes hold the point"""
    ops = [{"code": "A", "accuracy": 22, "box": (8.0, 68.0, 36.0, 98.0)},
           {"code": "B", "accuracy": 3, "box": (18.0, 73.0, 19.0, 74.0)},
           {"code": "C", "accuracy": None, "box": (18.0, 73.0, 19.0, 74.0)}]
    assert pick_operation(18.5204, 73.8567, ops) == "B"
    assert pick_operation(26.1445, 91.7362, ops) == "A"


def test_edges_and_unknown_accuracy():
    """counts a point on a box's edge as inside, and falls back to an operation of unknown accuracy only when alone"""
    assert pick_operation(20.52, 88.01, KALIANPUR_1937) == "EPSG:1155"
    lonely = [{"code": "C", "accuracy": None, "box": (18.0, 73.0, 19.0, 74.0)}]
    assert pick_operation(18.5, 73.5, lonely) == "C"


def test_both_limits_must_hold():
    """needs the latitude and the longitude inside the box, not one of them"""
    assert pick_operation(22.5726, 73.8567, KALIANPUR_1937) is None
    assert pick_operation(30.0, 90.0, KALIANPUR_1937) is None

A hint

Filter first: fits = [op for op in operations if south <= lat <= north and west <= lon <= east], reading the four limits from op["box"]. Then min(…, key=lambda op: op["accuracy"]) over the ones whose accuracy is not None.

Then test yourself on the five ideas, from a wrong label to a doubled step, starting with the one your answers leave weakest:

Interview questions

Warm-up (fresher): What is the difference between set_crs and to_crs in GeoPandas?

set_crs attaches a label: it tells GeoPandas which CRS the existing numbers are in, and changes none of them. to_crs transforms: it computes new numbers for the same places in another CRS, and updates the label. Use set_crs when a layer has no CRS or a wrong one, and to_crs to move a correctly labelled layer into the CRS an analysis needs. Mixing them up is silent: labelling Pune’s degrees as UTM zone 43N reads 73.86 and 18.52 as metres and puts the point 2,081 km away near the equator, while to_crs(32643) gives the right easting of 379,320.82 m and northing of 2,048,144.62 m.

Two questions in the module’s bank take this further: clinics that appear in the ocean after a “fix”, and a pipeline whose machines disagree about Kalianpur 1937. Each lists what it tests and a first step:

Model answers open with any pass: ₹79 for 1 week ($2.99), one payment that never renews.

Choose a pass

Fresher to mid

After a colleague "fixed the projection" of a layer of 2,000 clinics in Maharashtra, every clinic appears in the Indian Ocean near the equator, about 2,000 km to the south-west. What most likely happened, how do you confirm it, and how do you repair the layer?

Model answer, Pro. Model answers open with any pass; prices are in the notice on this page. First step and follow-ups
Mid to senior

A pipeline converts old survey points in Kalianpur 1937 for districts of India to WGS 84, and two machines produce positions about 100 m apart for the same points. Why can this happen, and what should the pipeline do instead?

Model answer, Pro. First step and follow-ups

Flashcards · 7 cards

Recall it, then check

Say the answer to yourself before you look: recalling it is what makes it stick.

What is the difference between set_crs and to_crs in GeoPandas, or ST_SetSRID and ST_Transform in PostGIS?

set_crs and ST_SetSRID attach a label saying which CRS the existing numbers are in and change none of them. to_crs and ST_Transform compute new numbers for the same places in another CRS and update the label.
How do you spot, and repair, a layer of degrees labelled as UTM zone 43N?

A range check: eastings in a UTM zone lie between about 100,000 and 900,000 m, so an easting of 73.86 m cannot be one, and Pune read that way lands near the equator, 2,081 km off. Put back the label that fits the numbers, EPSG:4326, then transform.
Which axis order do pyproj, DuckDB, GeoPandas and PostGIS expect for EPSG:4326 coordinates?

A pyproj Transformer and DuckDB's ST_Transform follow EPSG's order, latitude first, unless always_xy is set. GeoPandas and PostGIS always work in x and y, that is longitude, then latitude.
What happens when longitude-first coordinates reach a transformation that expects latitude first?

73.8567 is taken as a latitude, and Pune lands off Norway, about 7,000 km away, with no error raised. Fix the order once where data enters the code, always_xy=True in pyproj or always_xy := true in DuckDB, and keep longitude first after that.
How should you pick among several coordinate operations between two datums, or when there is none?

Take the one whose area of use contains the data, prefer the better stated accuracy and record its code. When none applies, refuse, for example with allow_ballpark=False in pyproj, instead of accepting a ballpark operation that makes no datum shift and can be hundreds of metres off.
Does QGIS's on-the-fly reprojection change what its measurements use?

No. It redraws layers in the project's CRS only. area($geometry) stays planimetric in the layer's own CRS, so a layer in EPSG:4326 gives square degrees. Reproject the layers into one projected CRS and save them before any analysis.
Which is the risk: a round trip to UTM and back, or a datum shift applied twice?

The round trip lost at most 37.8 nanometres. The doubled shift moved Pune another 130 m and nothing reported it, because the first step's output is valid input for the second. Transform once and write the operation's code into the metadata.

Key takeaways

Labelling a CRS (set_crs, ST_SetSRID) changes no coordinate and transforming (to_crs, ST_Transform) changes every one; a range check catches degrees labelled as UTM, which land 2,081 km from Pune.
pyproj and DuckDB follow EPSG’s latitude-first order for EPSG:4326 unless told to use x, y, while GeoPandas and PostGIS work in longitude, latitude; a swapped pair lands 7,000 km away.
Two datums may have several operations or none: choose by area of use and accuracy, record the code, and refuse rather than accept a ballpark operation with no datum shift.
QGIS reprojects layers on the fly for display only; measure and process after saving the data in one projected CRS.
A round trip loses nanometres, while a transformation applied twice adds its whole shift again, 130 m at Pune.

Terms in this lesson

ballpark transformation: PROJ's fallback between two geographic CRSs when it knows no operation: it handles only axis order, units and the prime meridian and attempts no datum shift, so its error can reach a few hundred metres.
coordinate operation: A method with parameters that changes coordinates from one CRS to another: a conversion such as a map projection, or a transformation between datums, which has an area of use and a stated accuracy.
on-the-fly reprojection: Reprojecting layers in other CRSs into a project's CRS only for drawing them, as QGIS does, without changing the layers' data or the CRS they are measured in.

References

GeoPandas user guide: managing projections (GeoPandas developers)
pyproj: Transformer and TransformerGroup (pyproj contributors)
PostGIS: ST_Transform (PostGIS Project Steering Committee)
PostGIS: ST_MakePoint (PostGIS Project Steering Committee)
DuckDB spatial extension: functions (ST_Transform) (DuckDB Foundation)
QGIS User Guide: working with projections (QGIS Project)
QGIS User Guide: Reproject layer (native:reprojectlayer) (QGIS Project)
QGIS User Guide: expression functions (area, $area) (QGIS Project)
PROJ glossary: ballpark transformation (PROJ contributors (OSGeo))
EPSG:1155, Kalianpur 1937 to WGS 84 (1) (MapTiler (EPSG.io), from the EPSG dataset)
EPSG:1533, Kalianpur 1937 to WGS 84 (2) (MapTiler (EPSG.io), from the EPSG dataset)
GeographicLib: online UTM conversion (GeoConvert) (GeographicLib (C. F. F. Karney))

Report a problem with this lesson

In this module: Coordinate reference systems

Earth's shape: ellipsoids, geoids and datums Free lesson, Pro answers: 25 minutes
Map projections and their distortions Free lesson, Pro answers: 20 minutes
EPSG codes, WKT2 and PROJJSON in practice Free lesson, Pro answers: 20 minutes
Reprojecting data without silent errors Free lesson, Pro answers: 20 minutes
Measuring distance and area correctly Free lesson, Pro answers: 20 minutes

Your country

Reprojecting data without silent errors

What you will learn

Before you start

What you will learn

Before you start

Labelling a CRS is not transforming it

Every engine has an axis-order rule

Between two datums there may be several operations, or none

QGIS reprojects for display; analysis needs one CRS

Round trips are safe; doubled steps are not

Practise: check a layer before you reproject it

Catch a layer labelled with the wrong CRS

Results of the sample tests

Pick a datum transformation, or refuse

Results of the sample tests

Interview questions

After a colleague "fixed the projection" of a layer of 2,000 clinics in Maharashtra, every clinic appears in the Indian Ocean near the equator, about 2,000 km to the south-west. What most likely happened, how do you confirm it, and how do you repair the layer?

A pipeline converts old survey points in Kalianpur 1937 for districts of India to WGS 84, and two machines produce positions about 100 m apart for the same points. Why can this happen, and what should the pipeline do instead?

Key takeaways

References