Ideas for crs accessor to xarray about geoxarray

[ghiggi]

Additional notes:

• GDAL does not support grid_mapping defined as: "CRS: x y"
• GDAL does not support grid_mapping with 2 CRS: "CRS: x y CRSWGS84: lat lon"
• GDAL looks to the spatial_ref (crs_wkt) attribute in the CRS coordinate
• GDAL expects the crs_wkt to start with GEOGCS but pyproj provides the crs_wkt starting with GEOGCRS !!!

[djhoese]

Sorry it took so long for me to get to this @ghiggi. I know we talked a little bit on slack about this, but when I started reading the issue and realized how long it was I wasn't going to have time to give it a proper review. I'm hoping I can start that discussion now...

First, thanks for putting all these ideas together. It is clear you've run into problems with existing tools and tried to come up with interfaces that support as much as possible. These are, in general, not easy things to figure out and nail down.

My main issue with what you've suggested here is that it conflates/confuses/combines two or more separate topics. The CRS is specifically the coordinate system which has a datum and possibly some bounds associated with it, but it does not have an idea of pixel resolution or shape/size or even coordinates (at least not ones specific to pixels or an array of data). This is made even more confusing by pyresample's hack of generating cartopy CRS objects where the bounds of the CRS are set based on the bounds of the data (a bad decision in hindsight by me). You might be trying to handle that with your area_crs keyword argument in some of the conversion methods, but I'm not sure I understand what that kwarg does.

So I think I'd be hesitant to want to put this type of functionality in an accessor called .crs. You also refer to all of this stuff as being related to CF, but it isn't clear to me why CF matters for most of this functionality. For the stuff that is CF specific or has to be done differently to make CF happy, would this functionality be more appropriate for cf_xarray (https://cf-xarray.readthedocs.io/en/latest/)?

I think the parts of this functionality that are xarray-specific could go in geoxarray in the sense that they need to know about coordinates and the CRS and other potentially other DataArrays in a Dataset. Otherwise I would hope that the functionality would go in GDAL, or rasterio, or pyproj, and then come "out of the box" by accessing a pyproj CRS object. For the parts of CF that aren't supported by GDAL or other libraries, could lowering the CF standard being supported allow for more flexibility currently (CF 1.6 or 1.7 instead of 1.8)?

Side note: I think cartopy can now accept pyproj CRS objects since it now subclasses from them. Or at least they can be easily made. Similarly, I would hope that libraries like pyresample or rioxarray would not require a user to explicitly convert to a certain CRS/area-like object before being called, but would rather depend on accessors like those in geoxarray/cf_xarray/rioxarray to get the information they need.

CC @snowman2 who is the maintainer or is a contributor to at least 25% of the packages mentioned in this issue.

[rabernat]

I would definitely pay attention to what odc-geo is doing: https://github.com/opendatacube/odc-geo