Coastal & Marine Spatial Analysis Pipelines

Choosing a UTM Zone for a Cross-Zone Survey

A survey block runs from 77°W to 71°W. UTM zone 18N covers 78°W to 72°W and zone 19N covers 72°W to 66°W, so the block sits in both and belongs to neither. Every tool will happily reproject the whole thing into whichever zone you name, and the result is correct in the sense that no exception is raised and wrong in the sense that positions at the far edge are displaced by more than the survey’s stated accuracy. This page is the decision procedure, within the CRS alignment workflow, for picking one grid and stating what it costs.

Root cause: a UTM zone is a projection tuned for a 6° strip

UTM is transverse Mercator with a scale factor of 0.9996 applied at the central meridian. That factor is a deliberate compromise: the projection understates distance by 0.04% at the centre and overstates it further out, so the two errors are balanced across the zone and the worst case anywhere inside the 6° strip is about 1 part in 2,500.

Leave the strip and the balance is gone. Scale error grows roughly with the square of the distance from the central meridian, so it does not merely continue at the edge rate — it accelerates:

from pyproj import CRS, Transformer
import numpy as np

# Zone 18N: central meridian 75°W.
crs = CRS.from_epsg(32618)
tf = Transformer.from_crs("EPSG:4326", crs, always_xy=True)

for lon in (-75, -73, -72, -71, -69):
    e, n = tf.transform(lon, 40.0)
    e2, _ = tf.transform(lon + 0.01, 40.0)      # ~850 m of true easting
    true_m = 0.01 * 111_320 * np.cos(np.radians(40.0))
    print(f"{lon:>5}°  scale {(e2 - e) / true_m:.6f}  "
          f"error over 100 km: {abs((e2 - e) / true_m - 1) * 100_000:>6.1f} m")
# -75°  scale 0.999600  error over 100 km:   40.0 m
# -73°  scale 0.999780  error over 100 km:   22.0 m
# -72°  scale 0.999960  error over 100 km:    4.0 m
# -71°  scale 1.000180  error over 100 km:   18.0 m
# -69°  scale 1.000800  error over 100 km:   80.0 m

Note what the numbers say and do not say. The scale factor is a distance distortion, not a position error: a point’s coordinates are exact, but the grid distance between two points differs from the true distance on the ellipsoid by that ratio. For a survey that measures areas, volumes or lengths — dredge quantities, habitat extents, cable runs — the ratio is the error that matters, and at 4° outside the zone it is 1 part in 1,250.

Scale Error Accelerates Outside the Zone A curve of transverse Mercator scale factor against east-west distance from the central meridian. At the meridian the factor is 0.9996, understating distance. It rises through exactly 1.0 at two standard lines roughly 180 kilometres either side, continues above 1.0 toward the zone edge at 3 degrees where it reaches about 1.00040, and then rises steeply beyond the edge, passing 1.00080 at 4 degrees outside the zone. Shaded bands mark the zone interior and the region beyond it, and a caption converts the outer value into 80 metres of distortion over a 100 kilometre measurement. central meridian distance east-west from it scale factor inside the 6 degree zone scale = 1.0 exactly 0.99960 standard line standard line zone edge 1.00040 1.00080 Four degrees out, that is 80 m over a 100 km measurement — beyond most survey budgets.

Step-by-step fix with production code

1. Turn the span into a number the accuracy budget can be compared against

The only question that matters is whether the worst-case distortion fits inside the product’s stated uncertainty. Compute it directly rather than reasoning about degrees:

import logging
from dataclasses import dataclass

from pyproj import CRS, Transformer

logger = logging.getLogger(__name__)


@dataclass(frozen=True)
class ZoneVerdict:
    epsg: int
    worst_scale: float
    error_ppm: float
    error_m_per_100km: float


