Incremental commit on doc

doc/source/api.md

@@ -84,8 +84,8 @@ documentation.
     Raster.reproject
     Raster.polygonize
     Raster.proximity
-    Raster.value_at_coords
     Raster.interp_points
+    Raster.reduce_points
 ```
 
 ### Plotting
@@ -121,7 +121,7 @@ documentation.
     Raster.load
     Raster.save
     Raster.to_pointcloud
-    Raster.from_regular_pointcloud
+    Raster.from_pointcloud_regular
     Raster.to_rio_dataset
     Raster.to_xarray
 ```
@@ -135,7 +135,7 @@ documentation.
     Raster.xy2ij
     Raster.ij2xy
     Raster.coords
-    Raster.shift
+    Raster.translate
     Raster.outside_image
 ```
 

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
+```

doc/source/distance_ops.md

@@ -14,7 +14,7 @@ kernelspec:
 # 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 later analysis. To facilitate this type of operations, GeoUtils implements distance-specific functionalities 
 for both vectors and rasters.
 
 ```{code-cell} ipython3
@@ -57,7 +57,7 @@ vect = gu.Vector(gu.examples.get_path("everest_rgi_outlines"))
 ```{code-cell} ipython3
 # Compute proximity to vector outlines
 proximity = vect.proximity(rast)
-proximity.plot(cmap="viridis")
+proximity.plot(cmap="viridis", cbar_title="Proximity (m)")
 ```
 
 ## Buffering without overlap
@@ -71,6 +71,6 @@ between shapes, which is sometimes undesirable. Using Voronoi polygons, we provi
 # Compute buffer without overlap on glacier outlines
 vect_buff_nolap = vect.buffer_without_overlap(buffer_size=500)
 # Plot with color to see that the attributes are retained for every feature
-vect_buff_nolap.plot(ax="new", column="Area")
-vect.plot(ec="k", column="Area", alpha=0.5)
+vect.plot(ax="new", ec="k", column="Area", alpha=0.5, add_cbar=False)
+vect_buff_nolap.plot(column="Area", cbar_title="Glacier area (km)")
 ```

doc/source/feature_overview.md

@@ -25,6 +25,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`.

doc/source/georeferencing.md

@@ -50,39 +50,51 @@ vect = gu.Vector(gu.examples.get_path("exploradores_rgi_outlines"))
 ```
 
 ```{code-cell} ipython3
-# Print raster and vector info
+# 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 a 
+[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`. More information on the manipulation 
-of {class}`pyproj.crs.CRS` objects can be found in [PyProj's documentation](https://pyproj4.github.io/pyproj/stable/).
+{class}`pyproj.crs.CRS` object in {attr}`~geoutils.Raster.crs`.
 
 ```{code-cell} ipython3
-# Show CRS attribute of a raster and vector
+# Show CRS attribute of raster
 print(rast.crs)
-print(vect.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 a raster and vector
-print(rast.bounds)
-print(vect.bounds)
+# Show bounds attribute of raster
+rast.bounds
+```
+```{code-cell} ipython3
+# Show bounds attribute of vector
+vect.bounds
 ```
 
 ```{note}
@@ -103,15 +115,20 @@ 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`, and its shape by the data array shape.
-From it are derived the resolution {attr}`~geoutils.Raster.res`, the 2D raster shape 
-{attr}`~geoutils.Raster.shape` made up of its {attr}`~geoutils.Raster.height` and {attr}`~geoutils.Raster.width`.
+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
@@ -126,7 +143,8 @@ A largely overlooked aspect of a raster's georeferencing is the pixel interpreta
 [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)**, used often for DEMs. 
+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. 
 
@@ -140,11 +158,11 @@ interpretation during analysis**, especially for raster–vector–point interfa
 or re-gridding, and might also be a problem if defined differently when comparing two rasters.
 
 ```{important}
-By default, **pixel interpretation can induce a half-pixel shift during raster–point interfacing** 
+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**.
+but only **raises a warning for raster–raster operations** if interpretations differ.
 
-This behaviour can be modified at the package-level by using [GeoUtils' configuration parameters]() 
+This behaviour can be modified at the package-level by using GeoUtils' {ref}`config` 
 `shift_area_or_point` and `warns_area_or_point`.
 ```
 
@@ -172,7 +190,7 @@ The metric system returned can be either "universal" (zone of the Universal Tran
 
 ```{code-cell} ipython3
 # Get local metric CRS
-print(rast.get_metric_crs())
+rast.get_metric_crs()
 ```
 
 ### Re-set georeferencing metadata
@@ -183,7 +201,6 @@ To add?
 {func}`geoutils.Raster.set_transform`
 {func}`geoutils.Raster.set_area_or_point`
 
-
 ### Ccoordinates 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.

doc/source/geotransformations.md

@@ -43,7 +43,7 @@ For rasters, it can be useful to use {func}`~geoutils.Raster.reproject` in the s
 for instance when downsampling to a new resolution {attr}`~geoutils.Raster.res`.
 
 ```{tip}
-Due to the loss information when re-gridding, it is important to **minimize the number of reprojections during the 
+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**.
@@ -82,12 +82,13 @@ rast_reproj = rast.reproject(
 
 f, ax = plt.subplots(1, 2)
 ax[0].set_title("Before reprojection")
-rast.plot(ax=ax[0])
+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])
+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}
@@ -101,13 +102,16 @@ We can also simply pass another raster as reference to reproject to match the sa
 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)
 ```
 
-As shown above, if the rasters have different {ref}`Pixel interpretation<pixel-interpretation>`, GeoUtils 
-raises a warning to ensure this is intended. 
-This warning can be turned off at the package-level using [GeoUtils' global configuration]().
+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]
@@ -117,16 +121,17 @@ This warning can be turned off at the package-level using [GeoUtils' global conf
 
 f, ax = plt.subplots(1, 3)
 ax[0].set_title("Raster 1")
-rast.plot(ax=ax[0])
+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",)
+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])
+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
@@ -156,23 +161,25 @@ vect_clipped = vect_reproj.crop(rast, clip=True)
 
 f, ax = plt.subplots(1, 2)
 ax[0].set_title("Before clipping")
-rast.plot(ax=ax[0])
+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])
+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`.
 
-Shifting **modifies the georeferencing of the data by a horizontal offset** without modifying the underlying data,
+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
-# Shift the raster by a certain offset
+# Translate the raster by a certain offset
 rast_shift = rast.translate(xoff=1000, yoff=1000)
 ```
 
@@ -183,12 +190,14 @@ rast_shift = rast.translate(xoff=1000, yoff=1000)
 :  code_prompt_hide: "Hide the code for plotting the figure"
 
 f, ax = plt.subplots(1, 2)
-ax[0].set_title("Before shifting")
-rast.plot(ax=ax[0])
+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 shifting")
-rast_shift.plot(ax=ax[1])
+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

doc/source/index.md

@@ -35,9 +35,9 @@ GeoUtils is built on top of core geospatial packages (Rasterio, GeoPandas) and n
 (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**,
+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 **to be computationally scalable** through Dask. 
+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!
 
@@ -127,6 +127,7 @@ analysis_examples/index
 
 api
 cli
+config
 background
 ```
 

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}