Coastal & Marine Spatial Analysis Pipelines

Converting NetCDF to Cloud-Optimized GeoTIFF

Turning a model output or a gridded product into something a web map or a tile server can read a window of is a routine step, and rioxarray makes it look like one line. The line works, and it silently produces a transposed, upside-down, CRS-less raster often enough that the conversion deserves a procedure. This page is that procedure, within the NetCDF versus GeoTIFF routing topic.

Root cause: NetCDF describes axes, GeoTIFF assumes them

A CF NetCDF variable carries named dimensions in any order, with coordinate variables that say what each one means and which way it runs. A GeoTIFF carries a geotransform: an origin and a pixel size, with the implicit convention that rows increase southward from a north-west origin. The conversion has to turn the first description into the second, and four pieces of information have no default that is right more often than it is wrong.

import xarray as xr

ds = xr.open_dataset("hycom_ssh.nc")
print(ds["surf_el"].dims)          # ('time', 'lat', 'lon')  — three dims, not two
print(float(ds.lat[0]), float(ds.lat[-1]))   # 25.0 45.0 — ascending, i.e. south-up
print(ds["surf_el"].attrs.get("grid_mapping"))   # 'crs'  — the CRS is in another variable
print(ds["surf_el"].attrs.get("_FillValue"))     # -30000  — packed, not NaN

Every one of those is a decision the writer must be told about. Left alone, the time dimension becomes bands in an arbitrary order, the ascending latitude produces a vertically mirrored image, the CRS defaults to whatever the reader guesses, and the fill value becomes a real elevation of −30000.

Four Assumptions the GeoTIFF Makes and the NetCDF Does Not A four-row comparison table. Row one, dimension order: NetCDF names its dimensions in any order such as time, lat, lon, while GeoTIFF assumes band, row, column, and the symptom of a mismatch is a transposed image. Row two, latitude direction: NetCDF records whether latitude ascends or descends, GeoTIFF assumes rows increase southward from a north-west origin, and the symptom is a vertically mirrored image that still georeferences to the correct bounding box. Row three, coordinate reference system: NetCDF points at a separate grid_mapping variable, GeoTIFF expects the CRS in its own header, and the symptom is a raster that opens with no CRS or with an assumed one. Row four, fill value: NetCDF carries an underscore FillValue that may be a packed integer, GeoTIFF carries a nodata tag, and the symptom is a fill value read as a real measurement. property NetCDF states it GeoTIFF assumes symptom if unreconciled dimension order(time, lat, lon)(band, row, col)transposed image latitude directionascending or descendingnorth-up rowsmirrored vertically CRSa grid_mapping variablein its own headerno CRS, or a guess fill value_FillValue, often packeda nodata tagfill read as data Only the second symptom is visible without checking. A mirrored raster still georeferences to the right bbox, so it overlays a basemap convincingly while every feature is in the wrong place. The other three are silent until a downstream consumer computes something from them.

Step-by-step fix with production code

1. Reduce to a 2D or banded array, deliberately 

import logging

import rioxarray  # noqa: F401 — registers the .rio accessor
import xarray as xr

logger = logging.getLogger(__name__)


def select_slice(ds: xr.Dataset, var: str, **selectors) -> xr.DataArray:
    """One variable, reduced to (y, x) by explicit selection on every other dim."""
    da = ds[var]
    if selectors:
        da = da.sel(**selectors, method="nearest")
    extra = [d for d in da.dims if d not in ("y", "x", "lat", "lon", "latitude", "longitude")]
    if extra:
        raise ValueError(
            f"{var} still has non-spatial dimensions {extra}; select on them "
            "explicitly rather than letting the writer choose a band order")
    logger.info("selected %s -> dims %s shape %s", var, da.dims, da.shape)
    return da

Raising on leftover dimensions rather than squeezing them is the point. A four-dimensional variable silently written as bands produces a file whose band 7 means something nobody recorded.

2. Orient the grid, then attach the CRS 

