diff --git a/doc/source/vertical_ref.md b/doc/source/vertical_ref.md
index 96f965f8..76fdba5e 100644
--- a/doc/source/vertical_ref.md
+++ b/doc/source/vertical_ref.md
@@ -224,17 +224,22 @@ xDEM uses only the best available (i.e. best accuracy) transformation returned b
 # Open a DEM and set its CRS
 filename_dem = xdem.examples.get_path("longyearbyen_ref_dem")
 dem = xdem.DEM(filename_dem, vcrs="Ellipsoid")
-dem.to_vcrs("EGM96")
-dem.vcrs
+trans_dem = dem.to_vcrs("EGM96")
+trans_dem.vcrs
 ```
 
-The operation updates the DEM array **in-place**, shifting each pixel by the transformation at their coordinates:
+The operation returns a new {class}`~xdem.DEM` by default, but can also be done in-place. It vertically shifts
+each pixel value by the transformation at their coordinates:
 
 ```{code-cell} ipython3
 import numpy as np
 
+diff = trans_dem - dem
 # Mean difference after transformation (about 30 m in Svalbard)
-dem_orig = xdem.DEM(filename_dem)
-diff = dem_orig - dem
 np.nanmean(diff)
 ```
+
+```{code-cell} ipython3
+# Plot the elevation differences due to the vertical transformation
+diff.plot(cmap="RdYlBu", cbar_title="Elevation differences (m)")
+```
diff --git a/tests/test_dem.py b/tests/test_dem.py
index 578eaafd..536e5d3d 100644
--- a/tests/test_dem.py
+++ b/tests/test_dem.py
@@ -273,9 +273,20 @@ def test_to_vcrs(self) -> None:
         dem.set_vcrs(new_vcrs="Ellipsoid")
         ccrs_init = dem.ccrs
         median_before = np.nanmean(dem)
-        # Transform to EGM96 geoid
-        dem.to_vcrs(vcrs="EGM96")
-        median_after = np.nanmean(dem)
+        # Transform to EGM96 geoid not inplace (default)
+        trans_dem = dem.to_vcrs(vcrs="EGM96")
+
+        # The output should be a DEM, input shouldn't have changed
+        assert isinstance(trans_dem, DEM)
+        assert dem.raster_equal(dem_before_trans)
+
+        # Compare to inplace
+        should_be_none = dem.to_vcrs(vcrs="EGM96", inplace=True)
+        assert should_be_none is None
+        assert dem.raster_equal(trans_dem)
+
+        # Save the median of after
+        median_after = np.nanmean(trans_dem)
 
         # About 32 meters of difference in Svalbard between EGM96 geoid and ellipsoid
         assert median_after - median_before == pytest.approx(-32, rel=0.1)
@@ -337,10 +348,10 @@ def test_to_vcrs__grids(self, grid_shifts: dict[str, Any]) -> None:
         dem.set_vcrs("Ellipsoid")
 
         # Transform to the vertical CRS of the grid
-        dem.to_vcrs(grid_shifts["grid"])
+        trans_dem = dem.to_vcrs(grid_shifts["grid"])
 
         # Compare the elevation difference
-        z_diff = 100 - dem.data[0, 0]
+        z_diff = 100 - trans_dem.data[0, 0]
 
         # Check the shift is the one expect within 10%
         assert z_diff == pytest.approx(grid_shifts["shift"], rel=0.1)
diff --git a/xdem/dem.py b/xdem/dem.py
index 697859a2..c235b897 100644
--- a/xdem/dem.py
+++ b/xdem/dem.py
@@ -3,7 +3,7 @@
 
 import pathlib
 import warnings
-from typing import Any, Literal
+from typing import Any, Literal, overload
 
 import geopandas as gpd
 import numpy as np
@@ -254,6 +254,7 @@ def ccrs(self) -> CompoundCRS | CRS | None:
         else:
             return None
 
+    @overload
     def to_vcrs(
         self,
         vcrs: Literal["Ellipsoid", "EGM08", "EGM96"] | str | pathlib.Path | VerticalCRS | int,
@@ -263,15 +264,61 @@ def to_vcrs(
         | VerticalCRS
         | int
         | None = None,
+        *,
+        inplace: Literal[False] = False,
+    ) -> DEM:
+        ...
+
+    @overload
+    def to_vcrs(
+        self,
+        vcrs: Literal["Ellipsoid", "EGM08", "EGM96"] | str | pathlib.Path | VerticalCRS | int,
+        force_source_vcrs: Literal["Ellipsoid", "EGM08", "EGM96"]
+        | str
+        | pathlib.Path
+        | VerticalCRS
+        | int
+        | None = None,
+        *,
+        inplace: Literal[True],
     ) -> None:
+        ...
+
+    @overload
+    def to_vcrs(
+        self,
+        vcrs: Literal["Ellipsoid", "EGM08", "EGM96"] | str | pathlib.Path | VerticalCRS | int,
+        force_source_vcrs: Literal["Ellipsoid", "EGM08", "EGM96"]
+        | str
+        | pathlib.Path
+        | VerticalCRS
+        | int
+        | None = None,
+        *,
+        inplace: bool = False,
+    ) -> DEM | None:
+        ...
+
+    def to_vcrs(
+        self,
+        vcrs: Literal["Ellipsoid", "EGM08", "EGM96"] | str | pathlib.Path | VerticalCRS | int,
+        force_source_vcrs: Literal["Ellipsoid", "EGM08", "EGM96"]
+        | str
+        | pathlib.Path
+        | VerticalCRS
+        | int
+        | None = None,
+        inplace: bool = False,
+    ) -> DEM | None:
         """
         Convert the DEM to another vertical coordinate reference system.
 
         :param vcrs: Destination vertical CRS. Either as a name ("WGS84", "EGM08", "EGM96"),
             an EPSG code or pyproj.crs.VerticalCRS, or a path to a PROJ grid file (https://github.com/OSGeo/PROJ-data)
         :param force_source_vcrs: Force a source vertical CRS (uses metadata by default). Same formats as for `vcrs`.
+        :param inplace: Whether to return a new DEM (default) or the same DEM updated in-place.
 
-        :return:
+        :return: DEM with vertical reference transformed, or None.
         """
 
         if self.vcrs is None and force_source_vcrs is None:
@@ -307,12 +354,25 @@ def to_vcrs(
         zz = self.data
         xx, yy = self.coords()
         zz_trans = _transform_zz(crs_from=src_ccrs, crs_to=dst_ccrs, xx=xx, yy=yy, zz=zz)
+        new_data = zz_trans.astype(self.dtype)  # type: ignore
 
-        # Update DEM
-        self._data = zz_trans.astype(self.dtype)  # type: ignore
-
-        # Update vcrs (which will update ccrs if called)
-        self.set_vcrs(new_vcrs=vcrs)
+        # If inplace, update DEM and vcrs
+        if inplace:
+            self._data = new_data
+            self.set_vcrs(new_vcrs=vcrs)
+            return None
+        # Otherwise, return new DEM
+        else:
+            return DEM.from_array(
+                data=new_data,
+                transform=self.transform,
+                crs=self.crs,
+                nodata=self.nodata,
+                area_or_point=self.area_or_point,
+                tags=self.tags,
+                vcrs=vcrs,
+                cast_nodata=False,
+            )
 
     @copy_doc(terrain, remove_dem_res_params=True)
     def slope(