diff --git a/docs/history.rst b/docs/history.rst
index aa417509..1f04ec04 100644
--- a/docs/history.rst
+++ b/docs/history.rst
@@ -4,6 +4,10 @@ History
 Latest
 -------
 
+0.10.3
+------
+- BUG: Remove xarray crs attribute in rio.write_crs (issue #488)
+
 0.10.2
 -------
 - BUG: Lazy load colormap through _manager.acquire() in merge (issue #479)
diff --git a/rioxarray/_io.py b/rioxarray/_io.py
index 61199c4e..f4f2d336 100644
--- a/rioxarray/_io.py
+++ b/rioxarray/_io.py
@@ -937,8 +937,9 @@ def open_rasterio(
     # always (0, 0, 1) per definition (see
     # https://github.com/sgillies/affine)
     result.rio.write_transform(_rio_transform(riods), inplace=True)
-    if riods.crs:
-        result.rio.write_crs(riods.crs, inplace=True)
+    rio_crs = riods.crs or result.rio.crs
+    if rio_crs:
+        result.rio.write_crs(rio_crs, inplace=True)
     if has_gcps:
         result.rio.write_gcps(*riods.gcps, inplace=True)
 
diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py
index 70aceb83..e21c4e39 100644
--- a/rioxarray/rioxarray.py
+++ b/rioxarray/rioxarray.py
@@ -447,6 +447,9 @@ def write_crs(self, input_crs=None, grid_mapping_name=None, inplace=False):
             )
         data_obj.coords[grid_mapping_name].rio.set_attrs(grid_map_attrs, inplace=True)
 
+        # remove old crs if exists
+        data_obj.attrs.pop("crs", None)
+
         return data_obj.rio.write_grid_mapping(
             grid_mapping_name=grid_mapping_name, inplace=True
         )
diff --git a/test/integration/test_integration_merge.py b/test/integration/test_integration_merge.py
index 31cf450b..1f0a225a 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,
-            "crs": rds.attrs["crs"],
         }
         arrays = [
             rds.isel(x=slice(100), y=slice(100)),
@@ -140,7 +139,6 @@ def test_merge_arrays__res():
     with open_rasterio(dem_test, masked=True) as rds:
         rds.attrs = {
             "_FillValue": rds.rio.nodata,
-            "crs": rds.attrs["crs"],
         }
         arrays = [
             rds.isel(x=slice(100), y=slice(100)),
@@ -173,7 +171,6 @@ def test_merge_arrays__res():
     assert merged.rio.crs == rds.rio.crs
     assert_almost_equal(merged.attrs.pop("_FillValue"), rds.attrs.pop("_FillValue"))
     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)
diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py
index 6bc19831..31923326 100644
--- a/test/integration/test_integration_rioxarray.py
+++ b/test/integration/test_integration_rioxarray.py
@@ -892,18 +892,22 @@ def test_reproject__grid_mapping(modis_reproject):
     ) as mda, modis_reproject["open"](modis_reproject["compare"], **mask_args) as mdc:
 
         # remove 'crs' attribute and add grid mapping
+        # keep a copy of the crs before setting the spatial_ref to 0 as that will
+        # set `rio.crs` to None.
+        mda_rio_crs = mda.rio.crs
         mda.coords["spatial_ref"] = 0
         mda.coords["spatial_ref"].attrs["spatial_ref"] = CRS.from_user_input(
-            _get_attr(mda, "crs")
+            mda_rio_crs
         ).wkt
         _mod_attr(mda, "grid_mapping", val="spatial_ref")
-        _del_attr(mda, "crs")
+        # keep a copy of the crs before setting the spatial_ref to 0 as that will
+        # set `rio.crs` to None.
+        mdc_rio_crs = mdc.rio.crs
         mdc.coords["spatial_ref"] = 0
         mdc.coords["spatial_ref"].attrs["spatial_ref"] = CRS.from_user_input(
-            _get_attr(mdc, "crs")
+            mdc_rio_crs
         ).wkt
         _mod_attr(mdc, "grid_mapping", val="spatial_ref")
-        _del_attr(mdc, "crs")
 
         # reproject
         mds_repr = mda.rio.reproject(modis_reproject["to_proj"])
@@ -2025,6 +2029,13 @@ def test_get_crs_dataset():
     assert test_ds.rio.crs.to_epsg() == 4326
 
 
+def test_crs_is_removed():
+    test_ds = xarray.Dataset(attrs=dict(crs="+init=epsg:4326"))
+    test_ds = test_ds.rio.write_crs(4326)
+
+    assert "crs" not in test_ds.attrs
+
+
 def test_write_crs_cf():
     test_da = xarray.DataArray(1)
     test_da = test_da.rio.write_crs(4326)