Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

grid_mapping attrs lead to failure of cf writer #2604

Closed
zxdawn opened this issue Oct 19, 2023 · 2 comments
Closed

grid_mapping attrs lead to failure of cf writer #2604

zxdawn opened this issue Oct 19, 2023 · 2 comments

Comments

@zxdawn
Copy link
Member

zxdawn commented Oct 19, 2023

Describe the bug

cf writer works well for a DataArray with crs and geotransform written by rioxarray.
But, if we remove these spatial attrs and add the area attrs, cf writer failed.

To Reproduce

import numpy as np
from pyproj import CRS
import xarray as xr
import rioxarray
from satpy import Scene
from affine import Affine
from pyresample.geometry import AreaDefinition


def fake_data():
    test = xr.DataArray(data=np.random.rand(1157, 1235), dims=('y', 'x'))

    test.rio.write_transform(Affine(30.0, 0.0, 262056.03103060447, 0.0, -30.0, 2143592.6343544293), inplace=True)
    test.rio.write_crs(CRS.from_epsg(32643), inplace=True)

    return test

# write data using netcdf directly (version1)
test = fake_data()
test.to_netcdf('test.nc')

# add to Scene (version2)
scn = Scene()
scn['test'] = test
scn.save_datasets(writer="cf", filename='test.nc')

# remove spatial info and write area (version3)
scn = Scene()
# add area
area_def = AreaDefinition('id', 'description', 'proj_id',
                          CRS.from_epsg(32643), 1235, 1157, (-15.0, 1141.0, 1249.0, 15.0))
test.attrs['area'] = area_def
test = test.drop_vars(['spatial_ref'])
# add to Scene
scn['test'] = test
scn.save_datasets(writer="cf", filename='test.nc')

Expected behavior

cf writer works with deleted rioxarray info and manually added area attrs.

Actual results

Error output of version3 in code

File [~/miniconda3/envs/hyperch4/lib/python3.11/site-packages/xarray/coding/variables.py:181](https://file+.vscode-resource.vscode-cdn.net/Users/xinz/Documents/github/HyperCH4/notebooks/~/miniconda3/envs/hyperch4/lib/python3.11/site-packages/xarray/coding/variables.py:181), in safe_setitem(dest, key, value, name)
    179 if key in dest:
    180     var_str = f" on variable {name!r}" if name else ""
--> 181     raise ValueError(
    182         f"failed to prevent overwriting existing key {key} in attrs{var_str}. "
    183         "This is probably an encoding field used by xarray to describe "
    184         "how a variable is serialized. To proceed, remove this key from "
    185         "the variable's attributes manually."
    186     )
    187 dest[key] = value

ValueError: failed to prevent overwriting existing key grid_mapping in attrs. This is probably an encoding field used by xarray to describe how a variable is serialized. To proceed, remove this key from the variable's attributes manually.

Actually, I don't see the grid_mapping attrs in scn['test']:
image

And, I suppose the output netcdf file from version2 has wrong grid_mapping name:
image

Environment Info:

  • Satpy Version: 0.14.3.dev7585+gf018f1276
  • PyResample Version:1.27.1
@zxdawn
Copy link
Member Author

zxdawn commented Oct 20, 2023

Well. I checked ds['test'] at this line:

for group_name, ds in grouped_datasets.items():

It does have grid_mapping attrs!

Anyway, this issue can be fixed by manually passing exclude_attrs=['grid_mapping'] to save_datasets.

@zxdawn zxdawn closed this as completed Oct 20, 2023
@zxdawn
Copy link
Member Author

zxdawn commented Oct 21, 2023

Oh, I mixed rioxarray with pyresample in this case. Because pyresample AreDefinition can handle similar things and is more powerful, we can just use AreaDefinition:

import numpy as np
import xarray as xr
import rioxarray
from satpy import Scene
from pyresample.geometry import AreaDefinition

test = xr.DataArray(data=np.random.rand(1157, 1235), dims=('y', 'x'), name='test')

area_def = AreaDefinition.from_ul_corner(area_id='area_id', projection='EPSG:32643',
                                        shape=(1157, 1235),
                                        upper_left_extent=(262056.03103060447, 2143592.6343544293),
                                        resolution=(30, 30))

test.attrs['area'] = area_def

scn = Scene()
scn['test'] = test
scn.save_datasets(writer="cf", filename='test.nc')

The exported netcdf can be decoded successfully:

ds = xr.open_dataset('test.nc', decode_coords='all')

print(ds.rio.crs)
# CRS.from_epsg(32643)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant