Add CRS to the spatial dimensions.

datacube/api/core.py

@@ -443,16 +443,18 @@ def empty_func(measurement_):
 
         data_func = data_func or empty_func
 
-        result = xarray.Dataset(attrs={'crs': geobox.crs})
+        result = xarray.Dataset(attrs={'crs': str(geobox.crs)})
         for name, coord in coords.items():
             result[name] = coord
         for name, coord in geobox.coordinates.items():
-            result[name] = (name, coord.values, {'units': coord.units, 'resolution': coord.resolution})
+            result[name] = (name, coord.values, {'units': coord.units,
+                                                 'resolution': coord.resolution,
+                                                 'crs': result.crs})
 
         for measurement in measurements:
             data = data_func(measurement)
             attrs = measurement.dataarray_attrs()
-            attrs['crs'] = geobox.crs
+            attrs['crs'] = result.crs
             dims = tuple(coords.keys()) + tuple(geobox.dimensions)
             result[measurement.name] = (dims, data, attrs)
 

datacube/drivers/netcdf/_write.py

@@ -79,11 +79,14 @@ def write_dataset_to_netcdf(dataset, filename, global_attributes=None, variable_
     if not dataset.data_vars.keys():
         raise DatacubeException('Cannot save empty dataset to disk.')
 
-    if not hasattr(dataset, 'crs'):
-        raise DatacubeException('Dataset does not contain CRS, cannot write to NetCDF file.')
+    if dataset.geobox is None:
+        raise DatacubeException('Dataset geobox property is None, cannot write to NetCDF file.')
+
+    if dataset.geobox.crs is None:
+        raise DatacubeException('Dataset geobox.crs property is None, cannot write to NetCDF file.')
 
     nco = create_netcdf_storage_unit(filename,
-                                     dataset.crs,
+                                     dataset.geobox.crs,
                                      dataset.coords,
                                      dataset.data_vars,
                                      variable_params,

datacube/utils/xarray_geoextensions.py

@@ -25,8 +25,34 @@ def _norm_crs(crs):
         raise ValueError('Can not interpret {} as CRS'.format(type(crs)))
 
 
+def _get_crs(obj):
+    if not isinstance(obj, xarray.Dataset) and not isinstance(obj, xarray.DataArray):
+        raise ValueError('Can not get crs from {}'.format(type(obj)))
+
+    if isinstance(obj, xarray.Dataset):
+        if len(obj.data_vars) > 0:
+            data_array = next(iter(obj.data_vars.values()))
+        else:
+            # fall back option
+            return obj.attrs.get('crs', None)
+    else:
+        data_array = obj
+
+    # Assumption: spatial dimensions are always the last 2)
+    spatial_dims = data_array.dims[-2:]
+    crs_set = set(data_array[d].attrs.get('crs', None) for d in spatial_dims)
+    if len(crs_set) > 1:
+        raise ValueError('Spatial dimensions have different crs.')
+    elif len(crs_set) == 1 and None not in crs_set:
+        crs = data_array[data_array.dims[-1]].crs
+    else:
+        # fall back option
+        crs = data_array.attrs.get('crs', None) or obj.attrs.get('crs', None)
+    return crs
+
+
 def _xarray_affine(obj):
-    crs = _norm_crs(obj.crs)
+    crs = _norm_crs(_get_crs(obj))
     if crs is None:
         return None
 
@@ -51,7 +77,7 @@ def _xarray_extent(obj):
 
 
 def _xarray_geobox(obj):
-    crs = _norm_crs(obj.crs)
+    crs = _norm_crs(_get_crs(obj))
     if crs is None:
         return None
 

datacube_apps/stacker/fixer.py

@@ -181,8 +181,14 @@ def do_fixer_task(config, task):
     data['dataset'] = datasets_to_doc(unwrapped_datasets)
 
     try:
+        if data.geobox is None:
+            raise DatacubeException('Dataset geobox property is None, cannot write to NetCDF file.')
+
+        if data.geobox.crs is None:
+            raise DatacubeException('Dataset geobox.crs property is None, cannot write to NetCDF file.')
+
         nco = create_netcdf_storage_unit(temp_filename,
-                                         data.crs,
+                                         data.geobox.crs,
                                          data.coords,
                                          data.data_vars,
                                          variable_params,

datacube_apps/stacker/stacker.py

@@ -158,8 +158,14 @@ def do_stack_task(config, task):
     data['dataset'] = datasets_to_doc(unwrapped_datasets)
 
     try:
+        if data.geobox is None:
+            raise DatacubeException('Dataset geobox property is None, cannot write to NetCDF file.')
+
+        if data.geobox.crs is None:
+            raise DatacubeException('Dataset geobox.crs property is None, cannot write to NetCDF file.')
+
         nco = create_netcdf_storage_unit(temp_filename,
-                                         data.crs,
+                                         data.geobox.crs,
                                          data.coords,
                                          data.data_vars,
                                          variable_params,

tests/storage/test_netcdfwriter.py

@@ -284,4 +284,4 @@ def test_useful_error_on_write_empty_dataset(tmpnetcdf_filename):
     with pytest.raises(DatacubeException) as excinfo:
         ds = xr.Dataset(data_vars={'blue': (('time',), numpy.array([0, 1, 2]))})
         write_dataset_to_netcdf(ds, tmpnetcdf_filename)
-    assert 'CRS' in str(excinfo.value)
+    assert 'geobox' in str(excinfo.value)

tests/test_utils.py

@@ -259,6 +259,9 @@ def test_write_geotiff_str_crs(tmpdir, odc_style_xr_dataset):
         assert (written_data == odc_style_xr_dataset['B10']).all()
 
     del odc_style_xr_dataset.attrs['crs']
+    del odc_style_xr_dataset.B10.attrs['crs']
+    for dim in odc_style_xr_dataset.B10.dims:
+        del odc_style_xr_dataset[dim].attrs['crs']
     with pytest.raises(ValueError):
         write_geotiff(filename, odc_style_xr_dataset)
 

tests/test_utils_cog.py

@@ -147,6 +147,8 @@ def test_cog_no_crs(tmpdir, with_dask):
 
     xx, ds = gen_test_data(pp, dask=with_dask)
     del xx.attrs['crs']
+    for dim in xx.dims:
+        del xx[dim].attrs['crs']
 
     with pytest.raises(ValueError):
         write_cog(xx, ":mem:")

tests/test_xarray_extension.py

@@ -22,6 +22,9 @@ def test_xr_extension(odc_style_xr_dataset):
     assert (zz0, zz1) == (0, 0)
 
     odc_style_xr_dataset.attrs['crs'] = None
+    odc_style_xr_dataset.B10.attrs['crs'] = None
+    for dim in odc_style_xr_dataset.B10.dims:
+        odc_style_xr_dataset[dim].attrs['crs'] = None
     assert _xarray_affine(odc_style_xr_dataset) is None
     assert _xarray_geobox(odc_style_xr_dataset) is None
     assert _xarray_extent(odc_style_xr_dataset) is None