Repository: corteva/rioxarray

Issue:
rioxarray export to_raster and CRS being dropped - user error or bug?
hi there @snowman2 
I have a question about CRS persistence when you export a raster using rioxarray. It seems like when we begin to squeeze dimensions from the data, that the CRS gets dropped. I tried to reset the CRS before exporting to a geotiff but when i import the data the crs is not there as far as i can see.

Am i doing something incorrectly in the approach below? The CRS is that climate data global longitude values (no negatives). i'm not sure what the EPSG is for that one but we have the WKT definition from the original data. 

Any thoughts on what is incorrect? OR how i can ensure a CRS is set when exporting the data?

```
data_path = http://thredds.northwestknowledge.net:8080/thredds/dodsC/agg_macav2metdata_pr_BNU-ESM_r1i1p1_rcp45_2006_2099_CONUS_monthly.nc

# Open the data
with xr.open_dataset(full_file_path) as file_nc:
    climate_data_xr = file_nc

# CRS -- wkt format returned
climate_data_xr.rio.crs

# natural earth data - https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/cultural/ne_50m_admin_1_states_provinces_lakes.zip
# Open the natural earth layer with geopandas
states_gdf = gpd.read_file(states_path)
# subset the data
cali_aoi = states_gdf[states_gdf.name == "California"]

# Get lat min, max for cali subset & convert to 0-360
aoi_lat = [float(cali_aoi.total_bounds[1]), float(cali_aoi.total_bounds[3])]
aoi_lon = [float(cali_aoi.total_bounds[0]), float(cali_aoi.total_bounds[2])]
aoi_lon[0] = aoi_lon[0] + 360
aoi_lon[1] = aoi_lon[1] + 360

# i'm using region mask to mask the data

state_mask = regionmask.mask_3D_geopandas(cali_aoi,
                                          climate_data_xr.lon,
                                          climate_data_xr.lat)

# Slice the data by time and spatial extent
start_date = "2010-01-15"
end_date = "2010-02-15"

# Subset the data by time and aoi & mask the data
climate_data_subset = climate_data_xr[climate_variable].sel(
    time=slice(start_date, end_date),
    lon=slice(aoi_lon[0], aoi_lon[1]),
    lat=slice(aoi_lat[0], aoi_lat[1])).where(state_mask)

#  This drops the 4th dimension which is region
new_masked_data = masked_data[:, :, :, 0]
# Reset the CRS given it gets dropped when you squeeze the data. 
new_masked_data.rio.set_crs(climate_crs)
# Save as new geotiff 
new_masked_data.rio.to_raster('subset_climate_raster_data.tif')

# reopen the data
saved_raster = xr.open_rasterio('subset_climate_raster_data.tif')
# The crs is empty now
saved_raster.rio.crs

```


Discussion:

Patch:
diff --git a/docs/history.rst b/docs/history.rst
index 7937a3d1..ac181108 100644
--- a/docs/history.rst
+++ b/docs/history.rst
@@ -3,7 +3,7 @@ History
 
 Latest
 ------
-
+- BUG: Check all CRS are the same in the dataset in crs() method
 
 0.1.0
 ------
diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py
index acbc922c..648fc6ff 100644
--- a/rioxarray/rioxarray.py
+++ b/rioxarray/rioxarray.py
@@ -1749,14 +1749,24 @@ class RasterDataset(XRasterBase):
         self._crs = super().crs
         if self._crs is not None:
             return self._crs
+        # ensure all the CRS of the variables are the same
+        crs_list = []
         for var in self.vars:
-            crs = self._obj[var].rio.crs
-            if crs is not None:
-                self._crs = crs
-                break
-        else:
+            if self._obj[var].rio.crs is not None:
+                crs_list.append(self._obj[var].rio.crs)
+        try:
+            crs = crs_list[0]
+        except IndexError:
+            crs = None
+        if crs is None:
             self._crs = False
             return None
+        elif all(crs_i == crs for crs_i in crs_list):
+            self._crs = crs
+        else:
+            raise RioXarrayError(
+                "CRS in DataArrays differ in the Dataset: {}".format(crs_list)
+            )
         return self._crs
 
     def reproject(
@@ -2060,12 +2070,10 @@ class RasterDataset(XRasterBase):
         scales = []
         offsets = []
         nodatavals = []
-        crs_list = []
         for data_var in data_array[variable_dim].values:
             scales.append(self._obj[data_var].attrs.get("scale_factor", 1.0))
             offsets.append(self._obj[data_var].attrs.get("add_offset", 0.0))
             nodatavals.append(self._obj[data_var].rio.nodata)
-            crs_list.append(self._obj[data_var].rio.crs)
         data_array.attrs["scales"] = scales
         data_array.attrs["offsets"] = offsets
         nodata = nodatavals[0]
@@ -2079,14 +2087,8 @@ class RasterDataset(XRasterBase):
                 "All nodata values must be the same when exporting to raster. "
                 "Current values: {}".format(nodatavals)
             )
-        crs = crs_list[0]
-        if all(crs_i == crs for crs_i in crs_list):
-            data_array.rio.write_crs(crs, inplace=True)
-        else:
-            raise RioXarrayError(
-                "All CRS must be the same when exporting to raster. "
-                "Current values: {}".format(crs_list)
-            )
+        if self.crs is not None:
+            data_array.rio.write_crs(self.crs, inplace=True)
         # write it to a raster
         data_array.rio.to_raster(
             raster_path=raster_path,


Test patch:
diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py
index 4c7bc28f..8492b388 100644
--- a/test/integration/test_integration_rioxarray.py
+++ b/test/integration/test_integration_rioxarray.py
@@ -1251,7 +1251,7 @@ def test_to_raster__dataset__different_crs(tmpdir):
         rds.blue.attrs = attrs
         rds = rds.drop_vars("spatial_ref")
         with pytest.raises(
-            RioXarrayError, match="All CRS must be the same when exporting to raster."
+            RioXarrayError, match="CRS in DataArrays differ in the Dataset"
         ):
             rds.rio.to_raster(str(tmp_raster))
 
@@ -1486,6 +1486,20 @@ def test_crs_writer__missing():
         test_da.to_dataset(name="test").rio.write_crs()
 
 
+def test_crs__dataset__different_crs(tmpdir):
+    green = xarray.DataArray(
+        numpy.zeros((5, 5)),
+        dims=("y", "x"),
+        coords={"y": numpy.arange(1, 6), "x": numpy.arange(2, 7)},
+        attrs={"crs": "EPSG:4326"},
+    )
+    blue = green.copy(deep=True)
+    blue.attrs = {"crs": "EPSG:32722"}
+
+    with pytest.raises(RioXarrayError, match="CRS in DataArrays differ in the Dataset"):
+        xarray.Dataset({"green": green, "blue": blue}).rio.crs
+
+
 def test_clip_missing_crs():
     test_da = xarray.DataArray(
         numpy.zeros((5, 5)),