Re-organize documentation to list features independently of data objects

doc/source/api.md

@@ -47,6 +47,7 @@ documentation.
     Raster.crs
     Raster.transform
     Raster.nodata
+    Raster.area_or_point
 ```
 
 ### Derived attributes
@@ -120,6 +121,7 @@ documentation.
     Raster.load
     Raster.save
     Raster.to_pointcloud
+    Raster.from_regular_pointcloud
     Raster.to_rio_dataset
     Raster.to_xarray
 ```

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")
@@ -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": [
     #         {

doc/source/rasters_index.md → doc/source/data_object_index.md

@@ -1,12 +1,14 @@
-(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`.
 
 ```{toctree}
 :maxdepth: 2
 
 raster_class
+vector_class
+pointcloud_class
 mask_class
 satimg_class
 ```

doc/source/distance_ops.md

@@ -0,0 +1,76 @@
+---
+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)
+proximity.plot(cmap="viridis")
+```
+
+## Buffering without overlap
+
+Buffering consists in expanding 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 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)
+```

doc/source/georeferencing.md

@@ -0,0 +1,214 @@
+---
+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 and vector info
+rast.info()
+vect.info()
+```
+
+### Coordinate reference systems
+
+[Coordinate reference systems (CRSs)](https://en.wikipedia.org/wiki/Spatial_reference_system), sometimes also called a 
+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/).
+
+```{code-cell} ipython3
+# Show CRS attribute of a raster and vector
+print(rast.crs)
+print(vect.crs)
+```
+
+```{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 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)
+```
+
+```{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)
+# 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`.
+
+```{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)**, used often for 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 can induce a half-pixel shift during raster–point interfacing** 
+(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**.
+
+This behaviour can be modified at the package-level by using [GeoUtils' configuration parameters]() 
+`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
+print(rast.get_metric_crs())
+```
+
+### Re-set georeferencing metadata
+
+To add?
+
+{func}`geoutils.Vector.set_crs`
+{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.
+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)
+```