xarray NetCDF with groups to GeoTIFF

Question

I am trying to convert NetCDF files, available from this link, to raster or geotiff files using Python 3. After reading the a NetCDF file using the package netcdf4, it appears to have groups: CO2, CO2_uncertainty, and O2.

import netCDF4 as nc
xds = nc.Dataset('GCP_Global_1959.nc')
xds

I then tried to open the file with xarray and had the following output where I can see the coordinates and time slices. Although, the coordinates seem to be out of bounds.

import xarray as xr
xds = xr.open_dataset('GCP_Global_1959.nc')

enter image description here

Finally, I tried to add a group to the open_dataset call but coordinates are missing.

xds = xr.open_dataset('GCP_Global_1959.nc', group='CO2')

enter image description here

I've tried the solution from Using GDAL, NetCDF4 to GeoTIFF wrongly rotates the output GeoTIFF file 90 degrees but it doesn't seem to work.

How can I convert this NetCDF file to a raster?

Answer 1

You were close already! :)

import rioxarray
import rasterio
import xarray as xr
# Open the coordinates as a xr.Dataset.
ds_coords = xr.open_dataset('GCP_Global_1959.nc')
# Open the data inside the group as a xr.Dataset.
ds = xr.open_dataset('GCP_Global_1959.nc', group='CO2')
# Add the coordinates to the "data"-dataset.
ds = ds.assign_coords(ds_coords.coords)
# Assign a CRS.
ds = ds.rio.write_crs(4326)
# Rename some dimensions.
ds = ds.rename({"lat": "y", "lon": "x"})
# Create 1 geoTIFF per variable.
for var in ds.data_vars:
 ds[var].rio.to_raster(f"{var}.tif")

Check if it works:

test_ds = xr.open_dataset("OIL.tif", engine = "rasterio")
print(test_ds)
<xarray.Dataset>
Dimensions: (band: 12, x: 3600, y: 1800)
Coordinates:
 * band (band) int64 1 2 3 4 5 6 7 8 9 10 11 12
 * x (x) float64 -179.9 -179.8 -179.7 -179.6 ... 179.7 179.8 179.9
 * y (y) float64 -89.95 -89.85 -89.75 -89.65 ... 89.75 89.85 89.95
 spatial_ref int64 ...
Data variables:
 band_data (band, y, x) float32 ...