issue_comments

@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.

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 ```

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.

@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.