From a77116c6f379c75653ee74865a46d31d6ec1a0cd Mon Sep 17 00:00:00 2001
From: Timo Millenaar <timo.m@live.nl>
Date: Sun, 7 Sep 2025 11:42:10 +0200
Subject: [PATCH 01/10] Do intersection based on bounds instead of ids to allow
 for interection call on grids that don't aling

---
 gridkit-rs/src/tile.rs | 27 +++++++++++++++++++++++----
 1 file changed, 23 insertions(+), 4 deletions(-)

diff --git a/gridkit-rs/src/tile.rs b/gridkit-rs/src/tile.rs
index 2e95ea1ef..14070650e 100644
--- a/gridkit-rs/src/tile.rs
+++ b/gridkit-rs/src/tile.rs
@@ -164,6 +164,12 @@ impl TileTraits for Tile {
     }
 
     fn corners(&self) -> Array2<f64> {
+        // Returns corners in order:
+        // top-left, top-right, bottom-right, bottom-left
+        // 0 > 1
+        //     v
+        // 3 < 2
+        //
         let start_corner_x = self.start_id.0 as f64 * self.grid.dx() + self.grid.offset()[0];
         let start_corner_y = self.start_id.1 as f64 * self.grid.dy() + self.grid.offset()[1];
         let side_length_x = self.nx as f64 * self.grid.dx();
@@ -207,15 +213,28 @@ impl TileTraits for Tile {
     }
 
     fn intersects(&self, other: &Tile) -> bool {
+        let corners_self = self.corners();
+        let corners_other = other.corners();
+
+        let left_self = corners_self[Ix2(3,0)];
+        let bottom_self = corners_self[Ix2(3,1)];
+        let right_self = corners_self[Ix2(1,0)];
+        let top_self = corners_self[Ix2(1,1)];
+
+        let left_other = corners_other[Ix2(3,0)];
+        let bottom_other = corners_other[Ix2(3,1)];
+        let right_other = corners_other[Ix2(1,0)];
+        let top_other = corners_other[Ix2(1,1)];
         return !(
-                self.start_id.0 >= (other.start_id.0 + other.nx as i64)
-            || (self.start_id.0 + self.nx as i64) <= other.start_id.0
-            ||  self.start_id.1 >= (other.start_id.1 + other.ny as i64)
-            || (self.start_id.1 + self.ny as i64) <= other.start_id.1
+               left_self   >= right_other
+            || right_self  <= left_other
+            || bottom_self >= top_other
+            || top_self    <= bottom_other
         );
     }
 
     fn overlap(&self, other: &Tile) -> Result<Tile, String> {
+        // TODO: check whether grids align
         if !self.intersects(&other) {
             return Err("Tiles do not overlap".to_string());
         }

From 28ac4ddf74110f9ded52470e4e9d778f2ebe703d Mon Sep 17 00:00:00 2001
From: Timo Millenaar <timo.m@live.nl>
Date: Tue, 23 Sep 2025 23:56:33 +0200
Subject: [PATCH 02/10] Fix intersection on bounds for rotated grids

---
 gridkit-py/src/lib.rs                         |  8 ++--
 .../test_statistical_functions.py             |  5 ++-
 gridkit-rs/src/tile.rs                        | 37 +++++++++++++------
 3 files changed, 33 insertions(+), 17 deletions(-)

diff --git a/gridkit-py/src/lib.rs b/gridkit-py/src/lib.rs
index 61dd3824a..9ef0ec8bf 100644
--- a/gridkit-py/src/lib.rs
+++ b/gridkit-py/src/lib.rs
@@ -954,7 +954,7 @@ impl PyO3TriGrid {
     ) -> PyO3TriGrid {
         let target_loc: [f64; 2] = target_loc.into();
         let cell_element =
-            CellElement::from_string(&cell_element).expect("Unsupported cell_element {}");
+            CellElement::from_string(&cell_element).expect(&format!("Unsupported cell_element: '{}'", cell_element));
         let mut new_grid = self.clone();
         new_grid._grid.anchor_inplace(&target_loc, cell_element);
         new_grid
@@ -968,7 +968,7 @@ impl PyO3TriGrid {
     ) {
         let target_loc: [f64; 2] = target_loc.into();
         let cell_element =
-            CellElement::from_string(&cell_element).expect("Unsupported cell_element {}");
+            CellElement::from_string(&cell_element).expect(&format!("Unsupported cell_element: '{}'", cell_element));
         self._grid.anchor_inplace(&target_loc, cell_element);
     }
 
@@ -1145,7 +1145,7 @@ impl PyO3RectGrid {
     ) -> PyO3RectGrid {
         let target_loc: [f64; 2] = target_loc.into();
         let cell_element =
-            CellElement::from_string(&cell_element).expect("Unsupported cell_element {}");
+            CellElement::from_string(&cell_element).expect(&format!("Unsupported cell_element: '{}'", cell_element));
         let mut new_grid = self.clone();
         new_grid._grid.anchor_inplace(&target_loc, cell_element);
         new_grid
@@ -1159,7 +1159,7 @@ impl PyO3RectGrid {
     ) {
         let target_loc: [f64; 2] = target_loc.into();
         let cell_element =
-            CellElement::from_string(&cell_element).expect("Unsupported cell_element {}");
+            CellElement::from_string(&cell_element).expect(&format!("Unsupported cell_element: '{}'", cell_element));
         self._grid.anchor_inplace(&target_loc, cell_element);
     }
 
diff --git a/gridkit-py/tests/test_gridkit/test_statistical_functions.py b/gridkit-py/tests/test_gridkit/test_statistical_functions.py
index 918a737a0..1915b2352 100644
--- a/gridkit-py/tests/test_gridkit/test_statistical_functions.py
+++ b/gridkit-py/tests/test_gridkit/test_statistical_functions.py
@@ -75,8 +75,9 @@ def test_combine_tiles(grid):
         HexGrid(size=1, shape="flat"),
     ],
 )
-def test_count_data_tile(grid):
-    grid.rotation = 12
+@pytest.mark.parametrize("rotation", [-12, 363, 186])
+def test_count_data_tile(grid, rotation):
+    grid.rotation = rotation
 
     tile1 = Tile(grid, (-2, -3), 5, 7)
     tile2 = Tile(grid, (1, 0), 7, 5).to_data_tile_with_value(3.2)
diff --git a/gridkit-rs/src/tile.rs b/gridkit-rs/src/tile.rs
index 14070650e..ad56823dd 100644
--- a/gridkit-rs/src/tile.rs
+++ b/gridkit-rs/src/tile.rs
@@ -213,20 +213,35 @@ impl TileTraits for Tile {
     }
 
     fn intersects(&self, other: &Tile) -> bool {
+
+        // Note: It is easier to match tiles on their ids. This seems elegant
+        //       because we don't need to care about rotation. This was how it
+        //       was first implemented but it required self and other to be on
+        //       the same grid. This was not enforced so you would get weird
+        //       results if this was not upheld. Rather than forcing the grids
+        //       to be aligned we use the tile corner coordinates. This way
+        //       any tile can be checked for intersection with any other tile
+        //       purely based on location, regardless of the host grid. This
+        //       does have the downside that we need to account for rotated
+        //       grids where the corner that was originally the on the left
+        //       side might now be the rightmost corner.
         let corners_self = self.corners();
         let corners_other = other.corners();
 
-        let left_self = corners_self[Ix2(3,0)];
-        let bottom_self = corners_self[Ix2(3,1)];
-        let right_self = corners_self[Ix2(1,0)];
-        let top_self = corners_self[Ix2(1,1)];
-
-        let left_other = corners_other[Ix2(3,0)];
-        let bottom_other = corners_other[Ix2(3,1)];
-        let right_other = corners_other[Ix2(1,0)];
-        let top_other = corners_other[Ix2(1,1)];
-        return !(
-               left_self   >= right_other
+        // Get min and max values.
+        // Note: This is a thorn in the eye as .min() does not work
+        //       on ndarrays of dtype float because of nan and inf values.
+        let left_self = corners_self.slice(s![..,0]).iter().cloned().reduce(f64::min).unwrap();
+        let bottom_self = corners_self.slice(s![..,1]).iter().cloned().reduce(f64::min).unwrap();
+        let right_self = corners_self.slice(s![..,0]).iter().cloned().reduce(f64::max).unwrap();
+        let top_self = corners_self.slice(s![..,1]).iter().cloned().reduce(f64::max).unwrap();
+
+        let left_other = corners_other.slice(s![..,0]).iter().cloned().reduce(f64::min).unwrap();
+        let bottom_other = corners_other.slice(s![..,1]).iter().cloned().reduce(f64::min).unwrap();
+        let right_other = corners_other.slice(s![..,0]).iter().cloned().reduce(f64::max).unwrap();
+        let top_other = corners_other.slice(s![..,1]).iter().cloned().reduce(f64::max).unwrap();
+
+        return !(left_self   >= right_other
             || right_self  <= left_other
             || bottom_self >= top_other
             || top_self    <= bottom_other

From b79c4c27df06dc7e41acdf4d03ef53dba39887e9 Mon Sep 17 00:00:00 2001
From: Timo Millenaar <tmillenaar@gmail.com>
Date: Mon, 29 Sep 2025 21:43:51 +0200
Subject: [PATCH 03/10] Make rust tile hashable and equatable

---
 gridkit-rs/src/grid.rs      |  4 ++--
 gridkit-rs/src/hex_grid.rs  | 26 ++++++++++++++++++++++++++
 gridkit-rs/src/rect_grid.rs | 26 ++++++++++++++++++++++++++
 gridkit-rs/src/tile.rs      |  2 +-
 gridkit-rs/src/tri_grid.rs  | 26 ++++++++++++++++++++++++++
 5 files changed, 81 insertions(+), 3 deletions(-)

diff --git a/gridkit-rs/src/grid.rs b/gridkit-rs/src/grid.rs
index be68ceee1..dbf17372c 100644
--- a/gridkit-rs/src/grid.rs
+++ b/gridkit-rs/src/grid.rs
@@ -135,7 +135,7 @@ pub trait GridTraits {
     // ) -> Array1<f64>;
 }
 
-#[derive(Clone, Debug)]
+#[derive(Clone, Debug, PartialEq, Hash)]
 #[enum_delegate::implement(GridTraits)]
 pub enum Grid {
     TriGrid(TriGrid),
@@ -168,7 +168,7 @@ impl CellElement {
     }
 }
 
-#[derive(Clone, Debug)]
+#[derive(Clone, Debug, Hash, PartialEq)]
 pub enum Orientation {
     Flat,
     Pointy,
diff --git a/gridkit-rs/src/hex_grid.rs b/gridkit-rs/src/hex_grid.rs
index 311e00aaa..5cdcc721c 100644
--- a/gridkit-rs/src/hex_grid.rs
+++ b/gridkit-rs/src/hex_grid.rs
@@ -1,3 +1,5 @@
+use std::hash::{Hash, Hasher};
+
 use crate::grid::{GridTraits, Orientation};
 use crate::utils::*;
 use ndarray::Slice;
@@ -13,6 +15,30 @@ pub struct HexGrid {
     pub _rotation_matrix_inv: Array2<f64>,
 }
 
+impl PartialEq for HexGrid {
+    // Needs manual implementation, derive PartialEq does not work on floats because of NaN etc.
+    fn eq(&self, other: &Self) -> bool {
+        self.cellsize.to_bits() == other.cellsize.to_bits() &&
+        self.offset[0].to_bits() == other.offset[0].to_bits() &&
+        self.offset[1].to_bits() == other.offset[1].to_bits() &&
+        self.orientation == other.orientation &&
+        self._rotation.to_bits() == other._rotation.to_bits()
+    }
+}
+
+impl Eq for HexGrid {}
+
+impl Hash for HexGrid {
+    // Needs manual implementation, derive Hash does not work on floats because of NaN etc.
+    fn hash<H: Hasher>(&self, state: &mut H) {
+        self.cellsize.to_bits().hash(state);
+        self.offset[0].to_bits().hash(state);
+        self.offset[1].to_bits().hash(state);
+        self.orientation.hash(state);
+        self._rotation.to_bits().hash(state);
+    }
+}
+
 impl GridTraits for HexGrid {
     fn get_grid(&self) -> crate::Grid {
         crate::Grid::HexGrid(self.to_owned())
diff --git a/gridkit-rs/src/rect_grid.rs b/gridkit-rs/src/rect_grid.rs
index 43bbbe61c..3fb864e3b 100644
--- a/gridkit-rs/src/rect_grid.rs
+++ b/gridkit-rs/src/rect_grid.rs
@@ -1,3 +1,5 @@
+use std::hash::{Hash, Hasher};
+
 use crate::grid::GridTraits;
 use crate::tile::*;
 use crate::utils::*;
@@ -13,6 +15,30 @@ pub struct RectGrid {
     pub _rotation_matrix_inv: Array2<f64>,
 }
 
+impl PartialEq for RectGrid {
+    // Needs manual implementation, derive PartialEq does not work on floats because of NaN etc.
+    fn eq(&self, other: &Self) -> bool {
+        self._dx.to_bits() == other._dx.to_bits() &&
+        self._dy.to_bits() == other._dy.to_bits() &&
+        self.offset[0].to_bits() == other.offset[0].to_bits() &&
+        self.offset[1].to_bits() == other.offset[1].to_bits() &&
+        self._rotation.to_bits() == other._rotation.to_bits()
+    }
+}
+
+impl Eq for RectGrid {}
+
+impl Hash for RectGrid {
+    // Needs manual implementation, derive Hash does not work on floats because of NaN etc.
+    fn hash<H: Hasher>(&self, state: &mut H) {
+        self._dx.to_bits().hash(state);
+        self._dy.to_bits().hash(state);
+        self.offset[0].to_bits().hash(state);
+        self.offset[1].to_bits().hash(state);
+        self._rotation.to_bits().hash(state);
+    }
+}
+
 impl GridTraits for RectGrid {
     fn get_grid(&self) -> crate::Grid {
         crate::Grid::RectGrid(self.to_owned())
diff --git a/gridkit-rs/src/tile.rs b/gridkit-rs/src/tile.rs
index ad56823dd..750476f76 100644
--- a/gridkit-rs/src/tile.rs
+++ b/gridkit-rs/src/tile.rs
@@ -4,7 +4,7 @@ use std::{f32::MAX, f64, i64, u64};
 use crate::{data_tile::DataTile, grid::*, hex_grid::HexGrid};
 use ndarray::*;
 
-#[derive(Clone, Debug)]
+#[derive(Clone, Debug, PartialEq, Hash)]
 pub struct Tile {
     pub grid: Grid,
     pub start_id: (i64, i64),
diff --git a/gridkit-rs/src/tri_grid.rs b/gridkit-rs/src/tri_grid.rs
index e95caa7a4..2157c4959 100644
--- a/gridkit-rs/src/tri_grid.rs
+++ b/gridkit-rs/src/tri_grid.rs
@@ -1,4 +1,5 @@
 use std::fmt::Debug;
+use std::hash::{Hash, Hasher};
 
 use crate::grid::{GridTraits, Orientation};
 use crate::interpolate;
@@ -15,6 +16,31 @@ pub struct TriGrid {
     pub _rotation_matrix_inv: Array2<f64>,
 }
 
+impl PartialEq for TriGrid {
+    // Needs manual implementation, derive PartialEq does not work on floats because of NaN etc.
+    fn eq(&self, other: &Self) -> bool {
+        println!("COWBANGA!!!!");
+        self.cellsize.to_bits() == other.cellsize.to_bits() &&
+        self.offset[0].to_bits() == other.offset[0].to_bits() &&
+        self.offset[1].to_bits() == other.offset[1].to_bits() &&
+        self.orientation == other.orientation &&
+        self._rotation.to_bits() == other._rotation.to_bits()
+    }
+}
+
+impl Eq for TriGrid {}
+
+impl Hash for TriGrid {
+    // Needs manual implementation, derive Hash does not work on floats because of NaN etc.
+    fn hash<H: Hasher>(&self, state: &mut H) {
+        self.cellsize.to_bits().hash(state);
+        self.offset[0].to_bits().hash(state);
+        self.offset[1].to_bits().hash(state);
+        self.orientation.hash(state);
+        self._rotation.to_bits().hash(state);
+    }
+}
+
 impl GridTraits for TriGrid {
     fn get_grid(&self) -> crate::Grid {
         crate::Grid::TriGrid(self.to_owned())

From 26022b99d86383995395d9619a5362eed4811086 Mon Sep 17 00:00:00 2001
From: Timo Millenaar <tmillenaar@gmail.com>
Date: Mon, 29 Sep 2025 21:51:59 +0200
Subject: [PATCH 04/10] Add equality check on Tile to PyO3Tile and PyO3DataTile

---
 gridkit-py/src/lib.rs | 56 ++++++++++++++++++++++++++++---------------
 1 file changed, 37 insertions(+), 19 deletions(-)

diff --git a/gridkit-py/src/lib.rs b/gridkit-py/src/lib.rs
index 9ef0ec8bf..eb8a7636d 100644
--- a/gridkit-py/src/lib.rs
+++ b/gridkit-py/src/lib.rs
@@ -474,7 +474,7 @@ impl PyO3DataTile {
     }
 }
 
-#[derive(Clone)]
+#[derive(Clone, PartialEq, Hash)]
 // #[enum_delegate::implement(GridTraits)]
 pub enum PyO3Grid {
     PyO3TriGrid(PyO3TriGrid),
@@ -483,7 +483,7 @@ pub enum PyO3Grid {
 }
 
 #[pyclass]
-#[derive(Clone)]
+#[derive(Clone, PartialEq, Hash)]
 struct PyO3Tile {
     // #[pyo3(get, set)]
     _grid: PyO3Grid,
@@ -803,6 +803,10 @@ impl PyO3Tile {
         }
     }
 
+    fn equal_to_tile<'py>(&self, _py: Python<'py>, other: &PyO3Tile) -> bool {
+        self._tile == other._tile
+    }
+
     fn corner_ids<'py>(&self, py: Python<'py>) -> &'py PyArray2<i64> {
         self._tile.corner_ids().into_pyarray(py)
     }
@@ -877,13 +881,19 @@ impl PyO3Tile {
             .grid_id_to_tile_id(&index, oob_value)
             .into_pyarray(py)
     }
+
+    fn __eq__<'py>(&self, py: Python<'py>, other: PyO3Tile) -> bool {
+        self._tile == other._tile
+    }
+
+    fn __ne__<'py>(&self, py: Python<'py>, other: PyO3Tile) -> bool {
+        self._tile != other._tile
+    }
 }
 
-#[derive(Clone)]
+#[derive(Clone, PartialEq, Hash)]
 #[pyclass]
 struct PyO3TriGrid {
-    cellsize: f64,
-    rotation: f64,
     _grid: tri_grid::TriGrid,
 }
 
@@ -902,8 +912,6 @@ impl PyO3TriGrid {
                 _grid.set_offset(offset.into());
                 _grid.set_rotation(rotation);
                 Ok(PyO3TriGrid {
-                    cellsize,
-                    rotation,
                     _grid,
                 })
             }
@@ -918,6 +926,10 @@ impl PyO3TriGrid {
         self._grid.offset.into()
     }
 
+    fn rotation(&self) -> f64 {
+        self._grid._rotation
+    }
+
     fn cell_height(&self) -> f64 {
         self._grid.cell_height()
     }
@@ -938,6 +950,10 @@ impl PyO3TriGrid {
         self._grid.dy()
     }
 
+    fn cellsize(&self) -> f64 {
+        self._grid.cellsize
+    }
+
     fn rotation_matrix<'py>(&self, py: Python<'py>) -> &'py PyArray2<f64> {
         &self._grid.rotation_matrix().clone().into_pyarray(py)
     }
@@ -1085,12 +1101,9 @@ impl PyO3TriGrid {
     }
 }
 
-#[derive(Clone)]
+#[derive(Clone, PartialEq, Hash)]
 #[pyclass]
 struct PyO3RectGrid {
-    dx: f64,
-    dy: f64,
-    rotation: f64,
     _grid: rect_grid::RectGrid,
 }
 
@@ -1102,9 +1115,6 @@ impl PyO3RectGrid {
         _grid.set_offset(offset.into());
         _grid.set_rotation(rotation);
         PyO3RectGrid {
-            dx,
-            dy,
-            rotation,
             _grid,
         }
     }
@@ -1129,6 +1139,10 @@ impl PyO3RectGrid {
         self._grid.offset.into()
     }
 
+    fn rotation(&self) -> f64 {
+        self._grid.rotation()
+    }
+
     fn rotation_matrix<'py>(&self, py: Python<'py>) -> &'py PyArray2<f64> {
         &self._grid.rotation_matrix().clone().into_pyarray(py)
     }
@@ -1201,11 +1215,9 @@ impl PyO3RectGrid {
     }
 }
 
-#[derive(Clone)]
+#[derive(Clone, PartialEq, Hash)]
 #[pyclass]
 struct PyO3HexGrid {
-    cellsize: f64,
-    rotation: f64,
     _grid: hex_grid::HexGrid,
 }
 
@@ -1224,8 +1236,6 @@ impl PyO3HexGrid {
                 _grid.set_offset(offset.into());
                 _grid.set_rotation(rotation);
                 Ok(PyO3HexGrid {
-                    cellsize,
-                    rotation,
                     _grid,
                 })
             }
@@ -1274,6 +1284,14 @@ impl PyO3HexGrid {
         self._grid.dy()
     }
 
+    fn rotation(&self) -> f64 {
+        self._grid.rotation()
+    }
+
+    fn cellsize(&self) -> f64 {
+        self._grid.cellsize
+    }
+
     fn rotation_matrix<'py>(&self, py: Python<'py>) -> &'py PyArray2<f64> {
         &self._grid.rotation_matrix().clone().into_pyarray(py)
     }

