Understanding NetCDF vs GeoTIFF for Marine Data

Format selection in automated coastal and marine spatial analysis pipelines sits within the broader Marine Spatial Data Fundamentals & Architecture domain and is not a stylistic preference — it is a deterministic routing decision that dictates memory topology, I/O throughput, CRS validation pathways, and long-term archival strategy. This page establishes the operational criteria for selecting between NetCDF and GeoTIFF, provides production-grade ingestion and transformation routines, and defines cloud-native pipeline boundaries for multi-dimensional oceanographic and coastal raster datasets. All workflows are engineered for reproducible execution in containerized, memory-constrained environments.

Reference Configuration and Specification Table #

The table below fixes the library versions, chunk dimensions, compression parameters, and accuracy thresholds used throughout this page. Pin these in your requirements.txt or environment.yml before running any ingestion routine.

---

Format Architecture and Routing Decision #

NetCDF (Network Common Data Form) and GeoTIFF solve fundamentally different spatial data problems and are not interchangeable. Understanding their binary architecture is the prerequisite for writing correct routing logic.

NetCDF is a self-describing, N-dimensional container typically backed by HDF5. Its design is optimised for temporal sequences, atmospheric and oceanographic model outputs, and CF-Compliant metadata hierarchies. ROMS and HYCOM model outputs, satellite-derived sea surface temperature (SST) stacks, and chlorophyll-a time-series almost always arrive in this format. The file carries its own dimension labels, coordinate variables, units, and calendar conventions internally.

GeoTIFF is a 2D/3D raster format optimised for spatial indexing, GDAL interoperability, and cloud-optimized tiling. Bathymetric grids exported from multibeam point-cloud filtering suites, habitat suitability rasters, and machine-learning training tiles arrive as GeoTIFF. CRS metadata is embedded in the file header via the GeoKeyDirectoryTag, and tiling is controlled by block layout parameters set at write time.

The diagram below maps the routing decision from raw dataset characteristics through to the two cloud-native output formats.

Memory-Access Topology #

The routing decision is ultimately a decision about how bytes reach RAM. NetCDF-4/HDF5 stores variables as N-dimensional chunks; a read of one time-slice pulls every chunk that the slice intersects, so a request that straddles chunk boundaries amplifies I/O. A Cloud-Optimized GeoTIFF stores 2D tiles plus pre-computed overviews, so a viewport read fetches only the intersecting tiles at the coarsest sufficient overview level. The diagram below contrasts the two access patterns — the source of both the read-amplification and overview-ordering failure modes described later.

Operational Routing Matrix #

Apply this matrix to any incoming dataset before writing a single line of ingestion code. Route by the first criterion that matches unambiguously; move down the table only if the upper criteria are inconclusive.

---

Memory-Constrained Python Implementation #

The routines below enforce memory-safe ingestion, explicit chunk alignment, and cloud-native output generation. They are designed for containerised execution where OOM termination must be prevented. All functions use strict type annotations, the logging module (never print), and raise exceptions on invalid state rather than silently returning None.

"""
netcdf_geotiff_routing.py
Production-grade ingestion and format-conversion routines for marine spatial pipelines.
Tested against xarray 2024.3, rasterio 1.3.10, rioxarray 0.15.5, GDAL 3.8.
"""
import logging
from pathlib import Path

import numpy as np
import rioxarray  # noqa: F401 — registers the .rio accessor on xarray objects
import xarray as xr
from rasterio.enums import Resampling

logger = logging.getLogger(__name__)


# ---------------------------------------------------------------------------
# NetCDF ingestion
# ---------------------------------------------------------------------------

