🔧 Fix Coastal & marine Raster & rasterio Coastal & Marine · step 16 of 20 Updated

Bathymetry depths are positive in one file and negative in another

Problem statement

Two bathymetry grids of the same area, and one has the sea at −40 and the other at +40. Combined without checking, the result is a seabed that rises where it should fall, a mean depth near zero, contours in the wrong place, and a hillshade lit from underneath.

Nothing raises. Both files are valid, both conventions are correct in their own field, and the combined array is full of plausible numbers. The only way to find it is to look, and the only way to stop it happening again is to normalise on read.

Quick answer

Normalise every grid to one convention as it is loaded, using the declared attribute where there is one and a probe where there is not:

import numpy as np, xarray as xr

def as_elevation(da, probe_lat=0.0, probe_lon=-150.0):
    """Return the array as elevation: positive up, sea negative."""
    declared = da.attrs.get("positive")
    if declared == "up":
        return da
    if declared == "down":
        return -da

    value = float(da.sel(lat=probe_lat, lon=probe_lon, method="nearest"))
    if value > 0:                      # deep Pacific is positive → it is depth
        return -da
    return da

z = as_elevation(xr.open_dataset("etopo_2022_60s.nc")["z"])
print(z.attrs.get("positive"), float(z.sel(lat=0, lon=-150, method="nearest")))
up -4462.0625

A point in the middle of the Pacific is about four and a half kilometres deep. If it comes back positive, the file stores depth; if negative, it stores elevation.

Triage of four symptoms of a mixed sign convention and the check that identifies each.
Four symptoms, one cause, and each is visible in a summary rather than in the data.

Step-by-step solution

1. Read the CF attribute first

positive: up means elevation, positive: down means depth. It is the authoritative statement and it is present in well-made files — ETOPO 2022 declares positive: up alongside units: meters and vert_crs_epsg: EPSG:3855.

2. Probe a known point when the attribute is missing

The centre of an ocean basin is unambiguous. (0, −150) is in the Pacific at −4,462 m; (0, −25) is in the Atlantic at −3,270 m; (−10, 80) is in the Indian Ocean at −5,343 m. One number settles the question.

3. Check the histogram, not the range

A grid holding land and sea is strongly bimodal. If most of the mass is positive and there is a small negative tail, it is elevation on a mostly-land tile or depth on a mostly-sea one — the range alone does not distinguish them, the shape does.

4. Look for the tell-tale symptoms

  • Mean depth near zero on a mixed mosaic.
  • An inverted hillshade — ridges look like valleys.
  • Contours in the wrong place, symmetric about the coast.
  • A sea-level mask that selects the land.

5. Normalise on read, not before write

Converting at the point of loading means every downstream function sees one convention, and the conversion is in one place. Converting files on disk leaves the question open for the next file that arrives.

6. Carry the convention into the output

Whatever you produce, write the CF attributes: positive, units and the vertical CRS. A file that does not declare its convention is the file that caused this.

7. Assert it in the pipeline

One assertion at the top of a processing function is enough, and it fails at the moment the wrong file arrives rather than three steps later.

Two panels showing the deep Pacific probe returning a negative value for elevation and a positive one for depth.
One probe point resolves the question when the attribute is missing.

Code examples

Example 1 — detect and report, before converting anything 

import numpy as np, xarray as xr

PROBES = {"pacific": (0.0, -150.0), "atlantic": (0.0, -25.0), "indian": (-10.0, 80.0)}

def sign_report(da):
    out = {"declared_positive": da.attrs.get("positive"),
           "units": da.attrs.get("units"),
           "vert_crs": da.attrs.get("vert_crs_name")}
    for name, (lat, lon) in PROBES.items():
        try:
            out[name] = float(da.sel(lat=lat, lon=lon, method="nearest"))
        except Exception:
            out[name] = None
    probes = [v for v in (out.get(k) for k in PROBES) if v is not None]
    if probes:
        out["inferred"] = "depth (positive down)" if np.mean(probes) > 0 else \
                          "elevation (positive up)"
    if out["declared_positive"] and out.get("inferred"):
        declared_up = out["declared_positive"] == "up"
        inferred_up = out["inferred"].startswith("elevation")
        if declared_up != inferred_up:
            out["conflict"] = "the declared attribute disagrees with the data"
    return out

print(sign_report(xr.open_dataset("etopo_2022_60s.nc")["z"]))

The conflict check is worth having: an array that was negated without its attribute being updated is a real and unpleasant case, and it is invisible to a probe alone.

Example 2 — assert before combining two grids 

import numpy as np, xarray as xr

def assert_same_convention(*arrays, probe=(0.0, -150.0)):
    signs = []
    for a in arrays:
        declared = a.attrs.get("positive")
        if declared in ("up", "down"):
            signs.append(declared)
            continue
        v = float(a.sel(lat=probe[0], lon=probe[1], method="nearest"))
        signs.append("down" if v > 0 else "up")
    if len(set(signs)) > 1:
        raise ValueError(f"mixed sign conventions: {signs}")
    return signs[0]

convention = assert_same_convention(grid_a, grid_b)
print("both grids are", convention)

Example 3 — a loader that normalises and records what it did 

import xarray as xr

def open_bathymetry(path, var="z", target="elevation"):
    ds = xr.open_dataset(path)
    a = ds[var]
    report = sign_report(a)
    is_depth = (report.get("declared_positive") == "down"
                or (report.get("declared_positive") is None
                    and report.get("inferred", "").startswith("depth")))

    if (target == "elevation") == is_depth:
        a = -a
        a.attrs["spatialworkflow_negated"] = "true"

    a.attrs["positive"] = "up" if target == "elevation" else "down"
    a.attrs.setdefault("units", "meters")
    a.attrs["spatialworkflow_source_convention"] = (
        "depth" if is_depth else "elevation")
    return a, report

Recording that the array was negated, in the attributes, is what stops the same file being negated twice by two different loaders in the same pipeline.

Explanation

Why both conventions exist and neither is going away

Hydrography measures depth below a datum because a mariner needs a positive number that must stay above the keel. Geodesy and geophysics measure height above a reference surface because that is continuous with the land and with the physics. Files move between the two communities constantly, and the CF positive attribute exists precisely because the format cannot assume either.

Why the range does not tell you

A grid running from −10,320 to 7,343 m could be elevation over land and sea, or depth over a basin with a trench and a seamount. A grid running from 0 to 5,000 could be depth in an ocean or elevation in a mountain range. The sign of a point you know settles it in one operation; the range never does.

Why an inverted hillshade is the most visible symptom

A hillshade computed from a negated surface lights the scene from the opposite side, and the human visual system reads the result as an inversion: ridges become valleys and craters become domes. It is the fastest visual check there is, and it is why anybody who has been caught once renders a hillshade before trusting a new bathymetry file.

Why normalising at the read boundary is the right place

Any other place leaves a window in which two conventions coexist. A single open_bathymetry used everywhere means the convention is established exactly once per file, the decision is recorded in the attributes, and every downstream function can state its assumption and rely on it.

Vertical steps for normalising a bathymetry grid: read the attribute, probe if absent, flag conflicts, negate and record, and write the convention on output.
The conflict check catches an array that was negated without its attribute being updated.

Edge cases or notes

  • A negated array with a stale attribute is worse than a missing attribute.
  • Some files store depth as positive with a _FillValue of 0. Check the mask too.
  • Nodata of −9999 looks like a trench under either convention.
  • Chart depths are positive by definition and referenced to a chart datum.
  • Do not negate in place twice. Record the conversion in the attributes.
  • Mixed-source mosaics need a per-tile check, not one for the mosaic.
  • positive applies to the coordinate too on a 3D ocean file.
  • Write the attribute on output. It is two lines and it ends the problem downstream.

FAQ

Is bathymetry positive or negative?

Either. Read the CF positive attribute: up means elevation with the sea negative, down means depth with the sea positive.

How do I check when the attribute is missing?

Probe a point you know is deep ocean — (0, −150) in the Pacific reads −4,462 m in ETOPO 2022. Positive means depth, negative means elevation.

Why does my mean depth come out near zero?

Because a mosaic mixes both conventions and the two halves cancel. Check each source tile rather than the mosaic.

Why does my hillshade look inverted?

Because the surface is negated relative to what the hillshade expects. It is the fastest visual check for this problem.

Should I convert the files on disk?

No. Normalise at the point of loading, so the decision is made once per file and recorded in the attributes.

What should I write in my own output?

positive, units and the vertical CRS, every time. A file that does not declare its convention is what caused the problem.