CF writer omits coordinate variable attributes if coordinates not named y, x

[gerritholl]

Describe the bug

I'm using the CF writer to write a dataset where the coordinates are named lon, lat rather than y, x. In the file written, the coordinate variables lon and lat are written without any data variable attributes and are not CF compliant.

To Reproduce

import datetime
import pyresample.geometry
import numpy as np
import xarray as xr
from satpy.writers.cf_writer import CFWriter
ar = pyresample.geometry.create_area_def(
                "test", 4326, area_extent=[0, 0, 10, 10], shape=(3, 4))
da = xr.DataArray(
                np.array([[[1, 2, 3, 4],
                           [5, 6, 7, 8],
                           [13, 14, 15, 16]],
                          [[0, 0.3, 0.6, 0.9],
                           [0.3, 0.6, 0.9, 1],
                           [0.9, 0.6, 0.3, 0]]]),
                #dims=("bands", "y", "x"),
                dims=("bands", "lon", "lat"),
                coords={"bands": ["L", "A"],
                        #"y": np.array([8.33333333, 5., 1.66666667]),
                        "lon": np.array([8.33333333, 5., 1.66666667]),
                        #"x": np.array([1.25, 3.75, 6.25, 8.75]),
                        "lat": np.array([1.25, 3.75, 6.25, 8.75]),
                        "crs": ar.crs},
                attrs={#"coordinates": "y x",
                       "coordinates": "lon lat",
                       "name": "test_dataset",
                       "platform_name": "FakeSat",
                       "mode": "LA",
                       "area": ar,
                       "start_time": datetime.datetime(2040, 4, 4, 4, 4, 40),
                       "end_time": datetime.datetime(2040, 4, 4, 5, 5, 50),
                       "standard_name": "image_ready"})

writer = CFWriter(default_config_filename="writers/cf.yaml")
writer.save_datasets([da], "/tmp/test2.nc", format="netCDF4",
                     include_lonlats=False)

Expected behavior

I expect a CF-confirm output file where all variables have variable attributes describing their contents. When coordinates are y and x, as is the default, this works:

netcdf test1 {
dimensions:
        bands = 2 ;
        y = 3 ;
        x = 4 ;
variables:
        int64 test ;
                test:crs_wkt = "GEOGCRS[\"WGS 84\",ENSEMBLE[\"World Geodetic System 1984 ensemble\",MEMBER[\"World Geodetic System 1984 (Transit)\"],MEMBER[\"World Geodetic System 1984 (G730)\"],MEMBER[\"World Geodetic System 1984 (G873)\"],MEMBER[\"World Geodetic System 1984 (G1150)\"],MEMBER[\"World Geodetic System 1984 (G1674)\"],MEMBER[\"World Geodetic System 1984 (G1762)\"],MEMBER[\"World Geodetic System 1984 (G2139)\"],MEMBER[\"World Geodetic System 1984 (G2296)\"],ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ENSEMBLEACCURACY[2.0]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],CS[ellipsoidal,2],AXIS[\"geodetic latitude (Lat)\",north,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433]],AXIS[\"geodetic longitude (Lon)\",east,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433]],USAGE[SCOPE[\"Horizontal component of 3D system.\"],AREA[\"World.\"],BBOX[-90,-180,90,180]],ID[\"EPSG\",4326]]" ;
                test:geographic_crs_name = "WGS 84" ;
                test:grid_mapping_name = "latitude_longitude" ;
                test:horizontal_datum_name = "World Geodetic System 1984 ensemble" ;
                test:inverse_flattening = 298.257223563 ;
                test:long_name = "test" ;
                test:longitude_of_prime_meridian = 0. ;
                test:prime_meridian_name = "Greenwich" ;
                test:reference_ellipsoid_name = "WGS 84" ;
                test:semi_major_axis = 6378137. ;
                test:semi_minor_axis = 6356752.31424518 ;
        string bands(bands) ;
        double y(y) ;
                y:standard_name = "latitude" ;
                y:units = "degrees_north" ;
        double x(x) ;
                x:standard_name = "longitude" ;
                x:units = "degrees_east" ;
        double test_dataset(bands, y, x) ;
                test_dataset:_FillValue = NaN ;
                test_dataset:end_time = "2040-04-04 05:05:50" ;
                test_dataset:grid_mapping = "test" ;
                test_dataset:mode = "LA" ;
                test_dataset:platform_name = "FakeSat" ;
                test_dataset:standard_name = "image_ready" ;
                test_dataset:start_time = "2040-04-04 04:04:40" ;

// global attributes:
                :history = "Created by pytroll/satpy on 2026-01-09 15:19:39.332566+00:00" ;
                :Conventions = "CF-1.7" ;
}

Actual results

As revealed by ncdump -h, the file has no attributes on the variables lon and lat:

$ ncdump -h /tmp/test2.nc
netcdf test2 {
dimensions:
        bands = 2 ;
        lon = 3 ;
        lat = 4 ;
variables:
        int64 test ;
                test:crs_wkt = "GEOGCRS[\"WGS 84\",ENSEMBLE[\"World Geodetic System 1984 ensemble\",MEMBER[\"World Geodetic System 1984 (Transit)\"],MEMBER[\"World Geodetic System 1984 (G730)\"],MEMBER[\"World Geodetic System 1984 (G873)\"],MEMBER[\"World Geodetic System 1984 (G1150)\"],MEMBER[\"World Geodetic System 1984 (G1674)\"],MEMBER[\"World Geodetic System 1984 (G1762)\"],MEMBER[\"World Geodetic System 1984 (G2139)\"],MEMBER[\"World Geodetic System 1984 (G2296)\"],ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ENSEMBLEACCURACY[2.0]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],CS[ellipsoidal,2],AXIS[\"geodetic latitude (Lat)\",north,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433]],AXIS[\"geodetic longitude (Lon)\",east,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433]],USAGE[SCOPE[\"Horizontal component of 3D system.\"],AREA[\"World.\"],BBOX[-90,-180,90,180]],ID[\"EPSG\",4326]]" ;
                test:geographic_crs_name = "WGS 84" ;
                test:grid_mapping_name = "latitude_longitude" ;
                test:horizontal_datum_name = "World Geodetic System 1984 ensemble" ;
                test:inverse_flattening = 298.257223563 ;
                test:long_name = "test" ;
                test:longitude_of_prime_meridian = 0. ;
                test:prime_meridian_name = "Greenwich" ;
                test:reference_ellipsoid_name = "WGS 84" ;
                test:semi_major_axis = 6378137. ;
                test:semi_minor_axis = 6356752.31424518 ;
        string bands(bands) ;
        double lon(lon) ;
        double lat(lat) ;
        double test_dataset(bands, lon, lat) ;
                test_dataset:_FillValue = NaN ;
                test_dataset:end_time = "2040-04-04 05:05:50" ;
                test_dataset:grid_mapping = "test" ;
                test_dataset:mode = "LA" ;
                test_dataset:platform_name = "FakeSat" ;
                test_dataset:standard_name = "image_ready" ;
                test_dataset:start_time = "2040-04-04 04:04:40" ;

// global attributes:
                :history = "Created by pytroll/satpy on 2026-01-09 15:00:09.854111+00:00" ;
                :Conventions = "CF-1.7" ;
}

Environment Info:

• OS: openSUSE Leap 15.6
• Satpy Version: main

Additional context

My aim is to create a file with coordinates "lon" and "lat", in an unprojected CRS (geographic CRS with lon/lat coordinates). Maybe there is a different way to do this.

The reason for this behaviour is in satpy.cf.coords.add_xy_coords_attrs, where y and x are hardcoded:

satpy/satpy/cf/coords.py

Lines 21 to 25 in 579e797

	def add_xy_coords_attrs(data_arr: xr.DataArray) -> xr.DataArray:
	"""Add relevant attributes to x, y coordinates."""
	# If there are no coords, return dataarray
	if not data_arr.coords.keys() & {"x", "y", "crs"}:
	return data_arr

The lower-level function _add_xy_geographic_coords_attrs does allow passing different names, but this doesn't seem to be reachable by public API, as far as I can tell:

satpy/satpy/cf/coords.py

Lines 93 to 101 in 579e797

	def _add_xy_geographic_coords_attrs(data_arr: xr.DataArray, x: str = "x", y: str = "y") -> xr.DataArray:
	"""Add relevant attributes to x, y coordinates of a geographic CRS."""
	if x in data_arr.coords:
	data_arr[x].attrs["standard_name"] = "longitude"
	data_arr[x].attrs["units"] = "degrees_east"
	if y in data_arr.coords:
	data_arr[y].attrs["standard_name"] = "latitude"
	data_arr[y].attrs["units"] = "degrees_north"
	return data_arr

[djhoese]

Where did your custom lon/lat coordinates in the DataArray come from? What happens if you remove those and do include_lonlats=True (the default)? How does that change the final output?

[gerritholl]

The lon/lat coordinates come from renaming the y/x coordinates:

da = da.rename({"y": "lon", "x": "lat"})

and

da.attrs["coordinates"] = "lon lat".

include_lonlats=True does something else. It does not add 1-D coordinates, but two 2-D data variables with lon/lat for every pixel. This is redundant information in case of unprojected coordinates (lat/lon grid).

I'm working on a custom in-house writer for client software that expects lon/lat coordinates.

What I need:

        double lon(lon) ;
                lon:units = "degrees_east" ;
                lon:long_name = "longitude" ;
                lon:standard_name = "longitude" ;
        double lat(lat) ;
                lat:units = "degrees_north" ;
                lat:long_name = "latitude" ;
                lat:standard_name = "latitude" ;

What include_lonlats=True does:

        double longitude(y, x) ;
                longitude:_FillValue = NaN ;
                longitude:name = "longitude" ;
                longitude:standard_name = "longitude" ;
                longitude:units = "degrees_east" ;
        double latitude(y, x) ;
                latitude:_FillValue = NaN ;
                latitude:name = "latitude" ;
                latitude:standard_name = "latitude" ;
                latitude:units = "degrees_north" ;

What I need is normally provided by:

        double y(y) ;
                y:standard_name = "latitude" ;
                y:units = "degrees_north" ;
        double x(x) ;
                x:standard_name = "longitude" ;
                x:units = "degrees_east" ;

But I need y to be called lat and x to be called lon.

This works only for a geographic CRS, but this is the only CRS the client understands, so this restriction is wanted.

[djhoese]

Hm interesting. Yeah some assumptions were probably made for simplicity and assumed the Satpy preference for "y" and "x" dimension naming. I'm not sure you'll get around that without either:

1. Modifying the writer.
2. Renaming the x/y variables and dimensions after creation, outside of Satpy.

[gerritholl]

Modifying the writer is essentially what I'm doing by writing my own writer, subclassing CFWriter. I need to do this in any case, even if CFWriter would grow the flexibility to write the y/x coordinates under user-supplied names, due to some specific oddities needed by the client application. It works when I redo the work in _add_xy_geographic_coords_attrs in my custom writer. It seems some functions have thought of the use case "what if your coordinates are not y/x", but not others. The CF writer is one of the few places, perhaps the only one, where Satpy's y/x preference becomes visible outside of Satpy, so maybe this is more of a feature request than a bug report.