def orient_north_up(da: xr.DataArray, y: str = "lat", x: str = "lon") -> xr.DataArray:
    """Transpose to (y, x) and flip so y descends, as a geotransform expects."""
    da = da.rename({y: "y", x: "x"}) if y in da.dims else da
    da = da.transpose("y", "x")
    if float(da.y[0]) < float(da.y[-1]):
        logger.info("y ascends; flipping so the origin is north-west")
        da = da.isel(y=slice(None, None, -1))
    return da


def attach_crs(da: xr.DataArray, ds: xr.Dataset) -> xr.DataArray:
    """Read the CF grid_mapping variable and write it as the raster CRS."""
    gm = da.attrs.get("grid_mapping")
    if gm and gm in ds.variables:
        wkt = ds[gm].attrs.get("crs_wkt") or ds[gm].attrs.get("spatial_ref")
        if wkt:
            return da.rio.write_crs(wkt)
        epsg = ds[gm].attrs.get("epsg_code")
        if epsg:
            return da.rio.write_crs(epsg)
    raise ValueError(
        "no usable CRS: the variable has no grid_mapping, or the mapping "
        "variable carries no crs_wkt/epsg_code. Do not default to EPSG:4326 — "
        "a model grid on a rotated pole will be silently wrong.")

The refusal to default to EPSG:4326 is deliberate. Most CF marine products are on a plain geographic grid and the default would usually be right, which is precisely what makes the exceptions — rotated-pole atmospheric grids, polar stereographic ice products — so damaging.

3. Write with the COG driver, and validate 

def write_cog(da: xr.DataArray, path: str, *, nodata: float, blocksize: int = 512) -> None:
    """Write a valid COG: tiled, with overviews, IFDs before the image data."""
    if da.rio.crs is None:
        raise ValueError("refusing to write a COG with no CRS")
    da = da.rio.write_nodata(nodata, encoded=True)
    da.rio.to_raster(
        path, driver="COG", compress="DEFLATE", predictor=2,
        blocksize=blocksize, overview_resampling="average",
        BIGTIFF="IF_SAFER",
    )
    logger.info("wrote %s (%s, %s)", path, da.rio.crs.to_string(), da.shape)

The COG driver rather than GTiff with manual overview building is what guarantees the header-first layout that makes range requests cheap — the property described in setting COG overviews and blocksize for S3. Building overviews after writing with GTiff produces a file that is tiled, has overviews, and is not a COG.

4. Carry the metadata that the format cannot hold

A CF NetCDF variable arrives with units, a standard name, a long name, a cell-methods string and a history. A GeoTIFF has none of those as first-class fields. Everything except the CRS and the nodata value has to be written into free-text tags or lost, and by default it is lost.

The practical minimum is a small, fixed set of tags written on every conversion, so a consumer opening the COG alone can tell what the numbers are:

CARRY = ("units", "standard_name", "long_name", "cell_methods", "source", "references")


def carry_cf_metadata(da: xr.DataArray, ds: xr.Dataset, extra: dict[str, str]) -> dict[str, str]:
    """CF attributes worth writing into the COG's metadata tags."""
    tags = {k: str(da.attrs[k]) for k in CARRY if k in da.attrs}
    for k in ("title", "institution", "Conventions"):
        if k in ds.attrs:
            tags[f"global_{k}"] = str(ds.attrs[k])
    if "time" in da.coords:
        tags["time"] = str(da.coords["time"].values)
    missing = [k for k in ("units", "standard_name") if k not in tags]
    if missing:
        raise ValueError(
            f"refusing to write a COG with no {', '.join(missing)}: a raster of "
            "unnamed numbers in unknown units is not a product")
    return tags | extra

Refusing to write without units and standard_name looks strict and is the cheapest possible defence against the most common downstream question, which is what the values mean. A CF NetCDF always has them; a COG converted from one has them only if somebody carried them across.

