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;
gdalwarpflipped 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).
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.
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.
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.
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"into_raster, orgdal_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=Trueandwindowed=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.