CF Coordinate Subsampling (section 8.3) and affine transformations

[ethanrd]

Topic for discussion

Hi all,

I'm looking at some of the concerns the GeoZarr community has with CRS handling in CF (see GeoZarr issue #20) and trying to understand if the affine transformation mentioned there can be handled by CF 8.3 "Lossy Compression by Coordinate Subsampling".

I believe section 8.3 could support a full affine transformation with the addition of a new interpolation method. I have been working on what that might look like and will post more on that here once I get a bit further.

In the meantime, I wanted share my attempt to represent an affine transformation that is limited to translation and scale (excluding shear and rotation) using the "linear" interpolation method. CF section 8.3 is pretty new to me so I wanted to ask those more familiar if this attempt looks good.

Here's a description of my limited-affine test dataset and the CDL I put together which I think is inline with CF section 8.3. It is an elevation dataset with equally spaced lat/lon:

lat(121) = [ 35.0, ..., 45.0]
lon(121) = [-110.0, ..., -100.0]

So the resolution in both lat and lon is 0.083333 = ( lat(120) - lat(0) ) / (121 - 1)) and I think the GDAL/GeoTIFF GeoTransform (affine transformation notation) would be

GT(6) = [ 35.0, 0.08333333, 0.0, -110.0, 0.0,  0.08333333 ]

And the equations would be:

lat(i, j) = GT(0) + GT(1)*i + GT(2)*j = 35.0 + (0.08333333 * i) + (0.0 * j)
lon(i, j) = GT(3) + GT(4)*i + GT(5)*j = -110.0 + (0.0 * i) + (0.08333333 * j)

which collapses to 1D lat and lon as above.

Here's the CDL snippet I came up with:

dimensions:
  lat = 121 ;
  lon = 121 ;
  lat_tp = 2 ;
  lon_tp = 2 ;

variables:
  float elev(lat, lon) ;
    elev:long_name = "elevation" ;
    elev:standard_name = "height" ;
    elev:units = "meters" ;
    elev:positive = "up" ;
    elev:coordinate_interpolation = "lat: lat_interpolation lon: lon_interpolation" ;

  int lat_interpolation ;
    lat_interpolation:interpolation_name = "linear" ;
    lat_interpolation:tie_point_mapping = "lat: lat_indices lat_tp" ;
    lat_interpolation:computational_precision = "32" ;

  int lon_interpolation ;
    lon_interpolation:interpolation_name = "linear" ;
    lon_interpolation:tie_point_mapping = "lon: lon_indices lon_tp" ;
    lon_interpolation:computational_precision = "32" ;


  float lat_tp(lat_tp) ;
  float lon_tp(lon_tp) ;

  int lat_indices(lat_tp) ;
  int lon_indices(lon_tp) ;

data:
  lat_indices = 0, 120 ;
  lon_indices = 0, 120 ;
  lat_tp = 35.0, 45.0 ;
  lon_tp = -110.0, -100.0 ;

Does that look right to those familiar with CF section 8.3?

[erget]

This sounds promising - I've reached out to Anders, who was the lead author for 8.3 - hope he's able to give some insights here! @mraspaud do you have an opinion?

[mraspaud]

I’m not familiar with affine transformations myself, sorry. However, when it comes to the subsampling, my understanding from the original algorithm is that the reconstruction/interpolation should really be done in 3d, so projected coordinates wouldn’t really make sense.
In CF though, it looks like there in a choice for which interpolation scheme to use is free(er), so I guess just applying the affine transformation on the subsampled coordinates makes sense.

[ethanrd]

Interesting. Thanks Martin @mraspaud

[davidhassell]

Hi Ethan,

The CDL is nearly right :) The name: elements of coordinate_interpolation attribute should be tie point coordinate variables, rather than (solely) dimensions. I tried a modified version of the CDL (see below) in cf-python, and it worked fine:

>>> import cf
>>> f = cf.read('example.cdl')[0]
>>> print(f)
Field: height (ncvar%elev)
--------------------------
Data            : height(latitude(121), longitude(121)) meters
Dimension coords: latitude(121) = [35.0, ..., 45.0] degrees_north
                : longitude(121) = [-110.0, ..., -100.0] degrees_east
>>> f.coordinate('Y').array
array([35.        , 35.08333333, 35.16666667, 35.25      , 35.33333333,
       35.41666667, 35.5       , 35.58333333, 35.66666667, 35.75      ,
       35.83333333, 35.91666667, 36.        , 36.08333333, 36.16666667,
       36.25      , 36.33333333, 36.41666667, 36.5       , 36.58333333,
       36.66666667, 36.75      , 36.83333333, 36.91666667, 37.        ,
       37.08333333, 37.16666667, 37.25      , 37.33333333, 37.41666667,
       37.5       , 37.58333333, 37.66666667, 37.75      , 37.83333333,
       37.91666667, 38.        , 38.08333333, 38.16666667, 38.25      ,
       38.33333333, 38.41666667, 38.5       , 38.58333333, 38.66666667,
       38.75      , 38.83333333, 38.91666667, 39.        , 39.08333333,
       39.16666667, 39.25      , 39.33333333, 39.41666667, 39.5       ,
       39.58333333, 39.66666667, 39.75      , 39.83333333, 39.91666667,
       40.        , 40.08333333, 40.16666667, 40.25      , 40.33333333,
       40.41666667, 40.5       , 40.58333333, 40.66666667, 40.75      ,
       40.83333333, 40.91666667, 41.        , 41.08333333, 41.16666667,
       41.25      , 41.33333333, 41.41666667, 41.5       , 41.58333333,
       41.66666667, 41.75      , 41.83333333, 41.91666667, 42.        ,
       42.08333333, 42.16666667, 42.25      , 42.33333333, 42.41666667,
       42.5       , 42.58333333, 42.66666667, 42.75      , 42.83333333,
       42.91666667, 43.        , 43.08333333, 43.16666667, 43.25      ,
       43.33333333, 43.41666667, 43.5       , 43.58333333, 43.66666667,
       43.75      , 43.83333333, 43.91666667, 44.        , 44.08333333,
       44.16666667, 44.25      , 44.33333333, 44.41666667, 44.5       ,
       44.58333333, 44.66666667, 44.75      , 44.83333333, 44.91666667,
       45.        ])
>>> f.coordinate('X').array 
array([-110.        , -109.91666667, -109.83333333, -109.75      ,
       -109.66666667, -109.58333333, -109.5       , -109.41666667,
       -109.33333333, -109.25      , -109.16666667, -109.08333333,
       -109.        , -108.91666667, -108.83333333, -108.75      ,
       -108.66666667, -108.58333333, -108.5       , -108.41666667,
       -108.33333333, -108.25      , -108.16666667, -108.08333333,
       -108.        , -107.91666667, -107.83333333, -107.75      ,
       -107.66666667, -107.58333333, -107.5       , -107.41666667,
       -107.33333333, -107.25      , -107.16666667, -107.08333333,
       -107.        , -106.91666667, -106.83333333, -106.75      ,
       -106.66666667, -106.58333333, -106.5       , -106.41666667,
       -106.33333333, -106.25      , -106.16666667, -106.08333333,
       -106.        , -105.91666667, -105.83333333, -105.75      ,
       -105.66666667, -105.58333333, -105.5       , -105.41666667,
       -105.33333333, -105.25      , -105.16666667, -105.08333333,
       -105.        , -104.91666667, -104.83333333, -104.75      ,
       -104.66666667, -104.58333333, -104.5       , -104.41666667,
       -104.33333333, -104.25      , -104.16666667, -104.08333333,
       -104.        , -103.91666667, -103.83333333, -103.75      ,
       -103.66666667, -103.58333333, -103.5       , -103.41666667,
       -103.33333333, -103.25      , -103.16666667, -103.08333333,
       -103.        , -102.91666667, -102.83333333, -102.75      ,
       -102.66666667, -102.58333333, -102.5       , -102.41666667,
       -102.33333333, -102.25      , -102.16666667, -102.08333333,
       -102.        , -101.91666667, -101.83333333, -101.75      ,
       -101.66666667, -101.58333333, -101.5       , -101.41666667,
       -101.33333333, -101.25      , -101.16666667, -101.08333333,
       -101.        , -100.91666667, -100.83333333, -100.75      ,
       -100.66666667, -100.58333333, -100.5       , -100.41666667,
       -100.33333333, -100.25      , -100.16666667, -100.08333333,
       -100.        ])

Modified CDL

(Note that it's easy to get confused (as I did) if the data variable dimensions have the same names as the tie point coordinate variables - I think it's legal, but for auxiliary coordinates is recommended against, so should no recommended for tie point coordinate variables, too ... but I need to think more about the consequences of that!)

