diff --git a/python/ribasim/ribasim/model.py b/python/ribasim/ribasim/model.py
index ccca20892..864d17116 100644
--- a/python/ribasim/ribasim/model.py
+++ b/python/ribasim/ribasim/model.py
@@ -161,6 +161,7 @@ def _write_toml(self, fn: Path) -> Path:
         return fn
 
     def _save(self, directory: DirectoryPath, input_dir: DirectoryPath):
+        # Set CRS of the tables to the CRS stored in the Model object
         self.set_crs(self.crs)
         db_path = directory / input_dir / "database.gpkg"
         db_path.parent.mkdir(parents=True, exist_ok=True)
@@ -182,14 +183,23 @@ def _save(self, directory: DirectoryPath, input_dir: DirectoryPath):
             sub._save(directory, input_dir)
 
     def set_crs(self, crs: str) -> None:
-        self.crs = crs
-        self.edge.df = self.edge.df.set_crs(crs)
+        self._apply_crs_function("set_crs", crs)
+
+    def to_crs(self, crs: str) -> None:
+        # Set CRS of the tables to the CRS stored in the Model object
+        self.set_crs(self.crs)
+        self._apply_crs_function("to_crs", crs)
+
+    def _apply_crs_function(self, function_name: str, crs: str) -> None:
+        """Apply `function_name`, with `crs` as the first and only argument to all spatial tables."""
+        self.edge.df = getattr(self.edge.df, function_name)(crs)
         for sub in self._nodes():
             if sub.node.df is not None:
-                sub.node.df = sub.node.df.set_crs(crs)
+                sub.node.df = getattr(sub.node.df, function_name)(crs)
             for table in sub._tables():
                 if isinstance(table, SpatialTableModel) and table.df is not None:
-                    table.df = table.df.set_crs(crs)
+                    table.df = getattr(table.df, function_name)(crs)
+        self.crs = crs
 
     def node_table(self) -> NodeTable:
         """Compute the full NodeTable from all node types."""
diff --git a/python/ribasim/tests/conftest.py b/python/ribasim/tests/conftest.py
index 2f0f3ae9f..30939e521 100644
--- a/python/ribasim/tests/conftest.py
+++ b/python/ribasim/tests/conftest.py
@@ -19,6 +19,11 @@ def basic_transient() -> ribasim.Model:
     return ribasim_testmodels.basic_transient_model()
 
 
+@pytest.fixture()
+def bucket() -> ribasim.Model:
+    return ribasim_testmodels.bucket_model()
+
+
 @pytest.fixture()
 def tabulated_rating_curve() -> ribasim.Model:
     return ribasim_testmodels.tabulated_rating_curve_model()
diff --git a/python/ribasim/tests/test_model.py b/python/ribasim/tests/test_model.py
index 001ab889f..d7311b63d 100644
--- a/python/ribasim/tests/test_model.py
+++ b/python/ribasim/tests/test_model.py
@@ -242,3 +242,14 @@ def test_xugrid(basic, tmp_path):
     uds.ugrid.to_netcdf(tmp_path / "ribasim.nc")
     uds = xugrid.open_dataset(tmp_path / "ribasim.nc")
     assert uds.attrs["Conventions"] == "CF-1.9 UGRID-1.0"
+
+
+def test_to_crs(bucket: Model):
+    model = bucket
+
+    # Reproject to World Geodetic System 1984
+    model.to_crs("EPSG:4326")
+
+    # Assert that the bucket is still at Deltares' headquarter
+    assert model.basin.node.df["geometry"].iloc[0].x == pytest.approx(4.38, abs=0.1)
+    assert model.basin.node.df["geometry"].iloc[0].y == pytest.approx(51.98, abs=0.1)
diff --git a/python/ribasim_testmodels/ribasim_testmodels/bucket.py b/python/ribasim_testmodels/ribasim_testmodels/bucket.py
index 66ffe2abb..87fb28e07 100644
--- a/python/ribasim_testmodels/ribasim_testmodels/bucket.py
+++ b/python/ribasim_testmodels/ribasim_testmodels/bucket.py
@@ -9,7 +9,7 @@
 
 
 def bucket_model() -> ribasim.Model:
-    """Bucket model with just a single basin."""
+    """Bucket model with just a single basin at Deltares' headquarter."""
 
     model = ribasim.Model(
         starttime="2020-01-01",
@@ -18,7 +18,7 @@ def bucket_model() -> ribasim.Model:
     )
 
     model.basin.add(
-        Node(1, Point(400.0, 200.0)),
+        Node(1, Point(85825.6, 444613.9)),
         [
             basin.Profile(
                 area=[1000.0, 1000.0],
@@ -38,7 +38,7 @@ def bucket_model() -> ribasim.Model:
 
 
 def leaky_bucket_model() -> ribasim.Model:
-    """Bucket model with dynamic forcing with missings."""
+    """Bucket model with dynamic forcing with missings at Deltares' headquarter."""
 
     model = ribasim.Model(
         starttime="2020-01-01",
@@ -47,7 +47,7 @@ def leaky_bucket_model() -> ribasim.Model:
     )
 
     model.basin.add(
-        Node(1, Point(400.0, 200.0)),
+        Node(1, Point(85825.6, 444613.9)),
         [
             basin.Profile(
                 area=[1000.0, 1000.0],