From e73d5c69adc4588566db9ec7997a52e61648d736 Mon Sep 17 00:00:00 2001
From: Timo Millenaar <tmillenaar@gmail.com>
Date: Mon, 29 Sep 2025 22:07:48 +0200
Subject: [PATCH 05/10] Add equality check to Grid in Python

---
 gridkit-py/gridkit/base_grid.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/gridkit-py/gridkit/base_grid.py b/gridkit-py/gridkit/base_grid.py
index 61f9c2ef0..8bb52f27c 100644
--- a/gridkit-py/gridkit/base_grid.py
+++ b/gridkit-py/gridkit/base_grid.py
@@ -751,3 +751,11 @@ def update(
         **kwargs,
     ):
         pass
+
+    def __eq__(self, other):
+        if isinstance(other, BaseGrid):
+            return self._grid == other._grid
+        return False
+
+    def __ne__(self, other):
+        return not self == other

From a4a7a29a85b33f179f368b6eb041a3a485261912 Mon Sep 17 00:00:00 2001
From: Timo Millenaar <tmillenaar@gmail.com>
Date: Tue, 30 Sep 2025 22:58:33 +0200
Subject: [PATCH 06/10] Rework how the tile and grid attributes are handled to
 prevent them getting out of sync. There were problems where for example
 tile.grid.rotation or any other grid property was modified but the
 tile._tile.grid was unchanged. Also, we cannot have changes to datatile.grid
 and datatile.tile affect the datatile because we would need to resample so we
 only allow getting a copy of the grid or tile there.