$ cat example.cdl
netcdf example {
dimensions:
  lat = 121 ;
  lon = 121 ;
  ntp_lat = 2 ;
  ntp_lon = 2 ;

variables:
  float elev(lat, lon) ;
    elev:long_name = "elevation" ;
    elev:standard_name = "height" ;
    elev:units = "meters" ;
    elev:positive = "up" ;
    elev:coordinate_interpolation = "lat: lat_interpolation lon: lon_interpolation" ;

  int lat_interpolation ;
    lat_interpolation:interpolation_name = "linear" ;
    lat_interpolation:tie_point_mapping = "lat: lat_indices ntp_lat" ;
    lat_interpolation:computational_precision = "32" ;

  int lon_interpolation ;
    lon_interpolation:interpolation_name = "linear" ;
    lon_interpolation:tie_point_mapping = "lon: lon_indices ntp_lon" ;
    lon_interpolation:computational_precision = "32" ;


  float lat(ntp_lat) ;
    lat:standard_name = "latitude" ;
    lat:units = "degrees_north" ;
  float lon(ntp_lon) ;
    lon:standard_name = "longitude" ;
    lon:units = "degrees_east" ;

  int lat_indices(ntp_lat) ;
  int lon_indices(ntp_lon) ;

data:
  lat_indices = 0, 120 ;
  lon_indices = 0, 120 ;
  lat = 35.0, 45.0 ;
  lon = -110.0, -100.0 ;
}

[ethanrd]

Thanks David! Ah, yes, referencing the tie point variables in the coordinate_interpolation attribute makes more sense.

Hmm. This does get a bit confusing. Though I guess you can think of the computed coordinate variables as an expanded version of the tie point variables. Then using the same variable/dimension names would match how one might name the (auxiliary) coordinate variables when not using subsampling of coordinates. (Or am I still confused?)

[benbovy]

Though I guess you can think of the computed coordinate variables as an expanded version of the tie point variables.

Yes that makes sense to me, e.g., implement a decoding step that turns lat(ntp_lat), lon(ntp_lon) into lat(lat), lon(lon) when opening the dataset (https://github.com/dcherian/rasterix/pull/4#discussion_r2039525245).

[AndersMS]

Hi all, good to see all the well-known names. And thank you to Daniel for reaching out.

For those, like me, who are not fully familiar with affine transformations, I found this short and easy to read explanation of affine transformations in both 2D and 3D: Affine Transformations

Is it correctly understood that the GeoZarr community wish to use the affine transformation as a lossy compression of coordinates, or would they only apply the affine transformations in situations where the coordinates could be reconstructed without loss?

[ethanrd]

Hmm. Good question. I'm not really sure how the affine transformations are determined. I will ask in the GeoZarr discussion if someone can comment here.

[benbovy]

Affine transformation is the most common way to represent geospatial raster positions (https://gdal.org/en/stable/user/raster_data_model.html#affine-geotransform), so I think that the goal is rather the latter (i.e., reconstruct the coordinates from the affine transformation parameters).

There are also various operations that can be expressed as input affine parameters -> output affine parameters (e.g., data slicing), which is often better than expressing these as input explicit coordinates -> output explicit coordinates.

[geospatial-jeff]

the goal is rather the latter (i.e., reconstruct the coordinates from the affine transformation parameters).

Yes this is the goal. Use an affine transform (6 numbers) to describe the coordinates in zarr. Instead of an explicit y and x variable that assigns an explicit latitude/longitude to each position. It doesn't necessarily matter whether or not this transformation is lossy or lossless, section 8.3 of the CF conventions even says this:

For some applications the coordinates of a data variable can require considerably more storage than the data itself. Space may be saved in the netCDF file by storing a subsample of the coordinates that describe the data. The uncompressed coordinate and auxiliary coordinate variables can be reconstituted by interpolation, from the subsampled coordinate values to the domain of the data (i.e. the target domain). This process will likely result in a loss in accuracy (as opposed to precision) in the uncompressed variables, due to rounding and approximation errors in the interpolation calculations, but it is assumed that these errors will be small enough to not be of concern to users of the uncompressed dataset.

[geospatial-jeff]

Coordinate subsampling (section 8.3) with linear interpolation is in fact a drop-in replacement for an affine transformation. By definition every linear transformation is affine, but not all affine transformations are linear. There are two key differences between CF section 8.3 and how GDAL (GeoZarr) thinks about affine transforms:

• CF supports other interpolation methods besides linear. Can anyone comment on how frequently these other interpolation methods appear in the wild? I imagine linear and quadratic are the two most commonly used interpolation methods for geographic datasets.
• CF supports "sub-interpolation" which allows the user to provide more accurate interpolation. Imagine a linear affine transformation applied over a curved surface, the farther you get away from the origin the more inaccurate your position becomes. Sub-interpolation lets the user break up the affine space across multiple transformations, each transformation applying to a subset of the overall space. There is no concept of "sub-interpolation" in GDAL.

[davidhassell]

Can anyone comment on how frequently these other interpolation methods appear in the wild?

I suspect that neither occurs very often, if at all, in the wild at this time. The original use case was for remote sensing products, but as far as I'm aware, CF-standardised coordinate subsampling is not yet being used operationally in this field (very happy to proved wrong, here!). It was always in our minds that there would be other applications, though. For instance, we're discussing using it for storing some types of HEALPix