Problem in reprojecting MSG SEVIRI data

[anikfal]

Reading a MSG SEVIRI datafile and trying to reproject it to "EPSG:4326":

from satpy import Scene
seviri_file = ["W_XX-EUMETSAT-Darmstadt,VIS+IR+HRV+IMAGERY,MSG1+SEVIRI_C_EUMG_20170709100010.nc"]
scn = Scene(reader="seviri_l1b_nc", filenames=seviri_file)
scn.load(['IR_108', 'IR_120'])
ir10 = scn['IR_120']
ir10.rio.reproject("EPSG:4326")

Output error:

Traceback (most recent call last):
File "/home/anikfal/extra_codes/skimage_codes/2.py", line 6, in
ir10.rio.reproject("EPSG:4326")
File "/home/anikfal/miniconda3/envs/mysat/lib/python3.10/site-packages/rioxarray/raster_array.py", line 433, in reproject
raise MissingCRS(
rioxarray.exceptions.MissingCRS: CRS not found. Please set the CRS with 'rio.write_crs()'. Data variable: ch10

It has a crs attribute, but it seems to be undefined:
ir10.crs

<xarray.DataArray 'crs' ()>
array(<Derived Projected CRS: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unk ...>
Name: unknown
Axis Info [cartesian]:

• E[east]: Easting (metre)
• N[north]: Northing (metre)
  Area of Use:
• undefined
  Coordinate Operation:
• name: unknown
• method: Geostationary Satellite (Sweep Y)
  Datum: unknown
• Ellipsoid: unknown
• Prime Meridian: Greenwich
  , dtype=object)
  Coordinates:
  crs object PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",E...

[djhoese]

Correct me if I'm wrong but I believe rioxarray uses DataArray coordinates, not attributes to hold the crs information. So let's be careful which we're referring to here. Could you please print out:

print(ir10.coords)

And if we're talking about attributes we should access it from them from the .attrs dictionary.

More importantly though, although we would love to advertise compatibility with rioxarray, it is not something we test or guarantee or advertise (unless I left something around in the documentation somewhere). Satpy uses a .attrs['area'] AreaDefinition object to define all geolocation for gridded data like SEVIRI. We did at one time add a crs coordinate along with a y and x coordinate, but I'm not sure we do that still. It is possible our crs coordinate, if it is added, is not done in the same exact way that rioxarray is doing.

Let's see if we can figure this out and reintroduce rioxarray compatibility.

[anikfal]

ir10.coords

Coordinates:
lat (y, x) float32 dask.array<chunksize=(1015, 1216), meta=np.ndarray>
lon (y, x) float32 dask.array<chunksize=(1015, 1216), meta=np.ndarray>
acq_time (y) datetime64[ns] 2017-07-09T10:11:16.114000 ... 2017-07-09T10...
crs object PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",...
* y (y) float64 4.424e+06 4.421e+06 4.418e+06 ... 1.385e+06 1.382e+06
* x (x) float64 -6.646e+05 -6.616e+05 ... 2.972e+06 2.975e+06

Actually my main concern is to rio.reproject_match an ERA5 data (u10 variable) with SEVIRI (ir10 variable). So raster calculations will be possible. It complains about crs:

u10_match_ir10 = u10.rio.reproject_match(ir10)

output errors:

---

MissingCRS Traceback (most recent call last)
Input In [2], in <cell line: 1>()
----> 1 var = u10.rio.reproject_match(ir10)

File ~/miniconda3/envs/mysat/lib/python3.10/site-packages/rioxarray/raster_array.py:554, in RasterArray.reproject_match(self, match_data_array, resampling, **reproject_kwargs)
520 def reproject_match(
521 self,
522 match_data_array: Union[xarray.DataArray, xarray.Dataset],
523 resampling: Resampling = Resampling.nearest,
524 **reproject_kwargs,
525 ) -> xarray.DataArray:
526 """
527 Reproject a DataArray object to match the resolution, projection,
528 and region of another DataArray.
(...)
552 match_data_array.
553 """
--> 554 reprojected_data_array = self.reproject(
555 match_data_array.rio.crs,
556 transform=match_data_array.rio.transform(recalc=True),
557 shape=match_data_array.rio.shape,
558 resampling=resampling,
559 **reproject_kwargs,
560 )
561 # hack to resolve: corteva/rioxarray#298
562 # may be resolved in the future by flexible indexes:
563 # pydata/xarray#4489 (comment)
564 x_attrs = reprojected_data_array[reprojected_data_array.rio.x_dim].attrs.copy()

File ~/miniconda3/envs/mysat/lib/python3.10/site-packages/rioxarray/raster_array.py:433, in RasterArray.reproject(self, dst_crs, resolution, shape, transform, resampling, nodata, **kwargs)
431 raise RioXarrayError("resolution cannot be used with shape or transform.")
432 if self.crs is None:
--> 433 raise MissingCRS(
434 "CRS not found. Please set the CRS with 'rio.write_crs()'."
435 f"{_get_data_var_message(self._obj)}"
436 )
437 gcps = self.get_gcps()
438 if gcps:

MissingCRS: CRS not found. Please set the CRS with 'rio.write_crs()'. Data variable: u10

[djhoese]

I think the simplest solution, for now, is to do as the error message suggests and do ir10.rio.write_crs(ir10.attrs['area'].crs, inplace=True). I would also accept a pull request to Satpy that modifies the CRS-related helper functions like add_crs_x_y_coords to check if rioxarray is installed and if so do write_crs in this same way.

Some things to note about details here: rioxarray uses rasterio's CRS object, but Satpy/pyresample use pyproj's CRS object. Luckily thanks to the efforts of Alan Snow (not mentioning him here to avoid unnecessary notifications) pyproj's CRS objects are accepted by rasterio/rioxarray. By doing write_crs you may be overwriting some coordinates that Satpy is setting, but luckily Satpy doesn't really use these...yet.

[anikfal]

Yeah, also u10 needs the same process:

u10.rio.write_crs("EPSG:4326", inplace=True)

Now u10.rio.reproject_match(ir10) works.