---
 gridkit-py/gridkit/base_grid.py               |   4 +-
 gridkit-py/gridkit/io.py                      |   4 +-
 gridkit-py/gridkit/tile.py                    | 223 ++++++++++++------
 gridkit-py/gridkit/tri_grid.py                |   6 +-
 gridkit-py/src/lib.rs                         |  12 +
 .../tests/test_gridkit/test_data_tile.py      |  16 +-
 gridkit-py/tests/test_gridkit/test_io.py      |  24 +-
 gridkit-py/tests/test_gridkit/test_tile.py    |  47 ++++
 8 files changed, 245 insertions(+), 91 deletions(-)

diff --git a/gridkit-py/gridkit/base_grid.py b/gridkit-py/gridkit/base_grid.py
index 8bb52f27c..ae6b4e131 100644
--- a/gridkit-py/gridkit/base_grid.py
+++ b/gridkit-py/gridkit/base_grid.py
@@ -753,8 +753,8 @@ def update(
         pass
 
     def __eq__(self, other):
-        if isinstance(other, BaseGrid):
-            return self._grid == other._grid
+        if type(self) == type(other):
+            return self._grid.equals(other._grid)
         return False
 
     def __ne__(self, other):
diff --git a/gridkit-py/gridkit/io.py b/gridkit-py/gridkit/io.py
index 8ee64b698..49ccbf169 100644
--- a/gridkit-py/gridkit/io.py
+++ b/gridkit-py/gridkit/io.py
@@ -90,7 +90,7 @@ def read_geotiff(*args, bands=1, **kwargs):
 
 def _write_data_tile_to_raster(data_tile, path):
     """Intended to be called only through :func:`.write_raster`"""
-    if data_tile.grid.rotation:
+    if data_tile.get_grid().rotation:
         raise ValueError(
             "Cannot write a data tile as raster if the grid has a rotaion."
         )
@@ -122,7 +122,7 @@ def _write_data_tile_to_raster(data_tile, path):
         width=data_tile.nx,
         count=1,  # nr bands
         dtype=data_tile.dtype,
-        crs=data_tile.grid.crs,
+        crs=data_tile.get_grid().crs,
         nodata=data_tile.nodata_value,
         transform=transform,
     ) as dst:
diff --git a/gridkit-py/gridkit/tile.py b/gridkit-py/gridkit/tile.py
index 4836b475f..6799fff2d 100644
--- a/gridkit-py/gridkit/tile.py
+++ b/gridkit-py/gridkit/tile.py
@@ -61,6 +61,10 @@ def __init__(
         nx: int,
         ny: int,
     ):
+        if not isinstance(grid, BaseGrid):
+            raise TypeError(
+                f"Unexpected type for 'grid', expected a TriGrid, RectGrid or HexGrid, got a: {type(grid)}"
+            )
 
         if not numpy.isclose(nx % 1, 0):
             raise ValueError(f"Expected an integer for 'nx', got: {nx}")
@@ -74,23 +78,15 @@ def __init__(
             raise ValueError(
                 "'start_id' must be a single pair of indices in the form (x,y), got: {start_id}"
             )
-        start_id = (
+        self.start_id = (
             tuple(start_id.index)
             if isinstance(start_id, GridIndex)
             else tuple(start_id)
         )
+        self.nx = int(nx)
+        self.ny = int(ny)
 
-        if isinstance(grid, TriGrid):
-            self._tile = PyO3Tile.from_tri_grid(grid._grid, start_id, nx, ny)
-        elif isinstance(grid, RectGrid):
-            self._tile = PyO3Tile.from_rect_grid(grid._grid, start_id, nx, ny)
-        elif isinstance(grid, HexGrid):
-            self._tile = PyO3Tile.from_hex_grid(grid._grid, start_id, nx, ny)
-        else:
-            raise TypeError(
-                f"Unexpected type for 'grid', expected a TriGrid, RectGrid or HexGrid, got a: {type(grid)}"
-            )
-        self.grid = grid.update()
+        self._grid = grid.update()
 
     @staticmethod
     def from_pyo3_tile(grid, pyo3_tile):
@@ -155,24 +151,48 @@ def to_data_tile_with_value(self, fill_value, nodata_value=None):
             raise AttributeError(f"Method {method_name} not found on tile")
 
         py03_data_tile = method(fill_value, nodata_value)
-        return DataTile.from_pyo3_data_tile(self.grid.update(), py03_data_tile)
+        return DataTile.from_pyo3_data_tile(self.get_grid().update(), py03_data_tile)
 
     @property
-    def start_id(self):
-        """The starting cell of the Tile.
-        The starting cell defines the bottom-left corner of the Tile if the associated grid is not rotated.
-        """
-        return GridIndex(self._tile.start_id)
+    def grid(self):
+        """ """
+        return self._grid
 
-    @property
-    def nx(self):
-        """The number of cells in x direction, starting from the ``start_id``"""
-        return self._tile.nx
+    @grid.setter
+    def grid(self, value):
+        self.update(grid=value)
+
+    def get_grid(self):
+        """Obtain a copy of the grid associated to this Tile"""
+        return self._grid.update()
 
     @property
-    def ny(self):
-        """The number of cells in y direction, starting from the ``start_id``"""
-        return self._tile.ny
+    def _tile(self):
+        """The PyO3Tile object.
+        Only meant for use internal to this class, access at risk of side-effects.
+
+        We create a new PyP3Tile every time. This is to avoid synchronisation issues.
+        Before we had a PyO3Tile created at Tile initialization which served as the
+        source of truth. This meant that parameters like 'nx' for example were properties
+        that returned the value of the PyO3Tile. There was however a problem when someone
+        modified a value on the grid property, because this was not reflected in the grid
+        on the PyO3Tile. This caused a divergence where tile.grid no longer reflected
+        the values used by the rust internals. Calling a function like 'corners' for
+        would then return unexpected values. To remedy this we use the properties on the
+        python object as the source of truth and create a new PyO3Tile whenever it is
+        requested.
+        """
+        grid = self.get_grid()
+        if isinstance(grid, TriGrid):
+            return PyO3Tile.from_tri_grid(grid._grid, self.start_id, self.nx, self.ny)
+        elif isinstance(grid, RectGrid):
+            return PyO3Tile.from_rect_grid(grid._grid, self.start_id, self.nx, self.ny)
+        elif isinstance(grid, HexGrid):
+            return PyO3Tile.from_hex_grid(grid._grid, self.start_id, self.nx, self.ny)
+
+        raise TypeError(
+            f"Unexpected type for 'grid', expected a TriGrid, RectGrid or HexGrid, got a: {type(grid)}"
+        )
 
     def corner_ids(self):
         """The ids at the corners of the Tile
@@ -243,16 +263,16 @@ def intersects(self, other):
     def centroid(self, index=None):
         if index is None:
             index = self.indices
-        return self.grid.centroid(index)
+        return self.get_grid().centroid(index)
 
     def overlap(self, other):
-        is_aligned, reason = self.grid.is_aligned_with(other.grid)
+        is_aligned, reason = self.get_grid().is_aligned_with(other.get_grid())
         if not is_aligned:
             raise AlignmentError(
                 f"Cannot find overlap of grids that are not aligned. Reason for misalignemnt: {reason}"
             )
         _tile = self._tile.overlap(other._tile)
-        return Tile.from_pyo3_tile(self.grid, _tile)
+        return Tile.from_pyo3_tile(self.get_grid(), _tile)
 
     @validate_index
     def to_shapely(self, index=None, as_multipolygon=None):
@@ -269,7 +289,7 @@ def to_shapely(self, index=None, as_multipolygon=None):
 
         if index is None:
             index = self.indices
-        return self.grid.to_shapely(index, as_multipolygon=as_multipolygon)
+        return self.get_grid().to_shapely(index, as_multipolygon=as_multipolygon)
 
     def tile_id_to_grid_id(self, tile_id=None, oob_value=numpy.iinfo(numpy.int64).max):
         """Convert the index referring to a cell from the tile reference frame to that of the grid.
@@ -425,6 +445,57 @@ def grid_id_to_tile_id(self, index=None, oob_value=numpy.iinfo(numpy.int64).max)
         #       is different when using arrays of integers compared to tuples.
         return tuple(self._tile.grid_id_to_tile_id(index, oob_value=oob_value).T)
 
+    def __eq__(self, other):
+        if isinstance(other, Tile):
+            return self._tile == other._tile
+        return False
+
+    def __ne__(self, other):
+        return not self == other
+
+    def update(
+        self,
+        grid=None,
+        start_id=None,
+        nx=None,
+        ny=None,
+    ):
+        """Return a modified copy of the original Tile.
+        The original Tile remains un-mutated.
+
+        Parameters
+        ----------
+        grid: :class:`.BaseGrid`
+            The :class:`.TriGrid`, :class:`.RectGrid` or :class:`.HexGrid` the Tile is associated with
+        start_id: Union[Tuple[int, int], GridIndex]
+            The starting cell of the Tile.
+            The starting cell defines the bottom-left corner of the Tile if the associated grid is not rotated.
+        nx: int
+            The number of cells in x direction, starting from the ``start_id``
+        ny: int
+            The number of cells in y direction, starting from the ``start_id``
+
+        Returns
+        -------
+        :class:`.Tile`
+            A modified copy of the current Tile
+        """
+        if grid is None:
+            grid = self.get_grid().update()
+        else:
+            if not isinstance(grid, BaseGrid):
+                raise TypeError(
+                    f"Expected a TriGrid, RectGrid or HexGrid, got a {type(value)}"
+                )
+        if start_id is None:
+            start_id = self.start_id
+        if nx is None:
+            nx = self.nx
+        if ny is None:
+            ny = self.ny
+
+        return Tile(grid, start_id, nx, ny)
+
 
 class DataTile(Tile):
 
@@ -450,7 +521,6 @@ def __init__(self, tile: Tile, data: numpy.ndarray, nodata_value=None):
                     raise TypeError(
                         f"Data type of the supplied array in argument 'data' (dtype: '{data.dtype}') did not match the supplied nodata value '{nodata_value}' of type '{type(nodata_value)}'"
                     ) from e
-        self.grid = tile.grid
 
         # Map numpy dtypes to method suffixes
         dtype_method_map = {
@@ -476,8 +546,10 @@ def __init__(self, tile: Tile, data: numpy.ndarray, nodata_value=None):
             raise AttributeError(f"Method {method_name} not found on tile")
         self._data_tile = method(data, nodata_value)
 
+        super(DataTile, self).__init__(tile.get_grid(), tile.start_id, tile.nx, tile.ny)
+
         # _tile is used by the Tile parent class
-        self._tile = self._data_tile.get_tile()
+        # self._tile = self._data_tile.get_tile()
 
     def from_bounds_as_rect(
         data,
@@ -700,29 +772,44 @@ def __setitem__(self, item, value):
 
         new_data = self.to_numpy()
         new_data[item] = value
-        if isinstance(self.grid, TriGrid):
-            tile = PyO3Tile.from_tri_grid(
-                self.grid._grid, tuple(self.start_id.index), self.nx, self.ny
-            )
-        elif isinstance(self.grid, RectGrid):
-            tile = PyO3Tile.from_rect_grid(
-                self.grid._grid, tuple(self.start_id.index), self.nx, self.ny
-            )
-        elif isinstance(self.grid, HexGrid):
-            tile = PyO3Tile.from_hex_grid(
-                self.grid._grid, tuple(self.start_id.index), self.nx, self.ny
-            )
+        grid = self.get_grid()
+        if isinstance(grid, TriGrid):
+            tile = PyO3Tile.from_tri_grid(grid._grid, self.start_id, self.nx, self.ny)
+        elif isinstance(grid, RectGrid):
+            tile = PyO3Tile.from_rect_grid(grid._grid, self.start_id, self.nx, self.ny)
+        elif isinstance(grid, HexGrid):
+            tile = PyO3Tile.from_hex_grid(grid._grid, self.start_id, self.nx, self.ny)
         else:
-            raise TypeError(f"Unrecognized grid type: {self.grid}")
+            raise TypeError(f"Unrecognized grid type: {type(self.get_grid())}")
         new_data_tile = self.to_data_tile(new_data, self.nodata_value)
         self._data_tile = new_data_tile._data_tile
 
+    @property
+    def grid(self):
+        raise RuntimeError(
+            """The 'grid' attribute on a DataTile should not be modified directly in order to avoid synchronization issues.
+            You can obtain a copy of the grid through the .get_grid() method.
+            """
+        )
+
+    @grid.setter
+    def grid(self, value):
+        raise RuntimeError(
+            """The 'grid' attribute on a DataTile should not be modified directly in order to avoid synchronization issues.
+            Consider calling .resample() to interpolate the values onto a differnt grid.
+            """
+        )
+
+    # @property
+    # def _tile(self):
+    #     return PyO3Tile()
+
     def update(
         self, data=None, grid=None, start_id=None, nx=None, ny=None, nodata_value=None
     ):
         # TODO: Make clear that update copies the data
         if grid is None:
-            grid = self.grid
+            grid = self.get_grid()
         if data is None:
             data = self.to_numpy()
         if start_id is None:
@@ -745,7 +832,7 @@ def update(
 
     def get_tile(self):
         """A Tile object with the same properties as this DataTile, but without the data attached."""
-        return Tile(self.grid, self.start_id, self.nx, self.ny)
+        return Tile(self.get_grid(), self.start_id, self.nx, self.ny)
 
     def corner_ids(self):
         """The ids at the corners of the Tile
@@ -839,7 +926,7 @@ def interpolate(
         :py:meth:`.BaseGrid.interp_from_points`
         """
         if method == "nearest":
-            new_ids = self.grid.cell_at_point(sample_points)
+            new_ids = self.get_grid().cell_at_point(sample_points)
             return self.value(new_ids, oob_value=self.nodata_value)
         elif method == "bilinear" or method == "linear":
             return_shape = sample_points.shape[:-1]
@@ -910,13 +997,13 @@ def resample(self, alignment_grid, method="nearest", **interp_kwargs):
             tile = alignment_grid
             alignment_grid = alignment_grid.grid
 
-        if self.grid.crs is None or alignment_grid.crs is None:
+        if self.get_grid().crs is None or alignment_grid.crs is None:
             warnings.warn(
                 "`crs` not set for one or both grids. Assuming both grids have an identical CRS."
             )
             different_crs = False
         else:
-            different_crs = not self.grid.crs.is_exact_same(alignment_grid.crs)
+            different_crs = not self.get_grid().crs.is_exact_same(alignment_grid.crs)
 
         if not tile_is_given:
             # make sure the bounds align with the grid
@@ -946,7 +1033,7 @@ def resample(self, alignment_grid, method="nearest", **interp_kwargs):
                 coords = numpy.hstack([top_xy, right_xy, bottom_xy, left_xy])
 
                 transformer = Transformer.from_crs(
-                    self.grid.crs, alignment_grid.crs, always_xy=True
+                    self.get_grid().crs, alignment_grid.crs, always_xy=True
                 )
                 corners_transformed = numpy.array(transformer.transform(*coords)).T
 
@@ -999,7 +1086,7 @@ def resample(self, alignment_grid, method="nearest", **interp_kwargs):
 
         if different_crs:
             transformer = Transformer.from_crs(
-                alignment_grid.crs, self.grid.crs, always_xy=True
+                alignment_grid.crs, self.get_grid().crs, always_xy=True
             )
             original_shape = new_points.shape
             raveled_new_points = new_points.reshape(-1, 2)
@@ -1051,7 +1138,7 @@ def to_crs(self, crs, resample_method="nearest"):
         :ref:`Example: coordinate transformations <example coordinate transformations>`
 
         """
-        new_inf_grid = self.grid.to_crs(crs)
+        new_inf_grid = self.get_grid().to_crs(crs)
         return self.resample(new_inf_grid, method=resample_method)
 
     def __add__(self, other):
@@ -1076,7 +1163,7 @@ def __add__(self, other):
             except:
                 raise TypeError(f"Cannot add DataTile and `{type(other)}`")
 
-        combined = DataTile.from_pyo3_data_tile(self.grid, _data_tile)
+        combined = DataTile.from_pyo3_data_tile(self.get_grid(), _data_tile)
         return combined
 
     def __radd__(self, other):
@@ -1100,7 +1187,7 @@ def __radd__(self, other):
             except:
                 raise TypeError(f"Cannot add DataTile and `{type(other)}`")
 
-        combined = DataTile.from_pyo3_data_tile(self.grid, _data_tile)
+        combined = DataTile.from_pyo3_data_tile(self.get_grid(), _data_tile)
         return combined
 
     def __sub__(self, other):
@@ -1123,7 +1210,7 @@ def __sub__(self, other):
                 _data_tile = self._data_tile._subtract_scalar(other_converted)
             except:
                 raise TypeError(f"Cannot subtract DataTile and `{type(other)}`")
-        combined = DataTile.from_pyo3_data_tile(self.grid, _data_tile)
+        combined = DataTile.from_pyo3_data_tile(self.get_grid(), _data_tile)
         return combined
 
     def __rsub__(self, other):
@@ -1146,7 +1233,7 @@ def __rsub__(self, other):
                 _data_tile = self._data_tile._subtract_scalar_reverse(other_converted)
             except:
                 raise TypeError(f"Cannot subtract DataTile and `{type(other)}`")
-        combined = DataTile.from_pyo3_data_tile(self.grid, _data_tile)
+        combined = DataTile.from_pyo3_data_tile(self.get_grid(), _data_tile)
         return combined
 
     def __mul__(self, other):
@@ -1169,7 +1256,7 @@ def __mul__(self, other):
                 _data_tile = self._data_tile._multiply_scalar(other_converted)
             except:
                 raise TypeError(f"Cannot multiply DataTile and `{type(other)}`")
-        combined = DataTile.from_pyo3_data_tile(self.grid, _data_tile)
+        combined = DataTile.from_pyo3_data_tile(self.get_grid(), _data_tile)
         return combined
 
     def __rmul__(self, other):
@@ -1192,7 +1279,7 @@ def __rmul__(self, other):
                 _data_tile = self._data_tile._multiply_scalar_reverse(other_converted)
             except:
                 raise TypeError(f"Cannot multiply DataTile and `{type(other)}`")
-        combined = DataTile.from_pyo3_data_tile(self.grid, _data_tile)
+        combined = DataTile.from_pyo3_data_tile(self.get_grid(), _data_tile)
         return combined
 
     def __truediv__(self, other):
@@ -1215,7 +1302,7 @@ def __truediv__(self, other):
                 _data_tile = self._data_tile._divide_scalar(other_converted)
             except:
                 raise TypeError(f"Cannot divide DataTile and `{type(other)}`")
-        combined = DataTile.from_pyo3_data_tile(self.grid, _data_tile)
+        combined = DataTile.from_pyo3_data_tile(self.get_grid(), _data_tile)
         return combined
 
     def __rtruediv__(self, other):
@@ -1238,7 +1325,7 @@ def __rtruediv__(self, other):
                 _data_tile = self._data_tile._divide_scalar_reverse(other_converted)
             except:
                 raise TypeError(f"Cannot divide DataTile and `{type(other)}`")
-        combined = DataTile.from_pyo3_data_tile(self.grid, _data_tile)
+        combined = DataTile.from_pyo3_data_tile(self.get_grid(), _data_tile)
         return combined
 
     def __pow__(self, other):
@@ -1254,7 +1341,7 @@ def __pow__(self, other):
                 _data_tile = self._data_tile._powf(other_converted)
             except:
                 raise TypeError(f"Cannot divide DataTile and `{type(other)}`")
-        combined = DataTile.from_pyo3_data_tile(self.grid, _data_tile)
+        combined = DataTile.from_pyo3_data_tile(self.get_grid(), _data_tile)
         return combined
 
     def __rpow__(self, other):
@@ -1268,7 +1355,7 @@ def __rpow__(self, other):
                 _data_tile = self._data_tile._powf_reverse(other_converted)
             except:
                 raise TypeError(f"Cannot divide DataTile and `{type(other)}`")
-        combined = DataTile.from_pyo3_data_tile(self.grid, _data_tile)
+        combined = DataTile.from_pyo3_data_tile(self.get_grid(), _data_tile)
         return combined
 
     def __eq__(self, other):
@@ -1390,7 +1477,7 @@ def combine_tiles(tiles):
                 f"Expected all Tile or DataTile objects but also got a: {type(tile)}"
             )
     pyo3_tile = tile_utils.combine_tiles(pyo3_tiles)
-    return Tile.from_pyo3_tile(tiles[0].grid, pyo3_tile)
+    return Tile.from_pyo3_tile(tiles[0].get_grid(), pyo3_tile)
 
 
 def count_tiles(tiles):
@@ -1426,11 +1513,13 @@ def count_tiles(tiles):
             )
     if pyo3_tiles:
         pyo3_data_tile = tile_utils.count_tiles(pyo3_tiles)
-        result_tiles_only = DataTile.from_pyo3_data_tile(tiles[0].grid, pyo3_data_tile)
+        result_tiles_only = DataTile.from_pyo3_data_tile(
+            tiles[0].get_grid(), pyo3_data_tile
+        )
     if pyo3_data_tiles:
         pyo3_data_tile = tile_utils.count_data_tiles(pyo3_data_tiles)
         result_data_tiles_only = DataTile.from_pyo3_data_tile(
-            tiles[0].grid, pyo3_data_tile
+            tiles[0].get_grid(), pyo3_data_tile
         )
 
     if pyo3_tiles and pyo3_data_tiles:
@@ -1467,8 +1556,8 @@ def sum_data_tiles(data_tiles: List):
     if len(data_tiles) == 0:
         raise ValueError("No data tiles were supplied")
 
-    grid = data_tiles[0].grid
-    if not all([t.grid.is_aligned_with(grid) for t in data_tiles]):
+    grid = data_tiles[0].get_grid()
+    if not all([t.get_grid().is_aligned_with(grid) for t in data_tiles]):
         raise AlignmentError(
             "Not all data tiles are on the same grid. Consider resampling them all to the same grid."
         )
diff --git a/gridkit-py/gridkit/tri_grid.py b/gridkit-py/gridkit/tri_grid.py
index 3027982d4..2184b9945 100644
--- a/gridkit-py/gridkit/tri_grid.py
+++ b/gridkit-py/gridkit/tri_grid.py
@@ -96,7 +96,11 @@ def __init__(
     @property
     def definition(self):
         return dict(
-            size=self.size, offset=self.offset, rotation=self.rotation, crs=self.crs
+            size=self.size,
+            orientation=self.orientation,
+            offset=self.offset,
+            rotation=self.rotation,
+            crs=self.crs,
         )
 
     @property
diff --git a/gridkit-py/src/lib.rs b/gridkit-py/src/lib.rs
index eb8a7636d..c333c41bc 100644
--- a/gridkit-py/src/lib.rs
+++ b/gridkit-py/src/lib.rs
@@ -1099,6 +1099,10 @@ impl PyO3TriGrid {
             )
             .into_pyarray(py)
     }
+
+    fn equals<'py>(&self, py: Python<'py>, other: PyO3TriGrid) -> bool {
+        self._grid == other._grid
+    }
 }
 
 #[derive(Clone, PartialEq, Hash)]
@@ -1213,6 +1217,10 @@ impl PyO3RectGrid {
             .cells_near_point(&point.as_array())
             .into_pyarray(py)
     }
+
+    fn equals<'py>(&self, py: Python<'py>, other: PyO3RectGrid) -> bool {
+        self._grid == other._grid
+    }
 }
 
 #[derive(Clone, PartialEq, Hash)]
@@ -1361,6 +1369,10 @@ impl PyO3HexGrid {
             .cells_near_point(&point.as_array())
             .into_pyarray(py)
     }
+
+    fn equals<'py>(&self, py: Python<'py>, other: PyO3HexGrid) -> bool {
+        self._grid == other._grid
+    }
 }
 
 #[pyfunction]
diff --git a/gridkit-py/tests/test_gridkit/test_data_tile.py b/gridkit-py/tests/test_gridkit/test_data_tile.py
index e9aea2257..6e68cd93f 100644
--- a/gridkit-py/tests/test_gridkit/test_data_tile.py
+++ b/gridkit-py/tests/test_gridkit/test_data_tile.py
@@ -98,7 +98,7 @@ def test_intersects(grid, tile_init, expected):
 
     with pytest.raises(TypeError):
         # Check if appropriate error is raised when non-tile is supplied
-        tile.intersects(tile.grid)
+        tile.intersects(tile.get_grid())
 
 
 @pytest.mark.parametrize(
@@ -588,7 +588,7 @@ def test_overlap(grid):
 
     def check(overlap_tile):
         assert isinstance(overlap_tile, Tile)
-        assert isinstance(overlap_tile.grid, type(grid))
+        assert isinstance(overlap_tile.get_grid(), type(grid))
         numpy.testing.assert_allclose(overlap_tile.start_id, tile2.start_id)
         assert overlap_tile.nx == 1
         assert overlap_tile.ny == 1
@@ -704,8 +704,8 @@ def test_resample(grid, interp_method):
     # For each nodata value make sure there is at least one neighbour that is outside of the original tile
     for id, value in zip(result.indices, result.to_numpy().ravel()):
         if result.is_nodata(value):
-            original_cells_near_nodata_cell = tile.grid.cells_near_point(
-                result.grid.centroid(id)
+            original_cells_near_nodata_cell = tile.get_grid().cells_near_point(
+                result.get_grid().centroid(id)
             )
             # check that at least one cell is outside the original tile to warrent the nodata value
             cells_otuside_orignal_tile = original_cells_near_nodata_cell.difference(
@@ -745,13 +745,13 @@ def test_from_bounds_as_rect(nodata_value):
     ]
     numpy.testing.assert_allclose(data_tile.nx, bounded_grid.width)
     numpy.testing.assert_allclose(data_tile.ny, bounded_grid.height)
-    numpy.testing.assert_allclose(data_tile.grid.dx, bounded_grid.dx)
-    numpy.testing.assert_allclose(data_tile.grid.dy, bounded_grid.dy)
-    numpy.testing.assert_allclose(data_tile.grid.offset, bounded_grid.offset)
+    numpy.testing.assert_allclose(data_tile.get_grid().dx, bounded_grid.dx)
+    numpy.testing.assert_allclose(data_tile.get_grid().dy, bounded_grid.dy)
+    numpy.testing.assert_allclose(data_tile.get_grid().offset, bounded_grid.offset)
     numpy.testing.assert_allclose(data_tile.to_numpy(), bounded_grid.data)
     numpy.testing.assert_allclose(bounded_grid.bounds, bounds)
     numpy.testing.assert_allclose(data_tile_bounds, bounds)
-    assert data_tile.grid.rotation == bounded_grid.rotation == 0
+    assert data_tile.get_grid().rotation == bounded_grid.rotation == 0
 
     if nodata_value is None:
         # Current nodata default for integer, likely to change in future
diff --git a/gridkit-py/tests/test_gridkit/test_io.py b/gridkit-py/tests/test_gridkit/test_io.py
index f0f136d08..938a2a568 100644
--- a/gridkit-py/tests/test_gridkit/test_io.py
+++ b/gridkit-py/tests/test_gridkit/test_io.py
@@ -31,10 +31,12 @@ def test_write_raster(tmp_path, basic_bounded_rect_grid):
 
     result = raster_to_data_tile(filepath)
     numpy.testing.assert_allclose(result, basic_bounded_rect_grid.data)
-    numpy.testing.assert_allclose(result.grid.offset, basic_bounded_rect_grid.offset)
-    assert result.grid.rotation == 0
-    numpy.testing.assert_allclose(result.grid.dx, basic_bounded_rect_grid.dx)
-    numpy.testing.assert_allclose(result.grid.dy, basic_bounded_rect_grid.dy)
+    numpy.testing.assert_allclose(
+        result.get_grid().offset, basic_bounded_rect_grid.offset
+    )
+    assert result.get_grid().rotation == 0
+    numpy.testing.assert_allclose(result.get_grid().dx, basic_bounded_rect_grid.dx)
+    numpy.testing.assert_allclose(result.get_grid().dy, basic_bounded_rect_grid.dy)
     b = basic_bounded_rect_grid.bounds
     corners = [
         [b[0], b[3]],
@@ -58,10 +60,10 @@ def test_write_raster_data_tile(tmp_path):
 
     result = raster_to_data_tile(filepath)
     numpy.testing.assert_allclose(result, data)
-    numpy.testing.assert_allclose(result.grid.offset, grid.offset)
-    assert result.grid.rotation == 0
-    numpy.testing.assert_allclose(result.grid.dx, grid.dx)
-    numpy.testing.assert_allclose(result.grid.dy, grid.dy)
+    numpy.testing.assert_allclose(result.get_grid().offset, grid.offset)
+    assert result.get_grid().rotation == 0
+    numpy.testing.assert_allclose(result.get_grid().dx, grid.dx)
+    numpy.testing.assert_allclose(result.get_grid().dy, grid.dy)
     numpy.testing.assert_allclose(result.start_id, tile.start_id)
 
 
@@ -76,6 +78,6 @@ def test_raster_to_data_tile():
         [4116200.0, 2575600.0],
     ]
     numpy.testing.assert_allclose(grid_tile.corners(), expected_corners)
-    assert grid_tile.grid.dx == 100.0
-    assert grid_tile.grid.dy == 100.0
-    assert grid_tile.grid.crs.to_epsg() == 3035
+    assert grid_tile.get_grid().dx == 100.0
+    assert grid_tile.get_grid().dy == 100.0
+    assert grid_tile.get_grid().crs.to_epsg() == 3035
diff --git a/gridkit-py/tests/test_gridkit/test_tile.py b/gridkit-py/tests/test_gridkit/test_tile.py
index 39b81abdf..85ee4824c 100644
--- a/gridkit-py/tests/test_gridkit/test_tile.py
+++ b/gridkit-py/tests/test_gridkit/test_tile.py
@@ -204,3 +204,50 @@ def test_tile_id_conversion_back_and_forth():
     numpy.testing.assert_allclose(ids.ravel(), ids_reverted)
     numpy.testing.assert_allclose(data_tile.value(tile.indices), data)
     numpy.testing.assert_allclose(data_tile.value(data_tile.indices), data)
+
+
+@pytest.mark.parametrize(
+    "grid",
+    [
+        TriGrid(size=1, orientation="flat"),
+        TriGrid(size=1, orientation="pointy"),
+        RectGrid(size=1),
+        HexGrid(size=1, orientation="flat"),
+        HexGrid(size=1, orientation="pointy"),
+    ],
+)
+@pytest.mark.parametrize(
+    "other_grid",
+    [
+        TriGrid(size=1, orientation="flat"),
+        TriGrid(size=1, orientation="pointy"),
+        RectGrid(size=1),
+        HexGrid(size=1, orientation="flat"),
+        HexGrid(size=1, orientation="pointy"),
+    ],
+)
+def test_equality(grid, other_grid):
+    tile = Tile(grid, (-1, -2), 3, 5)
+    other_tile = Tile(other_grid, (-1, -2), 3, 5)
+    if isinstance(other_tile.grid, type(tile.grid)) and (
+        getattr(grid, "orientation", None) == getattr(other_grid, "orientation", None)
+    ):
+        assert tile == other_tile
+    else:
+        assert not tile == other_tile
+        assert tile != other_tile
+
+    different_tiles = (
+        # Now test for same grid but different tile parameters
+        Tile(grid.update(), (-1, -2), 5, 3),
+        Tile(grid.update(), (-1, -2), 4, 5),
+        Tile(grid.update(), (-1, -2), 3, 4),
+        Tile(grid.update(), (1, 2), 3, 4),
+        # Now test for same grid type but different grid parameters
+        tile.update(grid.update(offset=(-0.1, 0.3))),
+        tile.update(grid.update(rotation=-182)),
+        tile.update(grid.update(size=123)),
+    )
+    for different_tile in different_tiles:
+        assert not tile == different_tile
+        assert tile != different_tile

From 9f5f57473e8ad7031fa485cc763bda3ebe98d27d Mon Sep 17 00:00:00 2001
From: Timo Millenaar <timo.m@live.nl>
Date: Thu, 18 Dec 2025 15:24:05 +0100
Subject: [PATCH 07/10] Add intersection_corners and intersection_bounds to
 Tile

---
 gridkit-rs/src/grid.rs |  2 +-
 gridkit-rs/src/tile.rs | 51 +++++++++++++++++++++++++++++++++++++++++-
 2 files changed, 51 insertions(+), 2 deletions(-)

diff --git a/gridkit-rs/src/grid.rs b/gridkit-rs/src/grid.rs
index dbf17372c..002dfe4d4 100644
--- a/gridkit-rs/src/grid.rs
+++ b/gridkit-rs/src/grid.rs
@@ -135,7 +135,7 @@ pub trait GridTraits {
     // ) -> Array1<f64>;
 }
 
-#[derive(Clone, Debug, PartialEq, Hash)]
+#[derive(Clone, Debug, Eq, PartialEq, Hash)]
 #[enum_delegate::implement(GridTraits)]
 pub enum Grid {
     TriGrid(TriGrid),
diff --git a/gridkit-rs/src/tile.rs b/gridkit-rs/src/tile.rs
index 750476f76..02b75f9b8 100644
--- a/gridkit-rs/src/tile.rs
+++ b/gridkit-rs/src/tile.rs
@@ -4,7 +4,7 @@ use std::{f32::MAX, f64, i64, u64};
 use crate::{data_tile::DataTile, grid::*, hex_grid::HexGrid};
 use ndarray::*;
 
-#[derive(Clone, Debug, PartialEq, Hash)]
+#[derive(Clone, Debug, Eq, PartialEq, Hash)]
 pub struct Tile {
     pub grid: Grid,
     pub start_id: (i64, i64),
@@ -52,6 +52,14 @@ pub trait TileTraits {
         self.get_tile().intersects(other)
     }
 
+    fn intersection_bounds(&self, other: &Tile) -> Option<(f64, f64, f64, f64)> {
+        self.get_tile().intersection_bounds(other)
+    }
+
+    fn intersection_corners(&self, other: &Tile) -> Option<Array2<f64>> {
+        self.get_tile().intersection_corners(other)
+    }
+
     fn overlap(&self, other: &Tile) -> Result<Tile, String> {
         self.get_tile().overlap(other)
     }
@@ -248,6 +256,47 @@ impl TileTraits for Tile {
         );
     }
 
+    fn intersection_bounds(&self, other: &Tile) -> Option<(f64, f64, f64, f64)> {
+        let corners_self = self.corners();
+        let corners_other = other.corners();
+
+        let left_self = corners_self.slice(s![.., 0]).iter().cloned().reduce(f64::min).unwrap();
+        let bottom_self = corners_self.slice(s![.., 1]).iter().cloned().reduce(f64::min).unwrap();
+        let right_self = corners_self.slice(s![.., 0]).iter().cloned().reduce(f64::max).unwrap();
+        let top_self = corners_self.slice(s![.., 1]).iter().cloned().reduce(f64::max).unwrap();
+
+        let left_other = corners_other.slice(s![.., 0]).iter().cloned().reduce(f64::min).unwrap();
+        let bottom_other = corners_other.slice(s![.., 1]).iter().cloned().reduce(f64::min).unwrap();
+        let right_other = corners_other.slice(s![.., 0]).iter().cloned().reduce(f64::max).unwrap();
+        let top_other = corners_other.slice(s![.., 1]).iter().cloned().reduce(f64::max).unwrap();
+
+        // Compute intersection bounds
+        let minx = left_self.max(left_other);
+        let miny = bottom_self.max(bottom_other);
+        let maxx = right_self.min(right_other);
+        let maxy = top_self.min(top_other);
+
+        // Check if there is an intersection
+        if minx < maxx && miny < maxy {
+            Some((minx, miny, maxx, maxy))
+        } else {
+            None
+        }
+    }
+
+    fn intersection_corners(&self, other: &Tile) -> Option<Array2<f64>> {
+        if let Some((minx, miny, maxx, maxy)) = self.intersection_bounds(other) {
+            Some(arr2(&[
+                [minx, maxy],
+                [maxx, maxy],
+                [maxx, miny],
+                [minx, miny],
+            ]))
+        } else {
+            None
+        }
+    }
+
     fn overlap(&self, other: &Tile) -> Result<Tile, String> {
         // TODO: check whether grids align
         if !self.intersects(&other) {

From 4c2fc4ce2dff6474723e171b85191d01ee3230a7 Mon Sep 17 00:00:00 2001
From: Timo Millenaar <tmillenaar@gmail.com>
Date: Sun, 18 Jan 2026 22:32:40 +0100
Subject: [PATCH 08/10] Fix .grid call that escaped replacing to .get_grid()

---
 .../cell_centric_operations/downslope_path.py      |  4 ++--
 gridkit-py/examples/patterns/game_of_life.py       |  8 ++++----
 .../coordinate_transformations.py                  |  2 +-
 gridkit-py/examples/raster_operations/ndvi.py      |  2 +-
 .../examples/raster_operations/resampling.py       | 14 +++++++++-----
 gridkit-py/gridkit/tile.py                         |  2 +-
 6 files changed, 18 insertions(+), 14 deletions(-)

diff --git a/gridkit-py/examples/cell_centric_operations/downslope_path.py b/gridkit-py/examples/cell_centric_operations/downslope_path.py
index fae71f2b7..a6f6810b8 100644
--- a/gridkit-py/examples/cell_centric_operations/downslope_path.py
+++ b/gridkit-py/examples/cell_centric_operations/downslope_path.py
@@ -58,7 +58,7 @@
 def downslope_path(cell_id):
     visited_cells.append(cell_id)
     # add new neigbours and their values to their respective lists
-    neighbours = dem.grid.neighbours(cell_id, connect_corners=True)
+    neighbours = dem.get_grid().neighbours(cell_id, connect_corners=True)
     for cell, value in zip(neighbours, dem.value(neighbours)):
         cell = tuple(cell.index)
         if not cell in visited_cells and not cell in all_neighbours:
@@ -80,7 +80,7 @@ def downslope_path(cell_id):
 
 # %%
 # Determine the starting cell and call the function
-start_cell_id = dem.grid.cell_at_point((29330, 167848))
+start_cell_id = dem.get_grid().cell_at_point((29330, 167848))
 downslope_path(tuple(start_cell_id.index))
 
 
diff --git a/gridkit-py/examples/patterns/game_of_life.py b/gridkit-py/examples/patterns/game_of_life.py
index cf95d9c46..449a64886 100644
--- a/gridkit-py/examples/patterns/game_of_life.py
+++ b/gridkit-py/examples/patterns/game_of_life.py
@@ -13,10 +13,10 @@
 
 The rules as taken from https://en.wikipedia.org/wiki/Conway%27s_Game_of_Life:
 
- 1. Any live cell with fewer than two live neighbours dies, as if by underpopulation.
- 2. Any live cell with two or three live neighbours lives on to the next generation.
- 3. Any live cell with more than three live neighbours dies, as if by overpopulation.
- 4. Any dead cell with exactly three live neighbours becomes a live cell, as if by reproduction.
+1. Any live cell with fewer than two live neighbours dies, as if by underpopulation.
+2. Any live cell with two or three live neighbours lives on to the next generation.
+3. Any live cell with more than three live neighbours dies, as if by overpopulation.
+4. Any dead cell with exactly three live neighbours becomes a live cell, as if by reproduction.
 
 Initial conditions
 ------------------
diff --git a/gridkit-py/examples/raster_operations/coordinate_transformations.py b/gridkit-py/examples/raster_operations/coordinate_transformations.py
index a927c49c0..064d8c255 100644
--- a/gridkit-py/examples/raster_operations/coordinate_transformations.py
+++ b/gridkit-py/examples/raster_operations/coordinate_transformations.py
@@ -56,7 +56,7 @@
 
 utm_epsg = 32629
 transformer = Transformer.from_crs(
-    data_tile_wgs84.grid.crs, CRS.from_user_input(utm_epsg), always_xy=True
+    data_tile_wgs84.get_grid().crs, CRS.from_user_input(utm_epsg), always_xy=True
 )  # UTM zone 28N
 
 wgs84_centroids_utm = transformer.transform(
diff --git a/gridkit-py/examples/raster_operations/ndvi.py b/gridkit-py/examples/raster_operations/ndvi.py
index 2c14aec5d..87a92d204 100644
--- a/gridkit-py/examples/raster_operations/ndvi.py
+++ b/gridkit-py/examples/raster_operations/ndvi.py
@@ -61,6 +61,6 @@
 fig.colorbar(im_ndvi, ax=ax, fraction=0.022, pad=0.01)
 ax.set_xlabel("lon")
 ax.set_ylabel("lat")
-ax.set_title(f"NDVI of scene in the alps \n EPSG:{ndvi.grid.crs.to_epsg()}")
+ax.set_title(f"NDVI of scene in the alps \n EPSG:{ndvi.get_grid().crs.to_epsg()}")
 
 plt.show()
diff --git a/gridkit-py/examples/raster_operations/resampling.py b/gridkit-py/examples/raster_operations/resampling.py
index 06f314c38..28815949e 100644
--- a/gridkit-py/examples/raster_operations/resampling.py
+++ b/gridkit-py/examples/raster_operations/resampling.py
@@ -25,7 +25,7 @@
 dem = raster_to_data_tile(
     "../../tests/data/alps_dem.tiff", bounds=(28300, 167300, 28700, 167700)
 )
-print(f"Original resolution in (dx, dy): ({dem.grid.dx}, {dem.grid.dy})")
+print(f"Original resolution in (dx, dy): ({dem.get_grid().dx}, {dem.get_grid().dy})")
 
 
 # %%
@@ -49,8 +49,10 @@
 # Let's define a rectangular grid with a cell size of 10x10 degrees.
 # I am calling it rectdem for later on I will define hexagonal ones as well.
 # To make sure it worked we can print out the cellsize after resampling
-rectdem = dem.resample(RectGrid(dx=10, dy=10, crs=dem.grid.crs), method="linear")
-print(f"Downsampled resolution in (dx, dy): ({rectdem.grid.dx}, {rectdem.grid.dy})")
+rectdem = dem.resample(RectGrid(dx=10, dy=10, crs=dem.get_grid().crs), method="linear")
+print(
+    f"Downsampled resolution in (dx, dy): ({rectdem.get_grid().dx}, {rectdem.get_grid().dy})"
+)
 
 # %%
 #
@@ -71,10 +73,12 @@
 # for a fair visual comparison.
 #
 hexdem_flat = dem.resample(
-    HexGrid(area=rectdem.grid.area, shape="flat", crs=dem.grid.crs), method="linear"
+    HexGrid(area=rectdem.get_grid().area, shape="flat", crs=dem.get_grid().crs),
+    method="linear",
 )
 hexdem_pointy = dem.resample(
-    HexGrid(area=rectdem.grid.area, shape="pointy", crs=dem.grid.crs), method="linear"
+    HexGrid(area=rectdem.get_grid().area, shape="pointy", crs=dem.get_grid().crs),
+    method="linear",
 )
 
 
diff --git a/gridkit-py/gridkit/tile.py b/gridkit-py/gridkit/tile.py
index 6799fff2d..f064c337e 100644
--- a/gridkit-py/gridkit/tile.py
+++ b/gridkit-py/gridkit/tile.py
@@ -995,7 +995,7 @@ def resample(self, alignment_grid, method="nearest", **interp_kwargs):
         tile_is_given = isinstance(alignment_grid, Tile)
         if tile_is_given:
             tile = alignment_grid
-            alignment_grid = alignment_grid.grid
+            alignment_grid = alignment_grid.get_grid()
 
         if self.get_grid().crs is None or alignment_grid.crs is None:
             warnings.warn(

From 1a475b37b1b0d6cdad6bc0a6ca6c95d24a192fd7 Mon Sep 17 00:00:00 2001
From: Timo Millenaar <tmillenaar@gmail.com>
Date: Sun, 18 Jan 2026 23:47:43 +0100
Subject: [PATCH 09/10] Fix aggregate_dask.py now that Dask needed the Grid
 object to be serializibale

---
 gridkit-py/examples/patterns/game_of_life.py       |  6 +++---
 .../examples/vector_synergy/aggregate_dask.py      | 14 ++++++++++----
 2 files changed, 13 insertions(+), 7 deletions(-)

diff --git a/gridkit-py/examples/patterns/game_of_life.py b/gridkit-py/examples/patterns/game_of_life.py
index 449a64886..5d0b44894 100644
--- a/gridkit-py/examples/patterns/game_of_life.py
+++ b/gridkit-py/examples/patterns/game_of_life.py
@@ -7,9 +7,9 @@
 Cellular automata are a small set of rules relating to cells on a grid that can create complex patterns when applied repeatedly.
 There is a near infinite veriety of rules, but far and away the most well known set of rules is that of the so called Conway's game of life.
 The setup is as follows:
- - The board is a 2d rectangular grid
- - The 8 neighbouring cells are considered, so this includes the diagonal neighbours
- - A cell has two possbile states: alive or dead
+- The board is a 2d rectangular grid
+- The 8 neighbouring cells are considered, so this includes the diagonal neighbours
+- A cell has two possbile states: alive or dead
 
 The rules as taken from https://en.wikipedia.org/wiki/Conway%27s_Game_of_Life:
 
diff --git a/gridkit-py/examples/vector_synergy/aggregate_dask.py b/gridkit-py/examples/vector_synergy/aggregate_dask.py
index 157229fc9..44753871b 100644
--- a/gridkit-py/examples/vector_synergy/aggregate_dask.py
+++ b/gridkit-py/examples/vector_synergy/aggregate_dask.py
@@ -88,14 +88,16 @@
 grid = HexGrid(size=1, shape="flat")
 
 
-def index_1d_at_point(df):
+def index_1d_at_point(df, **grid_params):
     """Find the cell for each point in the dataframe and return the 1D index.
     This operation is meant to map to each DataFrame partition
     """
+    # Note: Define 'grid' here based on primitives so Dask does not have to worry about serializing the Grid object.
+    grid = HexGrid(**grid_params)
     return grid.cell_at_point(df[["pnt_x", "pnt_y"]]).index_1d
 
 
-df["cell_id"] = df.map_partitions(index_1d_at_point)
+df["cell_id"] = df.map_partitions(index_1d_at_point, **grid.definition)
 df["nr_points"] = 1
 print(df)
 
@@ -121,14 +123,18 @@ def index_1d_at_point(df):
 #     Also, for the sake of the example I will plot the data with matplotlib.
 
 
-def shapely_geom_from_index_1d(id_1d):
+def shapely_geom_from_index_1d(id_1d, **grid_params):
     """Generate the Shapely geometry for each 1D index.
     This operation is meant to map to each DataFrame partition
     """
+    # Note: Define 'grid' here based on primitives so Dask does not have to worry about serializing the Grid object.
+    grid = HexGrid(**grid_params)
     return numpy.array(grid.to_shapely(GridIndex.from_index_1d(id_1d)).geoms)
 
 
-polygons = occurrences.index.map_partitions(shapely_geom_from_index_1d)
+polygons = occurrences.index.map_partitions(
+    shapely_geom_from_index_1d, **grid.definition
+)
 geoms, points_per_cell = dask.compute(polygons, occurrences.nr_points)
 
 # %%

From eb0d3d0c782f6434766f0eabba54566910765838 Mon Sep 17 00:00:00 2001
From: Timo Millenaar <tmillenaar@gmail.com>
Date: Sun, 18 Jan 2026 23:51:22 +0100
Subject: [PATCH 10/10] black reformatting

---
 .pre-commit-config.yaml    | 2 +-
 gridkit-py/gridkit/tile.py | 6 ++----
 2 files changed, 3 insertions(+), 5 deletions(-)

diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index dc323712d..8f51ce6e4 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,6 +1,6 @@
 repos:
   - repo: https://github.com/psf/black
-    rev: 25.1.0 # Replace by any tag/version: https://github.com/psf/black/tags
+    rev: 26.1.0 # Replace by any tag/version: https://github.com/psf/black/tags
     hooks:
       - id: black
         language_version: python3
diff --git a/gridkit-py/gridkit/tile.py b/gridkit-py/gridkit/tile.py
index f064c337e..e3595e73e 100644
--- a/gridkit-py/gridkit/tile.py
+++ b/gridkit-py/gridkit/tile.py
@@ -386,16 +386,14 @@ def tile_id_to_grid_id(self, tile_id=None, oob_value=numpy.iinfo(numpy.int64).ma
         else:
             tile_id = numpy.array(tile_id, dtype=int)
         if not tile_id.shape[0] == 2:
-            raise ValueError(
-                f"""
+            raise ValueError(f"""
                 Expected the first dimension of tile_id to have a length of two (for x,y). Got indices in shape: {tile_id.shape}.
                 This is different than the grid ids where we expect xy to be the last dimesnion.
                 This is in an effort to make indexing the data as a numpy array using data[tuple(numpy_ids)] the same
                 as querying values from datatile.value using grid ids botained through tile_id_to_grid_id(numpy_ids).
                 Note that not only does this function need the first demension to be of size two,
                 but numpy starts with the y-coordinate, so we expect: [[y0,y1,y2], [x0,x1,x2]].
-            """
-            )
+            """)
         if tile_id.ndim == 1:
             tile_id = tile_id[numpy.newaxis]
         else: