diff --git a/geoutils/interface/raster_point.py b/geoutils/interface/raster_point.py
index 8d817bd3..e7cc2d18 100644
--- a/geoutils/interface/raster_point.py
+++ b/geoutils/interface/raster_point.py
@@ -100,7 +100,7 @@ def _raster_to_pointcloud(
     as_array: bool,
     random_state: int | np.random.Generator | None,
     force_pixel_offset: Literal["center", "ul", "ur", "ll", "lr"],
-) -> NDArrayNum | gpd.GeoDataFrame:
+) -> NDArrayNum | gu.PointCloud:
     """
     Convert a raster to a point cloud. See Raster.to_pointcloud() for details.
     """
@@ -233,7 +233,7 @@ def _raster_to_pointcloud(
             geometry=gpd.points_from_xy(x_coords_2, y_coords_2),
             crs=source_raster.crs,
         )
-        return pc
+        return gu.PointCloud(pc, data_column=data_column_name)
     else:
         # Merge the coordinates and pixel data an array of N x K
         # This has the downside of converting all the data to the same data type
diff --git a/geoutils/pointcloud/pointcloud.py b/geoutils/pointcloud/pointcloud.py
index 7591a705..ea5451ea 100644
--- a/geoutils/pointcloud/pointcloud.py
+++ b/geoutils/pointcloud/pointcloud.py
@@ -98,7 +98,7 @@ class PointCloud(gu.Vector):  # type: ignore[misc]
     def __init__(
         self,
         filename_or_dataset: str | pathlib.Path | gpd.GeoDataFrame | gpd.GeoSeries | BaseGeometry,
-        data_column: str = "z",
+        data_column: str,
     ):
         """
         Instantiate a point cloud from either a data column name and a vector (filename, GeoPandas dataframe or series,
@@ -113,8 +113,9 @@ def __init__(
         self._crs: CRS | None = None
         self._bounds: BoundingBox
         self._data_column: str
+        self._data: np.ndarray
         self._nb_points: int
-        self._columns: pd.Index
+        self._all_columns: pd.Index
 
         # If PointCloud is passed, simply point back to PointCloud
         if isinstance(filename_or_dataset, PointCloud):
@@ -133,7 +134,7 @@ def __init__(
                 self._name = fn
                 self._crs = crs
                 self._nb_points = nb_points
-                self._columns = columns
+                self._all_columns = columns
                 self._bounds = bounds
                 self._ds = None
             # Check on filename are done with Vector.__init__
@@ -190,6 +191,26 @@ def bounds(self) -> BoundingBox:
         else:
             return self._bounds
 
+    #####################################
+    # NEW METHODS SPECIFIC TO POINT CLOUD
+    #####################################
+
+    @property
+    def data(self) -> NDArrayNum:
+        """
+        Data of the point cloud.
+
+        Points to the data column of the geodataframe, equivalent to calling self.ds[self.data_column].
+        """
+        # Triggers the loading mechanism through self.ds
+        return self.ds[self.data_column]
+
+    @data.setter
+    def data(self, new_data: NDArrayNum) -> None:
+        """Set new data for the point cloud."""
+
+        self.ds[self.data_column] = new_data
+
     @property
     def all_columns(self) -> pd.Index:
         """Index of all columns of the point cloud, excluding the column of 2D point geometries."""
@@ -199,11 +220,7 @@ def all_columns(self) -> pd.Index:
             all_columns_nongeom = all_columns[all_columns != "geometry"]
             return all_columns_nongeom
         else:
-            return self._columns
-
-    #####################################
-    # NEW METHODS SPECIFIC TO POINT CLOUD
-    #####################################
+            return self._all_columns
 
     @property
     def data_column(self) -> str:
diff --git a/geoutils/raster/raster.py b/geoutils/raster/raster.py
index 46cb9427..3d75f9f6 100644
--- a/geoutils/raster/raster.py
+++ b/geoutils/raster/raster.py
@@ -3390,8 +3390,7 @@ def to_pointcloud(
         :returns: A point cloud, or array of the shape (N, 2 + count) where N is the sample count.
         """
 
-        return gu.PointCloud(
-            _raster_to_pointcloud(
+        return _raster_to_pointcloud(
                 source_raster=self,
                 data_column_name=data_column_name,
                 data_band=data_band,
@@ -3401,9 +3400,7 @@ def to_pointcloud(
                 skip_nodata=skip_nodata,
                 as_array=as_array,
                 random_state=random_state,
-                force_pixel_offset=force_pixel_offset,
-            )
-        )
+                force_pixel_offset=force_pixel_offset)
 
     @classmethod
     def from_pointcloud_regular(
diff --git a/tests/test_pointcloud/test_pointcloud.py b/tests/test_pointcloud/test_pointcloud.py
index dadca55b..a57bc9d0 100644
--- a/tests/test_pointcloud/test_pointcloud.py
+++ b/tests/test_pointcloud/test_pointcloud.py
@@ -28,10 +28,10 @@ class TestPointCloud:
         geometry=gpd.points_from_xy(x=arr_points2[:, 0], y=arr_points2[:, 1]),
         crs=4326,
     )
-    # 2/ LAS file
+    # 3/ LAS file
     fn_las = "/home/atom/ongoing/own/geoutils/test.laz"
 
-    # 3/ Non-point vector (for error raising)
+    # 4/ Non-point vector (for error raising)
     poly = Polygon([(5, 5), (6, 5), (6, 6), (5, 6)])
     gdf3 = gpd.GeoDataFrame({"geometry": [poly]}, crs="EPSG:4326")
 
@@ -60,7 +60,7 @@ def test_init__errors(self) -> None:
 
         # If vector is not only comprised of points
         with pytest.raises(ValueError, match="This vector file contains non-point geometries*"):
-            PointCloud(self.gdf3)
+            PointCloud(self.gdf3, data_column="z")
 
     def test_load(self) -> None:
         """