How to Convert NetCDF to GeoTIFF in Python

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).

Flow for converting a NetCDF variable to GeoTIFF: squeeze extra dimensions, convert longitude, sort latitude north-up, set spatial dimensions and CRS, choose the data type, write and sample back.
Every step fixes something the NetCDF file allowed and a GeoTIFF reader will not expect.

Step-by-step solution

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.

Bar chart of GeoTIFF file size for the same OISST field written as packed int16 and float32, with and without deflate compression.
Compression mattered more than data type; the predictor made no difference for this field.

Code examples

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.

Explanation

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.

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.

Two panels comparing the GeoTIFF bounds and New York sample for the NCEP grid written at 0 to 360 and after converting to โˆ’180 to 180.
The 0โ€“360 file is a perfectly valid raster of the wrong longitudes; only a sample at a known place shows it.

Edge cases or notes

  • 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.

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.