Repository: Ouranosinc/xclim

Issue:
crs coordinate added during the `subset_shape` call
* xclim version: 0.15.0
* Python version: 3.8
* Operating System: ubuntu

### Description

This happened in the finch test suite: https://github.com/bird-house/finch/pull/107/commits/7f5ffff1ed422df8ddd12286f5990cf3fbe389c2

When I call `xclim.subset.subset_shape`, a 'crs' coordinate is added to the dataset: https://github.com/Ouranosinc/xclim/blob/a838560caa6b03678dbd1a25799f5fb8b562f9f4/xclim/subset.py#L560

However, when writing the dataset to disk, I get this error: 

```
File ".../python3.8/site-packages/xarray/backends/api.py", line
205, in check_attr
    raise TypeError(
TypeError: Invalid value for attr: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs lon_wrap=180 +type=crs must be a number, a string, an ndarray or a list/tuple of numbers/strings for serialization to
netCDF files
```

For now, I just dropped this coordinate in finch...

Discussion:
Zeitsperre: For more information on relevant fields: http://cfconventions.org/wkt-proj-4.html
coveralls: ## Pull Request Test Coverage Report for [Build 1718](https://coveralls.io/builds/29456563)

* **4** of **4**   **(100.0%)**  changed or added relevant lines in **1** file are covered.
* **6** unchanged lines in **1** file lost coverage.
* Overall coverage decreased (**-0.2%**) to **89.064%**

---


