Repository: scverse/spatialdata

Issue:
`rasterize_bins()` problem with coordinate transformations and with polygons.
I have noticed that `rasterize_bins()` computes the affine matrices on the intrinsic data against the row/col data, ignoring an eventual coordinate transformation. Also it fails when the bins are represented a polygons/multipolygons. I am working on a fix on this, which will solve both problems by calling `get_centroids()`.

- [x] I need to check that performance don't degrade. An ad-hoc solution would be preferred in that case.

CC @quentinblampey 

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/scverse/spatialdata/pull/651?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=scverse) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 89.96%. Comparing base [(`a3e62b6`)](https://app.codecov.io/gh/scverse/spatialdata/commit/a3e62b6a832ac09844707586f42103af081b1dfa?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=scverse) to head [(`36aa2ad`)](https://app.codecov.io/gh/scverse/spatialdata/commit/36aa2add544ca8e65b2a45380a9da5ed0c2134e6?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=scverse).

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main     #651      +/-   ##
==========================================
+ Coverage   89.94%   89.96%   +0.01%     
==========================================
  Files          44       44              
  Lines        6755     6768      +13     
==========================================
+ Hits         6076     6089      +13     
  Misses        679      679              
```

| [Files](https://app.codecov.io/gh/scverse/spatialdata/pull/651?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=scverse) | Coverage Δ | |
|---|---|---|
| [src/spatialdata/\_core/operations/rasterize\_bins.py](https://app.codecov.io/gh/scverse/spatialdata/pull/651?src=pr&el=tree&filepath=src%2Fspatialdata%2F_core%2Foperations%2Frasterize_bins.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=scverse#diff-c3JjL3NwYXRpYWxkYXRhL19jb3JlL29wZXJhdGlvbnMvcmFzdGVyaXplX2JpbnMucHk=) | `89.77% <100.00%> (+0.61%)` | :arrow_up: |
| [src/spatialdata/models/models.py](https://app.codecov.io/gh/scverse/spatialdata/pull/651?src=pr&el=tree&filepath=src%2Fspatialdata%2Fmodels%2Fmodels.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=scverse#diff-c3JjL3NwYXRpYWxkYXRhL21vZGVscy9tb2RlbHMucHk=) | `87.88% <100.00%> (+0.21%)` | :arrow_up: |

</details>

Patch:
diff --git a/src/spatialdata/_core/operations/rasterize_bins.py b/src/spatialdata/_core/operations/rasterize_bins.py
index 2e26f5e..fcd58b3 100644
--- a/src/spatialdata/_core/operations/rasterize_bins.py
+++ b/src/spatialdata/_core/operations/rasterize_bins.py
@@ -9,6 +9,7 @@ from dask.dataframe import DataFrame as DaskDataFrame
 from geopandas import GeoDataFrame
 from numpy.random import default_rng
 from scipy.sparse import csc_matrix
+from shapely import MultiPolygon, Point, Polygon
 from skimage.transform import estimate_transform
 from xarray import DataArray
 
@@ -155,8 +156,13 @@ def rasterize_bins(
 
     src = np.stack([sub_table.obs[col_key] - min_col, sub_table.obs[row_key] - min_row], axis=1)
     if isinstance(sub_df, GeoDataFrame):
-        sub_x = sub_df.geometry.x.values
-        sub_y = sub_df.geometry.y.values
+        if isinstance(sub_df.iloc[0].geometry, Point):
+            sub_x = sub_df.geometry.x.values
+            sub_y = sub_df.geometry.y.values
+        else:
+            assert isinstance(sub_df.iloc[0].geometry, (Polygon, MultiPolygon))
+            sub_x = sub_df.centroid.x
+            sub_y = sub_df.centroid.y
     else:
         assert isinstance(sub_df, DaskDataFrame)
         sub_x = sub_df.x.compute().values
diff --git a/src/spatialdata/models/models.py b/src/spatialdata/models/models.py
index 4314ba9..39911eb 100644
--- a/src/spatialdata/models/models.py
+++ b/src/spatialdata/models/models.py
@@ -647,14 +647,28 @@ class PointsModel:
                 )
         ndim = len(coordinates)
         axes = [X, Y, Z][:ndim]
+        index_monotonically_increasing = data.index.is_monotonic_increasing
+        if not isinstance(index_monotonically_increasing, bool):
+            index_monotonically_increasing = index_monotonically_increasing.compute()
         if isinstance(data, pd.DataFrame):
+            if not index_monotonically_increasing:
+                warnings.warn(
+                    "The index of the dataframe is not monotonic increasing. It is recommended to sort the data to "
+                    "adjust the order of the index before calling .parse(); this will make the division known when"
+                    "dask.dataframe.from_pandas() is called.",
+                    UserWarning,
+                    stacklevel=2,
+                )
             table: DaskDataFrame = dd.from_pandas(  # type: ignore[attr-defined]
                 pd.DataFrame(data[[coordinates[ax] for ax in axes]].to_numpy(), columns=axes, index=data.index),
+                # we need to pass sort=True also when the index is sorted to ensure that the divisions are computed
+                sort=index_monotonically_increasing,
                 **kwargs,
             )
             if feature_key is not None:
                 feature_categ = dd.from_pandas(
                     data[feature_key].astype(str).astype("category"),
+                    sort=index_monotonically_increasing,
                     **kwargs,
                 )  # type: ignore[attr-defined]
                 table[feature_key] = feature_categ
@@ -674,7 +688,10 @@ class PointsModel:
         if Z not in axes and Z in data.columns:
             logger.info(f"Column `{Z}` in `data` will be ignored since the data is 2D.")
         for c in set(data.columns) - {feature_key, instance_key, *coordinates.values(), X, Y, Z}:
-            table[c] = data[c]
+            column = data[c]
+            if isinstance(column, dd.Series):
+                column = column.compute()
+            table[c] = dd.from_pandas(column, npartitions=table.npartitions, sort=index_monotonically_increasing)
 
         validated = cls._add_metadata_and_validate(
             table, feature_key=feature_key, instance_key=instance_key, transformations=transformations


Test patch:
diff --git a/tests/core/operations/test_rasterize_bins.py b/tests/core/operations/test_rasterize_bins.py
index 6a24725..10eab97 100644
--- a/tests/core/operations/test_rasterize_bins.py
+++ b/tests/core/operations/test_rasterize_bins.py
@@ -3,9 +3,11 @@ from __future__ import annotations
 import numpy as np
 import pytest
 from anndata import AnnData
+from geopandas import GeoDataFrame
 from numpy.random import default_rng
 from pandas import DataFrame
 from scipy.sparse import csr_matrix
+from shapely.geometry import Polygon
 from spatialdata._core.data_extent import are_extents_equal, get_extent
 from spatialdata._core.operations.rasterize_bins import rasterize_bins
 from spatialdata._core.spatialdata import SpatialData
@@ -25,7 +27,7 @@ def _get_bins_data(n: int) -> tuple[ArrayLike, ArrayLike, ArrayLike]:
     return np.dot(data, rotation.T), x, y
 
 
-@pytest.mark.parametrize("geometry", ["points", "shapes"])
+@pytest.mark.parametrize("geometry", ["points", "circles", "squares"])
 @pytest.mark.parametrize("value_key", [None, "instance_id", ["gene0", "gene1"]])
 def test_rasterize_bins(geometry: str, value_key: str | list[str] | None):
     n = 10
@@ -34,9 +36,16 @@ def test_rasterize_bins(geometry: str, value_key: str | list[str] | None):
 
     if geometry == "points":
         points = PointsModel.parse(data, transformations={"global": scale})
-    else:
-        assert geometry == "shapes"
+    elif geometry == "circles":
         points = ShapesModel.parse(data, geometry=0, radius=1, transformations={"global": scale})
+    else:
+        assert geometry == "squares"
+
+        gdf = GeoDataFrame(
+            data={"geometry": [Polygon([(x, y), (x + 1, y), (x + 1, y + 1), (x, y + 1), (x, y)]) for x, y in data]}
+        )
+
+        points = ShapesModel.parse(gdf, transformations={"global": scale})
 
     obs = DataFrame(
         data={"region": ["points"] * n * n, "instance_id": np.arange(n * n), "col_index": x, "row_index": y}
diff --git a/tests/models/test_models.py b/tests/models/test_models.py
index c3d992c..9ee214d 100644
--- a/tests/models/test_models.py
+++ b/tests/models/test_models.py
@@ -453,3 +453,14 @@ def test_force2d():
     assert_elements_are_identical(circles_3d, expected_circles_2d)
     assert_elements_are_identical(polygons_3d, expected_polygons_2d)
     assert_elements_are_identical(multipolygons_3d, expected_multipolygons_2d)
+
+
+def test_dask_points_unsorted_index(points):
+    element = points["points_0"]
+    new_order = RNG.permutation(len(element))
+    with pytest.warns(
+        UserWarning,
+        match=r"The index of the dataframe is not monotonic increasing\.",
+    ):
+        ordered = PointsModel.parse(element.compute().iloc[new_order, :])
+    assert np.all(ordered.index.compute().to_numpy() == new_order)