diff --git a/geoutils/interface/raster_point.py b/geoutils/interface/raster_point.py
index b9127f47..516c93b8 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 | gu.Vector:
+) -> NDArrayNum | gpd.GeoDataFrame:
     """
     Convert a raster to a point cloud. See Raster.to_pointcloud() for details.
     """
@@ -227,15 +227,13 @@ def _raster_to_pointcloud(
     )
 
     if not as_array:
-        points = gu.Vector(
-            gpd.GeoDataFrame(
+        pc = gpd.GeoDataFrame(
                 pixel_data.T,
                 columns=all_column_names,
                 geometry=gpd.points_from_xy(x_coords_2, y_coords_2),
                 crs=source_raster.crs,
             )
-        )
-        return points
+        return pc
     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 24407b99..1e320f4c 100644
--- a/geoutils/pointcloud/pointcloud.py
+++ b/geoutils/pointcloud/pointcloud.py
@@ -3,9 +3,10 @@
 
 import os.path
 import warnings
-from typing import Iterable, Literal
+from typing import Iterable, Literal, Any
 import pathlib
 
+import pandas as pd
 from pyproj import CRS
 import numpy as np
 import geopandas as gpd
@@ -14,7 +15,9 @@
 
 import geoutils as gu
 from geoutils.interface.gridding import _grid_pointcloud
-from geoutils._typing import Number
+from geoutils.raster.sampling import subsample_array
+from geoutils.interface.raster_point import _raster_to_pointcloud
+from geoutils._typing import Number, NDArrayNum, ArrayLike
 
 try:
     import laspy
@@ -23,29 +26,25 @@
     has_laspy = False
 
 
-def _load_laspy_data(filename: str, main_data_column: str, auxiliary_data_columns: list[str] | None = None) -> gpd.GeoDataFrame:
+def _load_laspy_data(filename: str, columns: list[str]) -> gpd.GeoDataFrame:
     """Load point cloud data from LAS/LAZ/COPC file."""
 
+    # Read file
     las = laspy.read(filename)
 
     # Get data from requested columns
-    if auxiliary_data_columns is not None:
-        data = [las[n] for n in auxiliary_data_columns]
-    else:
-        data = []
-    data.insert(0, las[main_data_column])
-    data = np.vstack(data).transpose()
+    data = np.vstack([las[n] for n in columns]).T
 
     # Build geodataframe
     gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=las.x, y=las.y,
                                                        crs=las.header.parse_crs(prefer_wkt=False)),
-                           data=data)
+                           data=data,
+                           columns=columns)
 
     return gdf
 
 
-
-def _load_laspy_metadata(filename: str, ) -> tuple[CRS, int, BoundingBox, list[str]]:
+def _load_laspy_metadata(filename: str, ) -> tuple[CRS, int, BoundingBox, pd.Index]:
     """Load point cloud metadata from LAS/LAZ/COPC file."""
 
     with laspy.open(filename) as f:
@@ -53,7 +52,7 @@ def _load_laspy_metadata(filename: str, ) -> tuple[CRS, int, BoundingBox, list[s
         crs = f.header.parse_crs(prefer_wkt=False)
         nb_points = f.header.point_count
         bounds = BoundingBox(left=f.header.x_min, right=f.header.x_max, bottom=f.header.y_min, top=f.header.y_max)
-        columns_names = list(f.header.point_format.dimension_names)
+        columns_names = pd.Index(list(f.header.point_format.dimension_names))
 
     return crs, nb_points, bounds, columns_names
 
@@ -74,7 +73,7 @@ class PointCloud(gu.Vector):
     """
     The georeferenced point cloud.
 
-    A point cloud is a vector of 2D point geometries associated to values from a main data column, optionally with
+    A point cloud is a vector of 2D point geometries associated to numeric values from a data column, and potentially
     auxiliary data columns.
 
      Main attributes:
@@ -94,13 +93,13 @@ class PointCloud(gu.Vector):
 
     def __init__(self,
                  filename_or_dataset: str | pathlib.Path | gpd.GeoDataFrame | gpd.GeoSeries | BaseGeometry,
-                 data_column: str | None = None):
+                 data_column: str | None = "z"):
         """
         Instantiate a point cloud from either a data column name and a vector (filename, GeoPandas dataframe or series,
         or a Shapely geometry), or only with a point cloud file type.
 
         :param filename_or_dataset: Path to vector file, or GeoPandas dataframe or series, or Shapely geometry.
-        :param data_column: Name of data column defining the point cloud.
+        :param data_column: Name of main data column defining the point cloud.
         """
 
         self._ds: gpd.GeoDataFrame | None = None
@@ -109,7 +108,7 @@ def __init__(self,
         self._bounds: BoundingBox | None = None
         self._data_column: str | None = None
         self._nb_points: int | None = None
-        self._columns: list[str] | None = None
+        self._columns: pd.Index | None = None
 
         # If PointCloud is passed, simply point back to PointCloud
         if isinstance(filename_or_dataset, PointCloud):
@@ -136,21 +135,29 @@ def __init__(self,
                     raise ValueError("This vector file contains non-point geometries, "
                                      "cannot be instantiated as a point cloud.")
 
-        if data_column not in self.columns:
-            raise ValueError(f"Data column {data_column} not found among columns, available columns "
-                             f"are: {', '.join(self.columns)}.")
-        self._data_column = data_column
+        # Set data column following user input
+        self.set_data_column(new_data_column=data_column)
+
+
+    # TODO: Could also move to Vector directly?
+    ##############################################
+    # OVERRIDDEN VECTOR METHODS TO SUPPORT LOADING
+    ##############################################
+
 
     @property
     def ds(self) -> gpd.GeoDataFrame:
+        """Geodataframe of the point cloud."""
+        # We need to override the Vector method to introduce the is_loaded dynamic for LAS files
         if not self.is_loaded:
             self.load()
         return self._ds  # type: ignore
 
     @ds.setter
     def ds(self, new_ds: gpd.GeoDataFrame | gpd.GeoSeries) -> None:
-        """Set a new geodataframe."""
-
+        """Set a new geodataframe for the point cloud."""
+        # We need to override the setter Vector method because we have overridden the property method
+        # (even if the code below is the same)
         if isinstance(new_ds, gpd.GeoDataFrame):
             self._ds = new_ds
         elif isinstance(new_ds, gpd.GeoSeries):
@@ -161,22 +168,51 @@ def ds(self, new_ds: gpd.GeoDataFrame | gpd.GeoSeries) -> None:
     @property
     def crs(self) -> CRS:
         """Coordinate reference system of the vector."""
-        if self._ds is not None:
-            return self.ds.crs
+        # Overridding method in Vector
+        if self.is_loaded:
+            return super().crs
         else:
             return self._crs
 
+    @property
     def bounds(self) -> BoundingBox:
-        if self._ds is not None:
+        # Overridding method in Vector
+        if self.is_loaded:
             return super().bounds
         else:
             return self._bounds
 
-    def load(self):
-        """Load point cloud from disk (only supported for LAS files)."""
+    @property
+    def all_columns(self) -> pd.Index:
+        """Index of all columns of the point cloud, excluding the column of 2D point geometries."""
+        # Overridding method in Vector
+        if self.is_loaded:
+            all_columns = super().columns
+            all_columns_nongeom = all_columns[all_columns != "geometry"]
+            return all_columns_nongeom
+        else:
+            return self._columns
 
-        ds = _load_laspy_data(filename=self.name, main_data_column=self.data_column)
-        self._ds = ds
+    #####################################
+    # NEW METHODS SPECIFIC TO POINT CLOUD
+    #####################################
+
+    @property
+    def data_column(self) -> str:
+        """Name of data column of the point cloud."""
+        return self._data_column
+
+    @data_column.setter
+    def data_column(self, new_data_column: str) -> None:
+        self.set_data_column(new_data_column=new_data_column)
+
+    def set_data_column(self, new_data_column: str) -> None:
+        """Set new column as point cloud data column."""
+
+        if new_data_column not in self.all_columns:
+            raise ValueError(f"Data column {new_data_column} not found among columns, available columns "
+                             f"are: {', '.join(self.all_columns)}.")
+        self._data_column = new_data_column
 
     @property
     def is_loaded(self) -> bool:
@@ -184,20 +220,39 @@ def is_loaded(self) -> bool:
         return self._ds is not None
 
     @property
-    def data_column(self) -> str:
-        """Name of main data column of the point cloud."""
-        return self._data_column
-
-    @property
-    def columns(self) -> list[str]:
-        """Name of auxiliary data columns of the point cloud."""
+    def nb_points(self) -> int:
+        """Number of points in the point cloud."""
+        # New method for point cloud
         if self.is_loaded:
-            return self.ds.columns
+            return len(self.ds)
         else:
-            return self._columns
+            return self._nb_points
+
+    def load(self, columns: Literal["all", "main"] | list[str] = "main"):
+        """
+        Load point cloud from disk (only supported for LAS files).
+
+        :param columns: Columns to load. Defaults to main data column only.
+        """
+
+        if self.is_loaded:
+            raise ValueError("Data are already loaded.")
+
+        if self.name is None:
+            raise AttributeError(
+                "Cannot load as filename is not set anymore. Did you manually update the filename attribute?"
+            )
+
+        if columns == "all":
+            columns = self.all_columns
+        elif columns == "main":
+            columns = [self.data_column]
+
+        ds = _load_laspy_data(filename=self.name, columns=columns)
+        self._ds = ds
 
     @classmethod
-    def from_array(cls, array: np.ndarray, crs: CRS, data_column: str | None = "z") -> PointCloud:
+    def from_array(cls, array: NDArrayNum, crs: CRS, data_column: str | None = "z") -> PointCloud:
         """Create point cloud from a 3 x N or N x 3 array of X coordinate, Y coordinates and Z values."""
 
         # Check shape
@@ -212,28 +267,55 @@ def from_array(cls, array: np.ndarray, crs: CRS, data_column: str | None = "z")
 
         # Build geodataframe
         gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=array_in[0, :], y=array_in[1, :], crs=crs),
-                               data={data_column: array[2, :]})
+                               data={data_column: array_in[2, :]})
 
         return cls(filename_or_dataset=gdf, data_column=data_column)
 
-
     @classmethod
-    def from_tuples(cls, tuples: Iterable[tuple[Number, Number, Number]], crs: CRS, data_column: str | None = None):
-        """Create point cloud from a N-size list of tuples (X coordinate, Y coordinate, Z value)."""
+    def from_tuples(cls, tuples_xyz: Iterable[tuple[Number, Number, Number]], crs: CRS, data_column: str | None = "z") -> PointCloud:
+        """Create point cloud from an iterable of 3-tuples (X coordinate, Y coordinate, Z value)."""
 
-        cls.from_array(np.array(zip(*tuples)), crs=crs, data_column=data_column)
+        return cls.from_array(np.array(tuples_xyz), crs=crs, data_column=data_column)
 
+    @classmethod
+    def from_xyz(cls, x: ArrayLike, y: ArrayLike, z: ArrayLike, crs: CRS, data_column: str | None = "z") -> PointCloud:
+        """Create point cloud from three 1D array-like coordinates for X/Y/Z."""
+
+        return cls.from_array(np.stack((x, y, z)), crs=crs, data_column=data_column)
 
-    def to_array(self):
+    def to_array(self) -> NDArrayNum:
         """Convert point cloud to a 3 x N array of X coordinates, Y coordinates and Z values."""
 
         return np.stack((self.geometry.x.values, self.geometry.y.values, self.ds[self.data_column].values), axis=0)
 
-    def to_tuples(self):
-        """Convert point cloud to a N-size list of tuples (X coordinate, Y coordinate, Z value)."""
+    def to_tuples(self) -> Iterable[tuple[Number, Number, Number]]:
+        """Convert point cloud to a list of 3-tuples (X coordinate, Y coordinate, Z value)."""
+
+        return list(zip(self.geometry.x.values, self.geometry.y.values, self.ds[self.data_column].values))
+
+    def to_xyz(self) -> tuple[NDArrayNum, NDArrayNum, NDArrayNum]:
+        """Convert point cloud to three 1D arrays of coordinates for X/Y/Z."""
 
         return self.geometry.x.values, self.geometry.y.values, self.ds[self.data_column].values
 
+    def pointcloud_equal(self, other: PointCloud, **kwargs: Any):
+        """
+        Check if two point clouds are equal.
+
+        This means that:
+        - The two vectors (geodataframes) are equal.
+        - The data column is the same for both point clouds.
+
+        Keyword arguments are passed to geopandas.assert_geodataframe_equal.
+        """
+
+        # Vector equality
+        vector_eq = super().vector_equal(other, **kwargs)
+        # Data column equality
+        data_column_eq = self.data_column == other.data_column
+
+        return all([vector_eq, data_column_eq])
+
     def grid(self,
              ref: gu.Raster | None,
              grid_coords: tuple[np.ndarray, np.ndarray] | None,
@@ -254,8 +336,13 @@ def grid(self,
 
         return gu.Raster.from_array(data=array, transform=transform, crs=self.crs, nodata=None)
 
-    # def subsample(self) -> PointCloud | NDArrayf:
-    #
+    def subsample(self, subsample: float | int, random_state: int | np.random.Generator | None = None) -> PointCloud:
+
+        subsample = subsample_array(array=self.ds[self.data_column].values, subsample=subsample,
+                                    return_indices=True, random_state=random_state)
+
+        return PointCloud(self.ds[subsample])
+
     # @classmethod
     # def from_raster(cls, raster: gu.Raster) -> PointCloud:
     #     """Create a point cloud from a raster. Equivalent with Raster.to_pointcloud."""
diff --git a/geoutils/raster/raster.py b/geoutils/raster/raster.py
index b3f17ad6..e9a9b8bd 100644
--- a/geoutils/raster/raster.py
+++ b/geoutils/raster/raster.py
@@ -28,6 +28,7 @@
 from rasterio.enums import Resampling
 from rasterio.plot import show as rshow
 
+import geoutils as gu
 from geoutils._config import config
 from geoutils._typing import (
     ArrayLike,
@@ -69,7 +70,6 @@
     decode_sensor_metadata,
     parse_and_convert_metadata_from_filename,
 )
-from geoutils.vector.vector import Vector
 
 # If python38 or above, Literal is builtin. Otherwise, use typing_extensions
 try:
@@ -590,7 +590,7 @@ def bounds(self) -> rio.coords.BoundingBox:
         return _bounds(transform=self.transform, shape=self.shape)
 
     @property
-    def footprint(self) -> Vector:
+    def footprint(self) -> gu.Vector:
         """Footprint of the raster."""
         return self.get_footprint_projected(self.crs)
 
@@ -2212,7 +2212,7 @@ def __array_function__(
     @overload
     def crop(
         self: RasterType,
-        crop_geom: RasterType | Vector | list[float] | tuple[float, ...],
+        crop_geom: RasterType | gu.Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: Literal[False] = False,
@@ -2221,7 +2221,7 @@ def crop(
     @overload
     def crop(
         self: RasterType,
-        crop_geom: RasterType | Vector | list[float] | tuple[float, ...],
+        crop_geom: RasterType | gu.Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: Literal[True],
@@ -2230,7 +2230,7 @@ def crop(
     @overload
     def crop(
         self: RasterType,
-        crop_geom: RasterType | Vector | list[float] | tuple[float, ...],
+        crop_geom: RasterType | gu.Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: bool = False,
@@ -2238,7 +2238,7 @@ def crop(
 
     def crop(
         self: RasterType,
-        crop_geom: RasterType | Vector | list[float] | tuple[float, ...],
+        crop_geom: RasterType | gu.Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: bool = False,
@@ -2670,7 +2670,7 @@ def get_bounds_projected(self, out_crs: CRS, densify_points: int = 5000) -> rio.
 
         return new_bounds
 
-    def get_footprint_projected(self, out_crs: CRS, densify_points: int = 5000) -> Vector:
+    def get_footprint_projected(self, out_crs: CRS, densify_points: int = 5000) -> gu.Vector:
         """
         Get raster footprint projected in a specified CRS.
 
@@ -2682,7 +2682,7 @@ def get_footprint_projected(self, out_crs: CRS, densify_points: int = 5000) -> V
          Reduce if time computation is really critical (ms) or increase if extent is not accurate enough.
         """
 
-        return Vector(
+        return gu.Vector(
             _get_footprint_projected(
                 bounds=self.bounds, in_crs=self.crs, out_crs=out_crs, densify_points=densify_points
             )
@@ -3315,7 +3315,7 @@ def to_pointcloud(
         as_array: Literal[True],
         random_state: int | np.random.Generator | None = None,
         force_pixel_offset: Literal["center", "ul", "ur", "ll", "lr"] = "ul",
-    ) -> Vector: ...
+    ) -> gu.Vector: ...
 
     @overload
     def to_pointcloud(
@@ -3330,7 +3330,7 @@ def to_pointcloud(
         as_array: bool = False,
         random_state: int | np.random.Generator | None = None,
         force_pixel_offset: Literal["center", "ul", "ur", "ll", "lr"] = "ul",
-    ) -> NDArrayNum | Vector: ...
+    ) -> NDArrayNum | gu.Vector: ...
 
     def to_pointcloud(
         self,
@@ -3343,7 +3343,7 @@ def to_pointcloud(
         as_array: bool = False,
         random_state: int | np.random.Generator | None = None,
         force_pixel_offset: Literal["center", "ul", "ur", "ll", "lr"] = "ul",
-    ) -> NDArrayNum | Vector:
+    ) -> NDArrayNum | gu.PointCloud:
         """
         Convert raster to point cloud.
 
@@ -3390,7 +3390,7 @@ def to_pointcloud(
         :returns: A point cloud, or array of the shape (N, 2 + count) where N is the sample count.
         """
 
-        return _raster_to_pointcloud(
+        return gu.PointCloud(_raster_to_pointcloud(
             source_raster=self,
             data_column_name=data_column_name,
             data_band=data_band,
@@ -3401,7 +3401,7 @@ def to_pointcloud(
             as_array=as_array,
             random_state=random_state,
             force_pixel_offset=force_pixel_offset,
-        )
+        ))
 
     @classmethod
     def from_pointcloud_regular(
@@ -3447,7 +3447,7 @@ def polygonize(
         self,
         target_values: Number | tuple[Number, Number] | list[Number] | NDArrayNum | Literal["all"] = "all",
         data_column_name: str = "id",
-    ) -> Vector:
+    ) -> gu.Vector:
         """
         Polygonize the raster into a vector.
 
@@ -3456,14 +3456,14 @@ def polygonize(
         :param data_column_name: Data column name to be associated with target values in the output vector
             (defaults to "id").
 
-        :returns: Vector containing the polygonized geometries associated to target values.
+        :returns: gu.Vector containing the polygonized geometries associated to target values.
         """
 
         return _polygonize(source_raster=self, target_values=target_values, data_column_name=data_column_name)
 
     def proximity(
         self,
-        vector: Vector | None = None,
+        vector: gu.Vector | None = None,
         target_values: list[float] | None = None,
         geometry_type: str = "boundary",
         in_or_out: Literal["in"] | Literal["out"] | Literal["both"] = "both",
@@ -3479,7 +3479,7 @@ def proximity(
         passing "geometry", or any lower dimensional geometry attribute such as "centroid", "envelope" or "convex_hull".
         See all geometry attributes in the Shapely documentation at https://shapely.readthedocs.io/.
 
-        :param vector: Vector for which to compute the proximity to geometry,
+        :param vector: gu.Vector for which to compute the proximity to geometry,
             if not provided computed on this raster target pixels.
         :param target_values: (Only with raster) List of target values to use for the proximity,
             defaults to all non-zero values.
@@ -3777,7 +3777,7 @@ def reproject(
     @overload
     def crop(
         self: Mask,
-        crop_geom: Mask | Vector | list[float] | tuple[float, ...],
+        crop_geom: Mask | gu.Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: Literal[False] = False,
@@ -3786,7 +3786,7 @@ def crop(
     @overload
     def crop(
         self: Mask,
-        crop_geom: Mask | Vector | list[float] | tuple[float, ...],
+        crop_geom: Mask | gu.Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: Literal[True],
@@ -3795,7 +3795,7 @@ def crop(
     @overload
     def crop(
         self: Mask,
-        crop_geom: Mask | Vector | list[float] | tuple[float, ...],
+        crop_geom: Mask | gu.Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: bool = False,
@@ -3803,7 +3803,7 @@ def crop(
 
     def crop(
         self: Mask,
-        crop_geom: Mask | Vector | list[float] | tuple[float, ...],
+        crop_geom: Mask | gu.Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: bool = False,
@@ -3828,7 +3828,7 @@ def polygonize(
         self,
         target_values: Number | tuple[Number, Number] | list[Number] | NDArrayNum | Literal["all"] = 1,
         data_column_name: str = "id",
-    ) -> Vector:
+    ) -> gu.Vector:
         # If target values is passed but does not correspond to 0 or 1, raise a warning
         if not isinstance(target_values, (int, np.integer, float, np.floating)) or target_values not in [0, 1]:
             warnings.warn("In-value converted to 1 for polygonizing boolean mask.")
@@ -3847,7 +3847,7 @@ def polygonize(
 
     def proximity(
         self,
-        vector: Vector | None = None,
+        vector: gu.Vector | None = None,
         target_values: list[float] | None = None,
         geometry_type: str = "boundary",
         in_or_out: Literal["in"] | Literal["out"] | Literal["both"] = "both",
diff --git a/tests/test_pointcloud/test_pointcloud.py b/tests/test_pointcloud/test_pointcloud.py
index ab609795..21365c82 100644
--- a/tests/test_pointcloud/test_pointcloud.py
+++ b/tests/test_pointcloud/test_pointcloud.py
@@ -2,8 +2,11 @@
 from __future__ import annotations
 
 import pytest
+from pyproj import CRS
+from rasterio.coords import BoundingBox
 import numpy as np
 import geopandas as gpd
+from shapely import Polygon
 from geopandas.testing import assert_geodataframe_equal
 
 from geoutils import PointCloud
@@ -13,41 +16,216 @@ class TestPointCloud:
 
     # 1/ Synthetic point cloud with no auxiliary column
     rng = np.random.default_rng(42)
-    points = rng.integers(low=1, high=1000, size=(100, 2)) + rng.normal(0, 0.15, size=(100, 2))
-    val1 = rng.normal(scale=3, size=100)
-    pc1 = gpd.GeoDataFrame(data={"b1": val1}, geometry=gpd.points_from_xy(x=points[:, 0], y=points[:, 1]), crs=4326)
-
+    arr_points = rng.integers(low=1, high=1000, size=(100, 3)) + rng.normal(0, 0.15, size=(100, 3))
+    gdf1 = gpd.GeoDataFrame(data={"b1": arr_points[:, 2]},
+                            geometry=gpd.points_from_xy(x=arr_points[:, 0], y=arr_points[:, 1]), crs=4326)
+
+    # 2/ Synthetic point cloud with auxiliary column
+    arr_points2 = rng.integers(low=1, high=1000, size=(100, 4)) + rng.normal(0, 0.15, size=(100, 4))
+    gdf2 = gpd.GeoDataFrame(data=arr_points[:, 2:], columns=["b1", "b2"],
+                            geometry=gpd.points_from_xy(x=arr_points[:, 0], y=arr_points[:, 1]), crs=4326)
     # 2/ LAS file
-    filename = "/home/atom/ongoing/own/geoutils/points.laz"
+    fn_las = "/home/atom/ongoing/own/geoutils/points.laz"
+
+    # 3/ Non-point vector (for error raising)
+    poly = Polygon([(5, 5), (6, 5), (6, 6), (5, 6)])
+    gdf3 = gpd.GeoDataFrame({"geometry": [poly]}, crs="EPSG:4326")
 
-    def test_init(self):
+    def test_init(self) -> None:
         """Test instantiation of a point cloud."""
 
         # 1/ For a single column point cloud
-        pc = PointCloud(self.pc1, data_column="b1")
+        pc = PointCloud(self.gdf1, data_column="b1")
 
         # Assert that both the dataframe and data column name are equal
         assert pc.data_column == "b1"
-        assert_geodataframe_equal(pc.ds, self.pc1)
+        assert_geodataframe_equal(pc.ds, self.gdf1)
 
-        # 2/ For a point cloud from filie
-        pc = PointCloud(self.filename, data_column="Z")
+        # 2/ For a point cloud from LAS/LAZ file
+        pc = PointCloud(self.fn_las, data_column="Z")
 
         assert pc.data_column == "Z"
         assert not pc.is_loaded
 
+    def test_init__errors(self) -> None:
+        """Test errors raised by point cloud instantiation."""
 
-    def test_data_column_name(self):
-        pass
+        # If the data column does not exist
+        with pytest.raises(ValueError, match="Data column column_that_does_not_exist not found*"):
+            PointCloud(self.gdf1, data_column="column_that_does_not_exist")
 
-    def test_from_array(self):
-        """Test building point cloud from array."""
+        # If vector is not only comprised of points
+        with pytest.raises(ValueError, match="This vector file contains non-point geometries*"):
+            PointCloud(self.gdf3)
+
+    def test_load(self) -> None:
+        """
+        Test loading of a point cloud (only possible with a LAS file).
+
+        This test also serves to test the overridden methods "crs", "bounds", "nb_points", "all_columns" in relation
+        to loading.
+        """
+
+        # 1/ Check unloaded and loaded attributes are all the same
+
+        pc = PointCloud(self.fn_las, data_column="Z")
+
+        # Check point cloud is not loaded, and fetch before metadata
+        assert not pc.is_loaded
+        before_crs = pc.crs
+        before_bounds = pc.bounds
+        before_nb_points = pc.nb_points
+        before_columns = pc.all_columns
+
+        # Load and fetch after metadata
+        pc.load(columns="all")
+        assert pc.is_loaded
+
+        after_crs = pc.crs
+        after_bounds = pc.bounds
+        after_nb_points = pc.nb_points
+        after_columns = pc.all_columns
+
+        # Check those are equal
+        assert before_crs == after_crs
+        assert before_bounds == after_bounds
+        assert before_nb_points == after_nb_points
+        assert all(before_columns == after_columns)
+
+        # 2/ Check default column argument
+        pc = PointCloud(self.fn_las, data_column="Z")
+        pc.load()
+
+        assert pc.all_columns == ["Z"]
+
+        # 3/ Check implicit loading when calling a function requiring .ds
+        pc = PointCloud(self.fn_las, data_column="Z")
+        assert not pc.is_loaded
+
+        pc2 = pc.buffer(distance=0.1)
+        assert pc.is_loaded
+
+    def test_load__errors(self) -> None:
+        """Test errors raised by loading."""
 
-    def test_from_tuples(self):
-        pass
+        pc = PointCloud(self.fn_las, data_column="Z")
+        pc.load()
 
-    def test_to_array(self):
-        pass
+        # Error if already loaded
+        with pytest.raises(ValueError, match="Data are already loaded."):
+            pc.load()
+
+        pc = PointCloud(self.fn_las, data_column="Z")
+        pc._name = None
+        with pytest.raises(AttributeError, match="Cannot load as filename is not set anymore.*"):
+            pc.load()
+
+    def test_data_column(self) -> None:
+        """Test the setting and getting of the main data column."""
+
+        # Assert column is set properly at instantiation
+        pc = PointCloud(self.gdf1, data_column="b1")
+        assert pc.data_column == "b1"
+
+        # And can be reset to another name if it exists
+        pc2 = PointCloud(self.gdf2, data_column="b1")
+        assert pc2.data_column == "b1"
+        # First syntax
+        pc2.data_column = "b2"
+        assert pc2.data_column == "b2"
+        # Equivalent syntax
+        pc2.set_data_column("b1")
+        assert pc2.data_column == "b1"
+
+    def test_data_column__errors(self) -> None:
+        """Test errors raised during setting of data column."""
+
+        pc = PointCloud(self.gdf1, data_column="b1")
+        # If the data column does not exist
+        with pytest.raises(ValueError, match="Data column column_that_does_not_exist not found*"):
+            pc.data_column = "column_that_does_not_exist"
+        # Equivalent syntax
+        with pytest.raises(ValueError, match="Data column column_that_does_not_exist not found*"):
+            pc.set_data_column("column_that_does_not_exist")
+
+    def test_pointcloud_equal(self) -> None:
+        """Test pointcloud equality."""
+
+
+
+    def test_from_array(self) -> None:
+        """Test building point cloud from array."""
 
-    def test_to_tuples(self):
-        pass
\ No newline at end of file
+        # Build from array and compare
+        pc1 = PointCloud(self.gdf1, data_column="b1")
+        pc_from_arr = PointCloud.from_array(array=self.arr_points, crs=4326, data_column="b1")
+        assert pc_from_arr.pointcloud_equal(pc1)
+
+        # Should be the same witht transposed array
+        pc_from_arr = PointCloud.from_array(array=self.arr_points.T, crs=4326, data_column="b1")
+        assert pc_from_arr.pointcloud_equal(pc1)
+
+    def test_from_array__errors(self):
+        """Test errors raised during creation with array."""
+
+        array = np.ones((4, 5))
+        with pytest.raises(ValueError, match="Array must be of shape 3xN or Nx3."):
+            PointCloud.from_array(array=array, crs=4326)
+
+    def test_from_tuples(self) -> None:
+        """Test building point cloud from list of tuples."""
+
+        pc1 = PointCloud(self.gdf1, data_column="b1")
+        tuples_xyz = list(zip(self.arr_points[:, 0], self.arr_points[:, 1], self.arr_points[:, 2]))
+        pc_from_tuples = PointCloud.from_tuples(tuples_xyz=tuples_xyz, crs=4326, data_column="b1")
+        assert pc_from_tuples.pointcloud_equal(pc1)
+
+    def test_from_xyz(self) -> None:
+        """Test building point cloud from xyz array-like."""
+
+        # Build from array and compare
+        pc1 = PointCloud(self.gdf1, data_column="b1")
+        pc_from_xyz = PointCloud.from_xyz(x=self.arr_points[:, 0], y=self.arr_points[:, 1], z=self.arr_points[:, 2],
+                                              crs=4326, data_column="b1")
+        assert pc_from_xyz.pointcloud_equal(pc1)
+
+        # Test with lists
+        pc_from_xyz = PointCloud.from_xyz(x=list(self.arr_points[:, 0]),
+                                          y=list(self.arr_points[:, 1]),
+                                          z=list(self.arr_points[:, 2]),
+                                        crs=4326,
+                                        data_column="b1")
+        assert pc_from_xyz.pointcloud_equal(pc1)
+
+        # Test with tuples
+        pc_from_xyz = PointCloud.from_xyz(x=tuple(self.arr_points[:, 0]),
+                                            y=tuple(self.arr_points[:, 1]),
+                                            z=tuple(self.arr_points[:, 2]),
+                                            crs=4326,
+                                            data_column="b1")
+        assert pc_from_xyz.pointcloud_equal(pc1)
+
+    def test_to_array(self) -> None:
+        """Test exporting point cloud to array."""
+
+        # Convert point cloud and compare
+        pc1 = PointCloud(self.gdf1, data_column="b1")
+        arr_from_pc = pc1.to_array()
+        assert np.array_equal(arr_from_pc, self.arr_points.T)
+
+    def test_to_tuples(self) -> None:
+        """Test exporting point cloud to tuples."""
+
+        # Convert point cloud and compare
+        pc1 = PointCloud(self.gdf1, data_column="b1")
+        tuples_xyz = list(zip(self.arr_points[:, 0], self.arr_points[:, 1], self.arr_points[:, 2]))
+        tuples_from_pc = pc1.to_tuples()
+        assert tuples_from_pc == tuples_xyz
+
+    def test_to_xyz(self) -> None:
+        """Test exporting point cloud to xyz arrays."""
+
+        # Convert point cloud and compare
+        pc1 = PointCloud(self.gdf1, data_column="b1")
+        xyz_from_pc = pc1.to_xyz()
+        assert np.array_equal(np.stack(xyz_from_pc), self.arr_points.T)
\ No newline at end of file