Repository: corteva/rioxarray

Issue:
Use CRS.from_cf when retrieving the CRS
Since the `pyproj.CRS` class is there, might be worthwhile to use the `from_cf` option when getting the CRS in case they didn't use `open_rasterio`.


Discussion:

Patch:
diff --git a/docs/history.rst b/docs/history.rst
index a82c1a91..55ceddbc 100644
--- a/docs/history.rst
+++ b/docs/history.rst
@@ -4,6 +4,7 @@ History
 Latest
 -------
 - rio.reproject: change input kwarg dst_affine_width_height -> shape & transform (#125)
+- ENH: Use pyproj.CRS to read/write CF parameters (issue #124)
 
 0.0.27
 ------
@@ -11,7 +12,7 @@ Latest
 - ENH: Added optional `shape` argument to `rio.reproject` (pull #116)
 - Fix ``RasterioDeprecationWarning`` (pull #117)
 - BUG: Make rio.shape order same as rasterio dataset shape (height, width) (pull #121)
-- Fix open_rasterio() for WarpedVRT with specified src_crs (pydata/xarray/pull/4104 & pull 120)
+- Fix open_rasterio() for WarpedVRT with specified src_crs (pydata/xarray/pull/4104 & pull #120)
 - BUG: Use internal reprojection as engine for resampling window in merge (pull #123)
 
 0.0.26
diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py
index 70d270c7..139555b1 100644
--- a/rioxarray/rioxarray.py
+++ b/rioxarray/rioxarray.py
@@ -16,6 +16,7 @@ import warnings
 from uuid import uuid4
 
 import numpy as np
+import pyproj
 import rasterio.warp
 import xarray
 from affine import Affine
@@ -325,10 +326,12 @@ class XRasterBase(object):
             # look in grid_mapping
             grid_mapping_coord = self._obj.attrs.get("grid_mapping", DEFAULT_GRID_MAP)
             try:
-                crs_wkt = self._obj.coords[grid_mapping_coord].attrs["spatial_ref"]
-            except KeyError:
-                crs_wkt = self._obj.coords[grid_mapping_coord].attrs["crs_wkt"]
-            self.set_crs(crs_wkt, inplace=True)
+                self.set_crs(
+                    pyproj.CRS.from_cf(self._obj.coords[grid_mapping_coord].attrs),
+                    inplace=True,
+                )
+            except pyproj.exceptions.CRSError:
+                pass
         except KeyError:
             try:
                 # look in attrs for 'crs'
@@ -423,11 +426,10 @@ class XRasterBase(object):
             )
         # add grid mapping coordinate
         data_obj.coords[grid_mapping_name] = xarray.Variable((), 0)
+        grid_map_attrs = pyproj.CRS.from_user_input(data_obj.rio.crs).to_cf()
+        # spatial_ref is for compatibility with GDAL
         crs_wkt = crs_to_wkt(data_obj.rio.crs)
-        grid_map_attrs = dict()
         grid_map_attrs["spatial_ref"] = crs_wkt
-        # http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html#appendix-grid-mappings
-        # http://desktop.arcgis.com/en/arcmap/10.3/manage-data/netcdf/spatial-reference-for-netcdf-data.htm
         grid_map_attrs["crs_wkt"] = crs_wkt
         data_obj.coords[grid_mapping_name].rio.set_attrs(grid_map_attrs, inplace=True)
 
diff --git a/setup.py b/setup.py
index 0760c85f..42d48474 100644
--- a/setup.py
+++ b/setup.py
@@ -24,7 +24,7 @@ def get_version():
 with open("README.rst") as readme_file:
     readme = readme_file.read()
 
-requirements = ["pillow", "rasterio", "scipy", "xarray", "pyproj>=2"]
+requirements = ["pillow", "rasterio", "scipy", "xarray", "pyproj>=2.2"]
 
 test_requirements = ["pytest>=3.6", "pytest-cov", "mock"]
 


Test patch:
diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py
index 28ada416..eccaf22c 100644
--- a/test/integration/test_integration_rioxarray.py
+++ b/test/integration/test_integration_rioxarray.py
@@ -1330,6 +1330,26 @@ def test_get_crs_dataset():
     assert test_ds.rio.crs.to_epsg() == 4326
 
 
+def test_write_crs_cf():
+    test_da = xarray.DataArray(1)
+    test_da = test_da.rio.write_crs(4326)
+    assert test_da.attrs["grid_mapping"] == "spatial_ref"
+    assert test_da.rio.crs.to_epsg() == 4326
+    assert "spatial_ref" in test_da.spatial_ref.attrs
+    assert "crs_wkt" in test_da.spatial_ref.attrs
+    assert test_da.spatial_ref.attrs["grid_mapping_name"] == "latitude_longitude"
+
+
+def test_read_crs_cf():
+    test_da = xarray.DataArray(1)
+    test_da = test_da.rio.write_crs(4326)
+    assert test_da.attrs["grid_mapping"] == "spatial_ref"
+    attrs = test_da.spatial_ref.attrs
+    attrs.pop("spatial_ref")
+    attrs.pop("crs_wkt")
+    assert test_da.rio.crs.to_epsg() == 4326
+
+
 def test_get_crs_dataset__nonstandard_grid_mapping():
     test_ds = xarray.Dataset()
     test_ds = test_ds.rio.write_crs(4326, grid_mapping_name="frank")