def ingest_oceanographic_netcdf(filepath: str | Path) -> xr.Dataset:
    """
    Open a multi-dimensional oceanographic NetCDF with Dask-backed lazy evaluation.

    Chunk sizes are set to time=1 so a single time-step occupies one Dask task;
    lat/lon chunks of 256 align well with typical HDF5 block layouts produced by
    ROMS, HYCOM, and CMEMS distribution files.  Adjust if `ncdump -sh` reveals
    a different native chunk layout — misaligned chunks cause read amplification.

    Parameters
    ----------
    filepath : path to a NetCDF-4 (or NetCDF-3) file.

    Returns
    -------
    xr.Dataset with Dask-backed arrays and decoded CF conventions.

    Raises
    ------
    ValueError  if required spatial dimensions are absent.
    FileNotFoundError  if filepath does not exist.
    """
    filepath = Path(filepath)
    if not filepath.exists():
        raise FileNotFoundError(f"NetCDF file not found: {filepath}")

    logger.info("Opening NetCDF: %s", filepath)
    ds = xr.open_dataset(
        filepath,
        chunks={"time": 1, "lat": 256, "lon": 256},
        engine="netcdf4",
        decode_cf=True,       # parse units, calendar, and _FillValue
        mask_and_scale=True,  # apply scale_factor / add_offset and mask _FillValue
    )

    required_dims = {"lat", "lon"}
    missing = required_dims - set(ds.dims)
    if missing:
        # Some ROMS/HYCOM files use 'latitude'/'longitude' or 'y'/'x' — detect and
        # remap rather than failing silently with wrong spatial extent.
        raise ValueError(
            f"Dataset at {filepath} is missing required spatial dimensions: {missing}. "
            "Check dimension names with ds.dims and remap before pipeline ingestion."
        )

    logger.info(
        "Loaded dataset — variables: %s, dims: %s",
        list(ds.data_vars),
        dict(ds.dims),
    )
    return ds


def validate_cf_compliance(ds: xr.Dataset) -> None:
    """
    Assert minimum CF convention compliance for pipeline safety.

    Checks for grid_mapping attribute, units on coordinate variables, and
    presence of _FillValue or missing_value on all data variables.

    Raises
    ------
    ValueError  on any CF violation that would corrupt downstream spatial ops.
    """
    for var in ds.data_vars:
        da = ds[var]
        if "_FillValue" not in da.attrs and "missing_value" not in da.attrs:
            logger.warning(
                "Variable '%s' has no _FillValue — nodata masking will be unreliable.", var
            )
        if "units" not in da.attrs:
            logger.warning("Variable '%s' carries no 'units' attribute.", var)

    coord_vars = [c for c in ds.coords if c in ("lat", "lon", "latitude", "longitude")]
    if not coord_vars:
        raise ValueError(
            "No recognised geographic coordinate variables found. "
            "Ensure 'lat'/'lon' or 'latitude'/'longitude' exist as coordinate variables."
        )
    logger.info("CF compliance check passed for variables: %s", list(ds.data_vars))


# ---------------------------------------------------------------------------
# CRS validation
# ---------------------------------------------------------------------------

def validate_crs_parity(ds1: xr.Dataset, ds2: xr.Dataset) -> None:
    """
    Enforce strict CRS equivalence before any spatial join or concatenation.

    Silent CRS mismatches are among the most common sources of corrupt spatial
    joins in cross-format pipelines — particularly when fusing NetCDF model
    outputs (typically EPSG:4326) with projected GeoTIFFs (UTM or Web Mercator).

    Raises
    ------
    ValueError  if CRS objects differ.
    RuntimeError  if either dataset has no CRS attached.
    """
    crs1 = ds1.rio.crs
    crs2 = ds2.rio.crs

    if crs1 is None or crs2 is None:
        raise RuntimeError(
            "One or both datasets have no attached CRS. "
            "Assign CRS with ds.rio.write_crs('EPSG:4326', inplace=True) before comparison."
        )
    if crs1 != crs2:
        raise ValueError(
            f"CRS mismatch detected: {crs1!r} vs {crs2!r}. "
            "Reproject the lower-priority dataset with ds.rio.reproject() before concatenation."
        )
    logger.info("CRS parity confirmed: %s", crs1)


# ---------------------------------------------------------------------------
# GeoTIFF → Cloud-Optimized GeoTIFF conversion
# ---------------------------------------------------------------------------