|  Files with Coverage Reduction | New Missed Lines | % |
| :-----|--------------|--: |
| [xclim/core/checks.py](https://coveralls.io/builds/29456563/source?filename=xclim%2Fcore%2Fchecks.py#L89) | 6 | 56.86% |
<!-- | **Total:** | **6** |  | -->

|  Totals | [![Coverage Status](https://coveralls.io/builds/29456563/badge)](https://coveralls.io/builds/29456563) |
| :-- | --: |
| Change from base [Build 1712](https://coveralls.io/builds/29385288): |  -0.2% |
| Covered Lines: | 2142 |
| Relevant Lines: | 2405 |

---
##### 💛  - [Coveralls](https://coveralls.io)

Zeitsperre: Never seen that error referring to .tox/bin/python before. I need to look into this more.
RondeauG: To add to this, I just had an error with CanRCM4. It uses a [0-360] longitude, which triggers `wrap_lons` in `subset_shape`

    if wrap_lons:
        ds_copy.coords["crs"] = 0
        ds_copy.coords["crs"].attrs = dict(spatial_ref=wgs84)`

The error:
`File "/scen3/rondeau/.conda/envs/pyrondeau/lib/python3.7/site-packages/xarray/core/dataset.py", line 1545, in to_netcdf
    invalid_netcdf=invalid_netcdf,`
`File "/scen3/rondeau/.conda/envs/pyrondeau/lib/python3.7/site-packages/xarray/backends/api.py", line 1040, in to_netcdf
    _validate_attrs(dataset)`
`File "/scen3/rondeau/.conda/envs/pyrondeau/lib/python3.7/site-packages/xarray/backends/api.py", line 219, in _validate_attrs
    check_attr(k, v)`
`File "/scen3/rondeau/.conda/envs/pyrondeau/lib/python3.7/site-packages/xarray/backends/api.py", line 209, in check_attr
    "files".format(value)`
`TypeError: Invalid value for attr: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs lon_wrap=180 +type=crs must be a number, a string, an ndarray or a list/tuple of numbers/strings for serialization to netCDF files`

I think this PR will address this issue?
Zeitsperre: @RondeauG Yes, this should fix that particular issue. It should be good to be merged in a few hours, after review.
Zeitsperre: Adding this here for future reference: http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/ch05s06.html
Zeitsperre: @aulemahal Does the most recent change work for you? If so, we're good to go.

aulemahal: It looks fine, but I'm not sure I understand. Is this grid_mapping attribute not needed when you only have lat-lon dims? And, should'nt we test if the variables as lat and lon instead of the number of dims?
Zeitsperre: Testing for lat and lon as dims makes much more sense. I'll update this.
Zeitsperre: @huard Tests now using an pytest fixture to create a temporary netcdf.
Zeitsperre: New issue: `xarray`@0.15.1 is breaking some tests. Will open an issue to investigate.

Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index c70e0593..d282d386 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -5,6 +5,7 @@ History
 0.16.x
 ------
 * Allow lazy units conversion
+* CRS definitions of projected DataSets are now written to file according to Climate and Forecast-convention standards
 
 0.15.x (2020-03-12)
 -------------------
diff --git a/environment.yml b/environment.yml
index d39d3c3e..37c16c7d 100644
--- a/environment.yml
+++ b/environment.yml
@@ -4,7 +4,7 @@ channels:
  - defaults
 dependencies:
  - numpy>=1.16
- - xarray>=0.14.1
+ - xarray>=0.14.1,<=0.15.0
  - scipy>=1.2
  - numba
  - pandas>=0.23
diff --git a/setup.cfg b/setup.cfg
index 61d5926d..56c1566e 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.15.0
+current_version = 0.15.1-beta
 commit = False
 tag = False
 parse = (?P<major>\d+)\.(?P<minor>\d+).(?P<patch>\d+)(\-(?P<release>[a-z]+))?
diff --git a/setup.py b/setup.py
index 382e2a9b..b326344a 100644
--- a/setup.py
+++ b/setup.py
@@ -12,7 +12,7 @@ URL = "https://github.com/Ouranosinc/xclim"
 AUTHOR = "Travis Logan"
 AUTHOR_EMAIL = "logan.travis@ouranos.ca"
 REQUIRES_PYTHON = ">=3.6.0"
-VERSION = "0.15.0"
+VERSION = "0.15.1-beta"
 LICENSE = "Apache Software License 2.0"
 
 with open("README.rst") as readme_file:
@@ -23,7 +23,7 @@ with open("HISTORY.rst") as history_file:
 
 requirements = [
     "numpy>=1.16",
-    "xarray>=0.14.1",
+    "xarray>=0.14.1,<=0.15.0",
     "scipy>=1.2",
     "numba",
     "pandas>=0.23",
diff --git a/xclim/__init__.py b/xclim/__init__.py
index 1f2fb30f..55cac2d8 100644
--- a/xclim/__init__.py
+++ b/xclim/__init__.py
@@ -8,4 +8,4 @@ from xclim.indicators import seaIce
 
 __author__ = """Travis Logan"""
 __email__ = "logan.travis@ouranos.ca"
-__version__ = "0.15.0"
+__version__ = "0.15.1-beta"
diff --git a/xclim/subset.py b/xclim/subset.py
index 0d39c020..778876c5 100644
--- a/xclim/subset.py
+++ b/xclim/subset.py
@@ -486,7 +486,7 @@ def subset_shape(
 
     if ds_copy.lon.size == 0 or ds_copy.lat.size == 0:
         raise ValueError(
-            "No gridcell centroids found within provided polygon bounding box. "
+            "No grid cell centroids found within provided polygon bounding box. "
             'Try using the "buffer" option to create an expanded area'
         )
 
@@ -530,7 +530,7 @@ def subset_shape(
             CRS(4326).to_proj4() != raster_crs
         ):
             warnings.warn(
-                "CRS definitions are not similar or both not using WGS84 datum. Caveat emptor.",
+                "CRS definitions are not similar or both not using WGS84 datum. Tread with caution.",
                 UserWarning,
                 stacklevel=3,
             )
@@ -541,8 +541,8 @@ def subset_shape(
 
     if np.all(mask_2d.isnull()):
         raise ValueError(
-            "No gridcell centroids found within provided polygon. "
-            'Try using the "buffer" option to create an expanded areas or verify polygon '
+            f"No grid cell centroids found within provided polygon bounds ({poly.bounds}). "
+            'Try using the "buffer" option to create an expanded areas or verify polygon.'
         )
 
     # loop through variables
@@ -555,10 +555,15 @@ def subset_shape(
         mask_2d = mask_2d.dropna(dim, how="all")
     ds_copy = ds_copy.sel({dim: mask_2d[dim] for dim in mask_2d.dims})
 
-    # Add a CRS definition as a coordinate for reference purposes
+    # Add a CRS definition using CF conventions and as a global attribute WKT for reference purposes
     if wrap_lons:
-        ds_copy.coords["crs"] = 0
-        ds_copy.coords["crs"].attrs = dict(spatial_ref=wgs84)
+        ds_copy.attrs["crs"] = wgs84.to_string()
+        ds_copy["crs"] = 1
+        ds_copy["crs"].attrs.update(wgs84.to_cf())
+
+        for v in ds_copy.variables:
+            if {"lat", "lon"}.issubset(set(ds_copy[v].dims)):
+                ds_copy[v].attrs["grid_mapping"] = "crs"
 
     if isinstance(ds, xarray.DataArray):
         return ds._from_temp_dataset(ds_copy)


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index ce5a1df3..d2f05e60 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,3 +1,5 @@
+from pathlib import Path
+
 import numpy as np
 import pandas as pd
 import pytest
@@ -6,6 +8,11 @@ import xarray as xr
 from xclim.core.calendar import calendars
 
 
+@pytest.fixture
+def tmp_netcdf_filename(tmpdir):
+    return Path(tmpdir).joinpath("testfile.nc")
+
+
 @pytest.fixture
 def tas_series():
     def _tas_series(values, start="7/1/2000"):
diff --git a/tests/test_subset.py b/tests/test_subset.py
index 3c877a05..f24a88d9 100644
--- a/tests/test_subset.py
+++ b/tests/test_subset.py
@@ -608,7 +608,7 @@ class TestSubsetShape:
                     sub[vari].sel(lon=lon1, lat=lat1), ds[vari].sel(lon=lon1, lat=lat1)
                 )
 
-    def test_wraps(self):
+    def test_wraps(self, tmp_netcdf_filename):
         ds = xr.open_dataset(self.nc_file)
 
         # Polygon crosses meridian, a warning should be raised
@@ -629,6 +629,14 @@ class TestSubsetShape:
             float(np.mean(subtas.isel(time=0))), 281.091553
         )
 
+        assert sub.crs.prime_meridian_name == "Greenwich"
+        assert sub.crs.grid_mapping_name == "latitude_longitude"
+
+        sub.to_netcdf(tmp_netcdf_filename)
+        assert tmp_netcdf_filename.exists()
+        with xr.open_dataset(filename_or_obj=tmp_netcdf_filename) as f:
+            assert {"tas", "crs"}.issubset(set(f.data_vars))
+
     def test_no_wraps(self):
         ds = xr.open_dataset(self.nc_file)
 
@@ -695,7 +703,7 @@ class TestSubsetShape:
             not in [q.message for q in record]
         )
 
-    def test_small_poly_buffer(self):
+    def test_small_poly_buffer(self, tmp_netcdf_filename):
         ds = xr.open_dataset(self.nc_file)
 
         with pytest.raises(ValueError):
@@ -709,6 +717,14 @@ class TestSubsetShape:
         assert len(sub.lon.values) == 3
         assert len(sub.lat.values) == 3
 
+        assert sub.crs.prime_meridian_name == "Greenwich"
+        assert sub.crs.grid_mapping_name == "latitude_longitude"
+
+        sub.to_netcdf(tmp_netcdf_filename)
+        assert tmp_netcdf_filename.exists()
+        with xr.open_dataset(filename_or_obj=tmp_netcdf_filename) as f:
+            assert {"tas", "crs"}.issubset(set(f.data_vars))
+
     def test_mask_multiregions(self):
         ds = xr.open_dataset(self.nc_file)
         regions = gpd.read_file(self.multi_regions_geojson)