The time tag deserves particular attention. Once a single timestep is extracted, nothing in the raster records which one, and a directory of forty COGs distinguished only by filename is exactly the arrangement in which two files get swapped. Writing the timestamp into the file’s own metadata makes the filename a convenience rather than the sole record.

What Survives the Conversion, and What You Have to Carry Two columns listing NetCDF attributes. The survives-by-default column holds the pixel values, the geotransform derived from the coordinate variables, the coordinate reference system once written, and the nodata value once set. The dropped-unless-carried column holds units, standard name, long name, cell methods, the time coordinate of the extracted slice, the institution and title global attributes, and the processing history. A note beneath observes that the dropped column contains everything needed to interpret the numbers, and that a directory of similarly named COGs with no embedded timestamp is the arrangement in which two files get swapped. survives by default pixel values geotransform, from the coordinate vars CRS — once you write it nodata — once you set it enough to place the raster not enough to read it dropped unless carried units · standard_name · long_name cell_methods the time coordinate of this slice institution · title · Conventions processing history everything that says what the numbers are A directory of forty similarly named COGs with no embedded timestamp is how two files get swapped. Six tags, written once in the writer, remove the whole class of question.

Verification and acceptance test

Round-tripping a known cell is the test that catches all four failure modes at once, because a transpose, a flip, a wrong CRS or a fill-as-data all move or change the value at a given coordinate.

import numpy as np
import pytest
import rasterio


def test_cog_round_trips_a_known_cell(tmp_path, ds):
    da = attach_crs(orient_north_up(select_slice(ds, "surf_el", time="2026-03-14")), ds)
    out = tmp_path / "ssh.tif"
    write_cog(da, str(out), nodata=float("nan"))

    lon, lat = -75.0941, 38.7870
    expected = float(da.sel(x=lon, y=lat, method="nearest"))

    with rasterio.open(out) as src:
        assert src.crs is not None, "COG has no CRS"
        assert src.profile["blockxsize"] == 512, "not tiled at the requested block size"
        assert src.overviews(1), "no overview pyramid — this is not a COG"
        got = float(next(src.sample([(lon, lat)]))[0])

    assert np.isclose(got, expected, equal_nan=True), (
        f"value at ({lon}, {lat}) is {got}, expected {expected} — "
        "the raster is transposed, flipped, or in the wrong CRS")

Sampling by coordinate rather than by pixel index is what makes this test meaningful. Comparing array[10, 20] before and after would pass on a flipped raster.

Why the Bounding Box Check Passes on a Flipped Raster Two panels showing the same raster over the same coastline outline. In the correct panel, a plume feature sits in the northern half of the frame and matches the coastline shape. In the flipped panel, the frame occupies exactly the same bounding box and its edges align with the coastline frame just as well, but the plume now sits in the southern half. Annotations note that a bounds comparison, an extent check and a visual overlay at low zoom all pass on the flipped raster, and that only sampling a known coordinate detects it. correct plume bbox matches · sample at (−75.09, 38.79) → 0.42 flipped plume bbox matches · sample at the same point → 0.03 The bounds are identical, the overlay looks plausible at low zoom, and the extent check passes. Only a coordinate sample separates them.

Edge cases and gotchas

  • Packed integers. A variable with scale_factor and add_offset is stored as an integer and must be unpacked before the nodata value means what you think. xr.open_dataset(..., mask_and_scale=True) is the default and does this, but decode_cf=False — used to inspect raw values — turns it off, and code written under that flag carries the raw fill value forward.
  • Longitude on 0–360. A model grid running 0 to 360 produces a COG whose western hemisphere is on the right-hand side. Roll the array and rewrite the coordinate before writing, and verify with a coordinate sample rather than a bbox check.
  • Time as bands is legal and rarely wanted. A COG can hold a time series as bands, but nothing in the format records what each band’s timestamp is beyond a free-text description. For anything more than a handful of steps, the Zarr route preserves the time coordinate properly.

Up: Understanding NetCDF vs GeoTIFF for Marine Data