def convert_to_cog(
    input_path: str | Path,
    output_path: str | Path,
    block_size: int = 512,
) -> Path:
    """
    Convert a raw GeoTIFF to Cloud-Optimized GeoTIFF with explicit tiling and
    ZSTD compression.

    Uses rioxarray for memory-windowed reads so the full raster never lands in
    RAM simultaneously.  Block size of 512 × 512 balances HTTP range-request
    efficiency (< 1 MB per tile at float32) with filesystem read performance.
    Use 256 for very high-resolution coastal lidar tiles to reduce over-fetch.

    The 'overviews="AUTO"' argument instructs rioxarray/GDAL to build embedded
    pyramid overviews — mandatory for COG compliance and web tile rendering.

    Parameters
    ----------
    input_path  : source GeoTIFF (any GDAL-readable raster).
    output_path : destination path; overwrites if present.
    block_size  : tile dimensions in pixels; must be a power of two.

    Returns
    -------
    Path to the written COG file.

    Raises
    ------
    FileNotFoundError  if input_path does not exist.
    ValueError  if block_size is not a positive power of two.
    """
    input_path = Path(input_path)
    output_path = Path(output_path)

    if not input_path.exists():
        raise FileNotFoundError(f"Input raster not found: {input_path}")
    if block_size <= 0 or (block_size & (block_size - 1)) != 0:
        raise ValueError(f"block_size must be a positive power of two; got {block_size}.")

    logger.info("Converting %s → COG at %s (block=%d)", input_path, output_path, block_size)

    rds = rioxarray.open_rasterio(
        input_path,
        chunks={"band": 1, "y": block_size, "x": block_size},
        lock=False,  # thread-safe for Dask scheduler
    )

    rds.rio.to_raster(
        output_path,
        tiled=True,
        blockxsize=block_size,
        blockysize=block_size,
        compress="zstd",
        driver="GTiff",
        overviews="AUTO",
        resampling=Resampling.average,
    )

    logger.info("COG written: %s (%.1f MB)", output_path, output_path.stat().st_size / 1e6)
    return output_path


# ---------------------------------------------------------------------------
# NetCDF → Zarr conversion for distributed compute
# ---------------------------------------------------------------------------

def export_to_zarr(
    ds: xr.Dataset,
    store_path: str | Path,
    consolidated: bool = True,
) -> None:
    """
    Write a Dask-backed xarray Dataset to a Zarr store for object-store deployment.

    Consolidated metadata writes a single .zmetadata file that eliminates the
    N-file metadata scan on first open — critical for S3/GCS latency.
    Call this after all variable subsetting and reprojection is complete;
    writing partial datasets and appending later risks chunk boundary corruption.

    Parameters
    ----------
    ds             : fully validated, CRS-assigned xr.Dataset.
    store_path     : local path or fsspec-compatible URI (e.g. s3://bucket/store).
    consolidated   : write consolidated metadata (default True).

    Raises
    ------
    ValueError  if ds contains no data variables.
    """
    if not ds.data_vars:
        raise ValueError("Dataset has no data variables; refusing to write empty Zarr store.")

    store_path = str(store_path)
    logger.info("Writing Zarr store: %s (consolidated=%s)", store_path, consolidated)

    ds.to_zarr(
        store_path,
        mode="w",
        consolidated=consolidated,
    )
    logger.info("Zarr store complete: %s", store_path)

---

Validation Gates and Quality Control #

Pipeline stage separation is only meaningful when each stage emits a validated artefact. Apply all four gates before handing data to a downstream consumer. These gates catch the same defect classes that surface during DEM interpolation for seafloor mapping and bathymetric artifact and noise removal, so failing fast here prevents corruption from propagating into terrain models.

Gate 1 — Coordinate Bounds Check #

Reject datasets whose bounding box falls outside the expected study area. A single corrupt lon axis encoding degrees-east as 0–360 rather than −180–180 will silently shift every spatial join by 360°.

def assert_geographic_bounds(
    ds: xr.Dataset,
    lat_range: tuple[float, float] = (-90.0, 90.0),
    lon_range: tuple[float, float] = (-180.0, 180.0),
) -> None:
    """Raise if coordinate bounds fall outside expected geographic range."""
    lat = ds["lat"].values
    lon = ds["lon"].values

    if lat.min() < lat_range[0] or lat.max() > lat_range[1]:
        raise ValueError(f"lat bounds {lat.min():.3f}–{lat.max():.3f} outside {lat_range}.")
    if lon.min() < lon_range[0] or lon.max() > lon_range[1]:
        raise ValueError(
            f"lon bounds {lon.min():.3f}–{lon.max():.3f} outside {lon_range}. "
            "Dataset may use 0–360 longitude convention — apply lon = lon % 360 - 360."
        )
    logger.info("Coordinate bounds: lat [%.3f, %.3f], lon [%.3f, %.3f]",
                lat.min(), lat.max(), lon.min(), lon.max())

Gate 2 — Nodata Coverage Check #

A variable where more than 90 % of values are masked typically indicates a unit mismatch, a wrong _FillValue sentinel, or a corrupt scale factor.

def assert_data_coverage(ds: xr.Dataset, min_valid_fraction: float = 0.10) -> None:
    """Raise if any data variable has fewer than min_valid_fraction valid values."""
    for var in ds.data_vars:
        da = ds[var]
        valid_fraction = float(da.notnull().mean())
        if valid_fraction < min_valid_fraction:
            raise ValueError(
                f"Variable '{var}' has only {valid_fraction:.1%} valid values "
                f"(threshold: {min_valid_fraction:.0%}). "
                "Check _FillValue, scale_factor, and add_offset in the source file."
            )
    logger.info("Data coverage check passed for all variables.")

Gate 3 — Spatial Resolution Consistency #

When fusing two rasters, verify pixel spacing matches within tolerance before any arithmetic. A 0.5-pixel offset, invisible at print scale, introduces systematic bias in ML feature extraction and in tidal datum transformation interpolation.

def assert_resolution_parity(
    da1: xr.DataArray,
    da2: xr.DataArray,
    tolerance: float = 1e-5,
) -> None:
    """Raise if spatial resolution differs between two DataArrays beyond tolerance."""
    res1 = da1.rio.resolution()
    res2 = da2.rio.resolution()
    if abs(res1[0] - res2[0]) > tolerance or abs(res1[1] - res2[1]) > tolerance:
        raise ValueError(
            f"Resolution mismatch: {res1} vs {res2}. "
            "Resample to common grid with da.rio.reproject_match() before combining."
        )
    logger.info("Resolution parity confirmed: %s", res1)

Gate 4 — CF Grid Mapping Verification #

For NetCDF files used in oceanographic modelling, the absence of a grid_mapping attribute on data variables causes silent datum assumptions downstream.

def assert_grid_mapping(ds: xr.Dataset) -> None:
    """
    Verify that at least one data variable carries a grid_mapping attribute
    referencing a coordinate variable that contains a CRS WKT or proj string.
    """
    for var in ds.data_vars:
        gm = ds[var].attrs.get("grid_mapping")
        if gm and gm in ds.coords:
            logger.info("Variable '%s' has valid grid_mapping: '%s'.", var, gm)
            return
    logger.warning(
        "No grid_mapping attribute found on any data variable. "
        "CRS will be inferred from coordinate variable names — "
        "assign explicitly with ds.rio.write_crs() to prevent datum errors."
    )

---

Common Failure Modes and Diagnosis #

1. Silent Longitude Convention Mismatch (0–360 vs −180–180) #

Symptom: Spatial joins return empty GeoDataFrames; NetCDF grid appears to the east of the study area in QGIS.

Root cause: HYCOM and some CMEMS products distribute longitude on the 0–360 range. GDAL and rioxarray default to −180–180. The arithmetic overlap is zero, so every spatial join silently returns no rows rather than raising an error.

Remediation:

if ds["lon"].max() > 180:
    ds = ds.assign_coords(lon=(ds["lon"] + 180) % 360 - 180)
    ds = ds.sortby("lon")
    logger.info("Normalised longitude to −180–180 convention.")

2. Chunk Boundary Misalignment Causing Read Amplification #

Symptom: Dask task graph shows 10× more data read than the array size warrants; memory spikes during .compute() calls.

Root cause: Specifying chunks={"lat": 256, "lon": 256} on a file whose native HDF5 chunk layout is (1, 90, 180) forces every partial chunk read to fetch the full native block. For a global 0.083° CMEMS SST product this amplifies reads by roughly 8×.

Diagnosis: Run ncdump -sh <file>.nc | grep "_ChunkSizes" to inspect native chunk dimensions, then set chunks to match or to integer multiples of the native size.

3. COG Overviews Written After Data — File Is Not COG-Compliant #

Symptom: rio cogeo validate reports Overview is not valid. Spatial index should be after it; tile servers cannot use range requests efficiently.

Root cause: Writing overviews after the main image data with the GTiff driver creates an invalid COG. The GDAL COG driver (available since GDAL 3.1) writes overviews before data in the correct interleave.

Remediation: Use driver="COG" when available:

rds.rio.to_raster(output_path, driver="COG", compress="zstd",
                  blockxsize=512, blockysize=512,
                  resampling=Resampling.average)

Verify with python -c "from rio_cogeo.cogeo import cog_validate; print(cog_validate('output.tif'))".

4. rioxarray Writes Wrong Nodata Sentinel to COG #

Symptom: ArcGIS and QGIS render valid ocean pixels as transparent; nodata extent appears larger than expected.

Root cause: rioxarray writes the rio.nodata value from the opened raster. If the source NetCDF used np.nan (float32 NaN) but the COG output is int16, the sentinel is silently truncated to 0, masking real zero-valued bathymetry.

Remediation: Explicitly set the nodata value before writing:

rds = rds.rio.write_nodata(-9999, encoded=True)
rds.rio.to_raster(output_path, driver="COG", compress="zstd")

---

Pipeline Integration and Downstream Handoff #

Format routing determines which downstream stages can consume a dataset and at what latency. The handoff contract between this stage and its consumers must be documented in a metadata manifest written alongside every output artefact.

NetCDF → Zarr path feeds distributed compute workers that need to iterate over time steps without loading full grids into memory. When integrating these outputs with the AIS NMEA sentence parser, maintain strict temporal alignment: the Zarr store’s time axis must carry ISO-8601 timestamps so vessel position lookups by time index are unambiguous. Zarr stores on S3 or GCS should set consolidated=True and publish the .zmetadata URL in the manifest so downstream consumers skip the N-file scan.

GeoTIFF → COG path feeds GIS visualisation tools, bathymetric DEM interpolation workflows, and ML tile extraction. COGs must carry embedded overviews at factors 2, 4, 8, 16, and 32. Web tile servers (GDAL VSI, TiTiler, GeoServer ImageMosaic) will use overviews automatically for lower zoom levels; without them, every request triggers a full-resolution read and degrades under load.

Metadata Manifest Schema #

Write a manifest.json alongside every output artefact:

{
  "artefact": "sst_2024_Q1_cog.tif",
  "format": "COG",
  "crs": "EPSG:4326",
  "resolution_deg": 0.083,
  "nodata": -9999,
  "source_file": "sst_2024_Q1.nc",
  "source_variable": "sst",
  "time_coverage_start": "2024-01-01T00:00:00Z",
  "time_coverage_end": "2024-03-31T23:59:59Z",
  "pipeline_version": "1.4.2",
  "validation_gates_passed": ["bounds", "coverage", "resolution", "cog_validate"]
}

Long-term archival of raw NetCDF files should use immutable object-store versioning. Do not binary-commit multi-gigabyte spatial assets to Git — use DVC or LakeFS for pointer-based version control. GeoTIFFs must be finalised as COGs before archival; storing raw unoptimised TIFFs wastes object-store egress bandwidth on every downstream read.

---

Related #

• CRS Alignment for Coastal GIS Projects — CRS validation and reprojection workflows that gate every cross-format pipeline join
• Tidal Datum Transformations in Python — vertical datum handling for coastal rasters exported from NetCDF model outputs
• Parsing AIS NMEA Sentences with Python — temporal alignment of vessel telemetry with NetCDF oceanographic grids
• DEM Interpolation Techniques for Seafloor Mapping — downstream consumer of COG-format bathymetric rasters

Up: Marine Spatial Data Fundamentals & Architecture