diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml
index 9bfec108..87b5436e 100644
--- a/.github/workflows/tests.yaml
+++ b/.github/workflows/tests.yaml
@@ -37,9 +37,9 @@ jobs:
             rasterio-version: '*'
             xarray-version: '*'
           - os: ubuntu-latest
-            python-version: 3.6
+            python-version: 3.7
             rasterio-version: 1.1
-            xarray-version: 0.16
+            xarray-version: '*'
     steps:
       - uses: actions/checkout@v2
 
diff --git a/docs/history.rst b/docs/history.rst
index f7475a45..38575cf1 100644
--- a/docs/history.rst
+++ b/docs/history.rst
@@ -3,6 +3,9 @@ History
 
 Latest
 ------
+- DEP: Python 3.6+ (issue #215)
+- DEP: xarray 0.17+ (needed for issue #282)
+- REF: Store `grid_mapping` in `encoding` instead of `attrs` (issue #282)
 
 0.3.2
 -----
diff --git a/rioxarray/__init__.py b/rioxarray/__init__.py
index 05ca5c36..2fe3e212 100644
--- a/rioxarray/__init__.py
+++ b/rioxarray/__init__.py
@@ -5,12 +5,7 @@
 
 import rioxarray.raster_array  # noqa
 import rioxarray.raster_dataset  # noqa
+from rioxarray._io import open_rasterio  # noqa
 from rioxarray._options import set_options  # noqa
 from rioxarray._show_versions import show_versions  # noqa
 from rioxarray._version import __version__  # noqa
-
-try:
-    # This requires xarray >= 0.12.3
-    from rioxarray._io import open_rasterio  # noqa
-except ImportError:
-    from xarray import open_rasterio  # noqa
diff --git a/rioxarray/_io.py b/rioxarray/_io.py
index db7333b4..255d4fb7 100644
--- a/rioxarray/_io.py
+++ b/rioxarray/_io.py
@@ -607,6 +607,8 @@ def _handle_encoding(result, mask_and_scale, masked, da_name):
     """
     Make sure encoding handled properly
     """
+    if "grid_mapping" in result.attrs:
+        variables.pop_to(result.attrs, result.encoding, "grid_mapping", name=da_name)
     if mask_and_scale:
         if "scale_factor" in result.attrs:
             variables.pop_to(
@@ -851,7 +853,6 @@ def open_rasterio(
 
     # handle encoding
     _handle_encoding(result, mask_and_scale, masked, da_name)
-
     # Affine transformation matrix (always available)
     # This describes coefficients mapping pixel coordinates to CRS
     # For serialization store as tuple of 6 floats, the last row being
@@ -865,11 +866,7 @@ def open_rasterio(
         result = _prepare_dask(result, riods, filename, chunks)
 
     # Make the file closeable
-    try:
-        # xarray 0.17 +
-        result.set_close(manager.close)
-    except AttributeError:
-        result._file_obj = manager
+    result.set_close(manager.close)
     result.rio._manager = manager
     # add file path to encoding
     result.encoding["source"] = riods.name
diff --git a/rioxarray/raster_array.py b/rioxarray/raster_array.py
index 287ad8f2..9f68d9a0 100644
--- a/rioxarray/raster_array.py
+++ b/rioxarray/raster_array.py
@@ -50,9 +50,6 @@ def _generate_attrs(src_data_array, dst_nodata):
     if src_data_array.rio.encoded_nodata is None and fill_value is not None:
         new_attrs["_FillValue"] = fill_value
 
-    # add raster spatial information
-    new_attrs["grid_mapping"] = src_data_array.rio.grid_mapping
-
     return new_attrs
 
 
@@ -75,6 +72,7 @@ def _add_attrs_proj(new_data_array, src_data_array):
     new_data_array.rio.set_attrs(new_attrs, inplace=True)
 
     # make sure projection added
+    new_data_array.rio.write_grid_mapping(src_data_array.rio.grid_mapping, inplace=True)
     new_data_array.rio.write_crs(src_data_array.rio.crs, inplace=True)
     new_data_array.rio.write_coordinate_system(inplace=True)
     new_data_array.rio.write_transform(inplace=True)
diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py
index 54a54c72..4341db2f 100644
--- a/rioxarray/rioxarray.py
+++ b/rioxarray/rioxarray.py
@@ -230,10 +230,11 @@ def grid_mapping(self):
         """
         str: The CF grid_mapping attribute. 'spatial_ref' is the default.
         """
-        try:
-            return self._obj.attrs["grid_mapping"]
-        except KeyError:
-            pass
+        grid_mapping = self._obj.encoding.get(
+            "grid_mapping", self._obj.attrs.get("grid_mapping")
+        )
+        if grid_mapping is not None:
+            return grid_mapping
         grid_mapping = DEFAULT_GRID_MAP
         # search the dataset for the grid mapping name
         if hasattr(self._obj, "data_vars"):
@@ -245,11 +246,12 @@ def grid_mapping(self):
                     self._obj[var].rio.y_dim
                 except DimensionError:
                     continue
-                try:
-                    grid_mapping = self._obj[var].attrs["grid_mapping"]
+                var_grid_mapping = self._obj[var].encoding.get(
+                    "grid_mapping", self._obj[var].attrs.get("grid_mapping")
+                )
+                if var_grid_mapping is not None:
+                    grid_mapping = var_grid_mapping
                     grid_mappings.add(grid_mapping)
-                except KeyError:
-                    pass
             if len(grid_mappings) > 1:
                 raise RioXarrayError("Multiple grid mappings exist.")
         return grid_mapping
@@ -279,12 +281,22 @@ def write_grid_mapping(self, grid_mapping_name=DEFAULT_GRID_MAP, inplace=False):
                 except DimensionError:
                     continue
 
-                data_obj[var].rio.update_attrs(
+                # remove grid_mapping from attributes if it exists
+                # and update the grid_mapping in encoding
+                new_attrs = dict(data_obj[var].attrs)
+                new_attrs.pop("grid_mapping", None)
+                data_obj[var].rio.update_encoding(
                     dict(grid_mapping=grid_mapping_name), inplace=True
-                ).rio.set_spatial_dims(x_dim=x_dim, y_dim=y_dim, inplace=True)
-        return data_obj.rio.update_attrs(
+                ).rio.update_attrs(new_attrs, inplace=True).rio.set_spatial_dims(
+                    x_dim=x_dim, y_dim=y_dim, inplace=True
+                )
+        # remove grid_mapping from attributes if it exists
+        # and update the grid_mapping in encoding
+        new_attrs = dict(data_obj.attrs)
+        new_attrs.pop("grid_mapping", None)
+        return data_obj.rio.update_encoding(
             dict(grid_mapping=grid_mapping_name), inplace=True
-        )
+        ).rio.update_attrs(new_attrs, inplace=True)
 
     def write_crs(self, input_crs=None, grid_mapping_name=None, inplace=False):
         """
@@ -591,6 +603,57 @@ def update_attrs(self, new_attrs, inplace=False):
         data_attrs.update(**new_attrs)
         return self.set_attrs(data_attrs, inplace=inplace)
 
+    def set_encoding(self, new_encoding, inplace=False):
+        """
+        Set the encoding of the dataset/dataarray and reset
+        rioxarray properties to re-search for them.
+
+        .. versionadded:: 0.4
+
+        Parameters
+        ----------
+        new_encoding: dict
+            A dictionary for encoding.
+        inplace: bool, optional
+            If True, it will write to the existing dataset. Default is False.
+
+        Returns
+        -------
+        :obj:`xarray.Dataset` | :obj:`xarray.DataArray`:
+            Modified dataset with new attributes.
+        """
+        data_obj = self._get_obj(inplace=inplace)
+        # set the attributes
+        data_obj.encoding = new_encoding
+        # reset rioxarray properties depending
+        # on attributes to be generated
+        data_obj.rio._nodata = None
+        data_obj.rio._crs = None
+        return data_obj
+
+    def update_encoding(self, new_encoding, inplace=False):
+        """
+        Update the encoding of the dataset/dataarray and reset
+        rioxarray properties to re-search for them.
+
+        .. versionadded:: 0.4
+
+        Parameters
+        ----------
+        new_encoding: dict
+            A dictionary with encoding values to update with.
+        inplace: bool, optional
+            If True, it will write to the existing dataset. Default is False.
+
+        Returns
+        -------
+        :obj:`xarray.Dataset` | :obj:`xarray.DataArray`:
+            Modified dataset with updated attributes.
+        """
+        data_encoding = dict(self._obj.encoding)
+        data_encoding.update(**new_encoding)
+        return self.set_encoding(data_encoding, inplace=inplace)
+
     def set_spatial_dims(self, x_dim, y_dim, inplace=True):
         """
         This sets the spatial dimensions of the dataset.
diff --git a/setup.py b/setup.py
index 1562c2e8..69f24ce9 100644
--- a/setup.py
+++ b/setup.py
@@ -23,7 +23,7 @@ def get_version():
 with open("README.rst") as readme_file:
     readme = readme_file.read()
 
-requirements = ["rasterio", "scipy", "xarray", "pyproj>=2.2"]
+requirements = ["rasterio", "scipy", "xarray>=0.17", "pyproj>=2.2"]
 
 test_requirements = ["pytest>=3.6", "pytest-cov", "dask"]
 doc_requirements = ["sphinx-click==1.1.0", "nbsphinx", "sphinx_rtd_theme"]
@@ -52,7 +52,6 @@ def get_version():
         "Natural Language :: English",
         "Topic :: Scientific/Engineering :: GIS",
         "Programming Language :: Python :: 3",
-        "Programming Language :: Python :: 3.6",
         "Programming Language :: Python :: 3.7",
         "Programming Language :: Python :: 3.8",
         "Programming Language :: Python :: 3.9",
@@ -71,5 +70,5 @@ def get_version():
     url="https://github.com/corteva/rioxarray",
     version=get_version(),
     zip_safe=False,
-    python_requires=">=3.6",
+    python_requires=">=3.7",
 )
diff --git a/test/conftest.py b/test/conftest.py
index d66d338f..3d48568f 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -23,6 +23,7 @@ def _assert_attrs_equal(input_xr, compare_xr, decimal_precision):
             attr != "_FillValue"
             and attr not in UNWANTED_RIO_ATTRS
             and attr != "creation_date"
+            and attr != "grid_mapping"
         ):
             try:
                 assert_almost_equal(
@@ -80,10 +81,6 @@ def _assert_xarrays_equal(
             "_FillValue", input_xarray.encoding.get("_FillValue")
         )
         assert_array_equal([input_fill_value], [compare_fill_value])
-        assert "grid_mapping" in compare_xarray.attrs
-        assert (
-            input_xarray[input_xarray.attrs["grid_mapping"]]
-            == compare_xarray[compare_xarray.attrs["grid_mapping"]]
-        )
+        assert input_xarray.rio.grid_mapping == compare_xarray.rio.grid_mapping
         for unwanted_attr in UNWANTED_RIO_ATTRS:
             assert unwanted_attr not in input_xarray.attrs
diff --git a/test/integration/test_integration__io.py b/test/integration/test_integration__io.py
index 7dd6757f..43c457fc 100644
--- a/test/integration/test_integration__io.py
+++ b/test/integration/test_integration__io.py
@@ -241,7 +241,6 @@ def test_open_group_load_attrs():
         assert sorted(attrs) == [
             "_FillValue",
             "add_offset",
-            "grid_mapping",
             "long_name",
             "scale_factor",
             "units",
@@ -249,7 +248,7 @@ def test_open_group_load_attrs():
         assert attrs["long_name"] == "500m Surface Reflectance Band 5 - first layer"
         assert attrs["units"] == "reflectance"
         assert attrs["_FillValue"] == -28672.0
-        assert attrs["grid_mapping"] == "spatial_ref"
+        assert rds["sur_refl_b05_1"].encoding["grid_mapping"] == "spatial_ref"
 
 
 def test_open_rasterio_mask_chunk_clip():
@@ -267,14 +266,13 @@ def test_open_rasterio_mask_chunk_clip():
         assert np.isnan(xdi.values).sum() == 52119
         test_encoding = dict(xdi.encoding)
         assert test_encoding.pop("source").endswith("small_dem_3m_merged.tif")
-        assert test_encoding == {"_FillValue": 0.0}
+        assert test_encoding == {"_FillValue": 0.0, "grid_mapping": "spatial_ref"}
         attrs = dict(xdi.attrs)
         assert_almost_equal(
             tuple(xdi.rio._cached_transform())[:6],
             (3.0, 0.0, 425047.68381405267, 0.0, -3.0, 4615780.040546387),
         )
         assert attrs == {
-            "grid_mapping": "spatial_ref",
             "add_offset": 0.0,
             "scale_factor": 1.0,
         }
@@ -305,7 +303,7 @@ def test_open_rasterio_mask_chunk_clip():
         _assert_xarrays_equal(clipped, comp_subset)
         test_encoding = dict(clipped.encoding)
         assert test_encoding.pop("source").endswith("small_dem_3m_merged.tif")
-        assert test_encoding == {"_FillValue": 0.0}
+        assert test_encoding == {"_FillValue": 0.0, "grid_mapping": "spatial_ref"}
 
         # test dataset
         clipped_ds = xdi.to_dataset(name="test_data").rio.clip(
@@ -315,7 +313,7 @@ def test_open_rasterio_mask_chunk_clip():
         _assert_xarrays_equal(clipped_ds, comp_subset_ds)
         test_encoding = dict(clipped.encoding)
         assert test_encoding.pop("source").endswith("small_dem_3m_merged.tif")
-        assert test_encoding == {"_FillValue": 0.0}
+        assert test_encoding == {"_FillValue": 0.0, "grid_mapping": "spatial_ref"}
 
 
 ##############################################################################
@@ -895,6 +893,7 @@ def test_mask_and_scale():
                 "scale_factor": 0.1,
                 "_FillValue": 32767.0,
                 "missing_value": 32767,
+                "grid_mapping": "crs",
             }
             attrs = rds.air_temperature.attrs
             assert attrs == {
@@ -902,7 +901,6 @@ def test_mask_and_scale():
                 "coordinate_system": "WGS84,EPSG:4326",
                 "description": "Daily Maximum Temperature",
                 "dimensions": "lon lat time",
-                "grid_mapping": "crs",
                 "long_name": "tmmx",
                 "standard_name": "tmmx",
                 "units": "K",
@@ -923,6 +921,7 @@ def test_no_mask_and_scale():
         assert test_encoding == {
             "_FillValue": 32767.0,
             "missing_value": 32767,
+            "grid_mapping": "crs",
         }
         attrs = rds.air_temperature.attrs
         assert attrs == {
@@ -932,7 +931,6 @@ def test_no_mask_and_scale():
             "coordinate_system": "WGS84,EPSG:4326",
             "description": "Daily Maximum Temperature",
             "dimensions": "lon lat time",
-            "grid_mapping": "crs",
             "long_name": "tmmx",
             "scale_factor": 0.1,
             "standard_name": "tmmx",
diff --git a/test/integration/test_integration_merge.py b/test/integration/test_integration_merge.py
index da10fd3f..b0e98015 100644
--- a/test/integration/test_integration_merge.py
+++ b/test/integration/test_integration_merge.py
@@ -15,7 +15,6 @@ def test_merge_arrays(squeeze):
     with open_rasterio(dem_test) as rds:
         rds.attrs = {
             "_FillValue": rds.rio.nodata,
-            "grid_mapping": "spatial_ref",
             "crs": rds.attrs["crs"],
         }
         arrays = [
@@ -52,6 +51,7 @@ def test_merge_arrays(squeeze):
     assert sorted(merged.coords) == ["band", "spatial_ref", "x", "y"]
     assert merged.rio.crs == rds.rio.crs
     assert merged.attrs == rds.attrs
+    assert merged.encoding["grid_mapping"] == "spatial_ref"
     assert_almost_equal(merged.sum(), 11368261)
 
 
@@ -95,16 +95,15 @@ def test_merge__different_crs(dataset):
     assert sorted(merged.coords) == ["band", "spatial_ref", "x", "y"]
     assert merged.rio.crs == rds.rio.crs
     if dataset:
-        assert merged.attrs == {"grid_mapping": "spatial_ref"}
         assert_almost_equal(merged[merged.rio.vars[0]].sum(), -131013894)
     else:
         assert merged.attrs == {
             "_FillValue": -28672,
             "add_offset": 0.0,
-            "grid_mapping": "spatial_ref",
             "scale_factor": 1.0,
         }
         assert_almost_equal(merged.sum(), -131013894)
+    assert merged.encoding["grid_mapping"] == "spatial_ref"
 
 
 def test_merge_arrays__res():
@@ -112,7 +111,6 @@ def test_merge_arrays__res():
     with open_rasterio(dem_test, masked=True) as rds:
         rds.attrs = {
             "_FillValue": rds.rio.nodata,
-            "grid_mapping": "spatial_ref",
             "crs": rds.attrs["crs"],
         }
         arrays = [
@@ -140,6 +138,7 @@ def test_merge_arrays__res():
     compare_attrs = dict(rds.attrs)
     compare_attrs.pop("crs")
     assert merged.attrs == compare_attrs
+    assert merged.encoding["grid_mapping"] == "spatial_ref"
     assert_almost_equal(nansum(merged), 13760565)
 
 
@@ -193,9 +192,8 @@ def test_merge_datasets():
     assert merged.coords["band"].values == [1]
     assert sorted(merged.coords) == ["band", "spatial_ref", "x", "y"]
     assert merged.rio.crs == rds.rio.crs
-    base_attrs = dict(rds.attrs)
-    base_attrs["grid_mapping"] = "spatial_ref"
-    assert merged.attrs == base_attrs
+    assert merged.attrs == rds.attrs
+    assert merged.encoding["grid_mapping"] == "spatial_ref"
     assert_almost_equal(merged[data_var].sum(), 4543446965182987)
 
 
@@ -239,7 +237,6 @@ def test_merge_datasets__res():
     assert merged.coords["band"].values == [1]
     assert sorted(merged.coords) == ["band", "spatial_ref", "x", "y"]
     assert merged.rio.crs == rds.rio.crs
-    base_attrs = dict(rds.attrs)
-    base_attrs["grid_mapping"] = "spatial_ref"
-    assert merged.attrs == base_attrs
+    assert merged.attrs == rds.attrs
+    assert merged.encoding["grid_mapping"] == "spatial_ref"
     assert_almost_equal(merged[data_var].sum(), 974566547463955)
diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py
index b0da9a32..5563ed0e 100644
--- a/test/integration/test_integration_rioxarray.py
+++ b/test/integration/test_integration_rioxarray.py
@@ -1485,7 +1485,7 @@ def test_crs_writer__array__copy():
     test_da.rio._crs = None
     assert test_da.rio.crs is None
     assert "crs" not in test_da.coords
-    assert out_da.attrs["grid_mapping"] == "crs"
+    assert out_da.encoding["grid_mapping"] == "crs"
 
 
 def test_crs_writer__array__inplace():
@@ -1501,7 +1501,7 @@ def test_crs_writer__array__inplace():
     assert out_da.coords["spatial_ref"] == test_da.coords["spatial_ref"]
     test_da.rio._crs = None
     assert test_da.rio.crs.to_epsg() == 4326
-    assert test_da.attrs["grid_mapping"] == "spatial_ref"
+    assert test_da.encoding["grid_mapping"] == "spatial_ref"
     assert out_da.attrs == test_da.attrs
     out_da.rio._crs = None
     assert out_da.rio.crs.to_epsg() == 4326
@@ -1520,7 +1520,7 @@ def test_crs_writer__dataset__copy():
     assert "spatial_ref" in out_da.coords["crs"].attrs
     out_da.test.rio._crs = None
     assert out_da.rio.crs.to_epsg() == 4326
-    assert out_da.test.attrs["grid_mapping"] == "crs"
+    assert out_da.test.encoding["grid_mapping"] == "crs"
     # make sure input did not change the dataset
     test_da.test.rio._crs = None
     test_da.rio._crs = None
@@ -1545,7 +1545,7 @@ def test_crs_writer__dataset__inplace():
     test_da.test.rio._crs = None
     test_da.rio._crs = None
     assert test_da.rio.crs.to_epsg() == 4326
-    assert out_da.test.attrs["grid_mapping"] == "spatial_ref"
+    assert out_da.test.encoding["grid_mapping"] == "spatial_ref"
     assert out_da.test.attrs == test_da.test.attrs
 
 
@@ -1648,14 +1648,14 @@ def test_crs_get_custom():
 def test_get_crs_dataset():
     test_ds = xarray.Dataset()
     test_ds = test_ds.rio.write_crs(4326)
-    assert test_ds.attrs["grid_mapping"] == "spatial_ref"
+    assert test_ds.encoding["grid_mapping"] == "spatial_ref"
     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.encoding["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
@@ -1666,7 +1666,7 @@ def test_write_crs_cf__disable_grid_mapping():
     test_da = xarray.DataArray(1)
     with rioxarray.set_options(export_grid_mapping=False):
         test_da = test_da.rio.write_crs(4326)
-    assert test_da.attrs["grid_mapping"] == "spatial_ref"
+    assert test_da.encoding["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
@@ -1687,7 +1687,7 @@ def test_write_crs__missing_geospatial_dims():
 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"
+    assert test_da.encoding["grid_mapping"] == "spatial_ref"
     attrs = test_da.spatial_ref.attrs
     attrs.pop("spatial_ref")
     attrs.pop("crs_wkt")
@@ -1697,7 +1697,7 @@ def test_read_crs_cf():
 def test_get_crs_dataset__nonstandard_grid_mapping():
     test_ds = xarray.Dataset()
     test_ds = test_ds.rio.write_crs(4326, grid_mapping_name="frank")
-    assert test_ds.attrs["grid_mapping"] == "frank"
+    assert test_ds.encoding["grid_mapping"] == "frank"
     assert test_ds.rio.crs.to_epsg() == 4326
 
 
@@ -1839,7 +1839,7 @@ def test_isel_window():
 def test_write_pyproj_crs_dataset():
     test_ds = xarray.Dataset()
     test_ds = test_ds.rio.write_crs(pCRS(4326))
-    assert test_ds.attrs["grid_mapping"] == "spatial_ref"
+    assert test_ds.encoding["grid_mapping"] == "spatial_ref"
     assert test_ds.rio.crs.to_epsg() == 4326
 
 
@@ -2112,12 +2112,12 @@ def test_write_transform():
     ds.rio.write_transform(test_affine, inplace=True)
     assert ds.spatial_ref.GeoTransform == "425047.0 3.0 0.0 4615780.0 0.0 -3.0"
     assert ds.rio._cached_transform() == test_affine
-    assert ds.grid_mapping == "spatial_ref"
+    assert ds.rio.grid_mapping == "spatial_ref"
     da = xarray.DataArray(1)
     da.rio.write_transform(test_affine, inplace=True)
     assert da.rio._cached_transform() == test_affine
     assert da.spatial_ref.GeoTransform == "425047.0 3.0 0.0 4615780.0 0.0 -3.0"
-    assert da.grid_mapping == "spatial_ref"
+    assert da.rio.grid_mapping == "spatial_ref"
 
 
 def test_missing_transform():