diff --git a/RELEASE_NOTES.rst b/RELEASE_NOTES.rst
index 06dde4df..1583d28c 100755
--- a/RELEASE_NOTES.rst
+++ b/RELEASE_NOTES.rst
@@ -17,8 +17,15 @@ Upcoming Release
    To use the features already you have to install the ``master`` branch, e.g. 
    ``pip install git+https://github.com/pypsa/atlite``.
 
+`v0.6.1 <https://github.com/PyPSA/atlite/releases/tag/v0.6.1>`__ (21st April 2026)
+=======================================================================================
+
+**Bug fixes**
+
+* Fix backwards compatibility of ``aggregate_matrix``.
+
 `v0.6.0 <https://github.com/PyPSA/atlite/releases/tag/v0.6.0>`__ (15th April 2026)
-=======================================================================================   
+=======================================================================================
 
 **Features**
 
diff --git a/atlite/aggregate.py b/atlite/aggregate.py
index 44464c2e..eb100260 100644
--- a/atlite/aggregate.py
+++ b/atlite/aggregate.py
@@ -5,14 +5,19 @@
 Functions for aggregating results.
 """
 
+import pandas as pd
 import scipy.sparse as sp
 import xarray as xr
 from dask.array.core import Array
 
+from atlite.utils import ensure_coords
+
 
 def aggregate_matrix(
-    da: xr.DataArray, matrix: sp.csr_matrix, coords: xr.Coordinates
+    da: xr.DataArray, matrix: sp.csr_matrix, index: xr.Coordinates | pd.Index
 ) -> xr.DataArray:
+    coords = ensure_coords(index)
+
     if isinstance(da.data, Array):
         da = da.stack(spatial=("y", "x"))
         da = da.chunk(dict(spatial=-1))
diff --git a/atlite/convert.py b/atlite/convert.py
index fab5628e..fb626e9b 100644
--- a/atlite/convert.py
+++ b/atlite/convert.py
@@ -237,7 +237,7 @@ def convert_and_aggregate(
         if isinstance(shapes, geoseries_like) and index is None:
             index = shapes.index
 
-        matrix = cutout.indicatormatrix(shapes, shapes_crs)
+        matrix = cutout.indicatormatrix(shapes, shapes_crs).tocsr()
 
     if layout is not None:
         assert isinstance(layout, xr.DataArray)
@@ -254,7 +254,7 @@ def convert_and_aggregate(
     coords = ensure_coords(pd.RangeIndex(matrix.shape[0]) if index is None else index)
     if len(coords.dims) > 1:
         raise ValueError(f"index must have a single dimension, not: {coords.dims}")
-    results = aggregate_matrix(da, matrix=matrix, coords=coords)
+    results = aggregate_matrix(da, matrix=matrix, index=coords)
 
     if per_unit or return_capacity:
         caps = matrix.sum(-1)
diff --git a/atlite/gis.py b/atlite/gis.py
index 6b7edb56..0d76eef8 100644
--- a/atlite/gis.py
+++ b/atlite/gis.py
@@ -479,7 +479,9 @@ def open_files(self):
                 assert isinstance(raster, rio.DatasetReader)
 
             # Check if the raster has a valid CRS
-            if not (raster.crs.is_geographic or raster.crs.is_projected):
+            if raster.crs is None or not (
+                raster.crs.is_geographic or raster.crs.is_projected
+            ):
                 if d["crs"]:
                     raster._crs = CRS(d["crs"])
                 else: