Weird GeoTransform encoding in NetCDF file causing issues with transform() #714

fbunt · 2023-11-30T21:47:36Z

fbunt
Nov 30, 2023

I'm working with some gridMET data and running into issues getting the affine transform for the dataset. The dataset looks like this:

ds = xr.open_dataset("fm1000_2011.nc", decode_coords="all")
print(ds)
<xarray.Dataset>
Dimensions:                    (lon: 1386, lat: 585, day: 365, crs: 1)
Coordinates:
  * lon                        (lon) float64 -124.8 -124.7 ... -67.1 -67.06
  * lat                        (lat) float64 49.4 49.36 49.32 ... 25.11 25.07
  * day                        (day) datetime64[ns] 2011-01-01 ... 2011-12-31
  * crs                        (crs) uint16 3
Data variables:
    dead_fuel_moisture_1000hr  (day, lat, lon) float32 ...
Attributes: (12/19)
    geospatial_bounds_crs:      EPSG:4326
    Conventions:                CF-1.6
    geospatial_bounds:          POLYGON((-124.7666666333333 49.40000000000000...
    geospatial_lat_min:         25.066666666666666
    geospatial_lat_max:         49.40000000000000
    geospatial_lon_min:         -124.7666666333333
    ...                         ...
    date:                       04 July 2019
    note1:                      The projection information for this file is: ...
    note2:                      Citation: Abatzoglou, J.T., 2013, Development...
    note3:                      Data in slices after last_permanent_slice (1-...
    note4:                      Data in slices after last_provisional_slice (...
    note5:                      Days correspond approximately to calendar day...

When I try to get the transform, this error occurs:

ds.rio.transform()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[703], line 1
----> 1 ds.rio.transform()

File ~/anaconda3/envs/rstools/lib/python3.10/site-packages/rioxarray/rioxarray.py:649, in XRasterBase.transform(self, recalc)
    636 def transform(self, recalc: bool = False) -> Affine:
    637     """
    638     Parameters
    639     ----------
   (...)
    647         The affine of the :obj:`xarray.Dataset` | :obj:`xarray.DataArray`
    648     """
--> 649     transform = self._cached_transform()
    650     if transform and (
    651         not transform.is_rectilinear or _affine_has_rotation(transform)
    652     ):
    653         if recalc:

File ~/anaconda3/envs/rstools/lib/python3.10/site-packages/rioxarray/rioxarray.py:579, in XRasterBase._cached_transform(self)
    575     transform = numpy.fromstring(
    576         self._obj.coords[self.grid_mapping].attrs["GeoTransform"], sep=" "
    577     )
    578     # Calling .tolist() to assure the arguments are Python float and JSON serializable
--> 579     return Affine.from_gdal(*transform.tolist())
    581 except KeyError:
    582     try:

TypeError: Affine.from_gdal() missing 1 required positional argument: 'e'

This is cause by the GeoTransform string encoded in the dataset only having 5 values:
'-124.7666666333333 0.041666666666666 0 49.400000000000000 -0.041666666666666'.
Inserting a 0 inbetween the last two values in the string produces a valid affine transform for the dataset.

Are 5-value gdal transform strings considered valid? This is the first time I've encountered one with 5 values. I'm wondering if I need to report this to the data maintainers (and add a special case in my code) or if I should file a bug in this repo.

Thanks!

ETA: It seems that the 5th coefficient is commonly 0. Is there a common or known convention of omitting it when it's 0?