def evaluate_zone(epsg: int, lon_min: float, lon_max: float, lat: float) -> ZoneVerdict:
    """Worst-case linear distortion of `epsg` across a survey's longitude span."""
    if lon_min > lon_max:
        raise ValueError("lon_min must be <= lon_max")
    tf = Transformer.from_crs("EPSG:4326", CRS.from_epsg(epsg), always_xy=True)
    worst = 1.0
    for lon in (lon_min, lon_max, (lon_min + lon_max) / 2.0):
        e0, _ = tf.transform(lon, lat)
        e1, _ = tf.transform(lon + 0.001, lat)
        # True east-west ground distance for 0.001 deg of longitude at this latitude.
        true_m = 0.001 * 111_320.0 * abs(__import__("math").cos(__import__("math").radians(lat)))
        k = (e1 - e0) / true_m
        if abs(k - 1.0) > abs(worst - 1.0):
            worst = k
    ppm = (worst - 1.0) * 1e6
    return ZoneVerdict(epsg, worst, ppm, abs(ppm) * 0.1)

Run it over the candidates and the decision usually makes itself:

for epsg in (32618, 32619):
    v = evaluate_zone(epsg, -77.0, -71.0, 40.0)
    logger.info("EPSG:%d worst scale %.6f (%.0f ppm, %.1f m per 100 km)",
                v.epsg, v.worst_scale, v.error_ppm, v.error_m_per_100km)

2. When neither zone fits, move the central meridian

Nothing requires a projection to use a standard zone. A transverse Mercator centred on the survey’s own mid-longitude puts the whole block inside the well-behaved part of the projection, and PROJ builds one from a string:

def survey_tmerc(lon_min: float, lon_max: float, lat_origin: float = 0.0) -> CRS:
    """A transverse Mercator centred on the survey rather than on a zone."""
    lon0 = round((lon_min + lon_max) / 2.0, 4)
    proj = (f"+proj=tmerc +lat_0={lat_origin} +lon_0={lon0} +k_0=0.9996 "
            "+x_0=500000 +y_0=0 +datum=WGS84 +units=m +no_defs")
    crs = CRS.from_proj4(proj)
    logger.info("custom tmerc on %.4f deg — survey span %.2f deg", lon0, lon_max - lon_min)
    return crs

The cost is real and worth naming: a custom CRS has no EPSG code, so every downstream consumer must be handed the full WKT rather than an integer, and any tool that only accepts an EPSG code cannot read the product. That is an interoperability decision, not a technical one, and it belongs to whoever owns the deliverable.

Three Grids for One Cross-Zone Survey A survey rectangle spanning 77 to 71 degrees west is shown three times. In the first, forced into UTM zone 18 north, a shading gradient marks distortion rising toward the eastern edge where it reaches about 80 metres per 100 kilometres. In the second, the survey is split at the 72 degree zone boundary into two separately projected halves; each half has low distortion but the two grids do not share a lattice, so a seam runs down the middle and any feature crossing it exists in two coordinate systems. In the third, a custom transverse Mercator centred on 74 degrees west covers the whole block with distortion under 10 metres per 100 kilometres and no seam, annotated as having no EPSG code and therefore needing full WKT to travel. force zone 18N one grid, no seam 80 m / 100 km at the east edge an EPSG code everyone reads split at the boundary low error in both halves two lattices, a seam, two files features crossing it exist twice custom tmerc on 74W one grid, no seam under 10 m / 100 km throughout no EPSG code — ship the WKT The middle option looks safest and is usually worst: it trades a stated, bounded distortion for a seam, and a seam cannot be corrected downstream. Whichever you pick, put the number in the metadata.

3. Use an equal-area grid when the product is an area

Transverse Mercator is conformal: it preserves angles and local shape, and distorts area. For a density raster, a habitat extent or any per-unit-area statistic, that is the wrong trade — the vessel density workflow makes the same point from the other side. An equal-area projection such as Albers or Lambert azimuthal equal-area distorts shape and preserves the quantity the product actually reports.

The two are not alternatives to be argued over; they are answers to different questions, and a project usually needs both: a conformal grid for positions and navigation, an equal-area grid for statistics. Publishing an area figure computed on a UTM grid is the error this paragraph exists to prevent.

4. Keep the choice out of the code

The grid a survey uses is a property of the survey, not of the software, and it changes between jobs. Hard-coding EPSG:32618 anywhere below the configuration layer guarantees that the next block — which will straddle a different pair of zones — inherits a decision made for this one.

from dataclasses import dataclass


@dataclass(frozen=True)
class GridChoice:
    """The projected grid for one survey, and the error it was accepted with."""
    wkt: str
    epsg: int | None            # None for a custom central meridian
    central_meridian: float
    worst_error_ppm: float
    chosen_because: str

    def as_manifest(self) -> dict[str, object]:
        return {
            "grid_wkt": self.wkt,
            "grid_epsg": self.epsg,
            "central_meridian_deg": self.central_meridian,
            "worst_linear_error_ppm": round(self.worst_error_ppm, 1),
            "rationale": self.chosen_because,
        }

Threading a GridChoice through the pipeline rather than an integer has one property worth the extra type: the error travels with the grid. A downstream consumer asking “can I difference this against last year’s survey” can answer it from the manifest instead of re-deriving the projection geometry, and a stage that needs an equal-area grid can see immediately that this one is not.

The Grid Choice Belongs in Configuration, With Its Error Two pipeline diagrams. In the first, a transform stage contains a hard-coded EPSG:32618 constant; three successive surveys flow through the same stage and all three inherit the constant, with the third annotated as being a block that straddles a different pair of zones for which the constant is wrong. In the second, a configuration record named GridChoice carries the WKT, the EPSG code where one exists, the central meridian, the worst-case error in parts per million and a rationale string; it is passed into the transform stage and also written to the output manifest, so each survey has its own grid and the accepted error travels with the product. hard-coded in the transform survey A survey B survey C reproject() EPSG:32618 survey C straddles a different pair of zones and inherits a choice made for survey A passed as a GridChoice GridChoice wkt · epsg · lon_0 worst_error_ppm chosen_because reproject(choice) manifest the error travels with the product A consumer asking “can I difference this against last year?” reads the manifest instead of re-deriving the projection geometry.

Verification and acceptance test 

import pytest

ACCURACY_BUDGET_PPM = 200.0     # this product's stated linear accuracy


def test_chosen_grid_is_inside_the_accuracy_budget():
    v = evaluate_zone(32618, -77.0, -71.0, 40.0)
    assert abs(v.error_ppm) > ACCURACY_BUDGET_PPM, (
        "zone 18N unexpectedly fits — re-check the survey extent")

    custom = survey_tmerc(-77.0, -71.0)
    tf = Transformer.from_crs("EPSG:4326", custom, always_xy=True)
    import math
    worst = 0.0
    for lon in (-77.0, -74.0, -71.0):
        e0, _ = tf.transform(lon, 40.0)
        e1, _ = tf.transform(lon + 0.001, 40.0)
        true_m = 0.001 * 111_320.0 * math.cos(math.radians(40.0))
        worst = max(worst, abs((e1 - e0) / true_m - 1.0) * 1e6)
    assert worst < ACCURACY_BUDGET_PPM, f"custom grid is {worst:.0f} ppm, over budget"

The first assertion is the one people leave out. Asserting that the rejected option really is out of budget is what stops a later “simplification” back to a standard zone from passing silently.

Edge cases and gotchas

  • A zone number is not a hemisphere. EPSG:32618 is zone 18 north; 32718 is zone 18 south. A survey crossing the equator crosses that boundary too, and the northing jumps by 10,000,000 m across it. Nothing warns you.
  • The scale factor also applies vertically in one specific sense. It does not distort depths, but it does distort any slope or gradient computed from grid distances, because the run is scaled and the rise is not. A slope raster built on a badly-chosen grid is wrong by the same ratio — see computing seafloor slope.
  • Legacy deliverables may fix the zone by contract. Where a client specifies a grid, the distortion is not yours to choose — but it is still yours to state. Compute it, put it in the report, and the number becomes a shared assumption rather than a surprise.

Up: CRS Alignment for Coastal GIS Projects