issue_comments

Feel free to re-open if this becomes something of interest. Thanks!

I'm going to leave this decision to @shoyer and @jhamman !

If we go forward with this, we have a couple of possibilities:

1. add crs to a coordinate variable and keep the attribute too (not so nice because doubled)
2. add crs to a coordinate variable and deprecate the crs attribute

Same with _FillValue.

@fmaussion, I guess the main reason for the proposed change is to be able to read in the raster using open_rasterio() and be able to generate a netCDF using to_netcdf() that would be compatible with the rasterio model when using rasterio.open().

Also, since: "xarray.Dataset is an in-memory representation of a netCDF file", to me it makes more sense to store the data in a standard netCDF location that is consistent regardless of the backend. In doing so, other parts of xarray and libraries build off of xarray won't be required to find the same data in multiple locations.

Thanks for taking the time to consider this idea. I have used it quite a bit and it has made the datasets compatible across the various GIS tools I use (rasterio, GDAL, QGIS). It has definitely made my life easier and so I wanted to share what I have learned. But, it is just an idea. I really appreciate the work of the xarray developers and whatever y'all decide is fine with me.

@fmaussion I completely agree except now that all of this is being brought up I see why it may have been better to put the 'crs' in the coordinates of the DataArray returned by open_rasterio. Since two DataArrays in different projections are not and should not be considered "mergeable". But I can also see how this walks the line of special handling of a data format by trying to interpret things in a certain way, but that line is already pretty blurry in the case of reading rasterio compatible datasets.

Again, I'm yet to be convinced that this logic should live in xarray, regardless if netCDF or rasterio is used as a backend. The job of xarray is to read the rasterio files in a "xarray way" (e.g. in order to leverage dask) and exposing the file attributes to the user, ideally with semantics/attribute names which are as close as possible as the underlying data model (rasterio).

All the logic you describe could live in a dedicated library which would become a wrapper around xarray's open_* functions, and converting them to the format you describe.

This is the output of using xarray with a standard GOES-16 ABI L1B data file:

```python In [2]: import xarray as xr

In [3]: nc = xr.open_dataset('OR_ABI-L1b-RadF-M3C01_G16_s20181741200454_e20181741211221_c20181741211264.nc')

In [4]: nc.data_vars['Rad'] Out[4]: <xarray.DataArray 'Rad' (y: 10848, x: 10848)> array([[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]], dtype=float32) Coordinates: t datetime64[ns] ... * y (y) float32 0.151858 0.15183 0.151802 0.151774 0.151746 ... * x (x) float32 -0.151858 -0.15183 -0.151802 -0.151774 -0.151746 ... y_image float32 ... x_image float32 ... Attributes: long_name: ABI L1b Radiances standard_name: toa_outgoing_radiance_per_unit_wavelength sensor_band_bit_depth: 10 valid_range: [ 0 1022] units: W m-2 sr-1 um-1 resolution: y: 0.000028 rad x: 0.000028 rad grid_mapping: goes_imager_projection cell_methods: t: point area: point ancillary_variables: DQF

In [5]: nc.data_vars['goes_imager_projection'] Out[5]: <xarray.DataArray 'goes_imager_projection' ()> array(-2147483647, dtype=int32) Coordinates: t datetime64[ns] ... y_image float32 ... x_image float32 ... Attributes: long_name: GOES-R ABI fixed grid projection grid_mapping_name: geostationary perspective_point_height: 35786023.0 semi_major_axis: 6378137.0 semi_minor_axis: 6356752.31414 inverse_flattening: 298.2572221 latitude_of_projection_origin: 0.0 longitude_of_projection_origin: -75.0 sweep_angle_axis: x ```

That is strange. Was the variable put in data_vars and not in coords?

I wouldn't expect it to add crs if there wasn't a grid_mapping specified, but if it was then I would. In a simple test where I did xr.open_dataset('my_nc.nc') which has a grid_mapping attribute, xarray does nothing special to create a crs or other named coordinate variable referencing the associated grid_mapping variable.

I don't think so. The crs should already exist on the dataset beforehand and I think it may be a bit too much to expect it to add it to the dataset in that function as there are some that don't have it and not all xarray datasets are geospatially located.

@snowman2 This should mean that open_dataset should handle CRS specially too, right? Currently it doesn't seem to do anything special for the coordinate variable pointed to by grid_mapping.