diff --git a/doc/source/api.md b/doc/source/api.md
index d0966127..d1a0d999 100644
--- a/doc/source/api.md
+++ b/doc/source/api.md
@@ -47,6 +47,7 @@ documentation.
     Raster.crs
     Raster.transform
     Raster.nodata
+    Raster.area_or_point
 ```
 
 ### Derived attributes
@@ -83,8 +84,8 @@ documentation.
     Raster.reproject
     Raster.polygonize
     Raster.proximity
-    Raster.value_at_coords
     Raster.interp_points
+    Raster.reduce_points
 ```
 
 ### Plotting
@@ -120,6 +121,7 @@ documentation.
     Raster.load
     Raster.save
     Raster.to_pointcloud
+    Raster.from_pointcloud_regular
     Raster.to_rio_dataset
     Raster.to_xarray
 ```
@@ -133,7 +135,7 @@ documentation.
     Raster.xy2ij
     Raster.ij2xy
     Raster.coords
-    Raster.shift
+    Raster.translate
     Raster.outside_image
 ```
 
@@ -265,6 +267,7 @@ And reverse operations.
 .. autosummary::
     :toctree: gen_modules/
 
+    raster.load_multiple_rasters
     raster.stack_rasters
     raster.merge_rasters
 ```
diff --git a/doc/source/conf.py b/doc/source/conf.py
index c7d2d3e7..36ea5566 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -65,6 +65,7 @@
     "xdem": ("https://xdem.readthedocs.io/en/stable", None),
     "rioxarray": ("https://corteva.github.io/rioxarray/stable/", None),
     "pandas": ("https://pandas.pydata.org/docs/", None),
+    "scipy": ("https://docs.scipy.org/doc/scipy/", None),
 }
 
 example_path = os.path.join("../", "../", "examples")
@@ -110,7 +111,7 @@
 }
 
 # For matplotlib figures generate with sphinx plot: (suffix, dpi)
-plot_formats = [(".png", 400)]
+plot_formats = [(".png", 500)]
 
 # To avoid long path names in inheritance diagrams
 inheritance_alias = {
@@ -193,7 +194,8 @@ def setup(app):
         "⚠️ Our 0.1 release refactored several early-development functions for long-term stability, "
         'to update your code see <a href="https://github.com/GlacioHack/geoutils/releases/tag/v0.1.0">here</a>. ⚠️'
         "<br>Future changes will come with deprecation warnings! 🙂"
-    )
+    ),
+    "show_toc_level": 3
     # "logo_only": True,
     # "icon_links": [
     #         {
diff --git a/doc/source/config.md b/doc/source/config.md
new file mode 100644
index 00000000..1a702a9b
--- /dev/null
+++ b/doc/source/config.md
@@ -0,0 +1,37 @@
+---
+file_format: mystnb
+jupytext:
+  formats: md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+kernelspec:
+  display_name: geoutils-env
+  language: python
+  name: geoutils
+---
+(config)=
+# Configuration
+
+You can configure the default behaviour of GeoUtils at the package level for operations that depend on user preference
+(such as resampling method, or pixel interpretation).
+
+## Changing configuration during a session
+
+Using a global configuration setting ensures operations will always be performed consistently, even when used
+under-the-hood by higher-level methods (such as [Coregistration](https://xdem.readthedocs.io/en/stable/coregistration.html)),
+without having to rely on multiple keyword arguments to pass to subfunctions.
+
+```{code-cell}
+import geoutils as gu
+# Changing default behaviour for pixel interpretation for this session
+gu.config["shift_area_or_point"] = False
+```
+
+## Default configuration file
+
+Below is the full default configuration file, which is updated by changes in configuration during a session.
+
+```{literalinclude} ../../geoutils/config.ini
+:class: full-width
+```
diff --git a/doc/source/rasters_index.md b/doc/source/data_object_index.md
similarity index 73%
rename from doc/source/rasters_index.md
rename to doc/source/data_object_index.md
index 1becdad7..958cf079 100644
--- a/doc/source/rasters_index.md
+++ b/doc/source/data_object_index.md
@@ -1,5 +1,5 @@
-(rasters-index)=
-# Rasters
+(data-object-index)=
+# Geospatial data objects
 
 Prefer to learn by running examples? Explore our example galleries on {ref}`examples-io`, {ref}`examples-handling` and {ref}`examples-analysis`.
 
@@ -7,6 +7,8 @@ Prefer to learn by running examples? Explore our example galleries on {ref}`exam
 :maxdepth: 2
 
 raster_class
+vector_class
+pointcloud_class
 mask_class
 satimg_class
 ```
diff --git a/doc/source/distance_ops.md b/doc/source/distance_ops.md
new file mode 100644
index 00000000..bd52eef5
--- /dev/null
+++ b/doc/source/distance_ops.md
@@ -0,0 +1,99 @@
+---
+file_format: mystnb
+jupytext:
+  formats: md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+kernelspec:
+  display_name: geoutils-env
+  language: python
+  name: geoutils
+---
+(distance-ops)=
+# Distance operations
+
+Computing distance between sets of geospatial data or manipulating their shape based on distance is often important
+for later analysis. To facilitate this type of operations, GeoUtils implements distance-specific functionalities
+for both vectors and rasters.
+
+```{code-cell} ipython3
+:tags: [remove-cell]
+
+# To get a good resolution for displayed figures
+from matplotlib import pyplot
+pyplot.rcParams['figure.dpi'] = 600
+pyplot.rcParams['savefig.dpi'] = 600
+pyplot.rcParams['font.size'] = 9
+```
+
+```{tip}
+It is often important to compute distances in a metric CRS. For this, reproject (with
+{func}`~geoutils.Raster.reproject`) to a local metric CRS (that can be estimated with {func}`~geoutils.Raster.get_metric_crs`).
+```
+
+## Proximity
+
+Proximity corresponds to **the distance to the closest target geospatial data**, computed on each pixel of a raster's grid.
+The target geospatial data can be either a vector or a raster.
+
+{func}`geoutils.Raster.proximity` and {func}`geoutils.Vector.proximity`
+
+```{code-cell} ipython3
+:tags: [hide-cell]
+:mystnb:
+:  code_prompt_show: "Show the code for opening example files"
+:  code_prompt_hide: "Hide the code for opening example files"
+
+import matplotlib.pyplot as plt
+import geoutils as gu
+import numpy as np
+
+rast = gu.Raster(gu.examples.get_path("everest_landsat_b4"))
+rast.set_nodata(0)  # Annoying to have to do this here, should we update it in the example?
+vect = gu.Vector(gu.examples.get_path("everest_rgi_outlines"))
+```
+
+```{code-cell} ipython3
+# Compute proximity to vector outlines
+proximity = vect.proximity(rast)
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show the code for plotting the figure"
+:  code_prompt_hide: "Hide the code for plotting the figure"
+
+f, ax = plt.subplots(1, 2)
+ax[0].set_title("Raster and vector")
+rast.plot(ax=ax[0], cmap="gray", add_cbar=False)
+vect.plot(ref_crs=rast, ax=ax[0], ec="k", fc="none")
+ax[1].set_title("Proximity")
+proximity.plot(ax=ax[1], cmap="viridis", cbar_title="Distance to outlines (m)")
+_ = ax[1].set_yticklabels([])
+plt.tight_layout()
+```
+
+## Buffering without overlap
+
+Buffering consists in **expanding or collapsing vector geometries equally in all directions**. However, this can often lead to overlap
+between shapes, which is sometimes undesirable. Using Voronoi polygons, we provide a buffering method with overlap.
+
+{func}`geoutils.Vector.buffer_without_overlap`
+
+```{code-cell} ipython3
+# Compute buffer without overlap from vector exterior
+vect_buff_nolap = vect.buffer_without_overlap(buffer_size=500)
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show the code for plotting the figure"
+:  code_prompt_hide: "Hide the code for plotting the figure"
+
+# Plot with color to see that the attributes are retained for every feature
+vect.plot(ax="new", ec="k", column="Area", alpha=0.5, add_cbar=False)
+vect_buff_nolap.plot(column="Area", cbar_title="Buffer around initial features\ncolored by glacier area (km)")
+```
diff --git a/doc/source/feature_overview.md b/doc/source/feature_overview.md
index 52b655de..b53c9b68 100644
--- a/doc/source/feature_overview.md
+++ b/doc/source/feature_overview.md
@@ -16,8 +16,6 @@ kernelspec:
 
 The following presents a descriptive example show-casing all core features of GeoUtils.
 
-For more details, refer to the {ref}`core-index`, {ref}`rasters-index` or {ref}`vectors-index` pages.
-
 ```{tip}
 All pages of this documentation containing code cells can be **run interactively online without the need of setting up your own environment**. Simply click the top launch button!
 (MyBinder can be a bit capricious: you might have to be patient, or restart it after the build is done the first time 😅)
@@ -25,6 +23,16 @@ All pages of this documentation containing code cells can be **run interactively
 Alternatively, start your own notebook to test GeoUtils at [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/GlacioHack/geoutils/main).
 ```
 
+```{code-cell} ipython3
+:tags: [remove-cell]
+
+# To get a good resolution for displayed figures
+from matplotlib import pyplot
+pyplot.rcParams['figure.dpi'] = 600
+pyplot.rcParams['savefig.dpi'] = 600
+pyplot.rcParams['font.size'] = 9
+```
+
 ## The core {class}`~geoutils.Raster` and {class}`~geoutils.Vector` classes
 
 In GeoUtils, geospatial handling is object-based and revolves around {class}`~geoutils.Raster` and {class}`~geoutils.Vector`.
diff --git a/doc/source/georeferencing.md b/doc/source/georeferencing.md
new file mode 100644
index 00000000..73dbf764
--- /dev/null
+++ b/doc/source/georeferencing.md
@@ -0,0 +1,248 @@
+---
+file_format: mystnb
+jupytext:
+  formats: md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+kernelspec:
+  display_name: geoutils-env
+  language: python
+  name: geoutils
+---
+(georeferencing)=
+# Referencing
+
+Below, a summary of the **georeferencing attributes** of geospatial data objects and the **methods to manipulate these
+georeferencing attributes** in different projections, without any data transformation. For georeferenced transformations
+(such as reprojection, cropping), see {ref}`geotransformations`.
+
+```{code-cell} ipython3
+:tags: [remove-cell]
+
+# To get a good resolution for displayed figures
+from matplotlib import pyplot
+pyplot.rcParams['figure.dpi'] = 600
+pyplot.rcParams['savefig.dpi'] = 600
+pyplot.rcParams['font.size'] = 9
+```
+
+## Attributes
+
+In GeoUtils, the **georeferencing syntax is consistent across all geospatial data objects**. Additionally, **data objects
+load only their metadata by default**, allowing quick operations on georeferencing without requiring the array data
+(for a {class}`~geoutils.Raster`) or geometry data (for a {class}`~geoutils.Vector`) to be present in memory.
+
+### Metadata summary
+
+To summarize all the metadata of a geospatial data object, including its georeferencing, {func}`~geoutils.Raster.info` can be used:
+
+
+```{code-cell} ipython3
+:tags: [hide-cell]
+:mystnb:
+:  code_prompt_show: "Show the code for opening example files"
+:  code_prompt_hide: "Hide the code for opening example files"
+
+import geoutils as gu
+rast = gu.Raster(gu.examples.get_path("exploradores_aster_dem"))
+vect = gu.Vector(gu.examples.get_path("exploradores_rgi_outlines"))
+```
+
+```{code-cell} ipython3
+# Print raster info
+rast.info()
+```
+
+```{code-cell} ipython3
+# Print vector info
+vect.info()
+```
+
+### Coordinate reference systems
+
+[Coordinate reference systems (CRSs)](https://en.wikipedia.org/wiki/Spatial_reference_system), sometimes also called
+spatial reference systems (SRSs), define the 2D projection of the geospatial data. They are stored as a
+{class}`pyproj.crs.CRS` object in {attr}`~geoutils.Raster.crs`.
+
+```{code-cell} ipython3
+# Show CRS attribute of raster
+print(rast.crs)
+```
+```{code-cell} ipython3
+# Show CRS attribute of vector as a WKT
+print(vect.crs.to_wkt())
+```
+
+More information on the manipulation of {class}`pyproj.crs.CRS` objects can be found in [PyProj's documentation](https://pyproj4.github.io/pyproj/stable/).
+
+```{note}
+3D CRSs for elevation data are only emerging, and not consistently defined in the metadata.
+The [vertical referencing functionalities of xDEM](https://xdem.readthedocs.io/en/stable/vertical_ref.html)
+can help define a 3D CRS.
+```
+
+(bounds)=
+### Bounds
+
+Bounds define the spatial extent of geospatial data, composed of the "left", "right", "bottom" and "top" coordinates.
+The {attr}`~geoutils.Raster.bounds` of a raster or a vector is a {class}`rasterio.coords.BoundingBox` object:
+
+```{code-cell} ipython3
+# Show bounds attribute of raster
+rast.bounds
+```
+```{code-cell} ipython3
+# Show bounds attribute of vector
+vect.bounds
+```
+
+```{note}
+To define {attr}`~geoutils.Raster.bounds` consistently between rasters and vectors, {attr}`~geoutils.Vector.bounds`
+ corresponds to {attr}`geopandas.GeoSeries.total_bounds` (total bounds of all geometry features) converted to a {class}`rasterio.coords.BoundingBox`.
+
+To reproduce the behaviour of {attr}`geopandas.GeoSeries.bounds` (per-feature bounds) with a
+{class}`~geoutils.Vector`, use {attr}`~geoutils.Vector.geom_bounds`.
+```
+
+### Footprints
+
+As reprojections between CRSs deform shapes, including extents, it is often better to consider a vectorized footprint
+to calculate intersections in different projections. The {class}`~geoutils.Raster.footprint` is a
+{class}`~geoutils.Vector` object with a single polygon geometry for which points have been densified, allowing
+reliable computation of extents between CRSs.
+
+```{code-cell} ipython3
+# Print raster footprint
+rast.get_footprint_projected(rast.crs)
+```
+```{code-cell} ipython3
+# Plot vector footprint
+vect.get_footprint_projected(vect.crs).plot()
+```
+
+### Grid (only for rasters)
+
+A raster's grid origin and resolution are defined by its geotransform attribute, {attr}`~geoutils.Raster.transform`.
+Comined with the 2D shape of the data array {attr}`~geoutils.Raster.shape` (and independently of the number of
+bands {attr}`~geoutils.Raster.bands`), these two attributes define the georeferenced grid of a raster.
+
+From it are derived the resolution {attr}`~geoutils.Raster.res`, and {attr}`~geoutils.Raster.height` and
+{attr}`~geoutils.Raster.width`, as well as the bounds detailed above in {ref}`bounds`.
+
+```{code-cell} ipython3
+# Get raster transform and shape
+print(rast.transform)
+print(rast.shape)
+```
+
+(pixel-interpretation)=
+### Pixel interpretation (only for rasters)
+
+A largely overlooked aspect of a raster's georeferencing is the pixel interpretation stored in the
+[AREA_OR_POINT metadata](https://gdal.org/user/raster_data_model.html#metadata).
+Pixels can be interpreted either as **"Area"** (the most common) where **the value represents a sampling over the region
+of the pixel (and typically refers to the upper-left corner coordinate)**, or as **"Point"**
+where **the value relates to a point sample (and typically refers to the center of the pixel)**, the latter often used
+for digital elevation models (DEMs).
+
+Pixel interpretation is stored as a string in the {attr}`geoutils.Raster.area_or_point` attribute.
+
+```{code-cell} ipython3
+# Get pixel interpretation of raster
+rast.area_or_point
+```
+
+Although this interpretation is not intended to influence georeferencing, it **can influence sub-pixel coordinate
+interpretation during analysis**, especially for raster–vector–point interfacing operations such as point interpolation,
+or re-gridding, and might also be a problem if defined differently when comparing two rasters.
+
+```{important}
+By default, **pixel interpretation induces a half-pixel shift during raster–point interfacing for a "Point" interpretation**
+(mirroring [GDAL's default ground-control point behaviour](https://trac.osgeo.org/gdal/wiki/rfc33_gtiff_pixelispoint)),
+but only **raises a warning for raster–raster operations** if interpretations differ.
+
+This behaviour can be modified at the package-level by using GeoUtils' {ref}`config`
+`shift_area_or_point` and `warns_area_or_point`.
+```
+
+## Manipulation
+
+Several functionalities are available to facilitate the manipulation of the georeferencing.
+
+### Getting projected bounds and footprints
+
+Retrieving projected bounds or footprints in any CRS is possible using directly {func}`~geoutils.Raster.get_bounds_projected`
+and {func}`~geoutils.Raster.get_footprint_projected`.
+
+```{code-cell} ipython3
+# Get footprint of larger buffered vector in polar stereo CRS (to show deformations)
+vect.buffer_metric(10**6).get_footprint_projected(3995).plot()
+```
+
+### Getting a metric CRS
+
+A local metric coordinate system can be estimated for both {class}`Rasters<geoutils.Raster>` and {class}`Vectors<geoutils.Vector>` through the
+{func}`~geoutils.Raster.get_metric_crs` function.
+
+The metric system returned can be either "universal" (zone of the Universal Transverse Mercator or Universal Polar Stereographic system), or "custom"
+(Mercator or Polar projection centered on the {class}`Raster<geoutils.Raster>` or {class}`Vector<geoutils.Vector>`).
+
+```{code-cell} ipython3
+# Get local metric CRS
+rast.get_metric_crs()
+```
+
+### Re-set georeferencing metadata
+
+The georeferencing metadata of an object can be re-set (overwritten) by setting the corresponding attribute such as {func}`geoutils.Vector.crs` or
+{func}`geoutils.Raster.transform`. When specific options might be useful during setting, a set function exists,
+such as for {func}`geoutils.Raster.set_area_or_point`.
+
+```{warning}
+Re-setting should only be used if the **data was erroneously defined and needs to be corrected in-place**.
+To create geospatial data from its attributes, use the construction functions such as {func}`~geoutils.Raster.from_array`.
+```
+
+```{code-cell} ipython3
+# Re-set CRS
+import pyproj
+rast.crs = pyproj.CRS(4326)
+rast.crs
+```
+
+```{code-cell} ipython3
+# Re-set pixel interpretation
+rast.set_area_or_point("Point")
+rast.area_or_point
+```
+
+
+### Coordinates to indexes (only for rasters)
+
+Raster grids are notoriously unintuitive to manipulate on their own due to the Y axis being inverted and stored as first axis.
+GeoUtils' features account for this under-the-hood when plotting, interpolating, gridding, or performing any other operation involving the raster coordinates.
+
+Three functions facilitate the manipulation of coordinates, while respecting any {ref}`Pixel interpretation<pixel-interpretation>`:
+
+1. {func}`~geoutils.Raster.xy2ij` to convert array indices to coordinates,
+2. {func}`~geoutils.Raster.ij2xy` to convert coordinates to array indices (reversible with {func}`~geoutils.Raster.xy2ij` for any pixel interpretation),
+3. {func}`~geoutils.Raster.coords` to directly obtain coordinates in order corresponding to the data array axes, possibly as a meshgrid.
+
+```{code-cell} ipython3
+# Get coordinates from row/columns indices
+x, y = rast.ij2xy(i=[0, 1], j=[2, 3])
+(x, y)
+```
+
+```{code-cell} ipython3
+# Get indices from coordinates
+i, j = rast.xy2ij(x=x, y=y)
+(i, j)
+```
+
+```{code-cell} ipython3
+:tags: [hide-output]
+# Get vector X/Y coordinates corresponding to data array
+rast.coords(grid=False)
+```
diff --git a/doc/source/geotransformations.md b/doc/source/geotransformations.md
new file mode 100644
index 00000000..ba08df84
--- /dev/null
+++ b/doc/source/geotransformations.md
@@ -0,0 +1,266 @@
+---
+file_format: mystnb
+jupytext:
+  formats: md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+kernelspec:
+  display_name: geoutils-env
+  language: python
+  name: geoutils
+---
+(geotransformations)=
+# Transformations
+
+In GeoUtils, **for all geospatial data objects, georeferenced transformations are exposed through the same functions**
+{func}`~geoutils.Raster.reproject`, {func}`~geoutils.Raster.crop` and {func}`~geoutils.Raster.shift`. Additionally,
+for convenience and consistency during analysis, most operations can be passed a {class}`~geoutils.Raster` or
+{class}`~geoutils.Vector` as a reference to match.
+In that case, no other argument is necessary. For more details, see {ref}`core-match-ref`.
+
+```{code-cell} ipython3
+:tags: [remove-cell]
+
+# To get a good resolution for displayed figures
+from matplotlib import pyplot
+pyplot.rcParams['figure.dpi'] = 600
+pyplot.rcParams['savefig.dpi'] = 600
+pyplot.rcParams['font.size'] = 9
+```
+
+## Reproject
+
+{func}`geoutils.Raster.reproject` or {func}`geoutils.Vector.reproject`.
+
+Reprojections **transform geospatial data from one CRS to another**.
+
+For vectors, the transformation of geometry points is exact. However, in the case of rasters, the projected points
+do not necessarily fall on a regular grid and require re-gridding by a 2D resampling algorithm, which results in a slight
+loss of information (value interpolation, propagation of nodata).
+
+For rasters, it can be useful to use {func}`~geoutils.Raster.reproject` in the same CRS simply for re-gridding,
+for instance when downsampling to a new resolution {attr}`~geoutils.Raster.res`.
+
+```{tip}
+Due to the loss of information when re-gridding, it is important to **minimize the number of reprojections during the
+analysis of rasters** (performing only one, if possible). For the same reason, when comparing vectors and rasters in
+different CRSs, it is usually **better to reproject the vector with no loss of information, which is the default
+behaviour of GeoUtils in raster–vector–point interfacing**.
+```
+
+```{code-cell} ipython3
+:tags: [hide-cell]
+:mystnb:
+:  code_prompt_show: "Show the code for opening example files"
+:  code_prompt_hide: "Hide the code for opening example files"
+
+import matplotlib.pyplot as plt
+import geoutils as gu
+rast = gu.Raster(gu.examples.get_path("everest_landsat_b4"))
+rast.set_nodata(0)  # Annoying to have to do this here, should we update it in the example?
+rast2 = gu.Raster(gu.examples.get_path("everest_landsat_b4_cropped"))
+vect = gu.Vector(gu.examples.get_path("everest_rgi_outlines"))
+```
+
+```{code-cell} ipython3
+# Reproject vector to CRS of raster by simply passing the raster
+vect_reproj = vect.reproject(rast)
+# Reproject raster to smaller bounds and different X/Y resolution
+rast_reproj = rast.reproject(
+    res=(rast.res[0] * 2, rast.res[1] / 2),
+    bounds={"left": rast.bounds.left, "bottom": rast.bounds.bottom,
+            "right": rast.bounds.left + 10000, "top": rast.bounds.bottom + 10000},
+    resampling="cubic")
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show the code for plotting the figure"
+:  code_prompt_hide: "Hide the code for plotting the figure"
+
+f, ax = plt.subplots(1, 2)
+ax[0].set_title("Before reprojection")
+rast.plot(ax=ax[0], cmap="gray", add_cbar=False)
+vect.plot(rast, ax=ax[0], ec="k", fc="none")
+ax[1].set_title("After reprojection")
+rast_reproj.plot(ax=ax[1], cmap="gray", add_cbar=False)
+vect_reproj.plot(ax=ax[1], ec="k", fc="none")
+_ = ax[1].set_yticklabels([])
+plt.tight_layout()
+```
+
+```{note}
+In GeoUtils, `"bilinear"` is the default resampling method. A simple {class}`str` matching the naming of a {class}`rasterio.enums.Resampling` method can be
+passed.
+
+Resampling methods are listed in **[the dedicated section of Rasterio's API](https://rasterio.readthedocs.io/en/latest/api/rasterio.enums.html#rasterio.enums.Resampling)**.
+```
+
+We can also simply pass another raster as reference to reproject to match the same CRS, and re-grid to the same bounds
+and resolution:
+
+```{code-cell} ipython3
+---
+mystnb:
+  output_stderr: warn
+---
+# Reproject vector to CRS of raster by simply passing the raster
+rast_reproj2 = rast.reproject(rast2)
+```
+
+GeoUtils raises a warning because the rasters have different {ref}`Pixel interpretation<pixel-interpretation>`,
+to ensure this is intended. This warning can be turned off at the package-level using GeoUtils' {ref}`config`.
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show the code for plotting the figure"
+:  code_prompt_hide: "Hide the code for plotting the figure"
+
+f, ax = plt.subplots(1, 3)
+ax[0].set_title("Raster 1")
+rast.plot(ax=ax[0], cmap="gray", add_cbar=False)
+vect.plot(rast, ax=ax[0], ec="k", fc="none")
+ax[1].set_title("Raster 2")
+rast2.plot(ax=ax[1], cmap="Reds", add_cbar=False)
+vect.plot(rast, ax=ax[1], ec="k", fc="none")
+ax[2].set_title("Match-ref\nreprojection")
+rast_reproj2.plot(ax=ax[2], cmap="gray", add_cbar=False)
+vect_reproj.plot(ax=ax[2], ec="k", fc="none")
+_ = ax[1].set_yticklabels([])
+_ = ax[2].set_yticklabels([])
+plt.tight_layout()
+```
+
+## Crop or pad
+
+{func}`geoutils.Raster.crop` or {func}`geoutils.Vector.crop`.
+
+Cropping **modifies the spatial bounds of the geospatial data in a rectangular extent**, by removing or adding data
+(in which case it corresponds to padding) without resampling.
+
+For rasters, cropping removes or adds pixels to the sides of the raster grid.
+
+For vectors, cropping removes some geometry features around the bounds, with three options possible:
+1. Removing all features **not intersecting** the cropping geometry,
+2. Removing all features **not contained** in the cropping geometry,
+3. Making all features **exactly clipped** to the cropping geometry (modifies the geometry data).
+
+```{code-cell} ipython3
+# Clip the vector to the raster
+vect_clipped = vect_reproj.crop(rast, clip=True)
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show the code for plotting the figure"
+:  code_prompt_hide: "Hide the code for plotting the figure"
+
+f, ax = plt.subplots(1, 2)
+ax[0].set_title("Before clipping")
+rast.plot(ax=ax[0], cmap="gray", add_cbar=False)
+vect_reproj.plot(ax=ax[0], ec="k", fc="none")
+ax[1].set_title("After clipping")
+rast.plot(ax=ax[1], cmap="gray", add_cbar=False)
+vect_clipped.plot(ax=ax[1], ec="k", fc="none")
+_ = ax[1].set_yticklabels([])
+plt.tight_layout()
+```
+
+## Translate
+
+{func}`geoutils.Raster.translate` or {func}`geoutils.Vector.translate`.
+
+Translations **modifies the georeferencing of the data by a horizontal offset** without modifying the underlying data,
+which is especially useful to align the data due to positioning errors.
+
+```{code-cell} ipython3
+# Translate the raster by a certain offset
+rast_shift = rast.translate(xoff=1000, yoff=1000)
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show the code for plotting the figure"
+:  code_prompt_hide: "Hide the code for plotting the figure"
+
+f, ax = plt.subplots(1, 2)
+ax[0].set_title("Before translation")
+rast.plot(ax=ax[0], cmap="gray", add_cbar=False)
+vect_clipped.plot(ax=ax[0], ec="k", fc="none")
+ax[1].set_title("After translation")
+rast_shift.plot(ax=ax[1], cmap="gray", add_cbar=False)
+vect_clipped.plot(ax=ax[1], ec="k", fc="none")
+_ = ax[1].set_yticklabels([])
+plt.tight_layout()
+```
+
+:::{admonition} See also
+:class: tip
+
+For 3D coregistration tailored to georeferenced elevation data, see [xDEM's coregistration module](https://xdem.readthedocs.io/en/stable/coregistration.html).
+:::
+
+## Merge
+
+{func}`geoutils.raster.merge_rasters()`
+
+Merge operations **join multiple geospatial data spatially, possibly with different georeferencing, into a single geospatial
+data object**.
+
+For rasters, the merging operation consists in combining all rasters into a single, larger raster. Pixels that overlap
+are combined by a reductor function (defaults to the mean). The output georeferenced grid (CRS, transform and shape) can
+be set to that of any reference raster (defaults to the extent that contains exactly all rasters).
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show the code for creating multiple raster pieces"
+:  code_prompt_hide: "Show the code for creating multiple raster pieces"
+
+# Get 4 cropped bits from initial rasters
+rast1 = rast.crop((rast.bounds.left + 1000, rast.bounds.bottom + 1000,
+                   rast.bounds.left + 3000, rast.bounds.bottom + 3000))
+rast2 = rast.crop((rast.bounds.left + 3000, rast.bounds.bottom + 1000,
+                   rast.bounds.left + 5000, rast.bounds.bottom + 3000))
+rast3 = rast.crop((rast.bounds.left + 1000, rast.bounds.bottom + 3000,
+                   rast.bounds.left + 3000, rast.bounds.bottom + 5000))
+rast4 = rast.crop((rast.bounds.left + 3000, rast.bounds.bottom + 3000,
+                   rast.bounds.left + 5000, rast.bounds.bottom + 5000))
+# Reproject some in other CRS, with other resolution
+#rast3 = rast3.reproject(crs=4326, res=rast.res[0] * 3)
+#rast4 = rast4.reproject(crs=32610, res=rast.res[0] / 3)
+```
+
+```{code-cell} ipython3
+---
+mystnb:
+  output_stderr: remove
+---
+# Merging all rasters, uses first raster's CRS, res, and the extent of all by default
+merged_rast = gu.raster.merge_rasters([rast1, rast2, rast3, rast4])
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show the code for plotting the figure"
+:  code_prompt_hide: "Hide the code for plotting the figure"
+
+f, ax = plt.subplots(1, 4)
+ax[0].set_title("Raster 1")
+rast1.plot(ax=ax[0], cmap="gray", add_cbar=False)
+ax[1].set_title("Raster 2")
+rast2.plot(ax=ax[1], cmap="gray", add_cbar=False)
+ax[2].set_title("Raster 3")
+rast3.plot(ax=ax[2], cmap="gray", add_cbar=False)
+ax[3].set_title("Raster 4")
+rast4.plot(ax=ax[3], cmap="gray", add_cbar=False)
+plt.tight_layout()
+
+merged_rast.plot(ax="new", cmap="gray", add_cbar=False)
+```
diff --git a/doc/source/index.md b/doc/source/index.md
index 229d13b4..d8fb8e12 100644
--- a/doc/source/index.md
+++ b/doc/source/index.md
@@ -31,6 +31,16 @@ GeoUtils ``v0.1`` is released, with most features drafted 3 years ago now finali
 
 ----------------
 
+GeoUtils is built on top of core geospatial packages (Rasterio, GeoPandas, PyProj) and numerical packages
+(NumPy, Xarray, SciPy) to provide **consistent higher-level functionalities at the interface of raster, vector and point
+cloud objects** (such as match-reference reprojection, point interpolation or gridding).
+
+It is **tailored to perform quantitative analysis that implicitly understands the intricacies of geospatial data**
+(nodata values, projection, pixel interpretation), through **an intuitive object-based API to foster accessibility**,
+and strives **to be computationally scalable** through Dask.
+
+If you are looking to **port your GDAL or QGIS workflow in Python**, GeoUtils is made for you!
+
 # Where to start?
 
 ::::{grid} 1 2 2 3
@@ -95,9 +105,11 @@ feature_overview
 :maxdepth: 2
 
 core_index
-rasters_index
-vectors_index
-proj_tools
+data_object_index
+georeferencing
+geotransformations
+raster_vector_point
+distance_ops
 ```
 
 ```{toctree}
@@ -115,6 +127,7 @@ analysis_examples/index
 
 api
 cli
+config
 background
 ```
 
diff --git a/doc/source/pointcloud_class.md b/doc/source/pointcloud_class.md
new file mode 100644
index 00000000..68b6b56d
--- /dev/null
+++ b/doc/source/pointcloud_class.md
@@ -0,0 +1,8 @@
+(point-cloud)=
+# The georeferenced point cloud ({class}`~geoutils.PointCloud`)
+
+Coming soon!
+
+In the meantime, most point cloud operations (for instance in {ref}`raster-vector-point`) return a
+{class}`~geoutils.Vector` with only point geometries and a specific `data_column_name` corresponding to the point
+cloud values.
diff --git a/doc/source/proj_tools.md b/doc/source/proj_tools.md
deleted file mode 100644
index d58162fe..00000000
--- a/doc/source/proj_tools.md
+++ /dev/null
@@ -1,67 +0,0 @@
----
-file_format: mystnb
-jupytext:
-  formats: md:myst
-  text_representation:
-    extension: .md
-    format_name: myst
-kernelspec:
-  display_name: geoutils-env
-  language: python
-  name: geoutils
----
-(proj-tools)=
-
-# Projection tools
-
-This section describes projection tools that are common to {class}`Rasters<geoutils.Raster>` and {class}`Vectors<geoutils.Vector>`, and facilitate
-geospatial analysis.
-
-## Get a metric coordinate system
-
-A local metric coordinate system can be estimated for both {class}`Rasters<geoutils.Raster>` and {class}`Vectors<geoutils.Vector>` through the
-{func}`~geoutils.Raster.get_metric_crs` function.
-
-The metric system returned can be either "universal" (zone of the Universal Transverse Mercator or Universal Polar Stereographic system), or "custom"
-(Mercator or Polar projection centered on the {class}`Raster<geoutils.Raster>` or {class}`Vector<geoutils.Vector>`).
-
-```{code-cell} ipython3
-import geoutils as gu
-
-# Initiate a raster from disk
-rast = gu.Raster(gu.examples.get_path("exploradores_aster_dem"))
-rast.info()
-
-# Estimate a universal metric CRS for the raster
-rast.get_metric_crs()
-```
-
-## Get projected bounds
-
-Projected bounds can be directly derived from both {class}`Rasters<geoutils.Raster>` and {class}`Vectors<geoutils.Vector>` through the
-{func}`~geoutils.Raster.get_bounds_projected` function.
-
-```{code-cell} ipython3
-# Get raster bounds in geographic CRS by passing its EPSG code
-rast.get_bounds_projected(4326)
-```
-
-```{important}
-When projecting to a new CRS, the footprint shape of the data is generally deformed. To account for this, use {func}`~geoutils.Raster.get_footprint_projected` described below.
-```
-
-## Get projected footprint
-
-A projected footprint can be derived from both {class}`Rasters<geoutils.Raster>` and {class}`Vectors<geoutils.Vector>` through the
-{func}`~geoutils.Raster.get_footprint_projected` function.
-
-For this, the original rectangular footprint polygon lines are densified to respect the deformation during reprojection.
-
-```{code-cell} ipython3
-# Get raster footprint in geographic CRS
-rast_footprint = rast.get_footprint_projected(4326)
-
-rast_footprint.plot()
-```
-
-This is for instance useful to check for intersection with other data.
diff --git a/doc/source/quick_start.md b/doc/source/quick_start.md
index dd120038..0ea9e5b6 100644
--- a/doc/source/quick_start.md
+++ b/doc/source/quick_start.md
@@ -18,6 +18,16 @@ A short code example using several end-user functionalities of GeoUtils. For a m
 have a look at the {ref}`feature-overview` page! Or, to find an example about
 a specific functionality, jump to {ref}`quick-gallery` right below.
 
+```{code-cell} ipython3
+:tags: [remove-cell]
+
+# To get a good resolution for displayed figures
+from matplotlib import pyplot
+pyplot.rcParams['figure.dpi'] = 600
+pyplot.rcParams['savefig.dpi'] = 600
+pyplot.rcParams['font.size'] = 9
+```
+
 ## Short example
 
 ```{note}
diff --git a/doc/source/raster_vector_point.md b/doc/source/raster_vector_point.md
new file mode 100644
index 00000000..8088676e
--- /dev/null
+++ b/doc/source/raster_vector_point.md
@@ -0,0 +1,287 @@
+---
+file_format: mystnb
+jupytext:
+  formats: md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+kernelspec:
+  display_name: geoutils-env
+  language: python
+  name: geoutils
+---
+(raster-vector-point)=
+# Raster–vector–point interface
+
+GeoUtils provides functionalities at the interface of rasters, vectors and point clouds, allowing to consistently perform
+operations such as mask creation or point interpolation **respecting both georeferencing and nodata values, as well
+as pixel interpretation for point interfacing**.
+
+```{code-cell} ipython3
+:tags: [remove-cell]
+
+# To get a good resolution for displayed figures
+from matplotlib import pyplot
+pyplot.rcParams['figure.dpi'] = 600
+pyplot.rcParams['savefig.dpi'] = 600
+pyplot.rcParams['font.size'] = 9
+```
+
+## Raster–vector operations
+
+### Rasterize
+
+{func}`geoutils.Vector.rasterize`
+
+Rasterization of a vector is **an operation that allows to translate some information of the vector data into a raster**, by
+setting the values raster pixels intersecting a vector geometry feature to that of an attribute of the vector
+associated to the geometry (e.g., feature ID, area or any other value), which is the geometry index by default.
+
+Rasterization generally implies some loss of information, as there is no exact way of representing a vector on a grid.
+Rather, the choice of which pixels are attributed a value depends on the amount of intersection with the vector
+geometries and so includes several options (percent of area intersected, all touched, etc).
+
+```{code-cell} ipython3
+:tags: [hide-cell]
+:mystnb:
+:  code_prompt_show: "Show the code for opening example files"
+:  code_prompt_hide: "Hide the code for opening example files"
+
+import matplotlib.pyplot as plt
+import geoutils as gu
+import numpy as np
+
+rast = gu.Raster(gu.examples.get_path("everest_landsat_b4"))
+rast.set_nodata(0)  # Annoying to have to do this here, should we update it in the example?
+vect = gu.Vector(gu.examples.get_path("everest_rgi_outlines"))
+```
+
+```{code-cell} ipython3
+# Rasterize the vector features based on their glacier ID number
+rasterized_vect = vect.rasterize(rast)
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show the code for plotting the figure"
+:  code_prompt_hide: "Hide the code for plotting the figure"
+
+f, ax = plt.subplots(1, 2)
+ax[0].set_title("Raster and vector")
+rast.plot(ax=ax[0], cmap="gray", add_cbar=False)
+vect.plot(ref_crs=rast, ax=ax[0], ec="k", fc="none")
+ax[1].set_title("Rasterized vector")
+rasterized_vect.plot(ax=ax[1], cmap="viridis", cbar_title="Feature index")
+_ = ax[1].set_yticklabels([])
+plt.tight_layout()
+```
+
+### Create mask
+
+{func}`geoutils.Vector.create_mask`
+
+Mask creation from a vector **is a rasterization of all vector features that only categorizes geometry intersection as a boolean mask**
+(if any feature falls in a given pixel or not), and is therefore independent of any vector attribute values.
+
+```{code-cell} ipython3
+# Create a boolean mask from all vector features
+mask = vect.create_mask(rast)
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show the code for plotting the figure"
+:  code_prompt_hide: "Hide the code for plotting the figure"
+
+f, ax = plt.subplots(1, 2)
+ax[0].set_title("Raster and vector")
+rast.plot(ax=ax[0], cmap="gray", add_cbar=False)
+vect.plot(ref_crs=rast, ax=ax[0], ec="k", fc="none")
+ax[1].set_title("Mask from vector")
+mask.plot(ax=ax[1], cbar_title="Intersects vector (1=yes, 0=no)")
+_ = ax[1].set_yticklabels([])
+plt.tight_layout()
+```
+
+It returns a {class}`~geoutils.Mask`, a georeferenced boolean raster (or optionally, a boolean NumPy array), which
+can both be used for indexing or index assignment of a raster.
+
+```{code-cell} ipython3
+# Mean of values in the mask
+np.mean(rast[mask])
+```
+
+### Polygonize
+
+{func}`geoutils.Raster.polygonize`
+
+Polygonization of a raster **consists of delimiting contiguous raster pixels with the same target values into vector polygon
+geometries**. By default, all raster values are used as targets. When using polygonize on a {class}`~geoutils.Mask`,
+the targets are implicitly the valid values of the mask.
+
+```{code-cell} ipython3
+# Mask 0 values
+rasterized_vect.set_mask(rasterized_vect == 0)
+# Polygonize all non-zero values
+vect_repolygonized = rasterized_vect.polygonize()
+
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show the code for plotting the figure"
+:  code_prompt_hide: "Hide the code for plotting the figure"
+
+f, ax = plt.subplots(1, 2)
+ax[0].set_title("Raster (vector\n rasterized above)")
+rasterized_vect.plot(ax=ax[0], cmap="viridis", cbar_title="Feature index")
+ax[1].set_title("Polygonized raster")
+vect_repolygonized.plot(ax=ax[1], column="id", fc="none", cbar_title="Feature index")
+_ = ax[1].set_yticklabels([])
+plt.tight_layout()
+```
+
+## Raster–point operations
+
+### Point interpolation
+
+{func}`geoutils.Raster.interp_points`
+
+Point interpolation of a raster **consists in estimating the values at exact point coordinates by 2D regular-grid
+interpolation** such as nearest neighbour, bilinear (default), cubic, etc.
+
+```{note}
+In order to support all types of resampling methods with nodata values while maintaining the robustness of results,
+GeoUtils implements **a modified version of {func}`scipy.interpolate.interpn` that propagates nodata
+values** in surrounding pixels of initial nodata values depending on the order of the resampling method:
+- Nearest or linear (order 0 or 1): up to 1 pixel,
+- Cubic (order 3): 2 pixels,
+- Quintic (order 5): 3 pixels.
+```
+
+```{code-cell} ipython3
+# Let's change raster type for a DEM, often requiring interpolation
+rast =  gu.Raster(gu.examples.get_path("exploradores_aster_dem"))
+
+# Get 50 random points to sample within the raster extent
+rng = np.random.default_rng(42)
+x_coords = rng.uniform(rast.bounds.left, rast.bounds.right, 50)
+y_coords = rng.uniform(rast.bounds.bottom, rast.bounds.top, 50)
+
+vals = rast.interp_points(points=(x_coords, y_coords))
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show the code for plotting the figure"
+:  code_prompt_hide: "Hide the code for plotting the figure"
+
+f, ax = plt.subplots(1, 2)
+ax[0].set_title("Raster")
+rast.plot(ax=ax[0], cmap="terrain", cbar_title="Elevation (m)")
+ax[1].set_title("Interpolated\npoint cloud")
+# (Update with release of PointCloud class)
+import geopandas as gpd
+pc = gu.Vector(gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=x_coords, y=y_coords), data={"b1": vals}, crs=rast.crs))
+pc.plot(column="b1", ax=ax[1], cmap="terrain", legend=True, marker="x", cbar_title="Elevation (m)")
+_ = ax[1].set_yticklabels([])
+plt.tight_layout()
+```
+
+### Reduction around point
+
+Point reduction of a raster is **the estimation of the values at point coordinates by applying a reductor function (e.g., mean,
+median) to pixels contained in a window centered on the point**. For a window smaller than the pixel size, the value of
+the closest pixel is returned.
+
+{func}`geoutils.Raster.reduce_points`
+
+```{code-cell} ipython3
+vals = rast.reduce_points((x_coords, y_coords), window=5, reducer_function=np.nanmedian)
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show the code for plotting the figure"
+:  code_prompt_hide: "Hide the code for plotting the figure"
+
+f, ax = plt.subplots(1, 2)
+ax[0].set_title("Raster")
+rast.plot(ax=ax[0], cmap="terrain", cbar_title="Elevation (m)")
+ax[1].set_title("Reduced\npoint cloud")
+# (Update with release of PointCloud class)
+import geopandas as gpd
+pc = gu.Vector(gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=x_coords, y=y_coords), data={"b1": vals}, crs=rast.crs))
+pc.plot(column="b1", ax=ax[1], cmap="terrain", legend=True, cbar_title="Elevation (m)")
+_ = ax[1].set_yticklabels([])
+plt.tight_layout()
+```
+
+
+### Raster to points
+
+{func}`geoutils.Raster.to_pointcloud`
+
+**A raster can be converted exactly into a point cloud**, which each pixel in the raster is associated to its pixel
+values to create a point cloud on a regular grid.
+
+```{code-cell} ipython3
+pc = rast.to_pointcloud(subsample=10000)
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show the code for plotting the figure"
+:  code_prompt_hide: "Hide the code for plotting the figure"
+
+f, ax = plt.subplots(1, 2)
+ax[0].set_title("Raster")
+rast.plot(ax=ax[0], cmap="terrain", cbar_title="Elevation (m)")
+ax[1].set_title("Regular subsampled\npoint cloud")
+pc.plot(column="b1", ax=ax[1], cmap="terrain", legend=True, cbar_title="Elevation (m)", markersize=2)
+_ = ax[1].set_yticklabels([])
+plt.tight_layout()
+```
+
+
+### Regular points to raster
+
+{func}`geoutils.Raster.from_pointcloud_regular`
+
+**If a point cloud is regularly spaced in X and Y coordinates, it can be converted exactly into a raster**. Otherwise,
+it must be re-gridded using {ref}`point-gridding` described below. For a regular point cloud, every point is associated to a
+pixel in the raster grid, and the values are set to the raster. The point cloud does not necessarily need to contain
+points for all grid coordinates, as pixels with no corresponding point are set to nodata values.
+
+```{code-cell} ipython3
+rast_from_pc = gu.Raster.from_pointcloud_regular(pc, transform=rast.transform, shape=rast.shape)
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+:mystnb:
+:  code_prompt_show: "Show the code for plotting the figure"
+:  code_prompt_hide: "Hide the code for plotting the figure"
+
+f, ax = plt.subplots(1, 2)
+ax[0].set_title("Regular subsampled\npoint cloud")
+pc.plot(column="b1", ax=ax[0], cmap="terrain", legend=True, cbar_title="Elevation (m)", markersize=2)
+ax[1].set_title("Raster from\npoint cloud")
+rast_from_pc.plot(ax=ax[1], cmap="terrain", cbar_title="Elevation (m)")
+_ = ax[1].set_yticklabels([])
+plt.tight_layout()
+```
+
+
+(point-gridding)=
+### Point gridding
+
+{func}`geoutils.PointCloud.grid`
+
+Coming with the next release on point clouds!
diff --git a/doc/source/vectors_index.md b/doc/source/vectors_index.md
deleted file mode 100644
index 735530fc..00000000
--- a/doc/source/vectors_index.md
+++ /dev/null
@@ -1,10 +0,0 @@
-(vectors-index)=
-# Vectors
-
-Prefer to learn by running examples? Explore our example galleries on {ref}`examples-io`, {ref}`examples-handling` and {ref}`examples-analysis`.
-
-```{toctree}
-:maxdepth: 2
-
-vector_class
-```
diff --git a/examples/analysis/geospatial/proximity_metric.py b/examples/analysis/geospatial/proximity_metric.py
index b7b2f02c..c73964be 100644
--- a/examples/analysis/geospatial/proximity_metric.py
+++ b/examples/analysis/geospatial/proximity_metric.py
@@ -15,7 +15,7 @@
 rast = gu.Raster(filename_rast)
 vect = gu.Vector(filename_vect)
 vect = vect[vect["RGIId"] == "RGI60-15.10055"]
-rast.crop(vect)
+rast = rast.crop(vect)
 
 # Plot the raster and vector
 rast.plot(cmap="Blues")
diff --git a/examples/analysis/point_extraction/interpolation.py b/examples/analysis/point_extraction/interpolation.py
index 0a89cec9..5682b717 100644
--- a/examples/analysis/point_extraction/interpolation.py
+++ b/examples/analysis/point_extraction/interpolation.py
@@ -12,7 +12,7 @@
 
 filename_rast = gu.examples.get_path("exploradores_aster_dem")
 rast = gu.Raster(filename_rast)
-rast.crop([rast.bounds.left, rast.bounds.bottom, rast.bounds.left + 2000, rast.bounds.bottom + 2000])
+rast = rast.crop([rast.bounds.left, rast.bounds.bottom, rast.bounds.left + 2000, rast.bounds.bottom + 2000])
 
 # Plot the raster
 rast.plot(cmap="terrain")
@@ -44,8 +44,8 @@
 # %%
 # Let's look and redefine our pixel interpretation into ``"Point"``. This will shift interpolation by half a pixel.
 
-print(rast.tags["AREA_OR_POINT"])
-rast.tags["AREA_OR_POINT"] = "Point"
+rast.area_or_point
+rast.area_or_point = "Point"
 
 # %%
 # We can interpolate again by shifting according to our interpretation, and changing the resampling algorithm (default to "linear").
diff --git a/examples/analysis/point_extraction/reduction.py b/examples/analysis/point_extraction/reduction.py
index d9f68d3e..533c2d5d 100644
--- a/examples/analysis/point_extraction/reduction.py
+++ b/examples/analysis/point_extraction/reduction.py
@@ -12,7 +12,7 @@
 
 filename_rast = gu.examples.get_path("exploradores_aster_dem")
 rast = gu.Raster(filename_rast)
-rast.crop([rast.bounds.left, rast.bounds.bottom, rast.bounds.left + 2000, rast.bounds.bottom + 2000])
+rast = rast.crop([rast.bounds.left, rast.bounds.bottom, rast.bounds.left + 2000, rast.bounds.bottom + 2000])
 
 # Plot the raster
 rast.plot(cmap="terrain")
diff --git a/examples/handling/georeferencing/crop_raster.py b/examples/handling/georeferencing/crop_raster.py
index bdb7cbb9..369b8380 100644
--- a/examples/handling/georeferencing/crop_raster.py
+++ b/examples/handling/georeferencing/crop_raster.py
@@ -30,7 +30,7 @@
 # **First option:** using the vector as a reference to match, we reproject the raster. We simply have to pass the :class:`~geoutils.Vector`
 # as single argument to :func:`~geoutils.Raster.crop`. See :ref:`core-match-ref` for more details.
 
-rast.crop(vect, inplace=True)
+rast = rast.crop(vect)
 
 # %%
 # Now the bounds should be the same as that of the vector (within the size of a pixel as the grid was not warped).
@@ -46,7 +46,7 @@
 # **Second option:** we can pass other ``crop_geom`` argument to :func:`~geoutils.Raster.crop`, including another :class:`~geoutils.Raster` or a
 # simple :class:`tuple` of bounds. For instance, we can re-crop the raster to be smaller than the vector.
 
-rast.crop((rast.bounds.left + 1000, rast.bounds.bottom, rast.bounds.right, rast.bounds.top - 500), inplace=True)
+rast = rast.crop((rast.bounds.left + 1000, rast.bounds.bottom, rast.bounds.right, rast.bounds.top - 500))
 
 rast.plot(ax="new", cmap="Purples")
 vect.plot(ref_crs=rast, fc="none", ec="k", lw=2)
diff --git a/examples/handling/georeferencing/crop_vector.py b/examples/handling/georeferencing/crop_vector.py
index f2a5df0b..65587bb7 100644
--- a/examples/handling/georeferencing/crop_vector.py
+++ b/examples/handling/georeferencing/crop_vector.py
@@ -24,7 +24,7 @@
 # **First option:** using the raster as a reference to match, we crop the vector. We simply have to pass the :class:`~geoutils.Raster` as single argument to
 # :func:`~geoutils.Vector.crop`. See :ref:`core-match-ref` for more details.
 
-vect.crop(rast, inplace=True)
+vect = vect.crop(rast)
 
 # %%
 # .. note::
@@ -38,7 +38,7 @@
 # The :func:`~geoutils.Vector.crop` keeps all features with geometries intersecting the extent to crop to. We can also force a clipping of the geometries
 # within the bounds using ``clip=True``.
 
-vect.crop(rast, clip=True)
+vect = vect.crop(rast, clip=True)
 rast.plot(ax="new", cmap="Greys_r", alpha=0.7)
 vect.plot(ref_crs=rast, fc="none", ec="tab:purple", lw=3)
 
@@ -47,9 +47,7 @@
 # simple :class:`tuple` of bounds.
 
 bounds = rast.get_bounds_projected(out_crs=vect.crs)
-vect.crop(
-    crop_geom=(bounds.left + 0.5 * (bounds.right - bounds.left), bounds.bottom, bounds.right, bounds.top), inplace=True
-)
+vect = vect.crop(crop_geom=(bounds.left + 0.5 * (bounds.right - bounds.left), bounds.bottom, bounds.right, bounds.top))
 
 rast.plot(ax="new", cmap="Greys_r", alpha=0.7)
 vect.plot(ref_crs=rast, fc="none", ec="tab:purple", lw=3)
diff --git a/examples/handling/interface/polygonize.py b/examples/handling/interface/polygonize.py
index c1e697fd..ea8d1280 100644
--- a/examples/handling/interface/polygonize.py
+++ b/examples/handling/interface/polygonize.py
@@ -12,7 +12,7 @@
 
 filename_rast = gu.examples.get_path("exploradores_aster_dem")
 rast = gu.Raster(filename_rast)
-rast.crop([rast.bounds.left, rast.bounds.bottom, rast.bounds.left + 5000, rast.bounds.bottom + 5000])
+rast = rast.crop([rast.bounds.left, rast.bounds.bottom, rast.bounds.left + 5000, rast.bounds.bottom + 5000])
 # %%
 # Let's plot the raster.
 rast.plot(cmap="terrain")
diff --git a/examples/handling/interface/topoints.py b/examples/handling/interface/topoints.py
index 56e8b071..fcf3df54 100644
--- a/examples/handling/interface/topoints.py
+++ b/examples/handling/interface/topoints.py
@@ -12,7 +12,7 @@
 
 filename_rast = gu.examples.get_path("exploradores_aster_dem")
 rast = gu.Raster(filename_rast)
-rast.crop([rast.bounds.left, rast.bounds.bottom, rast.bounds.left + 500, rast.bounds.bottom + 500])
+rast = rast.crop([rast.bounds.left, rast.bounds.bottom, rast.bounds.left + 500, rast.bounds.bottom + 500])
 
 # %%
 # Let's plot the raster.
diff --git a/examples/io/import_export/import_vector.py b/examples/io/import_export/import_vector.py
index 231779e5..6ac821c4 100644
--- a/examples/io/import_export/import_vector.py
+++ b/examples/io/import_export/import_vector.py
@@ -20,7 +20,7 @@
 # %%
 # We plot the vector.
 
-vect.plot(column="RGIId")
+vect.plot(column="RGIId", add_cbar=False)
 
 # %%
 # To export, the :class:`geopandas.GeoDataFrame` is always stored as an attribute as :class:`~geoutils.Vector` is composed from it. See :ref:`core-composition`.
diff --git a/geoutils/config.ini b/geoutils/config.ini
index 931d5a9f..6cd070dd 100644
--- a/geoutils/config.ini
+++ b/geoutils/config.ini
@@ -1,6 +1,9 @@
-# This script sets default global parameters for GeoUtils, which can be customized by the user
+# Default global parameters for GeoUtils
 
-# Sections are useless for now, but might become useful later on
 [raster]
+
+# Shift by half a pixel for "Point" pixel interpretation during raster-point operations
 shift_area_or_point = True
+
+# Raise a warning if two rasters have different pixel interpretation
 warn_area_or_point = True
diff --git a/geoutils/raster/raster.py b/geoutils/raster/raster.py
index c817c6dd..fd2e0e47 100644
--- a/geoutils/raster/raster.py
+++ b/geoutils/raster/raster.py
@@ -3286,7 +3286,7 @@ def plot(
         # Add colorbar
         if add_cbar:
             divider = make_axes_locatable(ax0)
-            cax = divider.append_axes("right", size="5%", pad=0.05)
+            cax = divider.append_axes("right", size="5%", pad="2%")
             norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
             cbar = matplotlib.colorbar.ColorbarBase(cax, cmap=cmap, norm=norm)
             cbar.solids.set_alpha(alpha)
@@ -3296,9 +3296,11 @@ def plot(
         else:
             cbar = None
 
+        plt.sca(ax0)
+
         # If returning axes
         if return_axes:
-            return ax, cax
+            return ax0, cax
         else:
             return None
 
@@ -4173,7 +4175,8 @@ def polygonize(
 
         # Mask all valid values
         elif target_values == "all":
-            bool_msk = (~self.data.mask).astype("uint8")
+            # Using getmaskarray is necessary in case .data.mask is nomask (False)
+            bool_msk = (~np.ma.getmaskarray(self.data)).astype("uint8")
 
         else:
             raise ValueError("in_value must be a number, a tuple or a sequence")
@@ -4242,7 +4245,15 @@ def proximity(
             distance_unit=distance_unit,
         )
 
-        return self.copy(new_array=proximity)
+        out_nodata = _default_nodata(proximity.dtype)
+        return self.from_array(
+            data=proximity,
+            transform=self.transform,
+            crs=self.crs,
+            nodata=out_nodata,
+            area_or_point=self.area_or_point,
+            tags=self.tags,
+        )
 
     @overload
     def subsample(
diff --git a/geoutils/vector.py b/geoutils/vector.py
index cbed8641..815249d5 100644
--- a/geoutils/vector.py
+++ b/geoutils/vector.py
@@ -165,7 +165,7 @@ def info(self, verbose: bool = True) -> str | None:
         """
         as_str = [  # 'Driver:             {} \n'.format(self.driver),
             f"Filename:           {self.name} \n",
-            f"Coordinate System:  EPSG:{self.ds.crs.to_epsg()}\n",
+            f"Coordinate system:  EPSG:{self.ds.crs.to_epsg()}\n",
             f"Extent:             {self.ds.total_bounds.tolist()} \n",
             f"Number of features: {len(self.ds)} \n",
             f"Attributes:         {self.ds.columns.tolist()}",
@@ -185,7 +185,7 @@ def plot(
         vmax: float | int | None = None,
         alpha: float | int | None = None,
         cbar_title: str | None = None,
-        add_cbar: bool = False,
+        add_cbar: bool = True,
         ax: matplotlib.axes.Axes | Literal["new"] | None = None,
         return_axes: bool = False,
         **kwargs: Any,
@@ -202,7 +202,7 @@ def plot(
         :param vmax: Colorbar maximum value. Default is data max.
         :param alpha: Transparency of raster and colorbar.
         :param cbar_title: Colorbar label. Default is None.
-        :param add_cbar: Set to True to display a colorbar. Default is True.
+        :param add_cbar: Set to True to display a colorbar. Default is True if a "column" argument is passed.
         :param ax: A figure ax to be used for plotting. If None, will plot on current axes. If "new",
             will create a new axis.
         :param return_axes: Whether to return axes.
@@ -228,6 +228,12 @@ def plot(
         else:
             raise ValueError("ax must be a matplotlib.axes.Axes instance, 'new' or None.")
 
+        # Set add_cbar depending on column argument
+        if "column" in kwargs.keys() and add_cbar:
+            add_cbar = True
+        else:
+            add_cbar = False
+
         # Update with this function's arguments
         if add_cbar:
             legend = True
@@ -236,29 +242,27 @@ def plot(
 
         if "legend" in list(kwargs.keys()):
             legend = kwargs.pop("legend")
-        else:
-            legend = False
 
         # Get colormap arguments that might have been passed in the keyword args
         if "legend_kwds" in list(kwargs.keys()) and legend:
             legend_kwds = kwargs.pop("legend_kwds")
-            if "label" in list(legend_kwds):
-                cbar_title = legend_kwds.pop("label")
+            if cbar_title is not None:
+                legend_kwds.update({"label": cbar_title})  # Pad updates depending on figsize during plot,
         else:
-            legend_kwds = None
+            if cbar_title is not None:
+                legend_kwds = {"label": cbar_title}
+            else:
+                legend_kwds = None
 
         # Add colorbar
         if add_cbar or cbar_title:
             divider = make_axes_locatable(ax0)
-            cax = divider.append_axes("right", size="5%", pad=0.05)
+            cax = divider.append_axes("right", size="5%", pad="2%")
             norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
             cbar = matplotlib.colorbar.ColorbarBase(
                 cax, cmap=cmap, norm=norm
             )  # , orientation="horizontal", ticklocation="top")
             cbar.solids.set_alpha(alpha)
-
-            if cbar_title is not None:
-                cbar.set_label(cbar_title)
         else:
             cax = None
             cbar = None
@@ -276,9 +280,13 @@ def plot(
             **kwargs,
         )
 
+        cax
+
+        plt.sca(ax0)
+
         # If returning axes
         if return_axes:
-            return ax, cax
+            return ax0, cax
         else:
             return None
 
@@ -1553,7 +1561,15 @@ def proximity(
             raster=raster, vector=self, geometry_type=geometry_type, in_or_out=in_or_out, distance_unit=distance_unit
         )
 
-        return raster.copy(new_array=proximity)
+        out_nodata = gu.raster.raster._default_nodata(proximity.dtype)
+        return gu.Raster.from_array(
+            data=proximity,
+            transform=raster.transform,
+            crs=raster.crs,
+            nodata=out_nodata,
+            area_or_point=raster.area_or_point,
+            tags=raster.tags,
+        )
 
     def buffer_metric(self, buffer_size: float) -> Vector:
         """
diff --git a/tests/test_vector.py b/tests/test_vector.py
index f8adf909..cabae434 100644
--- a/tests/test_vector.py
+++ b/tests/test_vector.py
@@ -82,7 +82,7 @@ def test_info(self) -> None:
         # Otherwise returns info
         output2 = v.info(verbose=False)
         assert isinstance(output2, str)
-        list_prints = ["Filename", "Coordinate System", "Extent", "Number of features", "Attributes"]
+        list_prints = ["Filename", "Coordinate system", "Extent", "Number of features", "Attributes"]
         assert all(p in output2 for p in list_prints)
 
     def test_query(self) -> None: