How to Convert NetCDF to GeoTIFF in Python

On this page

Problem statement

GIS software, web maps and most raster tools want a GeoTIFF, and rioxarray will write one from any xarray array with rio.to_raster. The file it writes is only as correct as the array it was given. A NetCDF grid often has longitudes from 0 to 360, latitude stored south to north, packed integer values, extra length-one dimensions and no CRS — and each of those survives into the GeoTIFF, where it produces wrong locations or wrong values rather than an error.

Measured on the NCEP/NCAR Reanalysis 1 monthly air temperature and one day of NOAA OISST v2.1 sea surface temperature, with rioxarray 0.23 and GDAL 3.12:

A 0–360 grid written as-is spans longitudes −1.25 to 358.75. Sampling New York at −74.01° fell outside it and returned 0.0; the converted file returned the correct 2.28 °C.
OISST stores latitude south to north, so its GeoTIFF was south-up: a positive pixel height and row 0 at −89.875°. GDAL's COG driver kept that orientation; gdalwarp flipped it.
The NetCDF packing came along. Written with its encoding intact, the GeoTIFF held 16-bit integers with a 0.01 scale tag, so a reader that ignores the tag sees 2036 instead of 20.36 °C.
The same field was 4.15 MB as uncompressed float32, 1.34 MB as float32 with deflate and 0.89 MB as packed int16 with deflate.

Quick answer

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

sst = xr.open_dataset("oisst-avhrr-v02r01.20240115.nc")["sst"].squeeze(drop=True)
sst = sst.assign_coords(lon=((sst.lon + 180) % 360) - 180).sortby("lon").sortby("lat", ascending=False)
sst = sst.rio.set_spatial_dims(x_dim="lon", y_dim="lat").rio.write_crs("EPSG:4326")
sst.encoding = {}                                             # write physical values, not packed integers
sst.rio.write_nodata(float("nan"), encoded=False).rio.to_raster("sst.tif", compress="deflate")

Five things happen before the write: drop length-one dimensions, convert longitude to −180–180, put north at the top, declare the spatial dimensions and CRS, and decide the data type. Then open the GeoTIFF and sample it at known points (Example 2).

Every step fixes something the NetCDF file allowed and a GeoTIFF reader will not expect.

What you get — from The data used on this page

{'air.mon.mean.nc': 29583703, 'oisst-avhrr-v02r01.20240115.nc': 1545085}

Reduce dimensions, name axes, fix longitude and north-up

1. Reduce the array to two or three dimensions

A GeoTIFF is bands × rows × columns. OISST's sst is (time, zlev, lat, lon) with the first two of length one; rioxarray refused it with TooManyDimensions: Only 2D and 3D data arrays supported. squeeze(drop=True) removes length-one dimensions. A real time dimension becomes bands: select one step for a single-band file, or keep (time, lat, lon) for one band per step.

2. Name the spatial dimensions — again after any reordering

rioxarray finds x and y from coordinate attributes or from conventional names. OISST's lat and lon lacked the attributes, and the write raised MissingSpatialDimensionError: y dimension not found. Call rio.set_spatial_dims(x_dim="lon", y_dim="lat"), or rename the dimensions to y and x. The setting did not survive sortby: sorting after setting it brought the same error back, so sort first and set the dimensions last.

3. Convert longitude to −180 to 180

rioxarray builds the GeoTIFF's transform from the coordinate values, so a 0–360 grid produces a file whose bounds run from −1.25 to 358.75. Nothing in GDAL wraps longitudes: the New York sample at −74.01 landed at row 20, column −30, outside the raster, and returned 0.0. After ((lon + 180) % 360) - 180 and a sort, the transform started at −181.25 and the same sample returned 2.277, identical to xarray.

4. Put north at the top

OISST's latitude runs from −89.875 to 89.875, and a GeoTIFF written in that order has a positive pixel height: row 0 is the southernmost row. GDAL-based readers, including rasterio, sample it correctly — the value at 200.125° E, 30.125° N matched xarray — and gdal_translate -of COG preserved it. Many other tools assume a negative pixel height and draw such a file upside down. sortby("lat", ascending=False) before writing gives the conventional north-up transform, with origin at 90° and pixel height −0.25.

5. Write the CRS

A NetCDF latitude–longitude grid has no CRS of its own. Without rio.write_crs("EPSG:4326") the GeoTIFF had a valid transform and crs None, and GIS software then has to guess. gdal_translate NETCDF:"file.nc":sst out.tif flipped the grid north-up automatically but also wrote no CRS, and kept 0–360 longitudes.

6. Choose the data type and compression

to_raster uses the variable's encoding, so a decoded OISST field is written back as 16-bit integers with a scale tag unless you clear it. Sizes for the same 720 × 1440 field:

option                                   size      dtype    scale tag
encoding kept (packed int16)             2.08 MB   int16    0.01
encoding kept, deflate                   0.89 MB   int16    0.01
float32, no compression                  4.15 MB   float32  1.0
float32, deflate                         1.34 MB   float32  1.0
float32, deflate + predictor 3           1.40 MB   float32  1.0

Packed integers are smaller and exact to 0.01 °C, but only software that applies the scale tag reads them correctly. Float32 with deflate is the safe default for sharing. The floating-point predictor did not help this field.

7. Sample the GeoTIFF back

Open the written file with rasterio, check its transform, bounds, CRS, data type and nodata, and sample a handful of points against the NetCDF values (Example 2). The first New York sample above is the kind of error only a check like this catches: the file opened, displayed and had the right shape.

Compression mattered more than data type; the predictor made no difference for this field.

A conversion function checked against its source

The data used on this page

The NCEP/NCAR Reanalysis 1 monthly mean near-surface air temperature from NOAA PSL (2.5°, January 1948 onwards, updated monthly), and one day of NOAA's OISST v2.1 sea-surface temperature (0.25°, 15 January 2024).

import os
import subprocess

PSL = "Datasets/ncep.reanalysis/Monthlies/surface/air.mon.mean.nc"
FILES = {
    "air.mon.mean.nc": [f"https://downloads.psl.noaa.gov/{PSL}", f"https://psl.noaa.gov/thredds/fileServer/{PSL}"],
    "oisst-avhrr-v02r01.20240115.nc": ["https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/"
                                       "v2.1/access/avhrr/202401/oisst-avhrr-v02r01.20240115.nc"],
}
for name, urls in FILES.items():
    for url in urls * 5:                                 # PSL drops connections: resume and retry
        if os.path.exists(name):
            break
        if subprocess.run(["curl", "-sf", "-C", "-", "--retry", "3", "-o", name + ".part", url]).returncode == 0:
            os.replace(name + ".part", name)
    if not os.path.exists(name):
        raise RuntimeError(f"could not download {name}")
print({name: os.path.getsize(name) for name in FILES})

{'air.mon.mean.nc': 29583703, 'oisst-avhrr-v02r01.20240115.nc': 1545085}

Example 1 — a conversion function that applies every fix

import numpy as np
import rioxarray  # noqa: F401
import xarray as xr


def to_geotiff(da, path, x="lon", y="lat", keep_packing=False, **creation_options):
    """Write a latitude-longitude DataArray as a north-up, -180..180, EPSG:4326 GeoTIFF."""
    extra = [d for d in da.dims if d not in (x, y) and da.sizes[d] == 1]
    da = da.squeeze(extra, drop=True)
    if da.ndim > 3:
        raise ValueError(f"{da.dims}: select or stack down to (band, {y}, {x}) first")
    if float(da[x].max()) > 180:
        da = da.assign_coords({x: ((da[x] + 180) % 360) - 180})
    da = da.sortby(x).sortby(y, ascending=False)
    da = da.rio.set_spatial_dims(x_dim=x, y_dim=y).rio.write_crs("EPSG:4326")
    if not keep_packing:
        da.encoding = {}
        da = da.rio.write_nodata(np.nan, encoded=False)
    da.rio.to_raster(path, **creation_options)
    return path


sst = xr.open_dataset("oisst-avhrr-v02r01.20240115.nc")["sst"]
air = xr.open_dataset("air.mon.mean.nc")["air"].sel(time="2024-01-01")
to_geotiff(sst, "oisst_20240115.tif", compress="deflate")
to_geotiff(air, "ncep_202401.tif", compress="deflate")

Example 2 — check the file against the source

import rasterio


def check_geotiff(path, source, points, x="lon", y="lat"):
    with rasterio.open(path) as src:
        t = src.transform
        print(f"{path}: {src.width} x {src.height}, pixel {t.a} x {t.e}, "
              f"bounds {tuple(round(v, 3) for v in src.bounds)}, crs {src.crs}, dtype {src.dtypes[0]}, nodata {src.nodata}")
        if float(source[x].max()) > 180:
            source = source.assign_coords({x: ((source[x] + 180) % 360) - 180}).sortby(x)
        for name, (lon, lat) in points.items():
            value = float(next(src.sample([(lon, lat)]))[0]) * src.scales[0] + src.offsets[0]
            expected = float(source.sel({x: lon, y: lat}, method="nearest"))
            print(f"  {name:14} GeoTIFF {value:8.3f}   NetCDF {expected:8.3f}")


ocean = {"North Pacific": (-159.875, 30.125), "North Atlantic": (-29.875, 40.125), "Gulf of Guinea": (0.125, 0.125)}
check_geotiff("oisst_20240115.tif", sst.squeeze(drop=True), ocean)
check_geotiff("ncep_202401.tif", air, {"New York": (-74.01, 40.71), "London": (-0.13, 51.51)})

oisst_20240115.tif: 1440 x 720, pixel 0.25 x -0.25, bounds (-180.0, -90.0, 180.0, 90.0), crs EPSG:4326, dtype float32, nodata nan
  North Pacific  GeoTIFF   20.360   NetCDF   20.360
  North Atlantic GeoTIFF   16.970   NetCDF   16.970
  Gulf of Guinea GeoTIFF   28.420   NetCDF   28.420
ncep_202401.tif: 144 x 73, pixel 2.5 x -2.5, bounds (-181.25, -91.25, 178.75, 91.25), crs EPSG:4326, dtype float32, nodata nan
  New York       GeoTIFF    2.277   NetCDF    2.277
  London         GeoTIFF    6.722   NetCDF    6.722

The check multiplies by the file's scale and adds its offset, so it works for packed and unpacked files alike, and it converts a 0–360 source grid to −180–180 before looking up the expected value. Converting only the point instead — lon % 360 — failed at London: −0.13° became 359.87°, and the nearest column on the 0–360 grid was 357.5° rather than 0°, so the check compared 6.722 °C in the GeoTIFF with 6.390 °C from the wrong NetCDF column.

Example 3 — what each option costs in size

import os

variants = {
    "encoding kept (packed int16)": dict(keep_packing=True),
    "encoding kept, deflate": dict(keep_packing=True, compress="deflate"),
    "float32": dict(),
    "float32, deflate": dict(compress="deflate"),
    "float32, deflate + predictor 3": dict(compress="deflate", predictor=3),
}
for i, (label, options) in enumerate(variants.items()):
    path = to_geotiff(sst, f"sst_variant_{i}.tif", **options)
    with rasterio.open(path) as src:
        print(f"{label:32} {os.path.getsize(path) / 1e6:5.2f} MB  {src.dtypes[0]:8} scale {src.scales[0]}")

encoding kept (packed int16)      2.08 MB  int16    scale 0.009999999776482582
encoding kept, deflate            0.89 MB  int16    scale 0.009999999776482582
float32                           4.15 MB  float32  scale 1.0
float32, deflate                  1.34 MB  float32  scale 1.0
float32, deflate + predictor 3    1.40 MB  float32  scale 1.0

The scale is stored as a 32-bit float, so 0.01 reads back as 0.009999999776482582.

Where the transform comes from

Why the transform comes from the coordinates

A GeoTIFF stores an affine transform: the position of the top-left corner and the size of a pixel. A NetCDF file stores coordinate values at cell centres. rioxarray derives the transform from the first coordinate and the spacing, which reproduces whatever the coordinates say — including a first longitude of 0, a first latitude of −89.875, or an axis that is not sorted. It does not interpret longitude as circular.

Why south-up files are legal and still a problem

The TIFF format allows a positive pixel height, and GDAL handles it throughout: reading, sampling, translating to a COG. Tools that assume the conventional negative height — some web viewers, image libraries and hand-written readers — do not, and either draw the file upside down or put values at the mirrored latitude. Writing north-up costs a sort and removes the question.

Why packing travels with the data

xarray keeps each variable's on-disk encoding so that a round trip to NetCDF reproduces the original file. rioxarray uses the same encoding when writing GeoTIFF, which is sensible for a round trip and surprising for a conversion. The scale factor becomes GDAL band metadata that rasterio exposes as scales, and many readers never apply it.

Why GDAL's direct route is not enough on its own

gdal_translate NETCDF:"file.nc":variable reads the NetCDF with GDAL's own driver, flips the rows north-up and keeps the packing tags. It does not convert 0–360 longitudes and wrote no CRS for OISST. It is quick for a look at a file; for a GeoTIFF other people will use, set every property deliberately.

The 0–360 file is a perfectly valid raster of the wrong longitudes; only a sample at a known place shows it.

Many bands, nodata, seams and half-cell shifts

Many time steps become many bands. Band descriptions are not set from the time coordinate automatically; record the dates in metadata or file names.
Cloud-optimised output is driver="COG" in to_raster, or gdal_translate -of COG; see writing a cloud-optimised GeoTIFF.
NaN as nodata works for float types; integer types need a sentinel such as −999.
Curvilinear or rotated grids have no single affine transform; regrid or warp them first.
The seam column at 180° appears once after conversion; a grid with both −180 and 180 columns has a duplicate.
A half-cell shift happens when a transform is built from the first coordinate as if it were a corner; see fixing a NetCDF map that is shifted or flipped.
Large grids can be written in blocks with tiled=True and windowed=True.

Internal links

Fixing a NetCDF map that is shifted, flipped or split at 180° — orientation problems in detail
0–360 or −180–180: longitude conventions in gridded data explained — converting longitude
Fixing NetCDF values that look wrong: scale_factor, _FillValue and units — packed values
How to clip a NetCDF grid to a polygon in Python — preparing a regional export
How to write a cloud-optimised GeoTIFF in Python — COG output
The raster data model explained: bands, dtype, NoData and the transform — what a GeoTIFF stores
Rasterio returns the wrong values: NoData, scaling and dtype fixes — reading scaled GeoTIFFs
NetCDF and gridded data explained: dimensions, variables and attributes — what the source holds

FAQ

How do I convert a NetCDF variable to GeoTIFF in Python?

Open it with xarray, squeeze extra dimensions, convert longitude to −180–180, sort latitude descending, call rio.set_spatial_dims and rio.write_crs("EPSG:4326"), then rio.to_raster. Sample the result against the NetCDF values.

Why is my GeoTIFF from NetCDF upside down?

The NetCDF stored latitude from south to north, so the GeoTIFF has a positive pixel height. GDAL reads it correctly, but other tools may not; sort latitude descending before writing.

Why does my GeoTIFF have longitudes from 0 to 360?

rioxarray writes the coordinates it is given. Convert with ((lon + 180) % 360) - 180 and sort before writing; otherwise points at negative longitudes fall outside the raster.

Why are the values in my GeoTIFF integers?

The NetCDF packing was kept. OISST values were written as int16 with a 0.01 scale tag, so 2036 means 20.36 °C. Clear the variable's encoding before writing to store physical values.

Why does rioxarray say the y dimension was not found?

The latitude and longitude coordinates lack the attributes it looks for. Call rio.set_spatial_dims(x_dim="lon", y_dim="lat") after any sorting, or rename the dimensions to y and x.

Can gdal_translate convert NetCDF to GeoTIFF directly?

Yes, with NETCDF:"file.nc":variable as the source. It wrote a north-up file but no CRS and kept 0–360 longitudes and packed values, so set those explicitly.