From a1f33100562998f741f27d9fe8661e00728db153 Mon Sep 17 00:00:00 2001 From: Romain Hugonnet Date: Thu, 5 Mar 2026 21:41:35 -0900 Subject: [PATCH 1/3] First draft of scalability documentation --- doc/source/_static/css/custom.css | 43 ++++ doc/source/cheatsheet_osgeo.md | 230 ++++++++++++++++++ doc/source/conf.py | 2 +- doc/source/core_index.md | 1 + doc/source/ecosystem.md | 49 ++++ doc/source/index.md | 15 +- .../{georeferencing.md => referencing.md} | 2 +- doc/source/scalability_concept.md | 58 +++++ doc/source/scalability_index.md | 13 + doc/source/scalability_logic.md | 4 + doc/source/scalability_support.md | 174 +++++++++++++ doc/source/scalability_usage.md | 148 +++++++++++ doc/source/summary.md | 216 ++++++++++++++++ ...otransformations.md => transformations.md} | 2 +- geoutils/raster/base.py | 5 +- geoutils/raster/xr_accessor.py | 4 +- 16 files changed, 958 insertions(+), 8 deletions(-) create mode 100644 doc/source/cheatsheet_osgeo.md create mode 100644 doc/source/ecosystem.md rename doc/source/{georeferencing.md => referencing.md} (99%) create mode 100644 doc/source/scalability_concept.md create mode 100644 doc/source/scalability_index.md create mode 100644 doc/source/scalability_logic.md create mode 100644 doc/source/scalability_support.md create mode 100644 doc/source/scalability_usage.md create mode 100644 doc/source/summary.md rename doc/source/{geotransformations.md => transformations.md} (99%) diff --git a/doc/source/_static/css/custom.css b/doc/source/_static/css/custom.css index 43a5b3e2d..947e6ee7f 100644 --- a/doc/source/_static/css/custom.css +++ b/doc/source/_static/css/custom.css @@ -6,3 +6,46 @@ div.cell details.hide > summary { div.cell details[open].above-input div.cell_input { border-top: None; } + +/* ------------------------------------------- + GeoUtils tables: subsection rows + ------------------------------------------- */ + +/* Put the band color on the row itself */ +table.tight-table tr:has(.gu-table-section) { + background-color: var(--pst-color-surface-secondary) !important; +} + +/* Make all cells transparent so the row color shows through */ +table.tight-table tr:has(.gu-table-section) > td { + background: transparent !important; + background-image: none !important; + box-shadow: none !important; + + /* slightly shorter row */ + padding: 0.45rem 0.6rem !important; + + /* remove vertical borders so cells look merged */ + border-left: 0 !important; + border-right: 0 !important; + + /* top + bottom edge */ + border-top: 1px solid var(--pst-color-border) !important; + border-bottom: 1px solid var(--pst-color-border) !important; +} + +/* left accent bar */ +table.tight-table tr:has(.gu-table-section) > td:first-child { + border-left: 4px solid var(--pst-color-primary) !important; +} + +/* hide placeholder cells but keep layout */ +table.tight-table tr:has(.gu-table-section) > td:not(:first-child) { + visibility: hidden; +} + +/* label styling */ +.gu-table-section { + font-weight: 700; + letter-spacing: 0.02em; +} \ No newline at end of file diff --git a/doc/source/cheatsheet_osgeo.md b/doc/source/cheatsheet_osgeo.md new file mode 100644 index 000000000..53e2363c6 --- /dev/null +++ b/doc/source/cheatsheet_osgeo.md @@ -0,0 +1,230 @@ +(cheatsheet-osgeo)= +# Cheatsheet: From GDAL/OGR + +This page helps users familiar with **GDAL/OGR** migrate their operations to the GeoUtils API. + +Regarding function names, GeoUtils exposes **an API almost entirely consistent with the recently overhauled GDAL CLI**. + +Note that GeoUtils is **object-oriented** (methods run on {class}`~geoutils.Raster`, {class}`~geoutils.Vector`, {class}`~geoutils.PointCloud`, or on {class}`~xarray.DataArray` and {class}`~geopandas.GeoDataFrame` through `rst`, `vct` and `pc` accessors), while GDAL/OGR utilities are typically **file-oriented** (read from disk, write to disk). + +We also provide a conversion table for operations specific to DEMs (e.g. slope, aspect, roughness indexes) that are supported through our sister-package [xDEM](https://xdem.readthedocs.io/en/stable/). + +## GDAL/OGR utilities + +```{list-table} GDAL/OGR ⟶ GeoUtils +:header-rows: 1 +:widths: 3 4 4 2 +:align: left +:class: tight-table + +* - Operation + - Old GDAL/OGR + - New GDAL CLI + - GeoUtils + +* - Metadata + - + - + - + +* - Footprint + - `gdalinfo`/`ogrinfo` + - `gdal raster/vector footprint` + - {attr}`~geoutils.Raster.footprint` + +* - Bounding box + - `gdalinfo`/`ogrinfo` + - `gdal raster/vector bbox` + - {attr}`~geoutils.Raster.bounds` + +* - Info summary + - `gdalinfo`/`ogrinfo` + - `gdal raster/vector info` + - {attr}`~geoutils.Raster.info` + +* - Raster ⟶ Raster + - + - + - + +* - Reproject/warp + - `gdalwarp` + - `gdal raster reproject` + - {meth}`~geoutils.Raster.reproject` + +* - Crop/clip + - `gdal_translate` / `gdalwarp` + - `gdal raster clip` + - {meth}`~geoutils.Raster.crop` / {meth}`~geoutils.Raster.icrop` + +* - Edit referencing + - `gdal_edit` / `gdalmove.py` + - `gdal raster edit` + - {meth}`~geoutils.Raster.set_crs`, {meth}`~geoutils.Raster.set_transform`, {meth}`~geoutils.Raster.set_nodata`, {meth}`~geoutils.Raster.translate` + +* - Convert file format + - `gdal_translate` + - `gdal raster convert` + - {meth}`~geoutils.Raster.to_file` + +* - Filter + - — + - `gdal raster neighbors` + - {meth}`~geoutils.Raster.filter` + +* - Proximity distance + - `gdal_proximity` + - `gdal raster proximity` + - {meth}`~geoutils.Raster.proximity` + +* - Raster calculator + - `gdal_calc.py` + - `gdal raster calc` + - NumPy array interface on {attr}`~geoutils.Raster.data`) + +* - Mosaic / merge rasters + - `gdal_merge.py` + - `gdalbuildvrt`+`gdal raster reproject` + - {func}`~geoutils.raster.merge_rasters` + +* - Stack rasters into multiband raster + - `gdal_stack.py` + - `gdal raster stack` + - {func}`~geoutils.raster.stack_rasters` + +* - Fill nodata gaps + - `gdal_fillnodata.py` + - `gdal raster fill-nodata` + - Not implemented (planned) + +* - Remove small raster regions + - `gdal_sieve.py` + - `gdal raster sieve` + - Not implemented (planned) + +* - Generate contours + - `gdal_contour` + - `gdal raster contour` + - Not implemented + +* - Raster ⟶ Point + - + - + - + +* - Interpolate at coordinates + - `gdallocationinfo` + - `gdal raster pixel-info` + - {meth}`~geoutils.Raster.to_pointcloud`/{meth}`~geoutils.Raster.interp_points` + +* - Raster ⟶ Vector + - + - + - + +* - Polygonize + - `gdal_polygonize.py` + - `gdal raster polygonize` + - {meth}`~geoutils.Raster.polygonize` + +* - Vector ⟶ Vector + - + - + - + +* - Reproject + - `ogr2ogr` + - `gdal vector reproject` + - {meth}`~geoutils.Vector.reproject` + +* - Crop/clip + - `ogr2ogr -clipsrc` + - `gdal vector clip` + - {meth}`~geoutils.Vector.crop` + +* - Translate + - `ogr2ogr` (SQL transform) + - — + - {meth}`~geoutils.Vector.translate` + +* - Copy + - `ogr2ogr` + - `gdal vector convert` + - {meth}`~geoutils.Vector.copy` + +* - Geometric operations + - `ogr2ogr` (specific) + - `gdal vector simplify/buffer/...` + - {meth}`~geoutils.Vector.simplify` / {meth}`~geoutils.Vector.buffer` / ... + +* - Vector ⟶ Raster + - + - + - + +* - Rasterize + - `gdal_rasterize` + - `gdal vector rasterize` + - {meth}`~geoutils.Vector.rasterize` + +* - Point ⟶ Raster + - + - + - + +* - Grid points + - `gdal_grid` + - `gdal vector grid` + - {meth}`~geoutils.PointCloud.grid` +``` + +## GDAL DEM utilities + +The table below maps common **GDAL DEM analysis utilities** (historically provided by `gdaldem`) to their equivalents in **xDEM**. + +```{list-table} +:header-rows: 1 +:widths: 2 3 3 2 +:align: left +:class: tight-table + +* - Operation + - Old GDAL + - New GDAL CLI + - xDEM + +* - Terrain attributes + - + - + - + +* - Slope + - `gdaldem slope` + - `gdal raster slope` + - {func}`~xdem.DEM.slope` + +* - Aspect + - `gdaldem aspect` + - `gdal raster aspect` + - {func}`~xdem.DEM.aspect` + +* - Hillshade + - `gdaldem hillshade` + - `gdal raster hillshade` + - {func}`~xdem.DEM.hillshade` + +* - Terrain Ruggedness Index + - `gdaldem TRI` + - `gdal raster TRI` + - {func}`~xdem.DEM.terrain_ruggedness_index` + +* - Topographic Position Index + - `gdaldem TPI` + - `gdal raster TPI` + - {func}`~xdem.DEM.topographic_position_index` + +* - Roughness + - `gdaldem roughness` + - `gdal raster roughness` + - {func}`~xdem.DEM.roughness` +``` \ No newline at end of file diff --git a/doc/source/conf.py b/doc/source/conf.py index 0c8bb9112..152953dfe 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -15,7 +15,7 @@ from sphinx_gallery.sorting import ExampleTitleSortKey, ExplicitOrder project = "GeoUtils" -copyright = "2025, GeoUtils Developers" +copyright = "2026, GeoUtils Developers" author = "GeoUtils Developers" diff --git a/doc/source/core_index.md b/doc/source/core_index.md index adc444533..888ed5b80 100644 --- a/doc/source/core_index.md +++ b/doc/source/core_index.md @@ -13,4 +13,5 @@ core_array_funcs core_lazy_load core_parsing core_inheritance +core_scalability ``` diff --git a/doc/source/ecosystem.md b/doc/source/ecosystem.md new file mode 100644 index 000000000..896c5b43e --- /dev/null +++ b/doc/source/ecosystem.md @@ -0,0 +1,49 @@ +(ecosystem)= +# Ecosystem + +GeoUtils integrates naturally with the broader **geospatial ecosystem**. +It extends commonly used tools and works alongside many other for geospatial data access, processing, and analysis. + +See the {ref}`accessors` page for details on GeoUtils' accessors. + +```{see-also} +**[xDEM](https://github.com/GlacioHack/xdem)** is the sister package of GeoUtils, focused on the **analysis of digital elevation models (DEMs) and elevation point clouds**, including terrain attributes, coregistration and uncertainty propagation. +``` + +## Related Python libraries + +Several Python libraries provide complementary functionality: + +- **[EOReader](https://github.com/sertit/eoreader)** — Unified access to satellite imagery products +- **[stackstac](https://github.com/gjoseph92/stackstac)** — Load STAC datasets as large raster stacks + +Some libraries focus on raster operations specifically: + +- **[Xarray-Spatial](https://github.com/makepath/xarray-spatial)** — Spatial analysis on Xarray arrays +- **[ODC-Geo](https://github.com/opendatacube/odc-geo)** — Geospatial extensions for Xarray +- **[GeoWombat](https://github.com/jgrss/geowombat)** — Earth observation processing workflows + +## Tools in other languages + +Many geospatial workflows also rely on tools in other languages. + +**R** + +- **[terra](https://rspatial.org/terra/)** — raster and vector spatial data processing +- **[sf](https://r-spatial.github.io/sf/)** — vector geospatial data + +**Julia** + +- **[ArchGDAL.jl](https://github.com/yeesian/ArchGDAL.jl)** — GDAL bindings for Julia +- **[Rasters.jl](https://github.com/rafaqz/Rasters.jl)** — raster data processing + +## Learning resources and data portals + +Several platforms provide geospatial datasets, cloud environments and learning resources: + +- **[STAC](https://stacspec.org)** — SpatioTemporal Asset Catalog standard for geospatial data discovery +- **[NASA EarthData](https://earthdata.nasa.gov/)** — NASA Earth observation datasets +- **[USGS EarthExplorer](https://earthexplorer.usgs.gov/)** — satellite and elevation data +- **[Copernicus Data Space](https://dataspace.copernicus.eu/)** — Sentinel satellite data +- **[Radiant Earth MLHub](https://mlhub.earth/)** — open Earth observation datasets +- **[CryoCloud](https://www.cryocloud.io/)** — open cloud platform for geospatial data analysis \ No newline at end of file diff --git a/doc/source/index.md b/doc/source/index.md index ce98077a5..63ae82f0c 100644 --- a/doc/source/index.md +++ b/doc/source/index.md @@ -98,6 +98,7 @@ about_geoutils how_to_install quick_start feature_overview +summary ``` ```{toctree} @@ -105,15 +106,25 @@ feature_overview :maxdepth: 2 core_index +scalability_index data_object_index -georeferencing -geotransformations +referencing +transformations raster_vector_point distance_ops stats filters ``` +```{toctree} +:caption: Resources +:maxdepth: 2 + +cheatsheet_osgeo +ecosystem +``` + + ```{toctree} :caption: Examples :maxdepth: 2 diff --git a/doc/source/georeferencing.md b/doc/source/referencing.md similarity index 99% rename from doc/source/georeferencing.md rename to doc/source/referencing.md index 88841b86a..9636c1bab 100644 --- a/doc/source/georeferencing.md +++ b/doc/source/referencing.md @@ -10,7 +10,7 @@ kernelspec: language: python name: geoutils --- -(georeferencing)= +(referencing)= # Referencing Below, a summary of the **georeferencing attributes** of geospatial data objects and the **methods to manipulate these diff --git a/doc/source/scalability_concept.md b/doc/source/scalability_concept.md new file mode 100644 index 000000000..bd0d34a45 --- /dev/null +++ b/doc/source/scalability_concept.md @@ -0,0 +1,58 @@ +(scalability-concept)= +# Concept definitions + +This section describes scalability concepts important to grasp to manipulate our objects. + +Scalable execution relies on three complementary mechanisms: + +- **Deferred I/O with implicit loading:** Operations that update data-related metadata without loading the underlying array or geometries, +- **Chunked execution:** Operations that process data tile-by-tile to limit memory usage, +- **Lazy execution:** Operations whose computation is deferred until explicitly requested. + +These mechanisms are often combined (e.g., Dask operations are always chunked **and** lazy) but are conceptually independent. + +Finally, one should note that the above concepts only apply to operations that interact with the underlying **data arrays or geometries** of GeoUtils objects. +Naturally, all **metadata operations** (e.g., accessing {attr}`~geoutils.Raster.crs`, {attr}`~geoutils.Raster.bounds`, or {meth}`~geoutils.Raster.info`) have no effect on the array, and do not trigger any loading. + +## Deferred I/O and implicit loading + +**Deferred input/output** refers to operations that modify only **internal I/O metadata**, avoid reading the data entirely and postponing loading. + +Typical examples include {meth}`~geoutils.Raster.crop`, {meth}`~geoutils.Raster.copy`, and {meth}`~geoutils.Raster.translate`, +which behave similarly as Xarray's {meth}`~xarray.DataArray.sel`, {meth}`~xarray.DataArray.copy`, or {meth}`~xarray.DataArray.assign_coords`. + +When using the Xarray `rst` accessor, this behavior follows the **native Xarray deferred I/O model**. The {class}`~geoutils.Raster` class implements the +same behavior so that both APIs have consistent semantics. + +This behaviour pairs intrinsically with **implicit loading:** When an object is opened, only metadata is loaded. +Accessing {attr}`~geoutils.Raster.data`, or calling operations that require the array, will **implicitly load the data into memory**. + +An important aspect of **deferred I/O** is that it works with both **in-memory** (NumPy) and **scalable backends** (Dask), allowing +to extract parts of large files without any chunked or lazy considerations. + +## Chunked execution + +**Chunked execution** refers to processing raster data **tile-by-tile** instead of loading the full array into memory. + +This enables **out-of-core execution**, allowing datasets larger than available RAM to be processed safely. + +In GeoUtils, chunked execution is implemented through two backends: + +- **Dask**, used through the Xarray `rst` accessor, +- **Multiprocessing**, used through the {class}`~geoutils.Raster` object. + +Both backends read and process raster chunks sequentially, keeping peak memory usage proportional to the chunk size rather than the full dataset size. +Chunked execution therefore allows GeoUtils to scale to large datasets while maintaining a **predictable memory footprint**. +For a list of expected memory usage per operation, see the {ref}`scalability-support` page. + +## Lazy execution + +Lazy execution refers to **deferring computation until results are explicitly requested**. + +In GeoUtils, lazy execution is available through the Xarray `rst` accessor with **Dask-backed arrays**. + +Operations build a **Dask computation graph** instead of executing immediately. The computation is triggered only when required, for example when calling +`compute()` or when writing results to disk. It is particularly useful when **chaining multiple raster operations**, because intermediate results do not need to be materialized or +written/read from disk (which costs extra I/O time, often much longer than compute time). + +Lazy execution always relies on **chunked execution**, but the reverse is not true: chunked processing can also run eagerly, as in the Multiprocessing backend. \ No newline at end of file diff --git a/doc/source/scalability_index.md b/doc/source/scalability_index.md new file mode 100644 index 000000000..585251018 --- /dev/null +++ b/doc/source/scalability_index.md @@ -0,0 +1,13 @@ +(scalability-index)= +# Scalability + +The following sections present how to use our scalability features with both **Dask** and **Multiprocessing**, with a full description of **supported operations** and the **logic behind chunked implementations** of each functionality. + +```{toctree} +:maxdepth: 2 + +scalability_usage +scalability_concept +scalability_support +scalability_logic +``` diff --git a/doc/source/scalability_logic.md b/doc/source/scalability_logic.md new file mode 100644 index 000000000..635051415 --- /dev/null +++ b/doc/source/scalability_logic.md @@ -0,0 +1,4 @@ +(scalability-logic)= +# Implementation strategies + +TODO LAST \ No newline at end of file diff --git a/doc/source/scalability_support.md b/doc/source/scalability_support.md new file mode 100644 index 000000000..aa5626fb6 --- /dev/null +++ b/doc/source/scalability_support.md @@ -0,0 +1,174 @@ +(scalability-support)= +# Supported operations + +GeoUtils supports **scalable execution for most of its raster methods**, including nearly all **raster–point** and **raster–vector** interface operations. Support for **point-cloud methods** is partially supported and under development, while **vector methods** may also gain scalable support in the future (lower priority). + +Chunked implementations can run through either **Dask** (via the Xarray and Pandas accessors) or **Multiprocessing** (via GeoUtils objects such as {class}`~geoutils.Raster`). + +Both object types (accessors or GeoUtils) expose the **exact same API**, and both chunked backends use the same internal logic, and all methods are tested to **yield identical output** as in-memory. + +## Table summary + +The table below summarizes the **scalability support** of GeoUtils operations with respect to their input and output behavior. + +If you are unfamiliar with **chunked and lazy execution** or **deferred I/O**, see the {ref}`scalability-concept` page. + +**Legend:** +- {bdg-success}`Chunked` — Processes in chunks, without loading (for input) and/or returning (for output) the full data. This is **also lazy (deferred execution) when using Dask**, but not when using Multiprocessing. +- {bdg-secondary}`In-memory` — Loads (for input) or returns (for output) full data in-memory. +- {bdg-primary}`Deferred I/O` — Deferred input/output by updating internal metadata (as Xarray's {meth}`~xarray.DataArray.isel`). + +```{list-table} +:name: Scalability summary +:widths: 1 1 1 2 +:header-rows: 1 +:align: center +:class: tight-table + +* - Method + - Input + - Output + - Memory usage (# chunks) + +* - Raster ⟶ Raster + - + - + - + +* - {meth}`~geoutils.Raster.reproject` + - {bdg-success}`Chunked` + - {bdg-success}`Chunked` + - ~4 (default), or ~downsampling² +* - {meth}`~geoutils.Raster.crop` / {meth}`~geoutils.Raster.icrop` + - {bdg-primary}`Deferred I/O` + - {bdg-primary}`Deferred I/O` + - 0 +* - {meth}`~geoutils.Raster.translate` + - {bdg-primary}`Deferred I/O` + - {bdg-primary}`Deferred I/O` + - 0 +* - {meth}`~geoutils.Raster.copy` + - {bdg-primary}`Deferred I/O` + - {bdg-primary}`Deferred I/O` + - 0 +* - {meth}`~geoutils.Raster.filter` + - {bdg-success}`Chunked` + - {bdg-success}`Chunked` + - ~2–3 (if small filter window) +* - {meth}`~geoutils.Raster.proximity` + - {bdg-secondary}`In-memory` + - {bdg-secondary}`In-memory` + - — + +* - Raster ⟶ Point + - + - + - + +* - {meth}`~geoutils.Raster.subsample` + - {bdg-success}`Chunked` + - {bdg-secondary}`In-memory` + - ~1 +* - {meth}`~geoutils.Raster.interp_points` + - {bdg-success}`Chunked` + - {bdg-secondary}`In-memory` + - ~1 +* - {meth}`~geoutils.Raster.reduce_points` + - {bdg-secondary}`In-memory` + - {bdg-secondary}`In-memory` + - — + +* - Raster ⟶ Vector + - + - + - + +* - {meth}`~geoutils.Raster.polygonize` + - {bdg-success}`Chunked` + - {bdg-secondary}`In-memory` + - ~1–2 + +* - Raster ⟶ Other + - + - + - + +* - {meth}`~geoutils.Raster.plot` + - {bdg-secondary}`In-memory` + - {bdg-secondary}`In-memory` + - — +* - {meth}`~geoutils.Raster.get_stats` + - {bdg-secondary}`In-memory` + - {bdg-secondary}`In-memory` + - — + +* - Point ⟶ Point + - + - + - + +* - {meth}`~geoutils.PointCloud.reproject` + - {bdg-secondary}`In-memory` + - {bdg-secondary}`In-memory` + - — +* - {meth}`~geoutils.PointCloud.translate` + - {bdg-secondary}`In-memory` + - {bdg-secondary}`In-memory` + - — +* - {meth}`~geoutils.PointCloud.crop` + - {bdg-secondary}`In-memory` + - {bdg-secondary}`In-memory` + - — + +* - Point ⟶ Raster + - + - + - + +* - {meth}`~geoutils.PointCloud.grid` + - {bdg-secondary}`In-memory` + - {bdg-secondary}`In-memory` + - — +* - {meth}`~geoutils.Raster.from_pointcloud_regular` + - {bdg-secondary}`In-memory` + - {bdg-secondary}`In-memory` + - — + +* - Vector ⟶ Raster + - + - + - + +* - {meth}`~geoutils.Vector.rasterize` + - {bdg-secondary}`In-memory` + - {bdg-success}`Chunked` + - ~1 +* - {meth}`~geoutils.Vector.create_mask` + - {bdg-secondary}`In-memory` + - {bdg-success}`Chunked` + - ~1 + +* - Vector ⟶ Point + - + - + - + +* - {meth}`~geoutils.Vector.create_mask` + - {bdg-secondary}`In-memory` + - {bdg-secondary}`In-memory` + - — + +* - Point ⟶ Other + - + - + - + +* - {meth}`~geoutils.PointCloud.get_stats` + - {bdg-secondary}`In-memory` + - {bdg-secondary}`In-memory` + - — +``` + +Note that nearly all **raster inputs/outputs** methods support {bdg-success}`Chunked`, while **point and vector inputs/outputs** are currently {bdg-secondary}`In-memory`, as often less limiting. + +For more insights into chunked implementation strategies and behaviour expected for each operation, see the {ref}`scalability-logic` page. \ No newline at end of file diff --git a/doc/source/scalability_usage.md b/doc/source/scalability_usage.md new file mode 100644 index 000000000..9f786ce25 --- /dev/null +++ b/doc/source/scalability_usage.md @@ -0,0 +1,148 @@ +--- +file_format: mystnb +jupytext: + formats: md:myst + text_representation: + extension: .md + format_name: myst +kernelspec: + display_name: geoutils-env + language: python + name: geoutils +--- +(scalability-usage)= +# Usage and good practices + +GeoUtils supports scalable execution for most of its **raster** and (soon) **point cloud** operations (**vector** support may be added in the future, but is usually less limiting). + +It relies on two execution backends: + +- **Dask**, through its `rst` Xarray accessor and `pc` Pandas accessor (**lazy** and **chunked** execution), +- **Multiprocessing**, through its {class}`~geoutils.Raster` and {class}`~geoutils.PointCloud` objects (**chunked** execution only) . + +Both backends mirror the **exact same object operations and chunked logic**, and yield **identical** results as in-memory operations. +**Lazy** refers to deferred execution using Dask, while **chunked** refers to processing raster data tile-by-tile to limit memory usage. +For details on scalability concepts, see the {ref}`scalability-concept` page. + +As a rule of thumb: + +- Use **Dask** to work on **Xarray and GeoPandas objects** through our accessors `rst` and `pc`, and if you want to chain several operations lazily. +- Use **Multiprocessing** to work with our {class}`~geoutils.Raster` and {class}`~geoutils.PointCloud` objects, and if you are fine with intermediate writing/reading between steps. +- Use standard **in-memory execution** to work efficiently on small rasters, which is possible even if those were loaded from larger rasters (use {class}`~geoutils.Raster.crop`). + +## Using Dask through accessors + +With Dask, raster operations are both **chunked** and **lazy**. +This behavior is enabled by opening a raster with the `chunks` argument, which returns an Xarray object backed by Dask arrays. + +```{code-cell} python +import geoutils as gu + +filename_rast = gu.examples.get_path("exploradores_aster_dem") + +ds = gu.open_raster(filename_rast, chunks={"x": 200, "y": 200}) + +ds +``` + +GeoUtils, through the `rst` accessor, automatically detects the **Dask** input and switches to a chunked implementation. + +```{code-cell} python +# Change output resolution +out_res = (ds.rst.res[0] * 2, ds.rst.res[1] / 2) + +# Reproject lazily and out-of-memory +ds_reproj = ds.rst.reproject( + res=out_res, + resampling="bilinear", +) + +ds_reproj +``` + +The resulting raster remains **lazy**. Computation only happens when explicitly requested with `compute()`. + +For a raster output, one typically wants to write to file lazily to avoid loading it in-memory: + +```{code-cell} python +# ds_reproj.rst.to_file("reproj_rast.tif", compute=True) +``` + +This triggers computation and writes the raster **chunk-by-chunk** to the output file. + +Or, the output can be chained with another operation. For example, we can extract a small random subsample of the reprojected raster: + +```{code-cell} python +# Subsample lazily and out-of-memory +sub_ds = ds_reproj.rst.subsample( + subsample=5000, +) + +sub_ds +``` + +The output array is again lazy, and in this case we can use `compute()` to return the in-memory NumPy array: +```{code-cell} python +sub_ds.compute() +``` + +Lazy execution is particularly useful to **chain several operations** without ever loading the full raster into memory or writing any file to disk. + +## Using Multiprocessing through GeoUtils objects + +Chunked execution can also be enabled on the {class}`~geoutils.Raster` object using Multiprocessing. +Multiprocessing performs **chunked execution** (only loads chunks in-memory), but is not **lazy** (it runs immediately). + +```{code-cell} python +rast = gu.Raster(filename_rast) + +rast +``` + +By passing a {class}`~geoutils.multiproc.MultiprocConfig` configuration to an operation, the out-of-memory behaviour is triggered. +The configuration requires a `chunk_size`, and can optionally define the output file and cluster to use (defaults to temporary instances for both). + +```{code-cell} python +from geoutils.multiproc import MultiprocConfig + +# Optional: specify output file (otherwise a temporary file is created) +mp_config = MultiprocConfig(chunk_size=200, outfile="reproj_rast.tif") + +# Reproject out-of-memory, reading and writing chunk-by-chunk +rast_reproj_mp = rast.reproject( + res=out_res, + resampling="bilinear", + mp_config=mp_config, +) + +rast_reproj_mp +``` + +If the output is a {class}`~geoutils.Raster`, it is written to disk out-of-memory, and the returned object is a {class}`~geoutils.Raster` of that file without data loaded. +This keeps syntax consistent with in-memory code, and allow to easily chain operations. + +For other output types, the Multiprocessing backends will load the result in-memory. + +```{code-cell} python +# Subsample out-of-memory and return loaded array +samp_rast_mp = rast_reproj_mp.subsample( + subsample=5000, +) + +samp_rast_mp +``` + +This backend is convenient when working directly with {class}`~geoutils.Raster` objects and performing **step-by-step processing**. + +## Good practices with chunked and lazy operations + +- If **memory** is the limitating factor for you, use a **single-threaded scheduler** through Dask (```dask.config.set(scheduler='single-threaded')```) or Multiprocessing (default cluster), +- If **speed** is the limiting factor for you, use **parallelized processes** through Dask (see [Dask scheduler configuration](https://docs.dask.org/en/stable/scheduler-overview.html#scheduler-overview)) or Multiprocessing (see our Cluster configuration), +- Choose chunk sizes large enough to reduce scheduling overhead, but **small enough to fit comfortably in memory**, +- Keep chunk sizes **consistent across operations** to avoid unnecessary rechunking, +- Insert **breakpoints** (for example by writing intermediate results to disk) to prevent building overly large Dask graphs. + +For more guidance on chunk sizing and performance, see the [Dask array best practices](https://docs.dask.org/en/stable/array-best-practices.html). + +Finally, note that currently, operations returning **point** or **vector** outputs are often **eager** and scalable execution applies mostly to the **raster input/output**. +The full description of supported methods is available on the {ref}`scalability-support` page. \ No newline at end of file diff --git a/doc/source/summary.md b/doc/source/summary.md new file mode 100644 index 000000000..921efecda --- /dev/null +++ b/doc/source/summary.md @@ -0,0 +1,216 @@ +(method-summary)= + +# Feature overview + +GeoUtils provides a unified API for manipulating **raster**, **vector**, and **point-cloud** data, with **scalable execution** for most raster operations. + +The **tables below** summarize the core operations of GeoUtils, their scalability and backends. + +If you are interested in converting from GDAL/OGR, see our {ref}`cheatsheet-osgeo` page. + +## Summary of methods and scalability + +Methods of GeoUtils are shared across object types and expose a **consistent API** for clarity (similarly as the recent [GDAL CLI overhaul](https://gdal.org/en/stable/programs/index.html)). They also support convenient inputs such as **match-reference arguments** (e.g., matching a grid for reprojection or rasterization, matching bounds for cropping, matching point coordinates for interpolation). See the {ref}`core-match-ref` page for details. + +Nearly all **raster operations** support **scalable execution** through [Dask](https://www.dask.org/) or Multiprocessing, allowing large datasets to be processed chunk-by-chunk without loading the full array into memory. +While the table below provide a scalability summary, details on exact **supported operations** relative to inputs/outputs are available on the {ref}`scalability-support` section. + +Some operations also support multiple computational **backends** (for example SciPy or Numba implementations for numerical routines). + +All methods are tested to ensure they produce **identical results** whether executed in-memory, using chunked processing, or through alternative computational backends. + +## Data operations + +We first describe GeoUtils' core **data operations**, which operate on underlying arrays or geometries and therefore benefit from **scalable execution**. + +**Legend:** **“/”** indicates methods **shared across object types**, while **“⟷”** indicates methods **interfacing between two object types**. + +```{list-table} Common API for data operations +:widths: 3 5 1 2 +:header-rows: 1 +:align: left +:class: tight-table + +* - Method + - Notes + - Scalable + - Backend + +* - Raster / Vector / Point + - + - + - + +* - {meth}`~geoutils.Raster.reproject()` + - Reproject to other CRS. Default tolerance parameters ensure chunk-invariance. + - ✅ + - Rasterio / PyProj + +* - {meth}`~geoutils.Raster.crop()` + - Crop to bounds. For vectors, can return geometries intersecting (untouched) or clipped. + - ✅ + - Rasterio / GeoPandas + +* - {meth}`~geoutils.Raster.translate()` + - Apply a grid shift to object. + - ✅ + - NumPy / GeoPandas + +* - {meth}`~geoutils.Raster.proximity()` + - Estimate proximity distance to target values or geometries. + - ❌ + - SciPy + +* - {meth}`~geoutils.Raster.plot()` + - Visualization helper. + - ❌ + - Matplotlib + +* - Raster / Point + - + - + - + +* - {meth}`~geoutils.Vector.create_mask()` + - Create boolean mask of a vector geometries over raster or point. + - ✅ + - Rasterio / GeoPandas + +* - {meth}`~geoutils.Raster.get_stats()` + - Compute statistics of valid values over a valid mask. + - ❌ + - NumPy / SciPy + +* - {meth}`~geoutils.Raster.subsample()` + - Randomly sample valid values. Chunk-invariant seed ensures reproducibility. + - ✅ + - NumPy + +* - {meth}`~geoutils.Raster.filter()` + - Filter over window. Fast vectorized logic with NaN support. + - ✅ + - SciPy + +* - Raster ⟷ Vector + - + - + - + +* - {meth}`~geoutils.Raster.polygonize()` + - Convert raster regions to vector polygons. Multiple chunked strategies for performance. + - ✅ + - Rasterio / GeoPandas + +* - {meth}`~geoutils.Vector.rasterize()` + - Burn vector geometries onto a raster grid. + - ✅ + - Rasterio + +* - Raster ⟷ Point + - + - + - + +* - {meth}`~geoutils.Raster.interp_points()` + - Interpolate raster at point locations. Fast regular-grid logic with added NaN propagation. + - ✅ + - SciPy + +* - {meth}`~geoutils.Raster.reduce_points()` + - Aggregate raster values around points. + - ❌ + - NumPy + +* - {meth}`~geoutils.PointCloud.grid()` + - Grid irregular points onto a raster grid. Multiple approaches with added NaN propagation. + - ❌ + - SciPy + +* - {meth}`~geoutils.Raster.from_pointcloud_regular()` + - Direct conversion when points lie on a regular grid. + - ❌ + - NumPy + +* - {meth}`~geoutils.Raster.to_pointcloud()` + - Conversion to point cloud. + - ❌ + - NumPy +``` + +## Metadata properties and operations + +In addition to data operations, GeoUtils exposes **metadata** properties and methods consistently across geospatial objects. +These operate only on metadata and therefore **do not load or modify underlying data arrays**. + +```{list-table} Common API from metadata operations +:widths: 3 7 +:header-rows: 1 +:align: left +:class: tight-table + +* - Attribute / Method + - Description + +* - Raster / Vector / Point + - + +* - {attr}`~geoutils.Raster.crs` + - Coordinate reference system (CRS) of object. + +* - {attr}`~geoutils.Raster.bounds` + - Bounding box of object. + +* - {attr}`~geoutils.Raster.footprint` + - Footprint polygon geometry of object. + +* - {attr}`~geoutils.Raster.is_loaded` + - Whether geospatial object is loaded in-memory. + +* - {attr}`~geoutils.Raster.name` + - Filename of object on disk, if it exists. + +* - {meth}`~geoutils.Raster.get_bounds_projected()` + - Bounds projected in other CRS. + +* - {meth}`~geoutils.Raster.get_footprint_projected()` + - Footprint polygon geometry in other CRS. + +* - {meth}`~geoutils.Raster.get_metric_crs()` + - Get metric CRS suitable for this object. + +* - {meth}`~geoutils.Raster.info()` + - Summary of attributes for geospatial object. + +* - Raster / Point + - + +* - {attr}`~geoutils.Raster.data` + - Data array (2D grid for raster, 1D for point cloud). + +* - {attr}`~geoutils.Raster.shape` + - Shape of data array. + +* - {attr}`~geoutils.Raster.is_mask` + - Whether object is a mask. Clarifies ambiguity of raster/point file types often not supporting boolean types. + +* - Raster + - + +* - {attr}`~geoutils.Raster.transform` + - Geotransform to map raster indices to spatial coordinates. + +* - {attr}`~geoutils.Raster.nodata` + - Nodata value used to represent missing data on disk. + +* - {attr}`~geoutils.Raster.area_or_point` + - Pixel interpretation of raster values, either center point or area average. + +* - Point + - + +* - {attr}`~geoutils.PointCloud.point_count` + - Number of points in the point cloud. +``` + + + diff --git a/doc/source/geotransformations.md b/doc/source/transformations.md similarity index 99% rename from doc/source/geotransformations.md rename to doc/source/transformations.md index 8127f3a77..294c87c1b 100644 --- a/doc/source/geotransformations.md +++ b/doc/source/transformations.md @@ -10,7 +10,7 @@ kernelspec: language: python name: geoutils --- -(geotransformations)= +(transformations)= # Transformations In GeoUtils, **for all geospatial data objects, georeferenced transformations are exposed through the same functions** diff --git a/geoutils/raster/base.py b/geoutils/raster/base.py index 688f0c3bf..b8fc059f3 100644 --- a/geoutils/raster/base.py +++ b/geoutils/raster/base.py @@ -1265,7 +1265,10 @@ def reproject( # If return copy is True (target georeferenced grid was the same as input) if return_copy: - return self + if self._is_xr: + return self._obj + else: + return self # To make MyPy happy without overload for _reproject (as it might re-structured soon anyway) # assert data is not None diff --git a/geoutils/raster/xr_accessor.py b/geoutils/raster/xr_accessor.py index 138c15c5a..243a4f08d 100644 --- a/geoutils/raster/xr_accessor.py +++ b/geoutils/raster/xr_accessor.py @@ -306,10 +306,10 @@ def to_geoutils(self) -> RasterBase: area_or_point=self.area_or_point, ) - def to_file(self, **kwargs: Any) -> None: + def to_file(self, *args: Any, **kwargs: Any) -> None: """ Write raster to file. Wrapper around rioxarray.to_raster(). """ - self._obj.rio.to_raster(**kwargs) + self._obj.rio.to_raster(*args, **kwargs) From dd1701f6db4d0593df3b11d9167e36886e50542b Mon Sep 17 00:00:00 2001 From: Romain Hugonnet Date: Mon, 9 Mar 2026 21:33:34 -0800 Subject: [PATCH 2/3] Incremental commit on doc --- doc/source/_templates/raster_method.rst | 8 + doc/source/api.md | 248 +++-- doc/source/api_raster.md | 17 + doc/source/api_rst.md | 18 + doc/source/cheatsheet_osgeo.md | 2 +- .../code/diagram_chunked_interp_points.py | 746 +++++++++++++++ doc/source/code/diagram_chunked_polygonize.py | 789 ++++++++++++++++ doc/source/code/diagram_chunked_rasterize.py | 692 ++++++++++++++ doc/source/code/diagram_chunked_reproject.py | 850 ++++++++++++++++++ doc/source/code/diagram_chunked_subsample.py | 749 +++++++++++++++ doc/source/conf.py | 1 + doc/source/core_index.md | 2 - doc/source/core_lazy_load.md | 86 -- doc/source/ecosystem.md | 2 +- doc/source/feature_overview.md | 589 +++++++----- doc/source/index.md | 13 +- doc/source/pointcloud_class.md | 2 +- doc/source/raster_class.md | 121 +-- doc/source/release_notes.md | 2 +- doc/source/scalability_concept.md | 109 ++- doc/source/scalability_logic.md | 125 ++- doc/source/scalability_usage.md | 23 +- doc/source/summary.md | 216 ----- doc/source/vector_class.md | 2 +- 24 files changed, 4692 insertions(+), 720 deletions(-) create mode 100644 doc/source/_templates/raster_method.rst create mode 100644 doc/source/api_raster.md create mode 100644 doc/source/api_rst.md create mode 100644 doc/source/code/diagram_chunked_interp_points.py create mode 100644 doc/source/code/diagram_chunked_polygonize.py create mode 100644 doc/source/code/diagram_chunked_rasterize.py create mode 100644 doc/source/code/diagram_chunked_reproject.py create mode 100644 doc/source/code/diagram_chunked_subsample.py delete mode 100644 doc/source/core_lazy_load.md delete mode 100644 doc/source/summary.md diff --git a/doc/source/_templates/raster_method.rst b/doc/source/_templates/raster_method.rst new file mode 100644 index 000000000..4ee5f4b67 --- /dev/null +++ b/doc/source/_templates/raster_method.rst @@ -0,0 +1,8 @@ +{% set name = objname.split('.')[-1] %} + +Raster.{{ name }}{% if objtype == "method" %}(){% endif %} or ds.rst.{{ name }}{% if objtype == "method" %}(){% endif %} +======================================================================= + +.. currentmodule:: geoutils + +.. auto{{ objtype }}:: {{ fullname }} diff --git a/doc/source/api.md b/doc/source/api.md index 8558cd0f8..3de6ffa15 100644 --- a/doc/source/api.md +++ b/doc/source/api.md @@ -2,226 +2,218 @@ # API reference This page provides a summary of GeoUtils’ API. -For more details and examples, refer to the relevant chapters in the main part of the -documentation. +For more details and examples, refer to the relevant chapters in the main part of the documentation. ```{eval-rst} .. currentmodule:: geoutils ``` -## Raster -```{eval-rst} -.. minigallery:: geoutils.Raster - :add-heading: -``` + +```{toctree} +:maxdepth: 1 +:hidden: -### Opening a file - -```{eval-rst} -.. autosummary:: - :toctree: gen_modules/ - - Raster - Raster.info +Raster +RasterAccessor ``` -### Create from an array +(raster-api)= +## Raster API -```{eval-rst} -.. autosummary:: - :toctree: gen_modules/ +GeoUtils exposes the raster API through two mirrored interfaces: +- {class}`~geoutils.Raster`, an interface operating directly on a GeoUtils object, +- A {class}`rst ` accessor extending {class}`xarray.DataArray` objects as rasters. - Raster.from_array -``` -(api-raster-attrs)= - -### Main attributes +Both expose the **same methods and attributes**. -```{eval-rst} -.. autosummary:: - :toctree: gen_modules/ +Only **file opening** and **scalable execution** differ between the two interfaces: +- **File opening:** {class}`~geoutils.Raster` objects are opened by instantiating the class, whereas {meth}`~geoutils.open_raster` is used for an {class}`xarray.DataArray` object, +- **Scalable execution:** the {class}`rst ` accessor supports **Dask**, while the {class}`~geoutils.Raster` supports **Multiprocessing** instead. - Raster.data - Raster.crs - Raster.transform - Raster.nodata - Raster.area_or_point -``` +### Opening a raster file -### Derived attributes +Use {meth}`~geoutils.open_raster` for an {class}`xarray.DataArray`, or instantiate for a {class}`~geoutils.Raster`. ```{eval-rst} .. autosummary:: :toctree: gen_modules/ - Raster.shape - Raster.height - Raster.width - Raster.count - Raster.bands - Raster.res - Raster.bounds - Raster.dtype + open_raster + Raster.__init__ ``` -### Other attributes +### Create raster from an array ```{eval-rst} .. autosummary:: :toctree: gen_modules/ - - Raster.is_mask - Raster.is_loaded - Raster.name - Raster.driver - Raster.tags + :template: raster_method.rst + + ~raster.base.RasterBase.from_array ``` -(api-geo-handle)= +(api-raster-attrs)= -### Geospatial handling methods +### Main attributes ```{eval-rst} .. autosummary:: :toctree: gen_modules/ - - Raster.crop - Raster.icrop - Raster.reproject - Raster.polygonize - Raster.proximity - Raster.interp_points - Raster.reduce_points - Raster.filter + :template: raster_method.rst + + ~raster.base.RasterBase.data + ~raster.base.RasterBase.crs + ~raster.base.RasterBase.transform + ~raster.base.RasterBase.nodata + ~raster.base.RasterBase.area_or_point ``` -### Plotting +### Derived attributes ```{eval-rst} .. autosummary:: :toctree: gen_modules/ + :template: raster_method.rst - Raster.plot + ~raster.base.RasterBase.shape + ~raster.base.RasterBase.height + ~raster.base.RasterBase.width + ~raster.base.RasterBase.count + ~raster.base.RasterBase.bands + ~raster.base.RasterBase.res + ~raster.base.RasterBase.bounds + ~raster.base.RasterBase.footprint + ~raster.base.RasterBase.dtype ``` -### Get statistics +### Other attributes ```{eval-rst} .. autosummary:: :toctree: gen_modules/ + :template: raster_method.rst - Raster.get_stats + ~raster.base.RasterBase.is_mask + ~raster.base.RasterBase.is_loaded + ~raster.base.RasterBase.name + ~raster.base.RasterBase.driver + ~raster.base.RasterBase.tags ``` -### Get or update data methods +(api-geo-handle)= + +### Geospatial operations ```{eval-rst} .. autosummary:: :toctree: gen_modules/ + :template: raster_method.rst - Raster.copy - Raster.astype - Raster.set_mask - Raster.set_nodata - Raster.get_nanarray - Raster.get_mask - Raster.subsample + ~raster.base.RasterBase.crop + ~raster.base.RasterBase.icrop + ~raster.base.RasterBase.reproject + ~raster.base.RasterBase.polygonize + ~raster.base.RasterBase.proximity + ~raster.base.RasterBase.interp_points + ~raster.base.RasterBase.reduce_points + ~raster.base.RasterBase.filter ``` -### I/O methods +### Plotting ```{eval-rst} .. autosummary:: :toctree: gen_modules/ + :template: raster_method.rst - Raster.load - Raster.to_file - Raster.to_pointcloud - Raster.from_pointcloud_regular - Raster.to_rio_dataset - Raster.to_xarray + ~raster.base.RasterBase.plot ``` -### Coordinate and extent methods +### Statistics ```{eval-rst} .. autosummary:: :toctree: gen_modules/ + :template: raster_method.rst - Raster.xy2ij - Raster.ij2xy - Raster.coords - Raster.translate - Raster.outside_image + ~raster.base.RasterBase.get_stats ``` -### Projection methods +### Data manipulation ```{eval-rst} .. autosummary:: :toctree: gen_modules/ + :template: raster_method.rst - Raster.get_metric_crs - Raster.get_bounds_projected - Raster.get_footprint_projected - Raster.intersection + ~raster.base.RasterBase.copy + ~raster.base.RasterBase.astype + ~raster.base.RasterBase.set_mask + ~raster.base.RasterBase.set_nodata + ~raster.base.RasterBase.get_nanarray + ~raster.base.RasterBase.get_mask + ~raster.base.RasterBase.subsample ``` -### Testing methods +### Loading, writing and converting ```{eval-rst} .. autosummary:: :toctree: gen_modules/ + :template: raster_method.rst - Raster.raster_equal - Raster.georeferenced_grid_equal + ~raster.base.RasterBase.load + ~raster.base.RasterBase.to_file + ~raster.base.RasterBase.to_pointcloud + ~raster.base.RasterBase.from_pointcloud_regular + ~raster.base.RasterBase.to_rio_dataset + ~raster.base.RasterBase.to_xarray ``` -### Arithmetic with other rasters, arrays or numbers +### Georeferencing utilities ```{eval-rst} .. autosummary:: :toctree: gen_modules/ + :template: raster_method.rst - Raster.__add__ - Raster.__sub__ - Raster.__neg__ - Raster.__mul__ - Raster.__truediv__ - Raster.__floordiv__ - Raster.__mod__ - Raster.__pow__ + ~raster.base.RasterBase.xy2ij + ~raster.base.RasterBase.ij2xy + ~raster.base.RasterBase.coords + ~raster.base.RasterBase.translate + ~raster.base.RasterBase.outside_image ``` -And reverse operations. - -### Logical operators casting to mask (boolean raster) +### Projection utilities ```{eval-rst} .. autosummary:: :toctree: gen_modules/ + :template: raster_method.rst - Raster.__eq__ - Raster.__ne__ - Raster.__lt__ - Raster.__le__ - Raster.__gt__ - Raster.__ge__ + ~raster.base.RasterBase.get_metric_crs + ~raster.base.RasterBase.get_bounds_projected + ~raster.base.RasterBase.get_footprint_projected + ~raster.base.RasterBase.intersection ``` -### Array interface with NumPy +### Testing utilities ```{eval-rst} .. autosummary:: :toctree: gen_modules/ + :template: raster_method.rst - Raster.__array_ufunc__ - Raster.__array_function__ + ~raster.base.RasterBase.raster_equal + ~raster.base.RasterBase.georeferenced_grid_equal ``` -## Multiple rasters +### Multiple rasters ```{eval-rst} .. autosummary:: @@ -232,21 +224,6 @@ And reverse operations. raster.merge_rasters ``` -[//]: # (## Multiprocessing) - -[//]: # () -[//]: # (```{eval-rst}) - -[//]: # (.. autosummary::) - -[//]: # ( :toctree: gen_modules/) - -[//]: # () -[//]: # ( raster.MultiprocConfig) - -[//]: # () -[//]: # (```) - ## Vector ```{eval-rst} @@ -626,3 +603,14 @@ documentation](https://shapely.readthedocs.io/en/stable/properties.html). PointCloud.pointcloud_equal PointCloud.georeferenced_coords_equal ``` + +## Multiprocessing configuration + +To performed **chunked execution** on GeoUtils objects, pass this Multiprocessing configuration to function that support it. + +```{eval-rst} +.. autosummary:: + :toctree: gen_modules/ + + multiproc.MultiprocConfig +``` \ No newline at end of file diff --git a/doc/source/api_raster.md b/doc/source/api_raster.md new file mode 100644 index 000000000..0947b3af5 --- /dev/null +++ b/doc/source/api_raster.md @@ -0,0 +1,17 @@ +(api-raster)= +# Raster + +{class}`~geoutils.Raster` is the GeoUtils object for rasters. +Its API is mirrored in the {class}`rst `, which exposes the same methods and attributes on {class}`xarray.DataArray` objects. + +**For the main raster API reference page, see {ref}`raster-api`.** +Below is the full class description specific to {class}`~geoutils.Raster`: + +```{eval-rst} +.. currentmodule:: geoutils + +.. autoclass:: Raster + :members: + :inherited-members: + :undoc-members: false + :show-inheritance: \ No newline at end of file diff --git a/doc/source/api_rst.md b/doc/source/api_rst.md new file mode 100644 index 000000000..595446126 --- /dev/null +++ b/doc/source/api_rst.md @@ -0,0 +1,18 @@ +(api-rst)= +# The `rst` accessor + +The {class}`rst ` accessor is the interface to GeoUtils raster methods on {class}`xarray.DataArray` objects. + +It mirrors the API of {class}`~geoutils.Raster`, which exposes the same methods and attributes on GeoUtils raster objects. + +**For the main raster API reference page, see {ref}`raster-api`.** +Below is the full class description specific to {class}`rst `: + +```{eval-rst} +.. currentmodule:: geoutils + +.. autoclass:: RasterAccessor + :members: + :inherited-members: + :undoc-members: false + :show-inheritance: \ No newline at end of file diff --git a/doc/source/cheatsheet_osgeo.md b/doc/source/cheatsheet_osgeo.md index 53e2363c6..af2b7b3a7 100644 --- a/doc/source/cheatsheet_osgeo.md +++ b/doc/source/cheatsheet_osgeo.md @@ -5,7 +5,7 @@ This page helps users familiar with **GDAL/OGR** migrate their operations to the Regarding function names, GeoUtils exposes **an API almost entirely consistent with the recently overhauled GDAL CLI**. -Note that GeoUtils is **object-oriented** (methods run on {class}`~geoutils.Raster`, {class}`~geoutils.Vector`, {class}`~geoutils.PointCloud`, or on {class}`~xarray.DataArray` and {class}`~geopandas.GeoDataFrame` through `rst`, `vct` and `pc` accessors), while GDAL/OGR utilities are typically **file-oriented** (read from disk, write to disk). +Note that GeoUtils is **object-oriented** (methods run on {class}`~geoutils.Raster`, {class}`~geoutils.Vector`, {class}`~geoutils.PointCloud`, or on {class}`~xarray.DataArray` and {class}`~geopandas.GeoDataFrame` through {class}`rst `, `vct` and `pc` accessors), while GDAL/OGR utilities are typically **file-oriented** (read from disk, write to disk). We also provide a conversion table for operations specific to DEMs (e.g. slope, aspect, roughness indexes) that are supported through our sister-package [xDEM](https://xdem.readthedocs.io/en/stable/). diff --git a/doc/source/code/diagram_chunked_interp_points.py b/doc/source/code/diagram_chunked_interp_points.py new file mode 100644 index 000000000..b483670ef --- /dev/null +++ b/doc/source/code/diagram_chunked_interp_points.py @@ -0,0 +1,746 @@ +"""Script to generate a diagram for chunked raster-to-point interpolation.""" +from __future__ import annotations + +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.collections import LineCollection +from matplotlib.lines import Line2D +from matplotlib.patches import FancyArrowPatch, Rectangle + + +# ----------------------------------------------------------------------------- +# Example data +# ----------------------------------------------------------------------------- + +# Raster grid +NROWS = 8 +NCOLS = 10 +CHUNK_ROWS = (4, 4) +CHUNK_COLS = (5, 5) + +# Interpolation method and implied overlap depth +METHOD = "cubic" +DEPTH = 4 + +SELECTED_CHUNK_LW = 2.8 +EXAMPLE_CHUNK_LW = 4.2 + +# Input point coordinates in raster data coordinates +POINT_X = np.array([0.9, 1.8, 3.7, 4.4, 1.4, 6.2, 8.4, 8.7, 2.6, 6.9]) +POINT_Y = np.array([6.7, 5.2, 6.0, 2.2, 1.4, 5.6, 2.3, 1.3, 1.8, 2.7]) + +# Conceptual chunking of the input point ribbon +POINT_CHUNKS = (4, 3, 3) +HIGHLIGHT_POINT_CHUNK = 1 + +point_chunk_starts = np.concatenate([[0], np.cumsum(POINT_CHUNKS)]) +POINT_IDS_HIGHLIGHT = np.arange( + point_chunk_starts[HIGHLIGHT_POINT_CHUNK], + point_chunk_starts[HIGHLIGHT_POINT_CHUNK + 1], +) + + +def point_to_raster_chunk(x: float, y: float) -> tuple[int, int]: + """Assign point to raster chunk based on containing chunk.""" + col_edges = np.concatenate([[0], np.cumsum(CHUNK_COLS)]) + row_edges_top = np.concatenate([[0], np.cumsum(CHUNK_ROWS)]) + + ix = np.searchsorted(col_edges[1:], x, side="right") + row_from_top = np.searchsorted(row_edges_top[1:], NROWS - y, side="right") + iy = row_from_top + return iy, ix + + +POINT_CHUNK_LOC = [point_to_raster_chunk(float(x), float(y)) for x, y in zip(POINT_X, POINT_Y, strict=True)] +HIGHLIGHTED_RASTER_CHUNKS = sorted({POINT_CHUNK_LOC[i] for i in POINT_IDS_HIGHLIGHT}) +highlighted_raster_chunk = POINT_CHUNK_LOC[POINT_IDS_HIGHLIGHT[0]] + +# Example interpolated values for the output ribbon +OUTPUT_VALUES = np.array([12.4, 15.0, 14.1, 9.8, 11.3, 10.7, 13.2, 8.9, 7.4, 10.1]) + + +# ----------------------------------------------------------------------------- +# Styling +# ----------------------------------------------------------------------------- + +NEUTRAL = "#333333" +PIXEL_GRID = "0.80" +CHUNK_COLOR = "#222222" +HIGHLIGHT_COLOR = "#F58518" +POINT_COLOR = "#6BAED6" +POINT_EDGE = "#4C78A8" +POINT_FADED = "#BFD7EA" +OVERLAP_FILL_ALPHA = 0.18 + + +# ----------------------------------------------------------------------------- +# Geometry helpers +# ----------------------------------------------------------------------------- + +def raster_chunk_bounds(chunk_loc: tuple[int, int]) -> tuple[float, float, float, float]: + """Return chunk bounds as (left, bottom, right, top).""" + iy, ix = chunk_loc + + row_starts = np.concatenate([[0], np.cumsum(CHUNK_ROWS)]) + col_starts = np.concatenate([[0], np.cumsum(CHUNK_COLS)]) + + top_row = row_starts[iy] + bottom_row = row_starts[iy + 1] + left_col = col_starts[ix] + right_col = col_starts[ix + 1] + + left = float(left_col) + right = float(right_col) + top = float(NROWS - top_row) + bottom = float(NROWS - bottom_row) + + return left, bottom, right, top + + +def expanded_chunk_bounds( + bounds: tuple[float, float, float, float], + *, + depth_px: float, +) -> tuple[float, float, float, float]: + """Expand chunk bounds by interpolation overlap depth.""" + left, bottom, right, top = bounds + return left - depth_px, bottom - depth_px, right + depth_px, top + depth_px + + +# ----------------------------------------------------------------------------- +# Plot helpers +# ----------------------------------------------------------------------------- + +def _setup_axis(ax: plt.Axes, *, xlim: tuple[float, float], ylim: tuple[float, float]) -> None: + """Common axis formatting.""" + ax.set_xlim(*xlim) + ax.set_ylim(*ylim) + ax.set_aspect("equal") + ax.set_xticks([]) + ax.set_yticks([]) + for spine in ax.spines.values(): + spine.set_visible(False) + + +def _draw_pixel_grid(ax: plt.Axes, *, nrows: int, ncols: int, extend: float = 1.0) -> None: + """Draw thin raster pixel grid, slightly extended around the raster.""" + segments: list[np.ndarray] = [] + + for x in range(ncols + 1): + segments.append(np.array([[x, -extend], [x, nrows + extend]], dtype=float)) + + for y in range(nrows + 1): + segments.append(np.array([[-extend, y], [ncols + extend, y]], dtype=float)) + + coll = LineCollection( + segments, + colors=PIXEL_GRID, + linewidths=0.8, + capstyle="round", + joinstyle="round", + zorder=1, + clip_on=False, + ) + ax.add_collection(coll) + + +def _draw_chunk_boundaries(ax: plt.Axes, *, extend: float = 1.0) -> None: + """Draw all raster chunk boundaries.""" + row_starts = np.concatenate([[0], np.cumsum(CHUNK_ROWS)]) + col_starts = np.concatenate([[0], np.cumsum(CHUNK_COLS)]) + + for x in col_starts: + ax.plot( + [x, x], + [-extend, NROWS + extend], + color=CHUNK_COLOR, + linewidth=2.3, + zorder=3, + solid_capstyle="round", + clip_on=False, + ) + + for y_idx in row_starts: + y = NROWS - y_idx + ax.plot( + [-extend, NCOLS + extend], + [y, y], + color=CHUNK_COLOR, + linewidth=2.3, + zorder=3, + solid_capstyle="round", + clip_on=False, + ) + + +def _draw_points_on_raster( + ax: plt.Axes, + *, + highlight_ids: np.ndarray | None = None, + show_ids: bool = False, +) -> None: + """Draw points in raster coordinates.""" + highlight_set = set([] if highlight_ids is None else highlight_ids.tolist()) + + for i, (x, y) in enumerate(zip(POINT_X, POINT_Y, strict=True)): + highlighted = i in highlight_set + ax.scatter( + x, + y, + s=42 if highlighted else 34, + facecolor=POINT_COLOR if highlighted else POINT_FADED, + edgecolor=POINT_EDGE, + linewidth=1.0, + alpha=0.95 if highlighted else 0.75, + zorder=6 if highlighted else 4, + ) + if show_ids: + ax.text( + x + 0.12, + y + 0.12, + f"{i}", + fontsize=8.8, + color=NEUTRAL, + ha="left", + va="bottom", + zorder=7, + ) + + +def _draw_highlight_chunk( + ax: plt.Axes, + bounds: tuple[float, float, float, float], + *, + color: str, + lw: float = 2.8, +) -> None: + """Draw highlighted chunk outline.""" + left, bottom, right, top = bounds + rect = Rectangle( + (left, bottom), + right - left, + top - bottom, + fill=False, + edgecolor=color, + linewidth=lw, + zorder=7, + ) + ax.add_patch(rect) + + +def _draw_overlap_chunk(ax: plt.Axes, bounds: tuple[float, float, float, float]) -> None: + """Draw overlap-expanded source chunk.""" + left, bottom, right, top = bounds + fill = Rectangle( + (left, bottom), + right - left, + top - bottom, + facecolor=HIGHLIGHT_COLOR, + edgecolor="none", + alpha=OVERLAP_FILL_ALPHA, + zorder=2, + clip_on=False, + ) + edge = Rectangle( + (left, bottom), + right - left, + top - bottom, + fill=False, + edgecolor=HIGHLIGHT_COLOR, + linewidth=2.6, + zorder=7, + clip_on=False, + ) + ax.add_patch(fill) + ax.add_patch(edge) + + +def _draw_point_ribbon( + ax: plt.Axes, + *, + highlight_chunk: int, + values: np.ndarray | None = None, + title: str | None = None, +) -> None: + """Draw a 1D point ribbon with chunk separators.""" + total = len(POINT_X) + y_center = 0.0 + box_h = 0.9 + + starts = np.concatenate([[0], np.cumsum(POINT_CHUNKS)]) + + rect = Rectangle( + (0, y_center - 0.5 * box_h), + total, + box_h, + fill=False, + edgecolor=CHUNK_COLOR, + linewidth=2.2, + zorder=3, + ) + ax.add_patch(rect) + + for s in starts[1:-1]: + ax.plot([s, s], [y_center - 0.5 * box_h, y_center + 0.5 * box_h], color=CHUNK_COLOR, lw=2.2, zorder=3) + + x0 = starts[highlight_chunk] + x1 = starts[highlight_chunk + 1] + + fill = Rectangle( + (x0, y_center - 0.5 * box_h), + x1 - x0, + box_h, + facecolor=HIGHLIGHT_COLOR, + edgecolor="none", + alpha=0.12, + zorder=1, + ) + edge = Rectangle( + (x0, y_center - 0.5 * box_h), + x1 - x0, + box_h, + fill=False, + edgecolor=HIGHLIGHT_COLOR, + linewidth=2.6, + zorder=4, + ) + ax.add_patch(fill) + ax.add_patch(edge) + + for i in range(total): + x = i + 0.5 + highlighted = x0 <= i < x1 + ax.scatter( + x, + y_center, + s=42 if highlighted else 34, + facecolor=POINT_COLOR if highlighted else POINT_FADED, + edgecolor=POINT_EDGE, + linewidth=1.0, + alpha=0.95 if highlighted else 0.75, + zorder=5, + ) + + label = f"{i}" if values is None else f"{values[i]:.1f}" + ax.text( + x, + y_center - 0.62, + label, + ha="center", + va="top", + fontsize=8.5, + color=NEUTRAL, + zorder=6, + ) + + if title is not None: + ax.text( + x0 + 0.5 * (x1 - x0), + y_center + 0.72, + title, + ha="center", + va="bottom", + fontsize=10, + color=HIGHLIGHT_COLOR, + fontweight="bold", + ) + + +def _add_workflow_arrow( + fig: plt.Figure, + start: tuple[float, float], + end: tuple[float, float], + *, + rad: float = 0.0, +) -> None: + """Add a figure-level workflow arrow.""" + arrow = FancyArrowPatch( + start, + end, + transform=fig.transFigure, + arrowstyle="-|>", + connectionstyle=f"arc3,rad={rad}", + linewidth=1.8, + color=NEUTRAL, + mutation_scale=18, + zorder=30, + ) + fig.add_artist(arrow) + +def _add_visual_legend(fig: plt.Figure) -> None: + """Draw a two-row horizontal legend with balanced spacing.""" + y1 = 0.05 + y2 = 0.018 + + centers = [0.1, 0.35, 0.6, 0.8] + + sym_w = 0.032 + txt_gap = 0.012 + + # ---------------- Row 1 ---------------- + # Pixel grid + cx = centers[0] + x = cx - 0.075 + fig.add_artist( + Line2D( + [x, x + sym_w], + [y1, y1], + transform=fig.transFigure, + color=PIXEL_GRID, + lw=0.8, + solid_capstyle="round", + ) + ) + fig.text(x + sym_w + txt_gap, y1, "Raster pixel grid", transform=fig.transFigure, va="center", fontsize=10) + + # Chunk boundary + cx = centers[1] + x = cx - 0.085 + fig.add_artist( + Line2D( + [x, x + sym_w], + [y1, y1], + transform=fig.transFigure, + color=CHUNK_COLOR, + lw=2.3, + solid_capstyle="round", + ) + ) + fig.text(x + sym_w + txt_gap, y1, "Raster chunk boundary", transform=fig.transFigure, va="center", fontsize=10) + + # All input points (light blue) + cx = centers[2] + x = cx - 0.080 + fig.add_artist( + Line2D( + [x + sym_w / 2], + [y1], + transform=fig.transFigure, + marker="o", + markersize=7, + markerfacecolor=POINT_FADED, + markeredgecolor=POINT_EDGE, + linestyle="None", + ) + ) + fig.text(x + sym_w + txt_gap, y1, "All input points", transform=fig.transFigure, va="center", fontsize=10) + + # Selected point chunk (darker blue) + cx = centers[3] + x = cx - 0.090 + fig.add_artist( + Line2D( + [x + sym_w / 2], + [y1], + transform=fig.transFigure, + marker="o", + markersize=7, + markerfacecolor=POINT_COLOR, + markeredgecolor=POINT_EDGE, + linestyle="None", + ) + ) + fig.text(x + sym_w + txt_gap, y1, "Selected point chunk", transform=fig.transFigure, va="center", fontsize=10) + + # ---------------- Row 2 ---------------- + # Selected raster chunks + cx = centers[1] + x = cx - 0.25 + rect1 = Rectangle( + (x, y2 - 0.010), + sym_w, + 0.020, + transform=fig.transFigure, + fill=False, + edgecolor=HIGHLIGHT_COLOR, + linewidth=SELECTED_CHUNK_LW, + ) + fig.add_artist(rect1) + fig.text(x + sym_w + txt_gap, y2, "Raster chunks used by this point chunk", transform=fig.transFigure, + va="center", fontsize=10) + + # Example raster chunk / overlap + cx = centers[3] + x = cx - 0.25 + rect_fill = Rectangle( + (x, y2 - 0.010), + sym_w, + 0.020, + transform=fig.transFigure, + facecolor=HIGHLIGHT_COLOR, + edgecolor="none", + alpha=0.15, + ) + rect_edge = Rectangle( + (x, y2 - 0.010), + sym_w, + 0.020, + transform=fig.transFigure, + fill=False, + edgecolor=HIGHLIGHT_COLOR, + linewidth=EXAMPLE_CHUNK_LW, + ) + fig.add_artist(rect_fill) + fig.add_artist(rect_edge) + fig.text(x + sym_w + txt_gap, y2, "Example chunk processed with overlap", transform=fig.transFigure, va="center", fontsize=10) + + +def _add_left_column_separator(fig: plt.Figure, ax_top_left: plt.Axes, ax_bottom_left: plt.Axes) -> None: + """Add a subtle horizontal separator in the left column.""" + pos_t = ax_top_left.get_position() + pos_b = ax_bottom_left.get_position() + + x0 = pos_t.x0 - 0.2 * (pos_t.x1 - pos_t.x0) + x1 = pos_t.x1 + 0.2 * (pos_t.x1 - pos_t.x0) + y = 0.52 * (pos_t.y0 + pos_b.y1) + + line = Line2D( + [x0, x1], + [y, y], + transform=fig.transFigure, + color="0.75", + lw=2.5, + alpha=0.8, + solid_capstyle="round", + zorder=50, + ) + fig.add_artist(line) + +def _draw_highlight_chunks( + ax: plt.Axes, + bounds_list: list[tuple[float, float, float, float]], + *, + color: str, + example_bounds: tuple[float, float, float, float] | None = None, + lw_selected: float = SELECTED_CHUNK_LW, + lw_example: float = EXAMPLE_CHUNK_LW, +) -> None: + """Draw selected raster chunks, with one optional example chunk thicker.""" + for bounds in bounds_list: + left, bottom, right, top = bounds + + is_example = example_bounds is not None and np.allclose(bounds, example_bounds) + lw = lw_example if is_example else lw_selected + + rect = Rectangle( + (left, bottom), + right - left, + top - bottom, + fill=False, + edgecolor=color, + linewidth=lw, + zorder=7, + joinstyle="round", + ) + ax.add_patch(rect) + +# ----------------------------------------------------------------------------- +# Main figure +# ----------------------------------------------------------------------------- + +def make_chunked_interp_points_diagram() -> tuple[plt.Figure, np.ndarray]: + """Build a 4-panel schematic for chunked raster-to-point interpolation.""" + fig = plt.figure(figsize=(8, 6)) + gs = fig.add_gridspec(2, 2, hspace=0.5, wspace=0.30) + + ax_ul = fig.add_subplot(gs[0, 0]) # Input points + ax_ur = fig.add_subplot(gs[0, 1]) # Assign to raster chunks + ax_bl = fig.add_subplot(gs[1, 0]) # Reordered outputs + ax_br = fig.add_subplot(gs[1, 1]) # Interpolate on overlap-expanded chunks + + axes = np.array([ax_ul, ax_ur, ax_bl, ax_br], dtype=object) + + chunk_bounds_list = [raster_chunk_bounds(ch) for ch in HIGHLIGHTED_RASTER_CHUNKS] + chunk_bounds = max(chunk_bounds_list, key=lambda b: (b[3], b[ + 2])) # Use upper-right selected chunk + overlap_bounds = expanded_chunk_bounds(chunk_bounds, depth_px=1.0) + + # ------------------------------------------------------------------------- + # Upper left: input points + # ------------------------------------------------------------------------- + _draw_point_ribbon(ax_ul, highlight_chunk=HIGHLIGHT_POINT_CHUNK, title="Point chunk") + ax_ul.text( + 0.5, + 1.02, + "Input points", + transform=ax_ul.transAxes, + ha="center", + va="bottom", + fontsize=13, + fontweight="bold", + color=NEUTRAL, + ) + ax_ul.text( + 0.5, + 0, + "Points are chunked along the 1D input sequence", + transform=ax_ul.transAxes, + ha="center", + va="top", + fontsize=10.2, + color=NEUTRAL, + ) + _setup_axis(ax_ul, xlim=(-0.4, len(POINT_X) + 0.4), ylim=(-1.0, 1.4)) + ax_ul.set_aspect("auto") + + # ------------------------------------------------------------------------- + # Upper right: assign to raster chunks + # ------------------------------------------------------------------------- + _draw_pixel_grid(ax_ur, nrows=NROWS, ncols=NCOLS, extend=1.0) + _draw_chunk_boundaries(ax_ur, extend=1.0) + _draw_points_on_raster(ax_ur, highlight_ids=POINT_IDS_HIGHLIGHT, show_ids=True) + _draw_highlight_chunks( + ax_ur, + chunk_bounds_list, + color=HIGHLIGHT_COLOR, + example_bounds=chunk_bounds, + ) + ax_ur.text( + 0.5, + 1.02, + "Assign points to raster chunks", + transform=ax_ur.transAxes, + ha="center", + va="bottom", + fontsize=13, + fontweight="bold", + color=NEUTRAL, + ) + ax_ur.text( + 0.5, + -0.05, + "Map point chunk coordinates to\nall intersecting raster chunks", + transform=ax_ur.transAxes, + ha="center", + va="top", + fontsize=10.2, + color=NEUTRAL, + linespacing=1.25, + ) + _setup_axis(ax_ur, xlim=(-1.2, NCOLS + 1.2), ylim=(-1.2, NROWS + 1.2)) + + # ------------------------------------------------------------------------- + # Bottom right: interpolate on overlap-expanded chunks + # ------------------------------------------------------------------------- + _draw_pixel_grid(ax_br, nrows=NROWS, ncols=NCOLS, extend=1.0) + _draw_chunk_boundaries(ax_br, extend=1.0) + _draw_overlap_chunk(ax_br, overlap_bounds) + _draw_highlight_chunk(ax_br, chunk_bounds, color=HIGHLIGHT_COLOR, lw=EXAMPLE_CHUNK_LW) + _draw_points_on_raster(ax_br, highlight_ids=POINT_IDS_HIGHLIGHT, show_ids=False) + + ax_br.text( + 0.5, + 1.02, + "Interpolate on overlapped chunks", + transform=ax_br.transAxes, + ha="center", + va="bottom", + fontsize=13, + fontweight="bold", + color=NEUTRAL, + ) + ax_br.text( + 0.5, + -0.05, + f"Loop on expanded raster chunks to interpolate\n" + f"(overlap size depends on resampling method)", + transform=ax_br.transAxes, + ha="center", + va="top", + fontsize=10.2, + color=NEUTRAL, + linespacing=1.25, + ) + _setup_axis(ax_br, xlim=(-1.2, NCOLS + 1.2), ylim=(-1.2, NROWS + 1.2)) + + # ------------------------------------------------------------------------- + # Bottom left: reordered outputs + # ------------------------------------------------------------------------- + _draw_point_ribbon( + ax_bl, + highlight_chunk=HIGHLIGHT_POINT_CHUNK, + values=OUTPUT_VALUES, + title="Reordered output values", + ) + ax_bl.text( + 0.5, + 1.02, + "Reorder outputs", + transform=ax_bl.transAxes, + ha="center", + va="bottom", + fontsize=13, + fontweight="bold", + color=NEUTRAL, + ) + ax_bl.text( + 0.5, + 0, + "Per-block interpolation results are concatenated and\nreordered to the original point sequence", + transform=ax_bl.transAxes, + ha="center", + va="top", + fontsize=10.2, + color=NEUTRAL, + linespacing=1.25, + ) + _setup_axis(ax_bl, xlim=(-0.4, len(POINT_X) + 0.4), ylim=(-1.0, 1.4)) + ax_bl.set_aspect("auto") + + # ------------------------------------------------------------------------- + # Global title + # ------------------------------------------------------------------------- + fig.text( + 0.5, + 0.99, + "Chunked raster interpolation at points", + ha="center", + va="top", + fontsize=15, + fontweight="semibold", + color="0.35", + ) + + # ------------------------------------------------------------------------- + # Workflow arrows + # ------------------------------------------------------------------------- + fig.canvas.draw() + + pos_ul = ax_ul.get_position() + pos_ur = ax_ur.get_position() + pos_br = ax_br.get_position() + pos_bl = ax_bl.get_position() + + # UL to UR + _add_workflow_arrow( + fig, + (pos_ul.x1 + 0.01, 0.5 * (pos_ul.y0 + pos_ul.y1)), + (pos_ur.x0 - 0.01, 0.5 * (pos_ur.y0 + pos_ur.y1)), + rad=0.0, + ) + + # UR to BR (curved outward to avoid the BR title) + _add_workflow_arrow( + fig, + (pos_ur.x1 - 0.25, 0.5 * (pos_ur.y0 + pos_ur.y1) - 0.02), + (pos_br.x1 - 0.25, 0.5 * (pos_br.y0 + pos_br.y1) + 0.02), + rad=0.45, + ) + + # BR to BL + _add_workflow_arrow( + fig, + (pos_br.x0 - 0.01, 0.5 * (pos_br.y0 + pos_br.y1)), + (pos_bl.x1 + 0.01, 0.5 * (pos_bl.y0 + pos_bl.y1)), + rad=0.0, + ) + + _add_left_column_separator(fig, ax_ul, ax_bl) + _add_visual_legend(fig) + + fig.subplots_adjust(left=0.05, right=0.98, top=0.90, bottom=0.16) + return fig, axes + + +fig, _ = make_chunked_interp_points_diagram() +plt.show() \ No newline at end of file diff --git a/doc/source/code/diagram_chunked_polygonize.py b/doc/source/code/diagram_chunked_polygonize.py new file mode 100644 index 000000000..3ce7583c4 --- /dev/null +++ b/doc/source/code/diagram_chunked_polygonize.py @@ -0,0 +1,789 @@ +"""Script to make diagram for chunked polygonize in documentation.""" +from __future__ import annotations + +from dataclasses import dataclass + +import geopandas as gpd +import matplotlib.pyplot as plt +import numpy as np +import rasterio as rio +from matplotlib.collections import LineCollection +from matplotlib.lines import Line2D +from matplotlib.patches import ConnectionPatch, FancyArrowPatch, Rectangle + +from geoutils.multiproc.chunked import ChunkedGeoGrid, GeoGrid +from geoutils.interface.vectorization import ( + _chunked_build_dst_geotiling, + _chunked_clip_gdf_to_bounds_polygonal, + _chunked_label_block_per_value, + _chunked_polygonize_block_labels, + _chunked_seam_pairs_from_strips, + _polygonize_base, +) + + +# ----------------------------------------------------------------------------- +# Example raster data and chunk layout +# ----------------------------------------------------------------------------- + +ARR = np.array( + [ + [1, 1, 1, 1, 1, 0], + [1, 1, 0, 0, 1, 0], + [0, 0, 2, 2, 2, 0], + [0, 3, 3, 3, 0, 0], + [0, 3, 0, 3, 0, 0], + ], + dtype=np.uint8, +) + +NROWS, NCOLS = ARR.shape +SPLIT_COL = 3 +HALO = 1 +CONNECTIVITY = 4 +FLOAT_TOL = 0.001 + +VALUE_COLORS: dict[int, str] = { + 0: "#FFFFFF", + 1: "#8ECAE6", + 2: "#B7E4C7", + 3: "#FFD166", +} + +NEUTRAL = "#333333" +PIXEL_GRID = "0.80" +CHUNK_COLOR = "#222222" +STITCH_COLOR = "#F58518" +POLYGON_COLOR = "#555555" + + +# ----------------------------------------------------------------------------- +# Small containers +# ----------------------------------------------------------------------------- + +@dataclass(frozen=True) +class BlockData: + """Per-block arrays and metadata used for the schematic.""" + block_id: int + block_id_dict: dict[str, int] + geogrid: GeoGrid + values: np.ndarray + mask: np.ndarray + labels: np.ndarray + gdf_labels: gpd.GeoDataFrame + gdf_geom: gpd.GeoDataFrame + + +# ----------------------------------------------------------------------------- +# Coordinate helpers +# ----------------------------------------------------------------------------- + +def cell_xy(i: int, j: int, nrows: int) -> tuple[float, float]: + """Bottom-left of a raster cell in plotting coordinates.""" + return float(j), float(nrows - 1 - i) + + +def cell_center(i: int, j: int, nrows: int) -> tuple[float, float]: + """Center of a raster cell in plotting coordinates.""" + x, y = cell_xy(i, j, nrows) + return x + 0.5, y + 0.5 + + +def centroid_from_cells(cells: list[tuple[int, int]], nrows: int) -> tuple[float, float]: + """Centroid from a list of global raster cells.""" + xs = [] + ys = [] + for i, j in cells: + xc, yc = cell_center(i, j, nrows) + xs.append(xc) + ys.append(yc) + return float(np.mean(xs)), float(np.mean(ys)) + + +# ----------------------------------------------------------------------------- +# Build block data using our package functions +# ----------------------------------------------------------------------------- + +def build_demo_blocks() -> tuple[ChunkedGeoGrid, list[GeoGrid], list[dict[str, int]], list[BlockData]]: + """Build a 2-chunk horizontal tiling and derive block products with package helpers.""" + transform = rio.transform.from_origin(0.0, float(NROWS), 1.0, 1.0) + crs = 4326 + chunks = ((NROWS,), (SPLIT_COL, NCOLS - SPLIT_COL)) + + tiling, block_geogrids, block_ids = _chunked_build_dst_geotiling( + shape=ARR.shape, + transform=transform, + crs=crs, + chunks=chunks, + ) + + blocks: list[BlockData] = [] + + for block_id, (b, gg) in enumerate(zip(block_ids, block_geogrids, strict=True)): + ys, ye, xs, xe = b["ys"], b["ye"], b["xs"], b["xe"] + values = ARR[ys:ye, xs:xe] + mask = values != 0 + + labels = _chunked_label_block_per_value( + values, + mask, + connectivity=CONNECTIVITY, + float_tol=FLOAT_TOL, + ) + + gdf_labels = _chunked_polygonize_block_labels( + labels=labels, + values=values, + mask=mask, + transform=gg.transform, + value_column="raster_value", + connectivity=CONNECTIVITY, + local_id_column="local_id", + float_tol=FLOAT_TOL, + ) + + ys_h = max(0, ys - HALO) + ye_h = min(NROWS, ye + HALO) + xs_h = max(0, xs - HALO) + xe_h = min(NCOLS, xe + HALO) + + values_h = ARR[ys_h:ye_h, xs_h:xe_h] + mask_h = values_h != 0 + transform_h = transform * rio.Affine.translation(xs_h, ys_h) + + gdf_geom = _polygonize_base( + values_h, + mask_h, + transform=transform_h, + crs=crs, + data_column_name="polygon_id", + value_column="raster_value", + connectivity=CONNECTIVITY, + float_tol=FLOAT_TOL, + ) + + bb_interior = gg.bounds + gdf_geom = _chunked_clip_gdf_to_bounds_polygonal( + gdf_geom, + bb_interior, + keep_border=True, + area_eps=0.0, + ) + + blocks.append( + BlockData( + block_id=block_id, + block_id_dict=b, + geogrid=gg, + values=values, + mask=mask, + labels=labels, + gdf_labels=gdf_labels, + gdf_geom=gdf_geom, + ) + ) + + return tiling, block_geogrids, block_ids, blocks + + +def label_cells_global(labels: np.ndarray, block: BlockData, local_id: int) -> list[tuple[int, int]]: + """Global raster cells belonging to one local component label.""" + ys0 = block.block_id_dict["ys"] + xs0 = block.block_id_dict["xs"] + ii, jj = np.where(labels == local_id) + return [(ys0 + int(i), xs0 + int(j)) for i, j in zip(ii, jj, strict=True)] + + +def build_seam_pairs(left: BlockData, right: BlockData) -> list[tuple[int, int, int, int]]: + """Use the package seam helper on the vertical seam between the two blocks.""" + left_lab = left.labels[:, -1:] + right_lab = right.labels[:, :1] + left_val = left.values[:, -1:] + right_val = right.values[:, :1] + left_mask = left.mask[:, -1:] + right_mask = right.mask[:, :1] + + pairs64 = _chunked_seam_pairs_from_strips( + left_lab, + right_lab, + left_val, + right_val, + left_mask, + right_mask, + left_block_id=left.block_id, + right_block_id=right.block_id, + connectivity=CONNECTIVITY, + axis="v", + float_tol=FLOAT_TOL, + ) + + out: list[tuple[int, int, int, int]] = [] + seen: set[tuple[int, int, int]] = set() + + for a, b in pairs64: + left_label = int(np.int64(a) & 0xFFFFFFFF) + right_label = int(np.int64(b) & 0xFFFFFFFF) + + seam_rows = np.where( + (left.labels[:, -1] == left_label) + & (right.labels[:, 0] == right_label) + & (left.values[:, -1] == right.values[:, 0]) + & left.mask[:, -1] + & right.mask[:, 0] + )[0] + if seam_rows.size == 0: + continue + + row = int(seam_rows[0]) + value = int(left.values[row, -1]) + key = (value, left_label, right_label) + if key not in seen: + seen.add(key) + out.append((value, left_label, right_label, row)) + + return out + + +# ----------------------------------------------------------------------------- +# Plot helpers +# ----------------------------------------------------------------------------- + +def draw_cells(ax: plt.Axes, arr: np.ndarray, *, alpha: float = 1.0) -> None: + """Draw colored raster cells.""" + for i in range(arr.shape[0]): + for j in range(arr.shape[1]): + x, y = cell_xy(i, j, arr.shape[0]) + rect = Rectangle( + (x, y), + 1.0, + 1.0, + facecolor=VALUE_COLORS[int(arr[i, j])], + edgecolor="none", + linewidth=0.0, + alpha=alpha, + zorder=1, + ) + ax.add_patch(rect) + + +def draw_values(ax: plt.Axes, arr: np.ndarray) -> None: + """Draw raster values in cell centers.""" + for i in range(arr.shape[0]): + for j in range(arr.shape[1]): + val = int(arr[i, j]) + if val == 0: + continue + xc, yc = cell_center(i, j, arr.shape[0]) + ax.text( + xc, + yc, + str(val), + ha="center", + va="center", + fontsize=11.5, + color=NEUTRAL, + zorder=5, + ) + + +def draw_pixel_grid(ax: plt.Axes) -> None: + """Thin pixel grid for the toy raster.""" + segments: list[np.ndarray] = [] + for x in range(NCOLS + 1): + segments.append(np.array([[x, 0], [x, NROWS]], dtype=float)) + for y in range(NROWS + 1): + segments.append(np.array([[0, y], [NCOLS, y]], dtype=float)) + + ax.add_collection( + LineCollection( + segments, + colors=PIXEL_GRID, + linewidths=0.8, + capstyle="round", + joinstyle="round", + zorder=2, + ) + ) + + +def draw_chunk_boundaries(ax: plt.Axes, geogrids: list[GeoGrid]) -> None: + """Draw only the internal chunk boundary (seam), not the outer perimeter.""" + ax.plot( + [SPLIT_COL, SPLIT_COL], + [0, NROWS], + color=CHUNK_COLOR, + linewidth=2.5, + zorder=6, + solid_capstyle="round", + ) + + +def draw_gdf_boundaries( + ax: plt.Axes, + gdf: gpd.GeoDataFrame, + *, + color: str, + linewidth: float, + zorder: int, +) -> None: + """Plot GeoDataFrame boundaries as linework.""" + if len(gdf) == 0: + return + + try: + boundary = gdf.geometry.union_all().boundary + except Exception: + boundary = gdf.geometry.unary_union.boundary + + geoms = getattr(boundary, "geoms", [boundary]) + segments: list[np.ndarray] = [] + + for geom in geoms: + if geom.is_empty: + continue + if geom.geom_type == "LineString": + segments.append(np.asarray(geom.coords)) + elif geom.geom_type == "MultiLineString": + for sub in geom.geoms: + segments.append(np.asarray(sub.coords)) + + if not segments: + return + + ax.add_collection( + LineCollection( + segments, + colors=color, + linewidths=linewidth, + capstyle="round", + joinstyle="round", + zorder=zorder, + ) + ) + + +def draw_local_labels(ax: plt.Axes, block: BlockData, *, prefix: str) -> None: + """Draw local label ids from actual per-block label rasters.""" + local_ids = sorted(int(v) for v in np.unique(block.labels) if v > 0) + for lid in local_ids: + cells = label_cells_global(block.labels, block, lid) + xc, yc = centroid_from_cells(cells, NROWS) + ax.text( + xc, + yc, + f"{prefix}{lid}", + ha="center", + va="center", + fontsize=9.5, + color=NEUTRAL, + bbox=dict( + boxstyle="round,pad=0.18", + fc="white", + ec=NEUTRAL, + lw=0.8, + alpha=0.95, + ), + zorder=12, + ) + + +def draw_union_links(ax: plt.Axes, seam_pairs: list[tuple[int, int, int, int]], left: BlockData, right: BlockData) -> None: + """Draw seam equivalence links for label_union.""" + for gid, (value, ll, rr, row) in enumerate(seam_pairs, start=1): + left_cells = label_cells_global(left.labels, left, ll) + right_cells = label_cells_global(right.labels, right, rr) + + x1, y1 = centroid_from_cells(left_cells, NROWS) + x2, y2 = centroid_from_cells(right_cells, NROWS) + + conn = ConnectionPatch( + (x1 + 0.35, y1), + (x2 - 0.35, y2), + coordsA="data", + coordsB="data", + axesA=ax, + axesB=ax, + arrowstyle="<->", + color=NEUTRAL, + linewidth=1.5, + mutation_scale=11, + zorder=11, + ) + ax.add_artist(conn) + + ym = NROWS - 1 - row + 0.5 + ax.text( + SPLIT_COL, + ym, + f"G{gid}", + ha="center", + va="center", + fontsize=8.8, + color=NEUTRAL, + bbox=dict( + boxstyle="round,pad=0.15", + fc=VALUE_COLORS[value], + ec=NEUTRAL, + lw=0.7, + alpha=0.95, + ), + zorder=13, + ) + + +def draw_stitch_links(ax: plt.Axes, seam_pairs: list[tuple[int, int, int, int]], left: BlockData, right: BlockData) -> None: + """Draw vector-stitch links across the seam.""" + for _, ll, rr, _row in seam_pairs: + left_rows = np.where(left.labels[:, -1] == ll)[0] + right_rows = np.where(right.labels[:, 0] == rr)[0] + if left_rows.size == 0 or right_rows.size == 0: + continue + + li = int(left_rows[0]) + ri = int(right_rows[0]) + + x1, y1 = cell_center(li, SPLIT_COL - 1, NROWS) + x2, y2 = cell_center(ri, SPLIT_COL, NROWS) + + conn = ConnectionPatch( + (x1 + 0.35, y1), + (x2 - 0.35, y2), + coordsA="data", + coordsB="data", + axesA=ax, + axesB=ax, + arrowstyle="-", + connectionstyle="arc3,rad=0.25", + color=STITCH_COLOR, + linewidth=2.4, + alpha=0.95, + zorder=12, + ) + ax.add_artist(conn) + + +def draw_halo_seam(ax: plt.Axes) -> None: + """Draw halo region along the chunk boundary, extending slightly outside the raster.""" + rect = Rectangle( + (SPLIT_COL - HALO, -0.1), + 2 * HALO, + NROWS + 0.2, + facecolor=STITCH_COLOR, + edgecolor="none", + alpha=0.4, + zorder=3, + clip_on=False, + ) + ax.add_patch(rect) + + +def setup_axis(ax: plt.Axes, *, halo_pad: bool = False) -> None: + """Common axis formatting.""" + extra = 1.1 if halo_pad else 0.0 + ax.set_xlim(-extra, NCOLS + extra) + ax.set_ylim(-extra, NROWS + extra) + ax.set_aspect("equal") + + ax.set_xticks([]) + ax.set_yticks([]) + + for spine in ax.spines.values(): + spine.set_visible(False) + + +def _add_visual_legend(ax: plt.Axes) -> None: + """ + Draw a compact multi-row visual legend inside the bottom-right of an axis. + """ + + # Anchor in axis coordinates + x0 = 1.6 + y0 = 0.74 + dy = 0.15 + line_len = 0.10 + text_dx = 0.15 + + text_kw = dict( + transform=ax.transAxes, + ha="left", + va="center", + fontsize=11, + color=NEUTRAL, + clip_on=False, + ) + + # Row 1: pixel grid + y = y0 + ax.plot( + [x0, x0 + line_len], + [y, y], + transform=ax.transAxes, + color=PIXEL_GRID, + lw=0.8, + solid_capstyle="round", + clip_on=False, + zorder=20, + ) + ax.text(x0 + text_dx, y, "Pixel grid", **text_kw) + + # Row 2: chunk boundary + y = y0 - dy + ax.plot( + [x0, x0 + line_len], + [y, y], + transform=ax.transAxes, + color=CHUNK_COLOR, + lw=2.2, + solid_capstyle="round", + clip_on=False, + zorder=20, + ) + ax.text(x0 + text_dx, y, "Chunk boundary", **text_kw) + + # Row 3: polygon outlines + y = y0 - 2 * dy + ax.plot( + [x0, x0 + line_len], + [y, y], + transform=ax.transAxes, + color=POLYGON_COLOR, + lw=1.6, + solid_capstyle="round", + clip_on=False, + zorder=20, + ) + ax.text(x0 + text_dx, y, "Polygon outlines", **text_kw) + + # Row 4: vector stitch + y = y0 - 3 * dy + arrow = FancyArrowPatch( + (x0, y), + (x0 + line_len, y), + transform=ax.transAxes, + arrowstyle="-", + connectionstyle="arc3,rad=0.4", + lw=2.4, + color=STITCH_COLOR, + clip_on=False, + zorder=20, + ) + ax.add_artist(arrow) + ax.text(x0 + text_dx, y, "Vector stitch", **text_kw) + + # Row 5: halo window + y = y0 - 4 * dy + rect = Rectangle( + (x0, y - 0.025), + line_len, + 0.05, + transform=ax.transAxes, + facecolor=STITCH_COLOR, + edgecolor="none", + alpha=0.4, + clip_on=False, + zorder=20, + ) + ax.add_patch(rect) + ax.text(x0 + text_dx, y, "Halo window", **text_kw) + +def _add_polygonize_title(fig: plt.Figure, strategy_axes: list[plt.Axes]) -> None: + """Add grouped title and underline above the three strategy panels.""" + pos_left = strategy_axes[0].get_position() + pos_right = strategy_axes[-1].get_position() + + x_center = 0.5 * (pos_left.x0 + pos_right.x1) + y_text = pos_left.y1 + 0.14 + y_line = pos_left.y1 + 0.137 + + fig.text( + x_center, + y_text, + "Chunked polygonization strategies", + ha="center", + va="bottom", + fontsize=15, + fontweight="semibold", + color="0.35", + ) + + line = Line2D( + [pos_left.x0, pos_right.x1], + [y_line, y_line], + transform=fig.transFigure, + color="0.6", + lw=1.0, + alpha=0.9, + solid_capstyle="round", + zorder=100, + ) + fig.add_artist(line) + + +# ----------------------------------------------------------------------------- +# Main figure +# ----------------------------------------------------------------------------- + +def make_chunked_polygonize_diagram() -> tuple[plt.Figure, np.ndarray]: + """Build a 2-row schematic diagram for the chunked polygonize strategies.""" + _tiling, block_geogrids, _block_ids, blocks = build_demo_blocks() + left, right = blocks + seam_pairs = build_seam_pairs(left, right) + + fig = plt.figure(figsize=(10, 6)) + gs = fig.add_gridspec( + 2, + 3, + height_ratios=[1.0, 1.35], + hspace=0.32, + wspace=0.18, + ) + + ax_input = fig.add_subplot(gs[0, 1]) + ax_union = fig.add_subplot(gs[1, 0]) + ax_stitch = fig.add_subplot(gs[1, 1]) + ax_geom = fig.add_subplot(gs[1, 2]) + + axes = np.array([ax_input, ax_union, ax_stitch, ax_geom], dtype=object) + + titles = ["Input raster", "label_union", "label_stitch", "geometry_stitch"] + subtitles = [ + "Zoom on raster pixels\nat the vertical boundary\nbetween two chunks", + "1. Label connected components per chunk\n2. Merge labels across chunk boundaries\n3. Polygonize merged " + "labels and dissolve", + "1. Label connected components per chunk\n2. Polygonize labels per chunk\n3. Stitch polygons across chunk boundaries", + "1. Polygonize halo-expanded chunks\n2. Clip polygons to chunk interior\n3. Stitch polygons across chunk boundaries", + ] + + # Top centered panel: input + draw_cells(ax_input, ARR) + draw_pixel_grid(ax_input) + draw_chunk_boundaries(ax_input, block_geogrids) + draw_values(ax_input, ARR) + ax_input.text( + 0.5, + 1.02, + titles[0], + transform=ax_input.transAxes, + ha="center", + va="bottom", + fontsize=14, + fontweight="bold", + color=NEUTRAL, + ) + ax_input.text( + -0.5, + 0.5, + subtitles[0], + transform=ax_input.transAxes, + ha="center", + va="center", + fontsize=10.5, + color=NEUTRAL, + linespacing=1.25, + ) + setup_axis(ax_input) + + # Bottom row: label_union + draw_cells(ax_union, ARR, alpha=0.92) + draw_pixel_grid(ax_union) + draw_chunk_boundaries(ax_union, block_geogrids) + draw_local_labels(ax_union, left, prefix="L") + draw_local_labels(ax_union, right, prefix="R") + draw_union_links(ax_union, seam_pairs, left, right) + ax_union.text( + 0.5, + 1.02, + titles[1], + transform=ax_union.transAxes, + ha="center", + va="bottom", + fontsize=13, + fontweight="bold", + color=NEUTRAL, + ) + ax_union.text( + 0.5, + -0.1, + subtitles[1], + transform=ax_union.transAxes, + ha="center", + va="top", + fontsize=10.2, + color=NEUTRAL, + linespacing=1.25, + ) + setup_axis(ax_union) + + # Bottom row: label_stitch + draw_cells(ax_stitch, ARR, alpha=0.92) + draw_pixel_grid(ax_stitch) + draw_chunk_boundaries(ax_stitch, block_geogrids) + draw_gdf_boundaries(ax_stitch, left.gdf_geom, color=POLYGON_COLOR, linewidth=1.4, zorder=9) + draw_gdf_boundaries(ax_stitch, right.gdf_geom, color=POLYGON_COLOR, linewidth=1.4, zorder=9) + draw_stitch_links(ax_stitch, seam_pairs, left, right) + ax_stitch.text( + 0.5, + 1.02, + titles[2], + transform=ax_stitch.transAxes, + ha="center", + va="bottom", + fontsize=13, + fontweight="bold", + color=NEUTRAL, + ) + ax_stitch.text( + 0.5, + -0.1, + subtitles[2], + transform=ax_stitch.transAxes, + ha="center", + va="top", + fontsize=10.2, + color=NEUTRAL, + linespacing=1.25, + ) + setup_axis(ax_stitch) + + # Bottom row: geometry_stitch + draw_cells(ax_geom, ARR, alpha=0.92) + draw_pixel_grid(ax_geom) + draw_chunk_boundaries(ax_geom, block_geogrids) + draw_halo_seam(ax_geom) + draw_gdf_boundaries(ax_geom, left.gdf_geom, color=POLYGON_COLOR, linewidth=1.4, zorder=9) + draw_gdf_boundaries(ax_geom, right.gdf_geom, color=POLYGON_COLOR, linewidth=1.4, zorder=9) + draw_stitch_links(ax_geom, seam_pairs, left, right) + ax_geom.text( + 0.5, + 1.02, + titles[3], + transform=ax_geom.transAxes, + ha="center", + va="bottom", + fontsize=13, + fontweight="bold", + color=NEUTRAL, + ) + ax_geom.text( + 0.5, + -0.1, + subtitles[3], + transform=ax_geom.transAxes, + ha="center", + va="top", + fontsize=10.2, + color=NEUTRAL, + linespacing=1.25, + ) + setup_axis(ax_geom) + + _add_visual_legend(ax_input) + _add_polygonize_title(fig, [ax_union, ax_stitch, ax_geom]) + + fig.subplots_adjust(left=0.01, right=0.99, top=0.95, bottom=0.15) + + return fig, axes + + +fig, _ = make_chunked_polygonize_diagram() +plt.show() diff --git a/doc/source/code/diagram_chunked_rasterize.py b/doc/source/code/diagram_chunked_rasterize.py new file mode 100644 index 000000000..6db4bd4c0 --- /dev/null +++ b/doc/source/code/diagram_chunked_rasterize.py @@ -0,0 +1,692 @@ +"""Script to make diagram for chunked rasterize in documentation.""" +from __future__ import annotations + +import matplotlib.pyplot as plt +import numpy as np +import rasterio as rio +from matplotlib.collections import LineCollection +from matplotlib.lines import Line2D +from matplotlib.patches import FancyArrowPatch, Polygon as MplPolygon, Rectangle, ConnectionPatch +from shapely.geometry import Polygon, box + +from geoutils.interface.rasterization import _rasterio_rasterize_burn + +# ----------------------------------------------------------------------------- +# Example raster grid, chunks, and vector geometries +# ----------------------------------------------------------------------------- + +NROWS = 8 +NCOLS = 10 + +CHUNK_ROWS = (4, 4) +CHUNK_COLS = (5, 5) + +# Highlight bottom-right chunk +HIGHLIGHT_CHUNK = (1, 1) + +NEUTRAL = "#333333" +PIXEL_GRID = "0.80" +CHUNK_COLOR = "#222222" +HIGHLIGHT_COLOR = "#F58518" +POLYGON_COLOR = "#6BAED6" +POLYGON_EDGE = "#4C78A8" +FADED_ALPHA = 0.22 + +# Three example polygons in raster coordinates +POLYGONS = [ + Polygon([(0.8, 6.8), (3.2, 7.4), (4.1, 5.7), (2.4, 4.8), (0.9, 5.6)]), + Polygon([(4.2, 5.9), (7.4, 6.7), (8.2, 4.5), (6.7, 3.2), (4.6, 4.0)]), + Polygon([(6.0, 2.5), (8.7, 2.9), (9.1, 0.8), (6.8, 0.6), (5.7, 1.4)]), +] +POLY_VALUES = np.array([1, 2, 3], dtype=np.uint8) + + +# ----------------------------------------------------------------------------- +# Geometry helpers +# ----------------------------------------------------------------------------- + +def _data_point_to_fig(fig: plt.Figure, ax: plt.Axes, x_data: float, y_data: float) -> tuple[float, float]: + """Convert a data-coordinate point to figure coordinates.""" + xy_disp = ax.transData.transform((x_data, y_data)) + xy_fig = fig.transFigure.inverted().transform(xy_disp) + return float(xy_fig[0]), float(xy_fig[1]) + +def chunk_bounds(chunk_location: tuple[int, int]) -> tuple[float, float, float, float]: + """Return chunk bounds as (left, bottom, right, top).""" + iy, ix = chunk_location + + row_starts = np.concatenate([[0], np.cumsum(CHUNK_ROWS)]) + col_starts = np.concatenate([[0], np.cumsum(CHUNK_COLS)]) + + top_row = row_starts[iy] + bottom_row = row_starts[iy + 1] + left_col = col_starts[ix] + right_col = col_starts[ix + 1] + + left = float(left_col) + right = float(right_col) + top = float(NROWS - top_row) + bottom = float(NROWS - bottom_row) + + return left, bottom, right, top + + +def polygon_to_patch(poly: Polygon, **kwargs) -> MplPolygon: + """Convert shapely polygon to matplotlib patch.""" + return MplPolygon(np.asarray(poly.exterior.coords), closed=True, **kwargs) + + +def geometry_intersects_chunk(poly: Polygon, bounds: tuple[float, float, float, float]) -> bool: + """Return whether polygon intersects chunk bbox.""" + left, bottom, right, top = bounds + qbox = box(left, bottom, right, top) + return poly.intersects(qbox) + + +def build_rasterized_chunk( + bounds: tuple[float, float, float, float], + polygons: list[Polygon], + values: np.ndarray, +) -> np.ndarray: + """Rasterize candidate polygons into the highlighted chunk using the package helper.""" + left, bottom, right, top = bounds + width = int(right - left) + height = int(top - bottom) + + transform = rio.transform.from_bounds( + west=left, + south=bottom, + east=right, + north=top, + width=width, + height=height, + ) + + geoms = np.asarray(polygons, dtype=object) + + return _rasterio_rasterize_burn( + geoms=geoms, + values=values, + default_value=None, + out_shape=(height, width), + transform=transform, + fill=0, + dtype=np.uint8, + all_touched=False, + ) + + +# ----------------------------------------------------------------------------- +# Plot helpers +# ----------------------------------------------------------------------------- + +def draw_pixel_grid(ax: plt.Axes, *, nrows: int, ncols: int, extend: float = 1.0) -> None: + """Draw thin pixel grid, extended slightly outside the plotted raster.""" + segments: list[np.ndarray] = [] + + # Vertical pixel lines + for x in range(ncols + 1): + segments.append(np.array([[x, -extend], [x, nrows + extend]], dtype=float)) + + # Horizontal pixel lines + for y in range(nrows + 1): + segments.append(np.array([[-extend, y], [ncols + extend, y]], dtype=float)) + + coll = LineCollection( + segments, + colors=PIXEL_GRID, + linewidths=0.8, + capstyle="round", + joinstyle="round", + zorder=1, + clip_on=False, + ) + ax.add_collection(coll) + +def draw_chunk_boundaries(ax: plt.Axes, *, extend: float = 1.0) -> None: + """Draw all chunk boundaries, extended slightly outside the plotted raster.""" + row_starts = np.concatenate([[0], np.cumsum(CHUNK_ROWS)]) + col_starts = np.concatenate([[0], np.cumsum(CHUNK_COLS)]) + + # Draw all vertical chunk boundaries, including outside ones + for x in col_starts: + ax.plot( + [x, x], + [-extend, NROWS + extend], + color=CHUNK_COLOR, + linewidth=2.3, + zorder=3, + solid_capstyle="round", + clip_on=False, + ) + + # Draw all horizontal chunk boundaries, including outside ones + for y_idx in row_starts: + y = NROWS - y_idx + ax.plot( + [-extend, NCOLS + extend], + [y, y], + color=CHUNK_COLOR, + linewidth=2.3, + zorder=3, + solid_capstyle="round", + clip_on=False, + ) + +def draw_highlight_chunk(ax: plt.Axes, bounds: tuple[float, float, float, float]) -> None: + """Draw highlighted chunk bbox.""" + left, bottom, right, top = bounds + rect = Rectangle( + (left, bottom), + right - left, + top - bottom, + fill=False, + edgecolor=HIGHLIGHT_COLOR, + linewidth=2.8, + zorder=5, + ) + ax.add_patch(rect) + +def draw_polygon_ids( + ax: plt.Axes, + polygons: list[Polygon], + ids: list[int], + *, + candidate_mask: list[bool] | None = None, +) -> None: + """Draw polygon IDs at their centroid.""" + for i, poly in enumerate(polygons): + + is_candidate = True if candidate_mask is None else candidate_mask[i] + + cx, cy = poly.centroid.coords[0] + + ax.text( + cx, + cy, + f"ID: {ids[i]}", + ha="center", + va="center", + fontsize=9, + color=NEUTRAL, + zorder=20, + alpha=1.0 if is_candidate else 0.7, + bbox=dict( + boxstyle="round,pad=0.2", + fc="white", + ec="none", + alpha=0.8, + ), + ) + +def draw_chunk_query_box(ax: plt.Axes, bounds: tuple[float, float, float, float]) -> None: + """Draw filled highlight for candidate-query chunk.""" + left, bottom, right, top = bounds + + rect_fill = Rectangle( + (left, bottom), + right - left, + top - bottom, + facecolor=HIGHLIGHT_COLOR, + edgecolor="none", + alpha=0.10, + zorder=2, + ) + rect_edge = Rectangle( + (left, bottom), + right - left, + top - bottom, + fill=False, + edgecolor=HIGHLIGHT_COLOR, + linewidth=2.8, + zorder=5, + ) + ax.add_patch(rect_fill) + ax.add_patch(rect_edge) + + +def draw_polygons( + ax: plt.Axes, + polygons: list[Polygon], + *, + candidate_mask: list[bool] | None = None, +) -> None: + """Draw input polygons, optionally fading non-candidates.""" + for i, poly in enumerate(polygons): + is_candidate = True if candidate_mask is None else candidate_mask[i] + alpha = 0.45 if is_candidate else FADED_ALPHA + edgecolor = POLYGON_EDGE if is_candidate else "#7A8FA6" + facecolor = POLYGON_COLOR if is_candidate else "#BFD7EA" + + patch = polygon_to_patch( + poly, + facecolor=facecolor, + edgecolor=edgecolor, + linewidth=1.5, + alpha=alpha, + zorder=4 if is_candidate else 2, + ) + ax.add_patch(patch) + + +def draw_chunk_result(ax: plt.Axes, arr: np.ndarray) -> None: + """Draw rasterized chunk result as a small chunk-local raster.""" + nrows, ncols = arr.shape + + value_colors = { + 0: "#FFFFFF", + 1: "#6BAED6", + 2: "#9ECAE1", + 3: "#C6DBEF", + } + + for i in range(nrows): + for j in range(ncols): + val = int(arr[i, j]) + y = nrows - 1 - i + x = j + rect = Rectangle( + (x, y), + 1.0, + 1.0, + facecolor=value_colors[val], + edgecolor="none", + linewidth=0.0, + zorder=1, + ) + ax.add_patch(rect) + + if val != 0: + ax.text( + x + 0.5, + y + 0.5, + str(val), + ha="center", + va="center", + fontsize=11, + color=NEUTRAL, + zorder=4, + ) + + segments: list[np.ndarray] = [] + for x in range(ncols + 1): + segments.append(np.array([[x, 0], [x, nrows]], dtype=float)) + for y in range(nrows + 1): + segments.append(np.array([[0, y], [ncols, y]], dtype=float)) + + ax.add_collection( + LineCollection( + segments, + colors=PIXEL_GRID, + linewidths=0.8, + capstyle="round", + joinstyle="round", + zorder=2, + ) + ) + + rect = Rectangle( + (0, 0), + ncols, + nrows, + fill=False, + edgecolor=HIGHLIGHT_COLOR, + linewidth=2.8, + zorder=5, + ) + ax.add_patch(rect) + + +def setup_axis(ax: plt.Axes, *, nrows: int, ncols: int, pad: float = 0.6) -> None: + """Common axis formatting.""" + ax.set_xlim(-pad, ncols + pad) + ax.set_ylim(-pad, nrows + pad) + ax.set_aspect("equal") + ax.set_xticks([]) + ax.set_yticks([]) + for spine in ax.spines.values(): + spine.set_visible(False) + + +def add_mapblocks_title(fig: plt.Figure, axes: list[plt.Axes]) -> None: + """Add grouped title over panels 2–3.""" + pos_left = axes[0].get_position() + pos_right = axes[-1].get_position() + + x_center = 0.5 * (pos_left.x0 + pos_right.x1) + y_text = pos_left.y1 + 0.195 + y_line = pos_left.y1 + 0.182 + + fig.text( + x_center, + y_text, + "Chunked rasterization", + ha="center", + va="bottom", + fontsize=15, + fontweight="semibold", + color="0.35", + ) + + line = Line2D( + [pos_left.x0, pos_right.x1], + [y_line, y_line], + transform=fig.transFigure, + color="0.6", + lw=1.0, + alpha=0.9, + solid_capstyle="round", + zorder=100, + ) + fig.add_artist(line) + + +def add_visual_legend(fig: plt.Figure) -> None: + """Draw horizontal bottom legend.""" + y = 0.06 + x = 0.10 + dx = 0.18 + + # Pixel grid + fig.add_artist( + Line2D( + [x, x + 0.03], + [y, y], + transform=fig.transFigure, + color=PIXEL_GRID, + lw=0.8, + solid_capstyle="round", + ) + ) + fig.text(x + 0.038, y, "Pixel grid", transform=fig.transFigure, va="center", fontsize=10) + + # Chunk boundary + x += dx + fig.add_artist( + Line2D( + [x, x + 0.03], + [y, y], + transform=fig.transFigure, + color=CHUNK_COLOR, + lw=2.3, + ) + ) + fig.text(x + 0.038, y, "Chunk grid", transform=fig.transFigure, va="center", fontsize=10) + + # Queried chunk + x += dx + rect = Rectangle( + (x, y - 0.010), + 0.03, + 0.020, + transform=fig.transFigure, + facecolor=HIGHLIGHT_COLOR, + edgecolor=HIGHLIGHT_COLOR, + linewidth=1.4, + alpha=0.15, + ) + fig.add_artist(rect) + fig.text(x + 0.038, y, "Queried chunk", transform=fig.transFigure, va="center", fontsize=10) + + # Candidate geometry + x += dx + poly = MplPolygon( + np.array( + [ + [x, y - 0.010], + [x + 0.010, y + 0.012], + [x + 0.022, y + 0.010], + [x + 0.028, y - 0.008], + [x + 0.012, y - 0.014], + ] + ), + closed=True, + transform=fig.transFigure, + facecolor=POLYGON_COLOR, + edgecolor=POLYGON_EDGE, + linewidth=1.0, + alpha=0.45, + ) + fig.add_artist(poly) + fig.text(x + 0.038, y, "Candidate geometry", transform=fig.transFigure, va="center", fontsize=10) + +def chunk_center(bounds: tuple[float, float, float, float]) -> tuple[float, float]: + """Return center of chunk bounds.""" + left, bottom, right, top = bounds + return 0.5 * (left + right), 0.5 * (bottom + top) + + +def add_between_panel_arrows( + fig: plt.Figure, + ax0: plt.Axes, + ax1: plt.Axes, + ax2: plt.Axes, + highlight_bounds: tuple[float, float, float, float], + chunk_result_shape: tuple[int, int], +) -> None: + """Add arrows anchored to meaningful data coordinates in each panel.""" + c01 = chunk_center(highlight_bounds) + c12 = chunk_center(highlight_bounds) + + # Center of chunk-local raster in panel 3 + nrows, ncols = chunk_result_shape + c2 = (0.5 * ncols, 0.5 * nrows) + + # Panel 1 to Panel 2 + arr01 = ConnectionPatch( + xyA=c01, + coordsA=ax0.transData, + xyB=c01, + coordsB=ax1.transData, + arrowstyle="-|>", + linewidth=2.2, + color=NEUTRAL, + mutation_scale=22, + shrinkA=8, + shrinkB=8, + connectionstyle="arc3,rad=0.0", + zorder=30, + clip_on=False, + ) + fig.add_artist(arr01) + + # Panel 2 to Panel 3 + arr12 = ConnectionPatch( + xyA=c12, + coordsA=ax1.transData, + xyB=c2, + coordsB=ax2.transData, + arrowstyle="-|>", + linewidth=2.2, + color=NEUTRAL, + mutation_scale=22, + shrinkA=8, + shrinkB=8, + connectionstyle="arc3,rad=0.0", + zorder=30, + clip_on=False, + ) + fig.add_artist(arr12) + +# ----------------------------------------------------------------------------- +# Main figure +# ----------------------------------------------------------------------------- + +def make_chunked_rasterize_diagram() -> tuple[plt.Figure, np.ndarray]: + """Build a 3-panel schematic for chunked rasterize.""" + highlight_bounds = chunk_bounds(HIGHLIGHT_CHUNK) + candidate_mask = [geometry_intersects_chunk(poly, highlight_bounds) for poly in POLYGONS] + candidate_polygons = [poly for poly, keep in zip(POLYGONS, candidate_mask, strict=True) if keep] + candidate_values = POLY_VALUES[np.array(candidate_mask, dtype=bool)] + chunk_result = build_rasterized_chunk(highlight_bounds, candidate_polygons, candidate_values) + + fig = plt.figure(figsize=(10, 4)) + gs = fig.add_gridspec(1, 3, wspace=0.28) + + ax0 = fig.add_subplot(gs[0, 0]) + ax1 = fig.add_subplot(gs[0, 1]) + ax2 = fig.add_subplot(gs[0, 2]) + + axes = np.array([ax0, ax1, ax2], dtype=object) + + titles = [ + "Output raster grid", + "Per-chunk geometry filtering", + "Per-chunk rasterization", + ] + subtitles = [ + "Chunked output raster overlaid with\nvector geometries to burn", + "Only geometries intersecting the highlighted\nchunk bounds are passed to this block", + "Each chunk rasterizes only its candidate\ngeometries, or exits early if none intersect", + ] + + # Panel 1 + draw_pixel_grid(ax0, nrows=NROWS, ncols=NCOLS) + draw_chunk_boundaries(ax0) + draw_polygons(ax0, POLYGONS) + draw_polygon_ids(ax0, POLYGONS, POLY_VALUES.tolist()) + draw_highlight_chunk(ax0, highlight_bounds) + ax0.text( + 0.5, + 1.05, + titles[0], + transform=ax0.transAxes, + ha="center", + va="bottom", + fontsize=13, + fontweight="bold", + color=NEUTRAL, + ) + ax0.text( + 0.5, + -0.07, + subtitles[0], + transform=ax0.transAxes, + ha="center", + va="top", + fontsize=10.2, + color=NEUTRAL, + linespacing=1.25, + ) + setup_axis(ax0, nrows=NROWS, ncols=NCOLS, pad=1.2) + + # Panel 2 + draw_pixel_grid(ax1, nrows=NROWS, ncols=NCOLS) + draw_chunk_boundaries(ax1) + draw_polygons(ax1, POLYGONS, candidate_mask=candidate_mask) + draw_polygon_ids(ax1, POLYGONS, POLY_VALUES.tolist(), candidate_mask=candidate_mask) + draw_chunk_query_box(ax1, highlight_bounds) + ax1.text( + 0.5, + 1.05, + titles[1], + transform=ax1.transAxes, + ha="center", + va="bottom", + fontsize=13, + fontweight="bold", + color=NEUTRAL, + ) + ax1.text( + 0.5, + -0.07, + subtitles[1], + transform=ax1.transAxes, + ha="center", + va="top", + fontsize=10.2, + color=NEUTRAL, + linespacing=1.25, + ) + setup_axis(ax1, nrows=NROWS, ncols=NCOLS, pad=1.2) + + # Panel 3 + draw_chunk_result(ax2, chunk_result) + ax2.text( + 0.5, + 1.05, + titles[2], + transform=ax2.transAxes, + ha="center", + va="bottom", + fontsize=13, + fontweight="bold", + color=NEUTRAL, + ) + ax2.text( + 0.5, + -0.07, + subtitles[2], + transform=ax2.transAxes, + ha="center", + va="top", + fontsize=10.2, + color=NEUTRAL, + linespacing=1.25, + ) + setup_axis(ax2, nrows=chunk_result.shape[0], ncols=chunk_result.shape[1], pad=0.6) + + + add_mapblocks_title(fig, [ax0, ax2]) + add_visual_legend(fig) + + fig.subplots_adjust(left=0.02, right=0.98, top=0.96, bottom=0.1) + + # Finalize transforms before placing arrows + fig.canvas.draw() + + # Horizontal chunk boundary in panels 1 and 2 + y_chunk_boundary = NROWS - CHUNK_ROWS[0] + + # For panel 3, align with the middle of the chunk-local raster + y_chunk_result = chunk_result.shape[0] / 2 + + pos0 = ax0.get_position() + pos1 = ax1.get_position() + pos2 = ax2.get_position() + + # Arrow 1: panel 1 to panel 2 + _, y0 = _data_point_to_fig(fig, ax0, NCOLS / 2, y_chunk_boundary) + _, y1 = _data_point_to_fig(fig, ax1, NCOLS / 2, y_chunk_boundary) + y01 = 0.5 * (y0 + y1) + + arrow01 = FancyArrowPatch( + (pos0.x1 + 0.015, y01), + (pos1.x0 - 0.015, y01), + transform=fig.transFigure, + arrowstyle="-|>", + connectionstyle="arc3,rad=0.0", + linewidth=2.0, + color=NEUTRAL, + mutation_scale=20, + shrinkA=0, + shrinkB=0, + zorder=30, + ) + fig.add_artist(arrow01) + + # Arrow 2: panel 2 to panel 3 + _, y1b = _data_point_to_fig(fig, ax1, NCOLS / 2, y_chunk_boundary) + _, y2 = _data_point_to_fig(fig, ax2, chunk_result.shape[1] / 2, y_chunk_result) + y12 = 0.5 * (y1b + y2) + + arrow12 = FancyArrowPatch( + (pos1.x1 + 0.015, y12), + (pos2.x0 - 0.015, y12), + transform=fig.transFigure, + arrowstyle="-|>", + connectionstyle="arc3,rad=0.0", + linewidth=2.0, + color=NEUTRAL, + mutation_scale=20, + shrinkA=0, + shrinkB=0, + zorder=30, + ) + fig.add_artist(arrow12) + + return fig, axes + + +fig, _ = make_chunked_rasterize_diagram() +plt.show() \ No newline at end of file diff --git a/doc/source/code/diagram_chunked_reproject.py b/doc/source/code/diagram_chunked_reproject.py new file mode 100644 index 000000000..e3c6eaf5d --- /dev/null +++ b/doc/source/code/diagram_chunked_reproject.py @@ -0,0 +1,850 @@ +"""Script to make diagram for chunked reproject in documentation.""" +from __future__ import annotations + +import geopandas as gpd +import matplotlib.pyplot as plt +from matplotlib.patches import Polygon as MplPolygon +from matplotlib.collections import PatchCollection +import numpy as np +import rasterio as rio +from shapely.geometry import Polygon, MultiPolygon +from matplotlib.patches import FancyArrowPatch +from matplotlib.collections import LineCollection +from shapely.geometry import LineString, MultiLineString +from geoutils.multiproc.chunked import GeoGrid, ChunkedGeoGrid, cached_cumsum + + +# ----------------------------------------------------------------------------- +# Plotting utilities +# ----------------------------------------------------------------------------- + +def _geom_to_lines(geom) -> list[np.ndarray]: + """Convert shapely line geometry to matplotlib line segments.""" + lines: list[np.ndarray] = [] + + if geom.is_empty: + return lines + + if isinstance(geom, LineString): + lines.append(np.asarray(geom.coords)) + elif isinstance(geom, MultiLineString): + for line in geom.geoms: + lines.extend(_geom_to_lines(line)) + + return lines + + +def _add_gdf_boundary( + ax: plt.Axes, + gdf: gpd.GeoDataFrame, + *, + color: str = "black", + linewidth: float = 1.0, + alpha: float = 1.0, + zorder: int = 1, +) -> None: + """Draw GeoDataFrame boundaries as linework, with each edge drawn only once.""" + boundary = gdf.geometry.boundary.union_all() + lines = _geom_to_lines(boundary) + + if len(lines) == 0: + return + + collection = LineCollection( + lines, + colors=color, + linewidths=linewidth, + alpha=alpha, + zorder=zorder, + capstyle="round", + joinstyle="round", + ) + ax.add_collection(collection) + + +def _geom_to_patches(geom) -> list[MplPolygon]: + """Convert shapely polygon / multipolygon to matplotlib patches.""" + patches: list[MplPolygon] = [] + + if geom.is_empty: + return patches + + if isinstance(geom, Polygon): + patches.append(MplPolygon(np.asarray(geom.exterior.coords), closed=True)) + elif isinstance(geom, MultiPolygon): + for poly in geom.geoms: + patches.extend(_geom_to_patches(poly)) + else: + pass + + return patches + + +def _add_gdf( + ax: plt.Axes, + gdf: gpd.GeoDataFrame, + *, + facecolor: str = "none", + edgecolor: str = "black", + linewidth: float = 1.0, + alpha: float = 1.0, + zorder: int = 1, +) -> None: + """Add all geometries of a GeoDataFrame as matplotlib patches.""" + patches: list[MplPolygon] = [] + for geom in gdf.geometry: + patches.extend(_geom_to_patches(geom)) + + if len(patches) == 0: + return + + collection = PatchCollection( + patches, + facecolor=facecolor, + edgecolor=edgecolor, + linewidth=linewidth, + alpha=alpha, + zorder=zorder, + capstyle="round", + joinstyle="round", + ) + ax.add_collection(collection) + + +def _set_ax_extent_from_gdfs(ax: plt.Axes, gdfs: list[gpd.GeoDataFrame], pad_frac: float = 0.08) -> None: + """Set equal aspect and padded extent from one or more GeoDataFrames.""" + bounds = np.array([gdf.total_bounds for gdf in gdfs], dtype=float) + xmin = np.min(bounds[:, 0]) + ymin = np.min(bounds[:, 1]) + xmax = np.max(bounds[:, 2]) + ymax = np.max(bounds[:, 3]) + + dx = xmax - xmin + dy = ymax - ymin + padx = max(dx * pad_frac, 1e-12) + pady = max(dy * pad_frac, 1e-12) + + ax.set_xlim(xmin - padx, xmax + padx) + ax.set_ylim(ymin - pady, ymax + pady) + ax.set_aspect("equal") + ax.axis("off") + + +def _add_rectilinear_grid( + ax: plt.Axes, + grid: GeoGrid, + *, + color: str = "0.75", + linewidth: float = 0.55, + alpha: float = 0.8, + zorder: int = 0, +) -> None: + """Draw the internal pixel grid of a rectilinear GeoGrid.""" + left, bottom, right, top = grid.bounds + dx, dy = grid.res + + xs = left + np.arange(grid.width + 1) * dx + ys = top - np.arange(grid.height + 1) * dy + + segments: list[np.ndarray] = [] + + for x in xs: + segments.append(np.array([[x, bottom], [x, top]], dtype=float)) + + for y in ys: + segments.append(np.array([[left, y], [right, y]], dtype=float)) + + collection = LineCollection( + segments, + colors=color, + linewidths=linewidth, + alpha=alpha, + zorder=zorder, + capstyle="round", + joinstyle="round", + ) + ax.add_collection(collection) + +def _add_chunk_size_arrows( + ax: plt.Axes, + chunk: GeoGrid, + *, + label_x: str, + label_y: str, + color: str = "#333333", + fontsize: float = 8.5, +) -> None: + """Add double-headed arrows showing the width/height of a chunk. + + The arrows are positioned using display coordinates so they look visually + consistent even when x/y data units differ strongly. + """ + left, bottom, right, top = chunk.bounds + + # Chunk corners in data coordinates + p_bl = np.array([left, bottom]) + p_br = np.array([right, bottom]) + p_tl = np.array([left, top]) + + # Transform to display coordinates + to_disp = ax.transData.transform + to_data = ax.transData.inverted().transform + + bl_d = to_disp(p_bl) + br_d = to_disp(p_br) + tl_d = to_disp(p_tl) + + # Visual chunk size in display units + width_d = br_d[0] - bl_d[0] + height_d = tl_d[1] - bl_d[1] + + # Offsets/padding in display units, so appearance is consistent + ref_d = min(abs(width_d), abs(height_d)) + off_d = 0.25 * ref_d + pad_d = 0.15 * ref_d + + h0_d = np.array([bl_d[0] - pad_d, bl_d[1] - off_d]) + h1_d = np.array([br_d[0] + pad_d, bl_d[1] - off_d]) + + h0 = to_data(h0_d) + h1 = to_data(h1_d) + + ax.add_patch( + FancyArrowPatch( + tuple(h0), + tuple(h1), + arrowstyle="<->", + mutation_scale=11, + linewidth=1.1, + color=color, + zorder=10, + clip_on=False, + ) + ) + + htxt_d = np.array([0.5 * (h0_d[0] + h1_d[0]), bl_d[1] - 1.35 * off_d]) + htxt = to_data(htxt_d) + ax.text( + htxt[0], + htxt[1], + label_x, + ha="center", + va="top", + fontsize=fontsize, + color=color, + ) + + # Vertical arrow left of chunk (in display coordinates) + v0_d = np.array([bl_d[0] - off_d, bl_d[1] - pad_d]) + v1_d = np.array([bl_d[0] - off_d, tl_d[1] + pad_d]) + + v0 = to_data(v0_d) + v1 = to_data(v1_d) + + ax.add_patch( + FancyArrowPatch( + tuple(v0), + tuple(v1), + arrowstyle="<->", + mutation_scale=11, + linewidth=1.1, + color=color, + zorder=10, + clip_on=False, + ) + ) + + vtxt_d = np.array([bl_d[0] - 1.35 * off_d, 0.5 * (v0_d[1] + v1_d[1])]) + vtxt = to_data(vtxt_d) + ax.text( + vtxt[0], + vtxt[1], + label_y, + ha="right", + va="center", + rotation=90, + fontsize=fontsize, + color=color, + ) + +def _add_line_legend( + ax: plt.Axes, + *, + thin_label: str = "Pixel grid", + thick_label: str = "Chunk grid", + thin_color: str = "0.75", + thick_color: str = "#6BAED6", + thick_color_2: str = "#F58518", + dest_color: str = "#F58518", +) -> None: + """Add a compact inline legend in axis coordinates.""" + + xt = 0.65 + x0 = xt - 0.18 + x1 = xt - 0.12 + x2 = xt - 0.06 + + y_top = -0.40 + y_bot = -0.58 + + # Thin pixel line + ax.plot( + [x0, x2], + [y_top, y_top], + transform=ax.transAxes, + color=thin_color, + lw=0.8, + solid_capstyle="round", + clip_on=False, + zorder=1, + ) + + ax.text( + xt, + y_top, + thin_label, + transform=ax.transAxes, + ha="left", + va="center", + fontsize=8.5, + color="#333333", + zorder=3, + ) + + # Thick source chunk line (blue) + ax.plot( + [x0, x1], + [y_bot, y_bot], + transform=ax.transAxes, + color=thick_color, + lw=1.8, + solid_capstyle="round", + clip_on=False, + zorder=1, + ) + + # Thick destination chunk line (orange) + ax.plot( + [x1 + 0.01, x2 + 0.01], + [y_bot, y_bot], + transform=ax.transAxes, + color=thick_color_2, + lw=1.8, + solid_capstyle="round", + clip_on=False, + zorder=1, + ) + + ax.text( + xt, + y_bot, + thick_label, + transform=ax.transAxes, + ha="left", + va="center", + fontsize=8.5, + color="#333333", + zorder=3, + ) + +def build_demo_destination_grid( + source_grid: GeoGrid, + dst_crs: rio.crs.CRS, + dst_shape: tuple[int, int], +) -> GeoGrid: + """ + Build a destination GeoGrid from projected source bounds. + """ + dst_bounds = source_grid.bounds_projected(dst_crs) + dst_transform = rio.transform.from_bounds( + west=dst_bounds.left, + south=dst_bounds.bottom, + east=dst_bounds.right, + north=dst_bounds.top, + width=dst_shape[1], + height=dst_shape[0], + ) + return GeoGrid(transform=dst_transform, shape=dst_shape, crs=dst_crs) + + +def plot_reprojection_chunk_diagram( + source_chunked_grid: ChunkedGeoGrid, + destination_chunked_grid: ChunkedGeoGrid, + highlight_destination_chunk: tuple[int, int] = (1, 1), + source_color: str = "#6BAED6", + intersect_color: str = "#4C78A8", + destination_color: str = "#F58518", + figsize: tuple[float, float] = (10.0, 4.8), +) -> tuple[plt.Figure, np.ndarray]: + """ + Plot a 2-panel diagram illustrating chunked reprojection. + + Left: Source raster with source chunks in source CRS. + Right: Same source chunks projected to destination CRS, with destination chunk grid overlaid. + """ + # ------------------------------------------------------------------------- + # Styling + # ------------------------------------------------------------------------- + source_fill_alpha = 0.07 + source_chunk_alpha = 0.22 + intersect_alpha = 0.58 + + grid_lw = 1.3 + highlight_lw = 2.6 + arrow_lw = 1.6 + pixel_grid_lw = 0.55 + + neutral_color = "#333333" + pixel_grid_color = "0.75" + + # ------------------------------------------------------------------------- + # Geometries + # ------------------------------------------------------------------------- + src_full = source_chunked_grid.grid.footprint + src_blocks = source_chunked_grid.get_block_footprints() + + dst_full = destination_chunked_grid.grid.footprint + dst_blocks = destination_chunked_grid.get_block_footprints() + + src_full_in_dst = source_chunked_grid.grid.footprint_projected(destination_chunked_grid.grid.crs) + src_blocks_in_dst = source_chunked_grid.get_block_footprints(crs=destination_chunked_grid.grid.crs) + + # ------------------------------------------------------------------------- + # Create figure + # ------------------------------------------------------------------------- + fig, axes = plt.subplots( + 1, + 2, + figsize=(8.0, 3.5), + gridspec_kw={"width_ratios": [0.75, 1]}, + ) + ax0, ax1 = axes + + # ------------------------------------------------------------------------- + # Left panel: source raster + # ------------------------------------------------------------------------- + _add_rectilinear_grid( + ax0, + source_chunked_grid.grid, + color=pixel_grid_color, + linewidth=pixel_grid_lw, + alpha=0.8, + zorder=0, + ) + + _add_gdf( + ax0, + src_full, + facecolor=source_color, + edgecolor="none", + linewidth=0, + alpha=source_fill_alpha, + zorder=1, + ) + + _add_gdf( + ax0, + src_blocks, + facecolor=source_color, + edgecolor="none", + linewidth=0, + alpha=source_chunk_alpha, + zorder=2, + ) + + _add_gdf_boundary( + ax0, + src_blocks, + color=source_color, + linewidth=grid_lw, + alpha=1, + zorder=3, + ) + + src_blocks_list = source_chunked_grid.get_blocks_as_geogrids() + ny, nx = source_chunked_grid.num_chunks + + bottom_left_chunk = src_blocks_list[(ny - 1) * nx] + _add_chunk_size_arrows( + ax0, + bottom_left_chunk, + label_x=f"{bottom_left_chunk.width} px", + label_y=f"{bottom_left_chunk.height} px", + color=neutral_color, + fontsize=8.5, + ) + + ax0.text( + 0.5, + 0.95, + "Source raster", + transform=ax0.transAxes, + ha="center", + va="bottom", + fontsize=12.5, + fontweight="bold", + color=source_color, + ) + + _set_ax_extent_from_gdfs(ax0, [src_full, src_blocks]) + + _add_line_legend( + ax1, + thin_label="Pixel grid", + thick_label="Chunk grid", + thin_color=pixel_grid_color, + thick_color=source_color, + thick_color_2=destination_color, + ) + + # ------------------------------------------------------------------------- + # Right panel: projected source + destination raster + # ------------------------------------------------------------------------- + _add_rectilinear_grid( + ax1, + destination_chunked_grid.grid, + color=pixel_grid_color, + linewidth=pixel_grid_lw, + alpha=1, + zorder=0, + ) + + _add_gdf( + ax1, + src_full_in_dst, + facecolor=source_color, + edgecolor="none", + linewidth=0, + alpha=source_fill_alpha, + zorder=1, + ) + + highlight_index = destination_chunked_grid.flat_block_index(highlight_destination_chunk) + dst_blocks_one = dst_blocks.iloc[[highlight_index]] + highlight_geom = dst_blocks_one.geometry.iloc[0] + + intersection_area = src_blocks_in_dst.geometry.intersection(highlight_geom).area + intersect_mask = intersection_area > 0 + + src_intersect = src_blocks_in_dst[intersect_mask] + src_non_intersect = src_blocks_in_dst[~intersect_mask] + + _add_gdf( + ax1, + src_non_intersect, + facecolor=source_color, + edgecolor="none", + linewidth=0, + alpha=source_chunk_alpha, + zorder=2, + ) + + _add_gdf( + ax1, + src_intersect, + facecolor=intersect_color, + edgecolor="none", + linewidth=0, + alpha=intersect_alpha, + zorder=3, + ) + + _add_gdf_boundary( + ax1, + src_blocks_in_dst, + color=source_color, + linewidth=grid_lw, + alpha=0.9, + zorder=4, + ) + + _add_gdf_boundary( + ax1, + dst_blocks, + color=destination_color, + linewidth=grid_lw, + alpha=1.0, + zorder=5, + ) + + _add_gdf( + ax1, + dst_blocks_one, + facecolor="none", + edgecolor=destination_color, + linewidth=highlight_lw, + alpha=1.0, + zorder=6, + ) + + dst_blocks_list = destination_chunked_grid.get_blocks_as_geogrids() + ny_d, nx_d = destination_chunked_grid.num_chunks + + bottom_left_dst_chunk = dst_blocks_list[(ny_d - 1) * nx_d] + _add_chunk_size_arrows( + ax1, + bottom_left_dst_chunk, + label_x=f"{bottom_left_dst_chunk.width} px", + label_y=f"{bottom_left_dst_chunk.height} px", + color=neutral_color, + fontsize=8.5, + ) + + _set_ax_extent_from_gdfs(ax1, [src_full_in_dst, src_blocks_in_dst, dst_full, dst_blocks]) + + # ------------------------------------------------------------------------- + # Annotation for highlighted destination chunk + # ------------------------------------------------------------------------- + n_intersections = len(src_intersect) + target = dst_blocks_one.geometry.iloc[0].centroid + tx, ty = target.x, target.y + + from matplotlib.offsetbox import AnnotationBbox, HPacker, VPacker, TextArea + + # Build colored text pieces + line1 = HPacker( + children=[ + TextArea("Destination chunk", textprops=dict(color=destination_color, fontsize=9)), + ], + align="center", + pad=0, + sep=2, + ) + + line2 = HPacker( + children=[ + TextArea("requires", textprops=dict(color="black", fontsize=9)), + ], + align="center", + pad=0, + sep=2, + ) + line3 = HPacker( + children=[ + TextArea(f"{n_intersections}", textprops=dict(color=intersect_color, fontsize=9, fontweight="bold")), + TextArea(" source chunks", textprops=dict(color=intersect_color, fontsize=9)), + ], + align="center", + pad=0, + sep=1, + ) + + label_box = VPacker( + children=[line1, line2, line3], + align="center", + pad=0, + sep=2, + ) + + ann = AnnotationBbox( + label_box, + (tx, ty), + xycoords="data", + xybox=(0.2, -0.5), + boxcoords="axes fraction", + arrowprops=dict( + arrowstyle="-|>", + color="black", + lw=arrow_lw, + shrinkA=0, + shrinkB=1, + mutation_scale=5, + connectionstyle="arc3,rad=0", + ), + bboxprops=dict( + boxstyle="round,pad=0.3,rounding_size=0.15", + fc="white", + ec="0.85", + lw=0.8, + alpha=0.96, + ), + ) + + ax1.add_artist(ann) + + # Ensure annotation is above the chunk polygons + ann.set_zorder(20) + + if ann.arrow_patch is not None: + ann.arrow_patch.set_zorder(21) + ann.arrow_patch.set_clip_on(False) + + if ann.patch is not None: + ann.patch.set_zorder(20) + + # ------------------------------------------------------------------------- + # Panel labels + # ------------------------------------------------------------------------- + ax1.text( + 0.5, + 1.01, + "Destination raster", + transform=ax1.transAxes, + ha="center", + va="bottom", + fontsize=12.5, + fontweight="bold", + color=destination_color, + ) + + ax1.text( + 0.5, + -0.02, + "Projected source raster", + transform=ax1.transAxes, + ha="center", + va="top", + fontsize=10, + color=source_color, + style="italic", + alpha=0.95, + ) + + # ------------------------------------------------------------------------- + # Grouped title + underline + reprojection arrow + # ------------------------------------------------------------------------- + src_crs = source_chunked_grid.grid.crs.to_epsg() + dst_crs = destination_chunked_grid.grid.crs.to_epsg() + + # Finalize layout before placing figure-level decorations + fig.canvas.draw() + + pos_left = ax0.get_position() + pos_right = ax1.get_position() + + x_center = 0.5 * (pos_left.x0 + pos_right.x1) + y_text = pos_left.y1 + 0.085 + y_line = pos_left.y1 + 0.077 + + # Grey grouped title + fig.text( + x_center, + y_text, + "Chunked reprojection", + ha="center", + va="bottom", + fontsize=13, + fontweight="semibold", + color="0.35", + ) + + # Underline spanning both panels + line = LineCollection( + [np.array([[pos_left.x0, y_line], [pos_right.x1, y_line]])], + colors="0.6", + linewidths=1.0, + alpha=0.9, + transform=fig.transFigure, + zorder=100, + ) + fig.add_artist(line) + + # Curved arrow and label between panels + label = f"EPSG:{src_crs} → EPSG:{dst_crs}\nResolution × 2" + + arrow = FancyArrowPatch( + (0.455, y_line - 0.16), + (0.545, y_line - 0.16), + transform=fig.transFigure, + arrowstyle="-|>", + connectionstyle="arc3,rad=-0.35", + linewidth=1.7, + color=neutral_color, + mutation_scale=20, + ) + fig.add_artist(arrow) + + fig.text( + 0.50, + y_line - 0.11, + label, + ha="center", + va="bottom", + fontsize=10, + color=neutral_color, + ) + + bbox_info = dict( + boxstyle="round,pad=0.35,rounding_size=0.2", + fc="white", + ec="0.85", + lw=0.8, + alpha=0.95, + ) + + ax1.text( + 0.7, + 1.65, + "Chunk size adapts " + "\nto resolution change\n" + "4×4 px → 2×2 px", + transform=ax1.transAxes, + ha="center", + va="top", + fontsize=9, + color="0.35", + bbox=bbox_info, + zorder=10, + ) + # ------------------------------------------------------------------------- + # Final layout polish + # ------------------------------------------------------------------------- + fig.subplots_adjust( + left=0, + right=1, + top=0.88, + bottom=0.03, + wspace=0.03, + ) + + return fig, axes + +# ----------------------------------------------------------------------------- +# Main figure +# ----------------------------------------------------------------------------- + +# Source grid in a projected CRS, with a fairly large extent to accentuate +# visible deformation after reprojection to geographic coordinates. +src_grid = GeoGrid( + transform=rio.transform.from_origin( + west=350_000, + north=8_300_000, + xsize=8_000, + ysize=8_000, + ), + shape=(20, 20), + crs=rio.crs.CRS.from_epsg(32633), # UTM 33N +) + +# Regular chunking for a cleaner schematic +# Use chunk sizes divisible by 2 so the destination chunking can be scaled +# consistently with the resolution change. +src_chunks = ( + (4, 4, 4, 4, 4), + (4, 4, 4, 4, 4), +) +src_chunked = ChunkedGeoGrid(grid=src_grid, chunks=src_chunks) + +# Destination grid: same extent, but 2x coarser resolution in both dimensions +# Since the source is 20x20 pixels, the destination becomes 10x10 pixels. +dst_grid = build_demo_destination_grid( + source_grid=src_grid, + dst_crs=rio.crs.CRS.from_epsg(4326), + dst_shape=(10, 10), +) + +# Scale chunk sizes with the resolution change so chunk footprints remain +# approximately comparable between source and destination. +# Source chunk size of 4x4 so we'll have a destination chunk size of 2x2 px +dst_chunks = ( + (2, 2, 2, 2, 2), + (2, 2, 2, 2, 2), +) +dst_chunked = ChunkedGeoGrid(grid=dst_grid, chunks=dst_chunks) + +fig, _ = plot_reprojection_chunk_diagram( + source_chunked_grid=src_chunked, + destination_chunked_grid=dst_chunked, + highlight_destination_chunk=(1, 1), +) +plt.show() diff --git a/doc/source/code/diagram_chunked_subsample.py b/doc/source/code/diagram_chunked_subsample.py new file mode 100644 index 000000000..56847a034 --- /dev/null +++ b/doc/source/code/diagram_chunked_subsample.py @@ -0,0 +1,749 @@ +"""Script to make diagram for chunked polygonize in documentation.""" +from __future__ import annotations + +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.collections import LineCollection +from matplotlib.patches import Rectangle +from geoutils.raster.array import get_mask_from_array +from geoutils.stats.sampling import ( + _get_subsample_size_from_user_input, + _splitmix64, + _subsample_numpy, +) + + +# ----------------------------------------------------------------------------- +# Example data +# ----------------------------------------------------------------------------- + +ARR = np.array( + [ + [1.2, 2.1, 2.8, 3.4, 4.0, np.nan, np.nan, 1.1], + [0.9, 1.5, 2.7, 3.1, 4.2, np.nan, np.nan, 0.6], + [0.8, 1.4, 2.2, 2.9, 3.9, 4.1, np.nan, 1.2], + [0.7, 1.0, np.nan, np.nan, 3.0, 3.4, 2.0, 1.5], + [1.6, 1.9, np.nan, np.nan, 2.6, 2.9, 1.4, 1.0], + [2.2, 1.1, 1.7, 4.2, 3.8, 3.1, 2.1, 1.7], + ], + dtype=np.float32, +) + +NROWS, NCOLS = ARR.shape +CHUNK_ROWS = (3, 3) +CHUNK_COLS = (4, 4) + +USER_SUBSAMPLE = 0.25 +SEED = 7 + +DISPLAY_K = 5 +EXAMPLE_CHUNK = (0, 1) # Top-right chunk + +# ----------------------------------------------------------------------------- +# Styling +# ----------------------------------------------------------------------------- + +NEUTRAL = "#333333" +PIXEL_GRID = "0.80" +CHUNK_COLOR = "#222222" +VALID_COLOR = "#6BAED6" +VALID_EDGE = "#4C78A8" +INVALID_COLOR = "#F7F7F7" +HIGHLIGHT_COLOR = "#F58518" +TOPK_COLOR = "#4C78A8" +SEQUENTIAL_COLOR = "#F58518" + + +# ----------------------------------------------------------------------------- +# Helpers for subsampling and valid values per chunk +# ----------------------------------------------------------------------------- + +def _first_valid_gids_in_chunk(arr: np.ndarray, chunk_loc: tuple[int, int], n: int = 5) -> np.ndarray: + """Return the first n valid global linear indices in one chunk, in flattened local order.""" + r0, r1, c0, c1 = _chunk_bounds(*chunk_loc) + block = arr[r0:r1, c0:c1] + mask = ~get_mask_from_array(block) + + flat_local = np.flatnonzero(mask.ravel())[:n] + if flat_local.size == 0: + return np.array([], dtype=np.int64) + + rr = flat_local // block.shape[1] + r0 + cc = flat_local % block.shape[1] + c0 + return np.ravel_multi_index((rr, cc), arr.shape).astype(np.int64) + +def _chunk_bounds(iy: int, ix: int) -> tuple[int, int, int, int]: + """Return raster chunk bounds as (row0, row1, col0, col1).""" + row_starts = np.concatenate([[0], np.cumsum(CHUNK_ROWS)]) + col_starts = np.concatenate([[0], np.cumsum(CHUNK_COLS)]) + return row_starts[iy], row_starts[iy + 1], col_starts[ix], col_starts[ix + 1] + + +def _chunk_valid_counts(arr: np.ndarray) -> list[int]: + """Valid counts in row-major chunk order.""" + counts: list[int] = [] + for iy in range(len(CHUNK_ROWS)): + for ix in range(len(CHUNK_COLS)): + r0, r1, c0, c1 = _chunk_bounds(iy, ix) + block = arr[r0:r1, c0:c1] + counts.append(int(np.count_nonzero(~get_mask_from_array(block)))) + return counts + + +def _global_valid_indices(arr: np.ndarray) -> np.ndarray: + """Global linear indices of valid pixels.""" + return np.flatnonzero(~get_mask_from_array(arr).ravel()) + + +def _topk_choice_with_keys(arr: np.ndarray, k: int, seed: int) -> tuple[np.ndarray, np.ndarray]: + """Return top-k gids and keys using the package SplitMix64 implementation.""" + valids = _global_valid_indices(arr).astype(np.uint64) + keys = _splitmix64(np.uint64(seed) ^ valids) + sel = np.argpartition(keys, k - 1)[:k] + sel = sel[np.lexsort((valids[sel], keys[sel]))] + return valids[sel].astype(np.int64), keys[sel] + + +# ----------------------------------------------------------------------------- +# Plot helpers +# ----------------------------------------------------------------------------- + +def _setup_axis(ax: plt.Axes, *, xlim: tuple[float, float], ylim: tuple[float, float], equal: bool = True) -> None: + """Common axis formatting.""" + ax.set_xlim(*xlim) + ax.set_ylim(*ylim) + ax.set_aspect("equal" if equal else "auto") + ax.set_xticks([]) + ax.set_yticks([]) + for spine in ax.spines.values(): + spine.set_visible(False) + + +def _draw_raster_cells(ax: plt.Axes, arr: np.ndarray) -> None: + """Draw valid / invalid raster cells.""" + mask = get_mask_from_array(arr) + for i in range(arr.shape[0]): + for j in range(arr.shape[1]): + x = j + y = arr.shape[0] - 1 - i + is_valid = not bool(mask[i, j]) + + rect = Rectangle( + (x, y), + 1.0, + 1.0, + facecolor=VALID_COLOR if is_valid else INVALID_COLOR, + edgecolor="none", + alpha=0.25 if is_valid else 1.0, + zorder=1, + ) + ax.add_patch(rect) + + if is_valid: + ax.text( + x + 0.5, + y + 0.5, + f"{arr[i, j]:.1f}", + ha="center", + va="center", + fontsize=8.5, + color=NEUTRAL, + zorder=3, + ) + + +def _draw_pixel_grid(ax: plt.Axes, *, nrows: int, ncols: int, extend: float = 0.0) -> None: + """Draw thin pixel grid.""" + segments: list[np.ndarray] = [] + + for x in range(ncols + 1): + segments.append(np.array([[x, -extend], [x, nrows + extend]], dtype=float)) + for y in range(nrows + 1): + segments.append(np.array([[-extend, y], [ncols + extend, y]], dtype=float)) + + ax.add_collection( + LineCollection( + segments, + colors=PIXEL_GRID, + linewidths=0.8, + capstyle="round", + joinstyle="round", + zorder=2, + clip_on=False, + ) + ) + + +def _draw_chunk_boundaries(ax: plt.Axes, *, nrows: int, ncols: int, extend: float = 0.0) -> None: + """Draw chunk boundaries.""" + row_starts = np.concatenate([[0], np.cumsum(CHUNK_ROWS)]) + col_starts = np.concatenate([[0], np.cumsum(CHUNK_COLS)]) + + for x in col_starts: + ax.plot( + [x, x], + [-extend, nrows + extend], + color=CHUNK_COLOR, + linewidth=2.3, + zorder=4, + solid_capstyle="round", + clip_on=False, + ) + + for y_idx in row_starts: + y = nrows - y_idx + ax.plot( + [-extend, ncols + extend], + [y, y], + color=CHUNK_COLOR, + linewidth=2.3, + zorder=4, + solid_capstyle="round", + clip_on=False, + ) + + +def _draw_example_count_chunk(ax: plt.Axes, chunk_loc: tuple[int, int], *, color: str = HIGHLIGHT_COLOR) -> None: + """Highlight one example chunk used for counting valid pixels.""" + r0, r1, c0, c1 = _chunk_bounds(*chunk_loc) + + left = c0 + right = c1 + top = NROWS - r0 + bottom = NROWS - r1 + + rect = Rectangle( + (left, bottom), + right - left, + top - bottom, + fill=False, + edgecolor=color, + linewidth=4.0, + zorder=7, + joinstyle="round", + ) + ax.add_patch(rect) + + +def _draw_chunk_count_labels( + ax: plt.Axes, + counts: list[int], + *, + example_chunk: tuple[int, int] | None = None, +) -> None: + """Draw per-chunk valid counts, with one example chunk highlighted.""" + k = 0 + for iy in range(len(CHUNK_ROWS)): + for ix in range(len(CHUNK_COLS)): + r0, r1, c0, c1 = _chunk_bounds(iy, ix) + x = 0.5 * (c0 + c1) + y = NROWS - 0.5 * (r0 + r1) + + is_example = example_chunk is not None and (iy, ix) == example_chunk + + ax.text( + x, + y - 1.1, + f"n = {counts[k]}", + ha="center", + va="center", + fontsize=10, + color=HIGHLIGHT_COLOR if is_example else NEUTRAL, + fontweight="bold" if is_example else None, + bbox=dict( + boxstyle="round,pad=0.20", + fc="white", + ec=HIGHLIGHT_COLOR if is_example else "0.85", + lw=1.0 if is_example else 0.8, + alpha=0.95, + ), + zorder=6, + ) + k += 1 + +def _draw_topk_selected(ax: plt.Axes, arr: np.ndarray, chosen_gids: np.ndarray) -> None: + """Draw topk-selected pixels as orange dots.""" + for gid in chosen_gids: + i, j = np.unravel_index(int(gid), arr.shape) + x = j + 0.5 + y = arr.shape[0] - 1 - i + 0.5 + + ax.scatter( + x, + y, + s=60, + facecolor=HIGHLIGHT_COLOR, + edgecolor="white", + linewidth=1.0, + zorder=7, + ) + + +def _add_title(fig: plt.Figure, title: str) -> None: + """Global grey title.""" + fig.text( + 0.5, + 0.99, + title, + ha="center", + va="top", + fontsize=15, + fontweight="semibold", + color="0.35", + ) + + +def _add_visual_legend(ax: plt.Axes) -> None: + """Draw a compact single-column legend in the upper-right free space of the top panel.""" + x_sym0 = 0.95 + x_sym1 = 1.0 + x_txt = 1.03 + + y0 = 0.68 + dy = 0.1 + + entries = [ + ("pixel_grid", "Pixel grid"), + ("chunk_boundary", "Chunk boundary"), + ("valid", "Valid pixel"), + ("invalid", "Invalid / nodata"), + ("example", "Illustrated example pixel"), + ("selected", "Selected sample pixel"), + ] + + for i, (kind, label) in enumerate(entries): + y = y0 - i * dy + + if kind == "pixel_grid": + ax.plot( + [x_sym0, x_sym1], + [y, y], + transform=ax.transAxes, + color=PIXEL_GRID, + lw=0.8, + solid_capstyle="round", + clip_on=False, + ) + + elif kind == "chunk_boundary": + ax.plot( + [x_sym0, x_sym1], + [y, y], + transform=ax.transAxes, + color=CHUNK_COLOR, + lw=2.3, + solid_capstyle="round", + clip_on=False, + ) + + elif kind == "valid": + rect = Rectangle( + (x_sym0, y - 0.018), + x_sym1 - x_sym0, + 0.036, + transform=ax.transAxes, + facecolor=VALID_COLOR, + edgecolor=VALID_EDGE, + linewidth=1.0, + alpha=0.25, + clip_on=False, + ) + ax.add_patch(rect) + + elif kind == "invalid": + rect = Rectangle( + (x_sym0, y - 0.018), + x_sym1 - x_sym0, + 0.036, + transform=ax.transAxes, + facecolor=INVALID_COLOR, + edgecolor="0.75", + linewidth=1.0, + clip_on=False, + ) + ax.add_patch(rect) + + elif kind == "example": + ax.scatter( + [0.5 * (x_sym0 + x_sym1)], + [y], + transform=ax.transAxes, + s=45, + facecolor=TOPK_COLOR, + edgecolor="white", + linewidth=1.0, + zorder=20, + clip_on=False, + ) + + elif kind == "selected": + ax.scatter( + [0.5 * (x_sym0 + x_sym1)], + [y], + transform=ax.transAxes, + s=45, + facecolor=HIGHLIGHT_COLOR, + edgecolor="white", + linewidth=1.0, + zorder=20, + clip_on=False, + ) + + ax.text( + x_txt, + y, + label, + transform=ax.transAxes, + ha="left", + va="center", + fontsize=9.5, + color=NEUTRAL, + ) + +def _valid_rank_map(arr: np.ndarray) -> dict[int, int]: + """Map global linear index -> rank in flattened valid-value order.""" + valids = _global_valid_indices(arr).astype(np.int64) + return {int(gid): i for i, gid in enumerate(valids)} + +def _draw_sequential_selected(ax: plt.Axes, arr: np.ndarray, chosen_gids: np.ndarray) -> None: + """Draw sequentially selected pixels as orange dots.""" + for gid in chosen_gids: + i, j = np.unravel_index(int(gid), arr.shape) + x = j + 0.5 + y = arr.shape[0] - 1 - i + 0.5 + + ax.scatter( + x, + y, + s=60, + facecolor=HIGHLIGHT_COLOR, + edgecolor="white", + linewidth=1.0, + zorder=7, + ) + +def _draw_topk_example_labels( + ax: plt.Axes, + arr: np.ndarray, + example_gids: np.ndarray, + *, + seed: int, + selected_gids: np.ndarray, +) -> None: + """Label example pixels with deterministic topk keys using log10(key).""" + gids_u64 = example_gids.astype(np.uint64) + keys = _splitmix64(np.uint64(seed) ^ gids_u64) + selected_set = set(selected_gids.astype(np.int64).tolist()) + + for gid, key in zip(example_gids, keys, strict=True): + i, j = np.unravel_index(int(gid), arr.shape) + x = j + 0.5 + y = arr.shape[0] - 1 - i + 0.5 + + display_val = int(key >> 56) + is_selected = int(gid) in selected_set + color = HIGHLIGHT_COLOR if is_selected else TOPK_COLOR + + ax.text( + x, + y + 0.48, + f"k={display_val}", + ha="center", + va="bottom", + fontsize=7.6, + color=color, + bbox=dict( + boxstyle="round,pad=0.15", + fc="white", + ec="0.85", + lw=0.8, + alpha=0.96, + ), + zorder=8, + ) + +def _sequential_chunk_order_rank_map(arr: np.ndarray) -> dict[int, int]: + """ + Map global linear index to rank in flattened valid order by chunk. + + Ordering is: + - chunks in row-major order + - within each chunk, valid pixels in flattened local order + """ + rank_map: dict[int, int] = {} + rank = 1 + + for iy in range(len(CHUNK_ROWS)): + for ix in range(len(CHUNK_COLS)): + r0, r1, c0, c1 = _chunk_bounds(iy, ix) + block = arr[r0:r1, c0:c1] + mask = ~get_mask_from_array(block) + + flat_local = np.flatnonzero(mask.ravel()) + ncols = block.shape[1] + + for flat in flat_local: + rr = r0 + flat // ncols + cc = c0 + flat % ncols + gid = int(np.ravel_multi_index((rr, cc), arr.shape)) + rank_map[gid] = rank + rank += 1 + + return rank_map + +def _draw_sequential_example_labels( + ax: plt.Axes, + arr: np.ndarray, + example_gids: np.ndarray, + *, + selected_gids: np.ndarray, +) -> None: + """Label example pixels with their rank in chunk-wise flattened valid order.""" + rank_map = _sequential_chunk_order_rank_map(arr) + selected_set = set(selected_gids.astype(np.int64).tolist()) + + for gid in example_gids: + i, j = np.unravel_index(int(gid), arr.shape) + x = j + 0.5 + y = arr.shape[0] - 1 - i + 0.5 + + is_selected = int(gid) in selected_set + color = HIGHLIGHT_COLOR if is_selected else TOPK_COLOR + + ax.text( + x, + y + 0.48, + f"n={rank_map[int(gid)]}", + ha="center", + va="bottom", + fontsize=8.0, + color=color, + bbox=dict( + boxstyle="round,pad=0.15", + fc="white", + ec="0.85", + lw=0.8, + alpha=0.96, + ), + zorder=8, + ) + +def _example_valid_gids_sparse(arr: np.ndarray) -> np.ndarray: + """ + Return example valid gids through the first row and into the middle of the + second row, keeping one valid pixel out of two for readability. + """ + valids = _global_valid_indices(arr).astype(np.int64) + + row_break_1 = CHUNK_ROWS[0] * NCOLS + row_break_2 = (CHUNK_ROWS[0] + CHUNK_ROWS[1] // 2) * NCOLS + + gids = valids[valids < row_break_2] + return gids[::2] + +def _draw_example_reference_dots( + ax: plt.Axes, + arr: np.ndarray, + example_gids: np.ndarray, + *, + selected_gids: np.ndarray, +) -> None: + """Draw blue dots for illustrated example pixels that are not selected.""" + selected_set = set(selected_gids.astype(np.int64).tolist()) + + for gid in example_gids: + if int(gid) in selected_set: + continue + + i, j = np.unravel_index(int(gid), arr.shape) + x = j + 0.5 + y = arr.shape[0] - 1 - i + 0.5 + + ax.scatter( + x, + y, + s=60, + facecolor=TOPK_COLOR, + edgecolor="white", + linewidth=1.0, + zorder=6.5, + ) + +def _lowest_key_example_gids(example_gids: np.ndarray, *, seed: int, n_low: int = 2) -> np.ndarray: + """Return gids of the n_low smallest keys among example gids.""" + gids_u64 = example_gids.astype(np.uint64) + keys = _splitmix64(np.uint64(seed) ^ gids_u64) + + m = min(n_low, len(example_gids)) + sel = np.argpartition(keys, m - 1)[:m] + return example_gids[sel].astype(np.int64) + +# ----------------------------------------------------------------------------- +# Main figure +# ----------------------------------------------------------------------------- + +def make_chunked_subsample_diagram() -> tuple[plt.Figure, np.ndarray]: + """Build a 3-panel schematic for chunked subsampling.""" + + EXAMPLE_GIDS = _example_valid_gids_sparse(ARR) + LOW_KEY_EXAMPLE_GIDS = _lowest_key_example_gids(EXAMPLE_GIDS, seed=SEED, n_low=2) + + fig = plt.figure(figsize=(8, 6)) + gs = fig.add_gridspec(2, 2, height_ratios=[1.0, 1.0], hspace=0.28, wspace=0.28) + + ax_top = fig.add_subplot(gs[0, :]) + ax_bl = fig.add_subplot(gs[1, 0]) + ax_br = fig.add_subplot(gs[1, 1]) + + axes = np.array([ax_top, ax_bl, ax_br], dtype=object) + + valid_counts = _chunk_valid_counts(ARR) + total_valids = int(np.count_nonzero(~get_mask_from_array(ARR))) + k = _get_subsample_size_from_user_input(USER_SUBSAMPLE, total_valids) + + # Use package subsampling function for both strategies + sequential_gids_rc = _subsample_numpy( + ARR, + subsample=USER_SUBSAMPLE, + return_indices=True, + random_state=SEED, + strategy="sequential", + ) + sequential_gids = np.ravel_multi_index(sequential_gids_rc, ARR.shape).astype(np.int64) + + topk_gids_rc = _subsample_numpy( + ARR, + subsample=USER_SUBSAMPLE, + return_indices=True, + random_state=SEED, + strategy="topk", + ) + topk_gids = np.ravel_multi_index(topk_gids_rc, ARR.shape).astype(np.int64) + + # Keys only for labeling the topk panel + topk_gids_for_keys, topk_keys = _topk_choice_with_keys(ARR, k, SEED) + + r0, r1, c0, c1 = _chunk_bounds(*EXAMPLE_CHUNK) + example_count = int(np.count_nonzero(~get_mask_from_array(ARR[r0:r1, c0:c1]))) + + # Top panel + _draw_raster_cells(ax_top, ARR) + _draw_pixel_grid(ax_top, nrows=NROWS, ncols=NCOLS) + _draw_chunk_boundaries(ax_top, nrows=NROWS, ncols=NCOLS) + _draw_chunk_count_labels(ax_top, valid_counts, example_chunk=EXAMPLE_CHUNK) + _draw_example_count_chunk(ax_top, EXAMPLE_CHUNK) + + ax_top.text( + 0.5, + 1.0, + "Count valid values per chunk", + transform=ax_top.transAxes, + ha="center", + va="bottom", + fontsize=13, + fontweight="bold", + color=NEUTRAL, + ) + + ax_top.text( + -0.36, + 0.50, + "We sum per-chunk\nvalid values to get\nrequested subsample size\nfrom total valid count\n" + "(as it may be a fraction\n or exceed total count)", + transform=ax_top.transAxes, + ha="center", + va="center", + fontsize=10.2, + color=NEUTRAL, + linespacing=1.3, + ) + + _add_visual_legend(ax_top) + _setup_axis(ax_top, xlim=(-0.5, NCOLS + 2.8), ylim=(-0.5, NROWS + 0.5)) + + # Bottom left: topk + _draw_raster_cells(ax_bl, ARR) + _draw_pixel_grid(ax_bl, nrows=NROWS, ncols=NCOLS) + _draw_chunk_boundaries(ax_bl, nrows=NROWS, ncols=NCOLS) + _draw_example_count_chunk(ax_bl, (0, 0), color="0.55") + + _draw_example_reference_dots(ax_bl, ARR, EXAMPLE_GIDS, selected_gids=topk_gids) + _draw_topk_selected(ax_bl, ARR, topk_gids) + _draw_topk_example_labels( + ax_bl, + ARR, + EXAMPLE_GIDS, + seed=SEED, + selected_gids=LOW_KEY_EXAMPLE_GIDS + ) + + ax_bl.text( + 0.5, + 1.0, + "topk (chunk-invariant)", + transform=ax_bl.transAxes, + ha="center", + va="bottom", + fontsize=13, + fontweight="bold", + color=NEUTRAL, + ) + ax_bl.text( + 0.5, + -0.02, + "Each valid pixel gets a deterministic uint64 key\nfrom its row/col indexes (rounded rank shown)\n" + "the k smallest keys are selected", + transform=ax_bl.transAxes, + ha="center", + va="top", + fontsize=10.0, + color=NEUTRAL, + linespacing=1.25, + ) + _setup_axis(ax_bl, xlim=(-0.5, NCOLS + 0.5), ylim=(-0.5, NROWS + 0.5)) + + # Bottom right: sequential + _draw_raster_cells(ax_br, ARR) + _draw_pixel_grid(ax_br, nrows=NROWS, ncols=NCOLS) + _draw_chunk_boundaries(ax_br, nrows=NROWS, ncols=NCOLS) + _draw_example_count_chunk(ax_br, (0, 0), color="0.55") + + _draw_example_reference_dots(ax_br, ARR, EXAMPLE_GIDS, selected_gids=sequential_gids) + _draw_sequential_selected(ax_br, ARR, sequential_gids) + _draw_sequential_example_labels(ax_br, ARR, EXAMPLE_GIDS, selected_gids=sequential_gids) + + ax_br.text( + 0.5, + 1.0, + "sequential (chunk-dependent)", + transform=ax_br.transAxes, + ha="center", + va="bottom", + fontsize=13, + fontweight="bold", + color=NEUTRAL, + ) + + ax_br.text( + 0.5, + -0.02, + "A random draw is applied to flattened valid indexes,\n" + "fast selection but depends on chunking", + transform=ax_br.transAxes, + ha="center", + va="top", + fontsize=10.0, + color=NEUTRAL, + linespacing=1.25, + ) + + _setup_axis(ax_br, xlim=(-0.5, NCOLS + 0.5), ylim=(-0.5, NROWS + 0.5)) + + _add_title(fig, "Chunked subsampling of valid raster values") + + fig.subplots_adjust(left=0.02, right=0.98, top=0.90, bottom=0.13) + return fig, axes + + +fig, _ = make_chunked_subsample_diagram() +plt.show() diff --git a/doc/source/conf.py b/doc/source/conf.py index 152953dfe..133b0ba7c 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -69,6 +69,7 @@ "rioxarray": ("https://corteva.github.io/rioxarray/stable/", None), "pandas": ("https://pandas.pydata.org/docs/", None), "scipy": ("https://docs.scipy.org/doc/scipy/", None), + "dask": ("https://docs.dask.org/en/stable/", None), } example_path = os.path.join("../", "../", "examples") diff --git a/doc/source/core_index.md b/doc/source/core_index.md index 888ed5b80..8b2532fa8 100644 --- a/doc/source/core_index.md +++ b/doc/source/core_index.md @@ -10,8 +10,6 @@ core_composition core_match_ref core_py_ops core_array_funcs -core_lazy_load core_parsing core_inheritance -core_scalability ``` diff --git a/doc/source/core_lazy_load.md b/doc/source/core_lazy_load.md deleted file mode 100644 index eaf9b3d69..000000000 --- a/doc/source/core_lazy_load.md +++ /dev/null @@ -1,86 +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 ---- -(core-lazy-load)= - -# Implicit lazy loading - -Lazy loading, also known as "call-by-need", is the delay in loading or evaluating a dataset. - -In GeoUtils, we implicitly load and pass only metadata until the data is actually needed, and are working to implement lazy analysis tools relying on other packages. - -## Lazy instantiation of {class}`Rasters` - -By default, GeoUtils instantiate a {class}`~geoutils.Raster` from an **on-disk** file without loading its {attr}`geoutils.Raster.data` array. It only loads its -metadata ({attr}`~geoutils.Raster.transform`, {attr}`~geoutils.Raster.crs`, {attr}`~geoutils.Raster.nodata` and derivatives, as well as -{attr}`~geoutils.Raster.name` and {attr}`~geoutils.Raster.driver`). - -```{code-cell} ipython3 - -import geoutils as gu - -# Instantiate a raster from a filename on disk -filename_rast = gu.examples.get_path("everest_landsat_b4") -rast = gu.Raster(filename_rast) - -# This raster is not loaded -rast -``` - -To load the data explicitly during instantiation opening, `load_data=True` can be passed to {class}`~geoutils.Raster`. Or the {func}`~geoutils.Raster.load` -method can be called after. The two are equivalent. - -```{code-cell} ipython3 -# Initiate another raster just for the purpose of loading -rast_to_load = gu.Raster(gu.examples.get_path("everest_landsat_b4")) -rast_to_load.load() - -# This raster is loaded -rast_to_load -``` - -## Lazy passing of georeferencing metadata - -Operations relying on georeferencing metadata of {class}`Rasters` or {class}`Vectors` are always done by respecting the -possible lazy loading of the objects. - -For instance, using any {class}`~geoutils.Raster` or {class}`~geoutils.Vector` as a match-reference for a geospatial operation (see {ref}`core-match-ref`) will -always conserve the lazy loading of that match-reference object. - -```{code-cell} ipython3 ---- -mystnb: - output_stderr: show ---- - -# Use a smaller Raster as reference to crop the initial one -smaller_rast = gu.Raster(gu.examples.get_path("everest_landsat_b4_cropped")) -rast.crop(smaller_rast) - -# The reference raster is not loaded -smaller_rast -``` - -## Optimized geospatial subsetting - -```{important} -These features are a work in progress, we aim to make GeoUtils more lazy-friendly through [Dask](https://docs.dask.org/en/stable/) in future versions of the -package! -``` - -Some georeferencing operations can be done without loading the entire array. Right now, relying directly on Rasterio, GeoUtils supports optimized subsetting -through the {func}`~geoutils.Raster.crop` method. - -```{code-cell} ipython3 -# The previously cropped Raster was loaded without accessing the entire array -rast -``` diff --git a/doc/source/ecosystem.md b/doc/source/ecosystem.md index 896c5b43e..1df98e0e8 100644 --- a/doc/source/ecosystem.md +++ b/doc/source/ecosystem.md @@ -6,7 +6,7 @@ It extends commonly used tools and works alongside many other for geospatial dat See the {ref}`accessors` page for details on GeoUtils' accessors. -```{see-also} +```{seealso} **[xDEM](https://github.com/GlacioHack/xdem)** is the sister package of GeoUtils, focused on the **analysis of digital elevation models (DEMs) and elevation point clouds**, including terrain attributes, coregistration and uncertainty propagation. ``` diff --git a/doc/source/feature_overview.md b/doc/source/feature_overview.md index 9f542c299..c17f365d1 100644 --- a/doc/source/feature_overview.md +++ b/doc/source/feature_overview.md @@ -1,272 +1,453 @@ ---- -file_format: mystnb -jupytext: - formats: md:myst - text_representation: - extension: .md - format_name: myst -kernelspec: - display_name: geoutils-env - language: python - name: geoutils ---- (feature-overview)= +# Feature and scalability overview -# Feature overview +GeoUtils provides a unified API for manipulating **raster**, **vector**, and **point-cloud** data, and provides **scalable CPU execution** for most raster operations through Dask and Multiprocessing. -The following presents a descriptive example show-casing all core features of GeoUtils. +As many of our numerical operations rely on **NumPy, SciPy or Numba**, those are planned to be linked to their **GPU** counterparts (**CuPy** and **Numba CUDA**) in the future. -```{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 😅) +The **{ref}`summary tables` directly below** lists the core features of GeoUtils, their scalability and available backends. +Further below, a series of **{ref}`illustrated examples`** demonstrate these features. -Alternatively, start your own notebook to test GeoUtils at [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/GlacioHack/geoutils/main). +```{seealso} +If you are interested in porting from GDAL/OGR, see our {ref}`cheatsheet-osgeo` page. +While tables below provide a scalability summary, the detailed **input/output behaviour of all operations** is available on the {ref}`scalability-support` page. ``` -```{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 +## Summary + +GeoUtils exposes a **consistent API across raster, vector and point-cloud objects** where possible (similar in spirit to the recent [GDAL CLI overhaul](https://gdal.org/en/stable/programs/index.html)). Many operations also support convenient **match-reference arguments** (e.g., matching a grid for reprojection or rasterization, bounds for cropping, or point coordinates for interpolation). See the {ref}`core-match-ref` page for details. + +At its core, GeoUtils provides two interchangeable ways to work with geospatial data, exposing **identical APIs**: + +- **Accessors** that extend existing data structures ({class}`rst ` for **rasters** with **Xarray**, `pc` and `vct` for **point clouds** and **vectors** with **GeoPandas**), +- **GeoUtils objects** {class}`~geoutils.Raster`, {class}`~geoutils.PointCloud`, {class}`~geoutils.Vector`. + +Nearly all **raster operations** support **scalable execution** using [Dask](https://www.dask.org/) or Multiprocessing, allowing large datasets to be processed **chunk-by-chunk without loading the full array into memory**. Support for **point-cloud operations** is partial and ongoing, while **vector operations** may gain scalable support in the future. + +Additionally, some numerical routines of GeoUtils provide multiple computational **backends** (e.g., SciPy or Numba implementations). + +All methods are tested to ensure they produce **identical results** whether executed **in-memory**, **chunked**, or with **different computational backends**. +(tables-overview)= +## Data operations + +We first describe GeoUtils' core **data operations**, which operate on underlying arrays or geometries and can therefore benefit from **scalable execution**. + +**Legend:** **“/”** indicates methods **shared across object types**, while **“⟷”** indicates methods **interfacing between two object types**. + +```{list-table} Common API for data operations +:widths: 3 5 1 2 +:header-rows: 1 +:align: left +:class: tight-table + +* - Method + - Notes + - Scalable + - Backend + +* - Raster / Vector / Point + - + - + - + +* - {meth}`~geoutils.Raster.reproject()` + - Reproject to other CRS. Default tolerance parameters ensure chunk-invariance. + - ✅ + - Rasterio / PyProj + +* - {meth}`~geoutils.Raster.crop()` + - Crop to bounds, either intersecting (untouched) allowing efficient I/O, or clipped (data modified). + - ✅ + - Rasterio / GeoPandas + +* - {meth}`~geoutils.Raster.translate()` + - Apply a grid shift to object. + - ✅ + - NumPy / GeoPandas + +* - {meth}`~geoutils.Raster.plot()` + - Visualization helper. + - ❌ + - Matplotlib + +* - Raster / Point + - + - + - + +* - {meth}`~geoutils.Vector.create_mask()` + - Create boolean mask of a vector geometries over raster or point. + - ✅ + - Rasterio / GeoPandas + +* - {meth}`~geoutils.Raster.get_stats()` + - Compute statistics of valid values over a valid mask. + - ❌ + - NumPy / SciPy + +* - {meth}`~geoutils.Raster.subsample()` + - Randomly sample valid values. Chunk-invariant seed ensures reproducibility. + - ✅ + - NumPy + +* - {meth}`~geoutils.Raster.filter()` + - Filter over window. Fast vectorized logic with NaN support. + - ✅ + - SciPy + +* - {meth}`~geoutils.Raster.proximity()` + - Estimate proximity distance to target values. + - ❌ + - SciPy + + +* - Raster ⟷ Vector + - + - + - + +* - {meth}`~geoutils.Raster.polygonize()` + - Convert raster regions to vector polygons. Multiple chunked strategies for performance. + - ✅ + - Rasterio / GeoPandas + +* - {meth}`~geoutils.Vector.rasterize()` + - Burn vector geometries onto a raster grid. + - ✅ + - Rasterio + +* - Raster ⟷ Point + - + - + - + +* - {meth}`~geoutils.Raster.interp_points()` + - Interpolate raster at point locations. Fast regular-grid logic with added NaN propagation. + - ✅ + - SciPy + +* - {meth}`~geoutils.Raster.reduce_points()` + - Aggregate raster values around points. + - ❌ + - NumPy + +* - {meth}`~geoutils.PointCloud.grid()` + - Grid irregular points onto a raster grid. Multiple approaches with added NaN propagation. + - ❌ + - SciPy + +* - {meth}`~geoutils.Raster.from_pointcloud_regular()` + - Direct conversion when points lie on a regular grid. + - ❌ + - NumPy + +* - {meth}`~geoutils.Raster.to_pointcloud()` + - Conversion to point cloud. + - ❌ + - NumPy ``` -## The core {class}`~geoutils.Raster`, {class}`~geoutils.Vector` and {class}`~geoutils.PointCloud` objects +## Metadata properties and operations -In GeoUtils, geospatial handling is object-based and revolves around {class}`~geoutils.Raster` and {class}`~geoutils.Vector`. -These link to either **in-memory** or **on-disk** datasets, opened by calling the object from a filepath for the latter. +In addition to data operations, GeoUtils exposes **metadata** properties and methods consistently across geospatial objects. +These operate only on metadata and therefore **do not load or modify underlying data arrays**. -```{code-cell} ipython3 -import geoutils as gu +```{list-table} Common API from metadata operations +:widths: 3 7 +:header-rows: 1 +:align: left +:class: tight-table -# Examples files: infrared band of Landsat and glacier outlines -filename_rast = gu.examples.get_path("everest_landsat_b4") -filename_vect = gu.examples.get_path("everest_rgi_outlines") +* - Attribute / Method + - Description -# Open files by calling Raster and Vector -rast = gu.Raster(filename_rast) -vect = gu.Vector(filename_vect) -``` +* - Raster / Vector / Point + - -A {class}`~geoutils.Raster` is an object with four main attributes: a {class}`numpy.ma.MaskedArray` as {attr}`~geoutils.Raster.data`, a -{class}`pyproj.crs.CRS` as {attr}`~geoutils.Raster.crs`, an [{class}`affine.Affine`](https://rasterio.readthedocs.io/en/stable/topics/migrating-to-v1.html#affine-affine-vs-gdal-style-geotransforms) -as {attr}`~geoutils.Raster.transform`, and a {class}`float` or {class}`int` as {attr}`~geoutils.Raster.nodata`. +* - {attr}`~geoutils.Raster.crs` + - Coordinate reference system (CRS) of object. +* - {attr}`~geoutils.Raster.bounds` + - Bounding box of object. -```{code-cell} ipython3 -# The opened raster -rast -``` +* - {attr}`~geoutils.Raster.footprint` + - Footprint polygon geometry of object. + +* - {attr}`~geoutils.Raster.is_loaded` + - Whether geospatial object is loaded in-memory. -```{important} -When a file exists on disk, {class}`~geoutils.Raster` is linked to a {class}`rasterio.io.DatasetReader` object for loading the metadata. The array will be -**loaded in-memory implicitly** when {attr}`~geoutils.Raster.data` is required by an operation. +* - {attr}`~geoutils.Raster.name` + - Filename of object on disk, if it exists. + +* - {meth}`~geoutils.Raster.get_bounds_projected()` + - Bounds projected in other CRS. -See {ref}`core-lazy-load` for more details. -``` +* - {meth}`~geoutils.Raster.get_footprint_projected()` + - Footprint polygon geometry in other CRS. + +* - {meth}`~geoutils.Raster.get_metric_crs()` + - Get metric CRS suitable for this object. + +* - {meth}`~geoutils.Raster.info()` + - Summary of attributes for geospatial object. + +* - Raster / Point + - + +* - {attr}`~geoutils.Raster.data` + - Data array (2D grid for raster, 1D for point cloud). + +* - {attr}`~geoutils.Raster.shape` + - Shape of data array. + +* - {attr}`~geoutils.Raster.is_mask` + - Whether object is a mask. Clarifies ambiguity of raster/point file types often not supporting boolean types. + +* - Raster + - + +* - {attr}`~geoutils.Raster.transform` + - Geotransform to map raster indices to spatial coordinates. -A {class}`~geoutils.Vector` is an object with a single main attribute: a {class}`~geopandas.GeoDataFrame` as {attr}`~geoutils.Vector.ds`, for which -most methods are wrapped directly into {class}`~geoutils.Vector`. +* - {attr}`~geoutils.Raster.nodata` + - Nodata value used to represent missing data on disk. + +* - {attr}`~geoutils.Raster.area_or_point` + - Pixel interpretation of raster values, either center point or area average. -```{code-cell} ipython3 -# The opened vector -vect +* - Point + - + +* - {attr}`~geoutils.PointCloud.point_count` + - Number of points in the point cloud. ``` -All other attributes are derivatives of those main attributes, or of the filename on disk. Attributes of {class}`~geoutils.Raster` and -{class}`~geoutils.Vector` update with geospatial operations on themselves. +[//]: # ((examples-overview)=) +[//]: # (## Examples) -## Handling and match-reference +[//]: # () +[//]: # (The following presents a descriptive example show-casing core features of GeoUtils.) -In GeoUtils, geospatial handling operations are based on class methods, such as {func}`~geoutils.Raster.crop` or {func}`~geoutils.Raster.reproject`. +[//]: # () +[//]: # (```{code-cell} ipython3) -For convenience and consistency, nearly all of these methods can be passed solely another {class}`~geoutils.Raster` or {class}`~geoutils.Vector` as a -**reference to match** during the operation. A **reference {class}`~geoutils.Vector`** enforces a matching of {attr}`~geoutils.Vector.bounds` and/or -{attr}`~geoutils.Vector.crs`, while a **reference {class}`~geoutils.Raster`** can also enforce a matching of {attr}`~geoutils.Raster.res`, depending on the nature of the operation. +[//]: # (:tags: [remove-cell]) +[//]: # () +[//]: # (# To get a good resolution for displayed figures) -```{code-cell} ipython3 -# Print initial bounds of the vector -print(vect.bounds) -# Crop vector to raster's extent, and add clipping option (otherwise keeps all intersecting features) -vect_cropped = vect.crop(rast, clip=True) -# Print bounds of cropped + clipped vector -print(vect_cropped.bounds) -``` +[//]: # (from matplotlib import pyplot) -```{margin} -1The names of geospatial handling methods is largely based on [GDAL and OGR](https://gdal.org/)'s, with the notable exception of {func}`~geoutils.Vector.reproject` that better applies to vectors than `warp`. -``` +[//]: # (pyplot.rcParams['figure.dpi'] = 600) -Additionally, in GeoUtils, **methods that apply to the same georeferencing attributes have consistent naming**1 across {class}`~geoutils.Raster` and {class}`~geoutils.Vector`. - -A {func}`~geoutils.Raster.reproject` involves a change in {attr}`~geoutils.Raster.crs` or {attr}`~geoutils.Raster.transform`, while a {func}`~geoutils.Raster.crop` only involves a change -in {attr}`~geoutils.Raster.bounds`. Using {func}`~geoutils.Raster.polygonize` allows to generate a {class}`~geoutils.Vector` from a {class}`~geoutils.Raster`, -and the other way around for {func}`~geoutils.Vector.rasterize`. - -```{list-table} - :widths: 30 30 30 - :header-rows: 1 - - * - **Class** - - {class}`~geoutils.Raster` - - {class}`~geoutils.Vector` - * - Warping/Reprojection - - {func}`~geoutils.Raster.reproject` - - {func}`~geoutils.Vector.reproject` - * - Cropping/Clipping - - {func}`~geoutils.Raster.crop` - - {func}`~geoutils.Vector.crop` - * - Rasterize/Polygonize - - {func}`~geoutils.Raster.polygonize` - - {func}`~geoutils.Vector.rasterize` -``` +[//]: # (pyplot.rcParams['savefig.dpi'] = 600) -All methods can also be passed any number of georeferencing arguments such as {attr}`~geoutils.Raster.shape` or {attr}`~geoutils.Raster.res`, and will -naturally deduce others from the input {class}`~geoutils.Raster` or {class}`~geoutils.Vector`, much as in [GDAL](https://gdal.org/)'s command line. +[//]: # (pyplot.rcParams['font.size'] = 9) +[//]: # (```) -## Higher-level analysis tools +[//]: # () +[//]: # (### Core objects and accessors) -GeoUtils also implements higher-level geospatial analysis tools for both {class}`Rasters` and {class}`Vectors`. For -example, one can compute the distance to a {class}`~geoutils.Vector` geometry, or to target pixels of a {class}`~geoutils.Raster`, using -{func}`~geoutils.Vector.proximity`. +[//]: # () +[//]: # (GeoUtils operates on **rasters, vectors and point clouds**.) -As with the geospatial handling functions previously listed, many analysis functions can take a {class}`~geoutils.Raster` or {class}`~geoutils.Vector` as a -**reference to utilize** during the operation. In the case of {func}`~geoutils.Vector.proximity`, passing a {class}`~geoutils.Raster` serves as a reference -for the georeferenced grid on which to compute the distances. +[//]: # () +[//]: # (Below, we demonstrate features through **accessors**: {class}`rst ` for **Xarray** rasters, `vct` for **GeoPandas** vectors, `pc` for **GeoPandas** point clouds.) -```{code-cell} ipython3 -# Compute proximity to vector on raster's grid -rast_proximity_to_vec = vect.proximity(rast) -``` +[//]: # (The exact same operations are available through the {class}`~geoutils.Raster`, {class}`~geoutils.Vector` and {class}`~geoutils.PointCloud` objects.) -```{note} -Right now, the array {attr}`~geoutils.Raster.data` of `rast` is still not loaded. Applying {func}`~geoutils.Raster.crop` does not yet require loading, -and `rast`'s metadata is sufficient to provide a georeferenced grid for {func}`~geoutils.Vector.proximity`. The array will only be loaded when necessary. -``` +[//]: # () +[//]: # (**Accessors** operations always return another **Xarray** or **GeoPandas** object, while **GeoUtils object** operations always return a **GeoUtils object**.) -## Quick plotting +[//]: # (However, you can pass **any object** as an argument to a function.) -To facilitate the analysis process, GeoUtils includes quick plotting tools that support multiple colorbars and implicitly add layers to the current axis. -Those are build on top of {func}`rasterio.plot.show` and {func}`geopandas.GeoDataFrame.plot`, and relay any argument passed. +[//]: # () +[//]: # (First, we open datasets with {meth}`~geoutils.open_raster` and {meth}`~geoutils.open_vector`, passing a `chunks` argument to trigger out-of-memory **Dask** behaviour for the raster.) -```{seealso} -GeoUtils' plotting tools only aim to smooth out the most common hassles when quickly plotting raster and vectors. +[//]: # () +[//]: # (```{code-cell} ipython3) -For advanced plotting tools to create "publication-quality" figures, see [Cartopy](https://scitools.org.uk/cartopy/docs/latest/) or -[GeoPlot](https://residentmario.github.io/geoplot/index.html). -``` +[//]: # (import geoutils as gu) -The plotting functionality is named {func}`~geoutils.Raster.plot` everywhere, for consistency. Here again, a {class}`~geoutils.Raster` or -{class}`~geoutils.Vector` can be passed as a **reference to match** to ensure all data is displayed on the same grid and projection. +[//]: # (import geopandas as gpd) -```{code-cell} ipython3 -# Plot proximity to vector -rast_proximity_to_vec = vect.proximity(rast) -rast_proximity_to_vec.plot(cbar_title="Distance to glacier outline") -vect.plot(rast_proximity_to_vec, fc="none") -``` +[//]: # () +[//]: # (# Example files: infrared band of Landsat and glacier outlines) -```{tip} -To quickly visualize a raster directly from a terminal, without opening a Python console/notebook, check out our tool `geoviewer.py` in the {ref}`cli` documentation. -``` +[//]: # (filename_rast = gu.examples.get_path("everest_landsat_b4")) -## Pythonic arithmetic and NumPy interface +[//]: # (filename_vect = gu.examples.get_path("everest_rgi_outlines")) -All {class}`~geoutils.Raster` objects support Python arithmetic ({func}`+`, {func}`-`, {func}`/`, {func}`//`, {func}`*`, -{func}`**`, {func}`%`) with any other {class}`~geoutils.Raster`, {class}`~numpy.ndarray` or -number. With another {class}`~geoutils.Raster`, the georeferencing must match, while only the shape with a {class}`~numpy.ndarray`. +[//]: # () +[//]: # (# Open datasets) -```{code-cell} ipython3 -# Add 1 to the raster array -rast += 1 -``` +[//]: # (ds = gu.open_raster(filename_rast, chunks={"x": 200, "y": 200})) -Additionally, the {class}`~geoutils.Raster` object possesses a NumPy masked-array interface that allows to apply to it any [NumPy universal function](https://numpy.org/doc/stable/reference/ufuncs.html) and -most other NumPy array functions, while logically casting {class}`dtype` and respecting {attr}`~geoutils.Raster.nodata` values. +[//]: # (gdf = gpd.read_file(filename_vect)) -```{code-cell} ipython3 -# Apply a normalization to the raster -import numpy as np -rast = (rast - np.min(rast)) / (np.max(rast) - np.min(rast)) -``` +[//]: # () +[//]: # (# Raster and vector objects) -## Casting to raster mask, indexing and overload +[//]: # (ds) -All {class}`~geoutils.Raster` objects also support Python logical comparison operators ({func}`==`, {func}` != `, {func}`>=`, {func}`>`, {func}`<=`, -{func}`<`), or more complex NumPy logical functions. Those operations automatically casts them into a raster mask, i.e. a boolean -{class}`~geoutils.Raster`. +[//]: # (gdf) -```{code-cell} ipython3 -# Get mask of an AOI: infrared index above 0.6, at least 200 m from glaciers -mask_aoi = np.logical_and(rast > 0.6, rast_proximity_to_vec > 200) -``` +[//]: # (```) -Raster masks can then be used for indexing a {class}`~geoutils.Raster`, which returns a {class}`~numpy.ma.MaskedArray` of indexed values. +[//]: # () +[//]: # (GeoUtils accessors expose spatial methods directly on these objects while keeping compatibility with the broader **NumPy, Xarray and GeoPandas ecosystem**.) -```{code-cell} ipython3 -# Index raster with mask to extract a 1-D array -values_aoi = rast[mask_aoi] -``` +[//]: # () +[//]: # (### Plotting) -Raster masks also have simplified {class}`~geoutils.Raster` methods due to their boolean {class}`dtype` rendering many arguments implicit. -For instance, using {func}`~geoutils.Raster.polygonize` with a raster mask is straightforward, to retrieve a {class}`~geoutils.Vector` of the area-of-interest: +[//]: # () +[//]: # (For plotting, GeoUtils includes lightweight plotting helpers to simplify common tasks such as overlaying rasters and vectors with multiple colorbars, while passing through standard plotting arguments.) -```{code-cell} ipython3 -# Polygonize areas where mask is True -vect_aoi = mask_aoi.polygonize() -``` +[//]: # () +[//]: # (```{code-cell} ipython3) -```{code-cell} ipython3 -# Plot result -rast.plot(cmap='Reds', cbar_title='Normalized infrared') -vect_aoi.plot(fc='none', ec='k', lw=0.75) -``` +[//]: # (# Plot raster and vector) -## Saving to file +[//]: # (ds.rst.plot(cbar_title="Distance to glacier outline")) -Finally, for saving a {class}`~geoutils.Raster` or {class}`~geoutils.Vector` to file, simply call the {func}`~geoutils.Raster.to_file` function. +[//]: # (gdf.plot(rast_proximity_to_vec, fc="none")) -```{code-cell} ipython3 -# Save our AOI vector -vect_aoi.to_file("myaoi.gpkg") -``` +[//]: # (```) -```{code-cell} ipython3 -:tags: [remove-cell] -import os -os.remove("myaoi.gpkg") -``` +[//]: # () +[//]: # (```{seealso}) -## Parsing sensor metadata +[//]: # (For publication-quality cartographic figures, see) -In our case, `rast` would be better opened using the ``parse_sensor_metadata`` argument of a {class}`~geoutils.Raster`, -which tentatively parses metadata recognized from the filename or auxiliary files. +[//]: # (Cartopy (https://scitools.org.uk/cartopy/docs/latest/) or) -```{code-cell} ipython3 -# Name of the image we used -import os -print(os.path.basename(filename_rast)) -``` +[//]: # (GeoPlot (https://residentmario.github.io/geoplot/index.html).) -```{code-cell} ipython3 -# Open while parsing metadata -rast = gu.Raster(filename_rast, parse_sensor_metadata=True, silent=False) -``` +[//]: # (```) + +[//]: # () +[//]: # (### Referencing, transformations and interfacing) + +[//]: # () +[//]: # (Most transformations and interfacing methods can accept another dataset as a **reference to match** during the operation (for example matching bounds, CRS or grid resolution). See the {ref}`core-match-ref` page for details.) + +[//]: # () +[//]: # (Here, we crop the vector to the raster extent.) + +[//]: # () +[//]: # (```{code-cell} ipython3) + +[//]: # (# Print initial bounds) + +[//]: # (print(gdf.vct.bounds)) + +[//]: # () +[//]: # (# Crop vector to raster extent) + +[//]: # (gdf_crop = gdf.vct.crop(ds, clip=True)) + +[//]: # () +[//]: # (print(gdf_crop.bounds)) + +[//]: # (```) + +[//]: # () +[//]: # (As GeoUtils uses **consistent method names across object types** whenever possible, the same {meth}`~geoutils.Raster.bounds` attribute ) + +[//]: # (and {meth}`~geoutils.Raster.crop` method would also work on a **raster** or **point cloud**.) + +[//]: # () +[//]: # (### Numerical operations and NumPy interface) + +[//]: # () +[//]: # (GeoUtils integrates naturally with the **scientific Python stack**.) + +[//]: # () +[//]: # (Objects, whether {class}`~geoutils.Raster` and {class}`~geoutils.PointCloud` or **Xarray/GeoPandas** data structures, can be used directly with **NumPy or SciPy**. ) + +[//]: # (Standard Python operations such as arithmetic, indexing or logical comparisons behave as expected.) + +[//]: # () +[//]: # (```{code-cell} ipython3) + +[//]: # (import numpy as np) + +[//]: # () +[//]: # (# Normalize raster values) + +[//]: # (ds = (ds - np.min(ds)) / (np.max(ds) - np.min(ds))) + +[//]: # (```) + +[//]: # () +[//]: # (### Masking and implicit behaviour) + +[//]: # () +[//]: # (Logical operations automatically create **boolean raster masks**:) + +[//]: # () +[//]: # (```{code-cell} ipython3) + +[//]: # (mask_aoi = np.logical_and(ds > 0.6, rast_proximity_to_vec > 200)) + +[//]: # (```) + +[//]: # () +[//]: # (These masks can be used for indexing, or to trigger implicit behaviour in certain spatial operations such as proximity calculation or polygonization where target pixels are clearly defined.) + +[//]: # () +[//]: # (```{code-cell} ipython3) + +[//]: # (vect_aoi = mask_aoi.rst.polygonize()) + +[//]: # (```) + +[//]: # () +[//]: # (Finally, GeoUtils contains its own {attr}`is_mask` attribute and option on file opening to manipulate masks conveniently under-the-hood (opening, reading, writing), as many geospatial file formats ) + +[//]: # (do not support **boolean types** and thus expect conversion to an integer type.) + +[//]: # () +[//]: # (### Higher-level analysis tools) + +[//]: # () +[//]: # (GeoUtils also provides higher-level spatial analysis tools.) + +[//]: # () +[//]: # (For example, one can compute the distance to vector geometries on a raster grid:) + +[//]: # () +[//]: # (```{code-cell} ipython3) + +[//]: # (# Compute proximity to vector on raster grid) + +[//]: # (rast_prox = vect.vct.proximity(rast)) + +[//]: # (```) + +[//]: # () +[//]: # (```{note}) + +[//]: # (The raster array may still not be loaded in memory.) + +[//]: # (GeoUtils loads data only when required by an operation.) + +[//]: # (```) + +[//]: # () +[//]: # (## Saving results) + +[//]: # () +[//]: # (GeoUtils objects and derived datasets can easily be written to disk.) + +[//]: # () +[//]: # (```{code-cell} ipython3) + +[//]: # (vect_aoi.to_file("myaoi.gpkg")) + +[//]: # (```) -```{admonition} Wrap-up -In a few lines, we: - - **easily handled georeferencing** operations on rasters and vectors, - - performed numerical calculations **inherently respecting invalid data**, - - **casted to a mask** implicitly from a logical operation on raster, and - - **vectorized a mask** without need for any additional metadata, simply using the nature of the mask object! -**Our result:** a vector of high infrared absorption indexes at least 200 meters away from glaciers -near Everest, which likely corresponds to **perennial snowfields**. -Otherwise, for more **hands-on** examples, explore GeoUtils' gallery of examples! -``` diff --git a/doc/source/index.md b/doc/source/index.md index 63ae82f0c..afc0a876c 100644 --- a/doc/source/index.md +++ b/doc/source/index.md @@ -21,7 +21,7 @@ title: GeoUtils :class: sd-fs-3 :child-align: center -GeoUtils is a Python package for **accessible** and **consistent** geospatial analysis. +GeoUtils is a Python package for **accessible**, **consistent** and **scalable** geospatial analysis. :::: ```{important} @@ -30,12 +30,14 @@ GeoUtils ``v0.2`` is released with more consistent point cloud support! We are w ``` 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 +(NumPy, SciPy, Numba) 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** (Dask support in development for future Xarray accessor). +It strives **to be computationally scalable** by adding **lazy and chunked implementations** to most **raster and point cloud** operations (Dask, Multiprocessing) +and **provides accessors to naturally extend existing Python data-structures** (Xarray, Pandas). + +GeoUtils is **tailored to perform quantitative analysis that implicitly understands the intricacies of geospatial data** +(nodata values, projection, pixel interpretation), through **an intuitive API to foster accessibility** (similar spirit as GDAL's new overhauled CLI). If you are looking to **port your GDAL or QGIS workflow in Python**, GeoUtils is made for you! @@ -98,7 +100,6 @@ about_geoutils how_to_install quick_start feature_overview -summary ``` ```{toctree} diff --git a/doc/source/pointcloud_class.md b/doc/source/pointcloud_class.md index 455041657..a6ebeee23 100644 --- a/doc/source/pointcloud_class.md +++ b/doc/source/pointcloud_class.md @@ -12,7 +12,7 @@ kernelspec: --- (point-cloud)= -# The georeferenced point cloud ({class}`~geoutils.PointCloud`) +# The georeferenced point cloud A point cloud represents 2D point geometries of georeferenced coordinates associated with a main 1D data array, and optionally auxiliary data. diff --git a/doc/source/raster_class.md b/doc/source/raster_class.md index bfd0149aa..c1032d496 100644 --- a/doc/source/raster_class.md +++ b/doc/source/raster_class.md @@ -12,30 +12,51 @@ kernelspec: --- (raster-class)= -# The georeferenced raster ({class}`~geoutils.Raster`) +# The georeferenced raster -Below, a summary of the {class}`~geoutils.Raster` object and its methods. +In GeoUtils, the georeferenced raster object is mirrored through two objects: + +- The Xarray {class}`rst ` accessor for a {class}`xarray.DataArray`, +- The {class}`~geoutils.Raster`. + +We recommend using **only one object type or the other**. While their behaviour is almost entirely similar, there are some differences that are summarized directly below. + +## Accessor {class}`rst ` versus {class}`~geoutils.Raster` + +The main differences are the following: + +1. A {class}`~geoutils.Raster` relies on {class}`~numpy.ma.MaskedArray` to **manipulate integer-type arrays while respecting {attr}`~geoutils.Raster.nodata`** values, while through Xarray the {class}`rst ` accessor only supports floating-type {class}`~numpy.ndarray` to propagate NaNs. We thus **enforce conversion of {class}`~xarray.DataArray` to floating-type on file opening**. + +2. A {class}`~geoutils.Raster` has **control over data-structure operations** (e.g. {func}`+`, {func}`-`, or NumPy interfacing), allowing to raise errors where appropriate +(e.g., if two rasters have same shape but different {attr}`~geoutils.Raster.crs`). For a {class}`~xarray.DataArray`, these operations **rely on user's rigour** but errors should be scarce as {attr}`~geoutils.Raster.shape` is checked during Xarray casting. + +3. A {class}`~geoutils.Raster` can have **ambiguous casting behaviour** for subclasses (e.g., a {class}`~xdem.DEM`) making maintenance difficult, while accessors' **data-structure-centered mechanism enables clearer interfacing**. + +4. A {class}`~geoutils.Raster` currently only supports **Multiprocessing** as scalable backend which is not lazy, while the {class}`rst ` accessor support **Dask** allowing lazy graph building. + +```{important} +To avoid accidentally propagating raw nodata values through operations, **{meth}`~geoutils.open_raster` forces conversion of the {class}`~xarray.DataArray` to floating-type to use NaNs**. +While this increases memory usage, it can be mitigated by using **Dask**. See the {ref}`scalability-index` section for more details. +``` (raster-obj-def)= ## Object definition and attributes -A {class}`~geoutils.Raster` contains **four main attributes**: +A **raster** has **four main attributes**: -1. a {class}`numpy.ma.MaskedArray` as {attr}`~geoutils.Raster.data`, of either {class}`~numpy.integer` or {class}`~numpy.floating` {class}`~numpy.dtype`, +1. a {class}`numpy.ma.MaskedArray` ({class}`~geoutils.Raster`) or {class}`np.ndarray` ({class}`rst ` accessor) as {attr}`~geoutils.Raster.data`, of either {class}`~numpy.integer` or {class}`~numpy.floating` {class}`~numpy.dtype` (forced to {class}`~numpy.floating` for {class}`rst ` accessor), 2. an [{class}`affine.Affine`](https://rasterio.readthedocs.io/en/stable/topics/migrating-to-v1.html#affine-affine-vs-gdal-style-geotransforms) as {attr}`~geoutils.Raster.transform`, 3. a {class}`pyproj.crs.CRS` as {attr}`~geoutils.Raster.crs`, and 4. a {class}`float` or {class}`int` as {attr}`~geoutils.Raster.nodata`. -For more details on {class}`~geoutils.Raster` class composition, see {ref}`core-composition`. - -A {class}`~geoutils.Raster` also contains many derivative attributes, with naming generally consistent with that of a {class}`rasterio.io.DatasetReader`. +A **raster** also contains many derivative attributes, with naming generally consistent with that of [GDAL's recently overhauled CLI](https://gdal.org/en/stable/programs/index.html) or Rasterio. A first category includes georeferencing attributes directly derived from {attr}`~geoutils.Raster.transform`, namely: {attr}`~geoutils.Raster.shape`, {attr}`~geoutils.Raster.height`, {attr}`~geoutils.Raster.width`, {attr}`~geoutils.Raster.res`, {attr}`~geoutils.Raster.bounds`. A second category concerns the attributes derived from the raster array shape and type: {attr}`~geoutils.Raster.count`, {attr}`~geoutils.Raster.bands` and -{attr}`~geoutils.Raster.dtype`. The two former refer to the number of bands loaded in a {class}`~geoutils.Raster`, and the band indexes. +{attr}`~geoutils.Raster.dtype`. The two former refer to the number of bands loaded in a **raster**, and the band indexes. The {attr}`~geoutils.Raster.is_mask` describes if the raster is a mask (i.e. a boolean raster), which overrides the behaviour of some methods to facilitate their manipulation, as boolean data types are not natively supported by raster filetypes or many operations despite their usefulness for analysis (see @@ -56,14 +77,14 @@ The {attr}`~geoutils.Raster.count` and {attr}`~geoutils.Raster.bands` attributes {attr}`~geoutils.Raster.bands_on_disk` only refers to the number of bands on the **on-disk** dataset, if it exists. For example, {attr}`~geoutils.Raster.count` and {attr}`~geoutils.Raster.count_on_disk` will differ when a single band is loaded from a -3-band **on-disk** file, by passing a single index to the `bands` argument in {class}`~geoutils.Raster` or {func}`~geoutils.Raster.load`. +3-band **on-disk** file, by passing a single index to the `bands` argument in **raster** or {func}`~geoutils.Raster.load`. ``` -The complete list of {class}`~geoutils.Raster` attributes with description is available in {ref}`dedicated sections of the API`. +The complete list of **raster** attributes with description is available in {ref}`dedicated sections of the API`. ## Open and save -A {class}`~geoutils.Raster` is opened by instantiating with either a {class}`str`, a {class}`pathlib.Path`, a {class}`rasterio.io.DatasetReader` or a +A **raster** is opened by instantiating with either a {class}`str`, a {class}`pathlib.Path`, a {class}`rasterio.io.DatasetReader` or a {class}`rasterio.io.MemoryFile`. @@ -76,7 +97,7 @@ rast = gu.Raster(filename_rast) rast ``` -Detailed information on the {class}`~geoutils.Raster` is printed using {func}`~geoutils.Raster.info`, along with basic statistics using `stats=True`: +Detailed information on the **raster** is printed using {func}`~geoutils.Raster.info`, along with basic statistics using `stats=True`: ```{code-cell} ipython3 # Print details of raster @@ -87,7 +108,7 @@ print(rast.info(stats=True)) Calling {class}`~geoutils.Raster.info()` with `stats=True` automatically loads the array in-memory, like any other operation calling {attr}`~geoutils.Raster.data`. ``` -A {class}`~geoutils.Raster` is saved to file by calling {func}`~geoutils.Raster.to_file` with a {class}`str` or a {class}`pathlib.Path`. +A **raster** is saved to file by calling {func}`~geoutils.Raster.to_file` with a {class}`str` or a {class}`pathlib.Path`. ```{code-cell} ipython3 # Save raster to disk @@ -101,7 +122,7 @@ os.remove("myraster.tif") ## Create from {class}`~numpy.ndarray` -A {class}`~geoutils.Raster` is created from an array by calling the class method {func}`~geoutils.Raster.from_array` and passing the +A **raster** is created from an array by calling the class method {func}`~geoutils.Raster.from_array` and passing the {ref}`four main attributes`. ```{code-cell} ipython3 @@ -133,7 +154,7 @@ matching the {attr}`~geoutils.Raster.nodata` value passed to {func}`~geoutils.Ra ## Get array -The array of a {class}`~geoutils.Raster` is available in {class}`~geoutils.Raster.data` as a {class}`~numpy.ma.MaskedArray`. +The array of a **raster** is available in {class}`~geoutils.Raster.data` as a {class}`~numpy.ma.MaskedArray`. ```{code-cell} ipython3 # Get raster's masked-array @@ -150,27 +171,27 @@ rast.get_nanarray() ```{important} Getting a {class}`~numpy.ndarray` filled with {class}`~numpy.nan` will automatically cast the {class}`dtype` to {class}`numpy.float32`. This -might result in larger memory usage than in the original {class}`~geoutils.Raster` (if of {class}`int` type). +might result in larger memory usage than in the original **raster** (if of {class}`int` type). -Thanks to the {ref}`core-array-funcs`, **NumPy functions applied directly to a {class}`~geoutils.Raster` will respect {class}`~geoutils.Raster.nodata` +Thanks to the {ref}`core-array-funcs`, **NumPy functions applied directly to a **raster** will respect {class}`~geoutils.Raster.nodata` values** as well as if computing with the {class}`~numpy.ma.MaskedArray` or an unmasked {class}`~numpy.ndarray` filled with {class}`~numpy.nan`. -Additionally, the {class}`~geoutils.Raster` will automatically cast between different {class}`dtype`, and possibly re-define missing +Additionally, the **raster** will automatically cast between different {class}`dtype`, and possibly re-define missing {class}`nodatas`. ``` ## Arithmetic -A {class}`~geoutils.Raster` can be applied any pythonic arithmetic operation ({func}`+`, {func}`-`, {func}`/`, {func}`//`, {func}`*`, -{func}`**`, {func}`%`) with another {class}`~geoutils.Raster`, {class}`~numpy.ndarray` or number. It will output one or two {class}`Rasters`. NumPy coercion rules apply for {class}`dtype`. +A **raster** can be applied any pythonic arithmetic operation ({func}`+`, {func}`-`, {func}`/`, {func}`//`, {func}`*`, +{func}`**`, {func}`%`) with another **raster**, {class}`~numpy.ndarray` or number. It will output one or two {class}`Rasters`. NumPy coercion rules apply for {class}`dtype`. ```{code-cell} ipython3 # Add 1 and divide raster by 2 (rast + 1)/2 ``` -A {class}`~geoutils.Raster` can also be applied any pythonic logical comparison operation ({func}`==`, {func}` != `, {func}`>=`, {func}`>`, {func}`<=`, -{func}`<`) with another {class}`~geoutils.Raster`, {class}`~numpy.ndarray` or number. It will cast to a raster mask, i.e. a boolean {class} +A **raster** can also be applied any pythonic logical comparison operation ({func}`==`, {func}` != `, {func}`>=`, {func}`>`, {func}`<=`, +{func}`<`) with another **raster**, {class}`~numpy.ndarray` or number. It will cast to a raster mask, i.e. a boolean {class} `~geoutils.Raster`. ```{code-cell} ipython3 @@ -182,15 +203,15 @@ See {ref}`core-py-ops` for more details. ## Array interface -A {class}`~geoutils.Raster` can be applied any NumPy universal functions and most mathematical, logical or masked-array functions with another -{class}`~geoutils.Raster`, {class}`~numpy.ndarray` or number. +A **raster** can be applied any NumPy universal functions and most mathematical, logical or masked-array functions with another +**raster**, {class}`~numpy.ndarray` or number. ```{code-cell} ipython3 # Compute the element-wise square-root np.sqrt(rast) ``` -Logical comparison functions will cast to a raster mask, i.e. a boolean {class}`~geoutils.Raster` (True or False). +Logical comparison functions will cast to a raster mask, i.e. a boolean **raster** (True or False). ```{code-cell} ipython3 # Is the raster close to another within tolerance? @@ -202,15 +223,15 @@ See {ref}`core-array-funcs` for more details. ## Reproject -Reprojecting a {class}`~geoutils.Raster` is done through the {func}`~geoutils.Raster.reproject` function, which enforces new {attr}`~geoutils.Raster.transform` +Reprojecting a **raster** is done through the {func}`~geoutils.Raster.reproject` function, which enforces new {attr}`~geoutils.Raster.transform` and/or {class}`~geoutils.Raster.crs`. ```{important} -As with all geospatial handling methods, the {func}`~geoutils.Raster.reproject` function can be passed a {class}`~geoutils.Raster` or +As with all geospatial handling methods, the {func}`~geoutils.Raster.reproject` function can be passed a **raster** or {class}`~geoutils.Vector` as a reference to match. In that case, no other argument is necessary. -A {class}`~geoutils.Raster` reference will enforce to match its {attr}`~geoutils.Raster.transform` and {class}`~geoutils.Raster.crs`. +A **raster** reference will enforce to match its {attr}`~geoutils.Raster.transform` and {class}`~geoutils.Raster.crs`. A {class}`~geoutils.Vector` reference will enforce to match its {attr}`~geoutils.Vector.bounds` and {class}`~geoutils.Vector.crs`. See {ref}`core-match-ref` for more details. @@ -249,7 +270,7 @@ Resampling methods are listed in **[the dedicated section of Rasterio's API](htt [//]: # (```{note}) -[//]: # (Reprojecting a {class}`~geoutils.Raster` can be done out-of-memory in multiprocessing by passing a) +[//]: # (Reprojecting a **raster** can be done out-of-memory in multiprocessing by passing a) [//]: # ({class}`~geoutils.raster.MultiprocConfig` parameter to the {func}`~geoutils.Raster.reproject` function.) @@ -284,12 +305,12 @@ Resampling methods are listed in **[the dedicated section of Rasterio's API](htt ## Crop -Cropping a {class}`~geoutils.Raster` is done through the {func}`~geoutils.Raster.crop` function, which enforces new {attr}`~geoutils.Raster.bounds`. +Cropping a **raster** is done through the {func}`~geoutils.Raster.crop` function, which enforces new {attr}`~geoutils.Raster.bounds`. Additionally, you can use the {func}`~geoutils.Raster.icrop` method to crop the raster using pixel coordinates instead of geographic bounds. Both cropping methods can be used before loading the raster's data into memory. This optimization can prevent loading unnecessary parts of the data, which is particularly useful when working with large rasters. ```{important} -As with all geospatial handling methods, the {func}`~geoutils.Raster.crop` function can be passed only a {class}`~geoutils.Raster` or {class}`~geoutils.Vector` +As with all geospatial handling methods, the {func}`~geoutils.Raster.crop` function can be passed only a **raster** or {class}`~geoutils.Vector` as a reference to match. In that case, no other argument is necessary. See {ref}`core-match-ref` for more details. @@ -301,14 +322,12 @@ The {func}`~geoutils.Raster.icrop` function accepts only a bounding box in pixel By default, {func}`~geoutils.Raster.crop` and {func}`~geoutils.Raster.icrop` return a new Raster unless the inplace parameter is set to True, in which case the cropping operation is performed directly on the original raster object. For more details, see the {ref}`specific section and function descriptions in the API`. -### Example for {func}`~geoutils.Raster.crop` ```{code-cell} ipython3 # Crop raster to smaller bounds rast_crop = rast.crop(bbox=(0.3, 0.3, 1, 1)) print(rast_crop.bounds) ``` -### Example for {func}`~geoutils.Raster.icrop` ```{code-cell} ipython3 # Crop raster using pixel coordinates rast_icrop = rast.icrop(bbox=(2, 2, 6, 6)) @@ -317,11 +336,11 @@ print(rast_icrop.bounds) ## Polygonize -Polygonizing a {class}`~geoutils.Raster` is done through the {func}`~geoutils.Raster.polygonize` function, which converts target pixels into a multi-polygon +Polygonizing a **raster** is done through the {func}`~geoutils.Raster.polygonize` function, which converts target pixels into a multi-polygon {class}`~geoutils.Vector`. ```{note} -For a boolean {class}`~geoutils.Raster`, {func}`~geoutils.Raster.polygonize` implicitly targets `True` values and thus does not require target pixels. +For a boolean **raster**, {func}`~geoutils.Raster.polygonize` implicitly targets `True` values and thus does not require target pixels. ``` ```{code-cell} ipython3 @@ -332,11 +351,11 @@ vect_lt_100 ## Proximity -Computing proximity from a {class}`~geoutils.Raster` is done through by the {func}`~geoutils.Raster.proximity` function, which computes the closest distance -to any target pixels in the {class}`~geoutils.Raster`. +Computing proximity from a **raster** is done through by the {func}`~geoutils.Raster.proximity` function, which computes the closest distance +to any target pixels in the **raster**. ```{note} -For a boolean {class}`~geoutils.Raster`, {func}`~geoutils.Raster.proximity` implicitly targets `True` values and thus does not require target pixels. +For a boolean **raster**, {func}`~geoutils.Raster.proximity` implicitly targets `True` values and thus does not require target pixels. ``` ```{code-cell} ipython3 @@ -355,9 +374,9 @@ prox_lt_100_from_vect ## Interpolate or reduce to point -Interpolating or extracting {class}`~geoutils.Raster` values at specific points can be done through: +Interpolating or extracting **raster** values at specific points can be done through: - the {func}`~geoutils.Raster.reduce_points` function, that applies a reductor function ({func}`numpy.ma.mean` by default) to a surrounding window for each coordinate, or -- the {func}`~geoutils.Raster.interp_points` function, that interpolates the {class}`~geoutils.Raster`'s regular grid to each coordinate using a resampling algorithm. +- the {func}`~geoutils.Raster.interp_points` function, that interpolates the **raster**'s regular grid to each coordinate using a resampling algorithm. ```{code-cell} ipython3 # Extract median value in a 3 x 3 pixel window @@ -376,7 +395,7 @@ Both {func}`~geoutils.Raster.reduce_points` and {func}`~geoutils.Raster.interp_p ## Export -A {class}`~geoutils.Raster` can be exported to different formats, to facilitate inter-compatibility with different packages and code versions. +A **raster** can be exported to different formats, to facilitate inter-compatibility with different packages and code versions. Those include exporting to: - a {class}`xarray.Dataset` with {class}`~geoutils.Raster.to_xarray`, @@ -418,9 +437,9 @@ pc_sub = rast.subsample(500) See {ref}`stats` for more details. (mask-type)= -# The georeferenced raster mask (boolean {class}`~geoutils.Raster`) +# The georeferenced raster mask (boolean **raster**) -A raster mask is a boolean {class}`~geoutils.Raster` (True or False). +A raster mask is a boolean **raster** (True or False). While boolean data types are typically not supported in raster filetypes or in-memory operations, they are incredibly useful for various logical and arithmetical operation in geospatial analysis, so GeoUtils facilitates their manipulation to support these operations natively and implicitly. @@ -430,12 +449,12 @@ Most raster file formats such a [GeoTIFFs](https://gdal.org/drivers/raster/gtiff on-disk**, and **most of Rasterio functionalities also do not support {class}`bool` {class}`dtype`**. To address this, during opening, saving and other geospatial handling operations, raster masks are automatically converted to and from {class}`numpy.uint8`. -The {class}`~geoutils.Raster.nodata` of a boolean {class}`~geoutils.Raster` can now be defined to save to a file, and defaults to `255`. +The {class}`~geoutils.Raster.nodata` of a boolean **raster** can now be defined to save to a file, and defaults to `255`. ``` ## Open, cast and save -A raster mask can be opened from a file through instantiation with {class}`~geoutils.Raster` with the argument `is_mask=True`. +A raster mask can be opened from a file through instantiation with **raster** with the argument `is_mask=True`. On opening, all data will be forced to a {class}`bool` {class}`numpy.dtype`. @@ -448,8 +467,8 @@ mask = gu.Raster(filename_mask, load_data=True, is_mask=True) mask ``` -Raster masks are automatically cast by a logical comparison operation performed on a {class}`~geoutils.Raster` with either another -{class}`~geoutils.Raster`, a {class}`~numpy.ndarray` or a number. +Raster masks are automatically cast by a logical comparison operation performed on a **raster** with either another +**raster**, a {class}`~numpy.ndarray` or a number. ```{code-cell} ipython3 # Instantiate a raster from disk @@ -492,7 +511,7 @@ mask Raster masks can also be created from a {class}`~geoutils.Vector` using {class}`~geoutils.Vector.create_mask`, which rasterizes all input geometries to a boolean array through {class}`~geoutils.Vector.rasterize`. -Georeferencing attributes to create the {class}`~geoutils.Raster` mask can also be passed individually, using `bounds`, `crs`, and `res`. +Georeferencing attributes to create the **raster** mask can also be passed individually, using `bounds`, `crs`, and `res`. ```{code-cell} ipython3 # Open a vector of glacier outlines @@ -522,7 +541,7 @@ Raster masks can be used for indexing and index assignment operations ({func}`[] {class}`Raster`. ```{important} -When indexing, a flattened {class}`~numpy.ma.MaskedArray` is returned with the indexed values of the {class}`~geoutils.Raster` **excluding those masked in its +When indexing, a flattened {class}`~numpy.ma.MaskedArray` is returned with the indexed values of the **raster** **excluding those masked in its {class}`~geoutils.Raster.data`'s {class}`~numpy.ma.MaskedArray`**. ``` @@ -545,8 +564,8 @@ The {func}`~geoutils.Raster.polygonize` function is one of those, implicitly app mask.polygonize() ``` -The {func}`~geoutils.Raster.proximity` function is another method of {class}`~geoutils.Raster` implicitly applying to the `True` values of the mask as -target pixels. It outputs a {class}`~geoutils.Raster` of the distances to the input mask. +The {func}`~geoutils.Raster.proximity` function is another method of **raster** implicitly applying to the `True` values of the mask as +target pixels. It outputs a **raster** of the distances to the input mask. ```{code-cell} ipython3 # Proximity to mask diff --git a/doc/source/release_notes.md b/doc/source/release_notes.md index 779ccc82b..5b6d481ff 100644 --- a/doc/source/release_notes.md +++ b/doc/source/release_notes.md @@ -54,7 +54,7 @@ Based on recent and ongoing progress, we envision the following roadmap. **Releases of 0.2, 0.3, 0.4, etc**, for the following planned (ongoing) additions: - The **addition of a point cloud `PointCloud` data object**, inherited from the `Vector` object alongside many features at the interface of point and raster, -- The **addition of a Xarray accessor `rst`** mirroring the `Raster` object, to work natively with Xarray objects and add support on out-of-memory Dask operations for most of GeoUtils' features, +- The **addition of a Xarray accessor {class}`rst `** mirroring the `Raster` object, to work natively with Xarray objects and add support on out-of-memory Dask operations for most of GeoUtils' features, - The **addition of a GeoPandas accessor `pc`** mirroring the `PointCloud` object, to work natively with GeoPandas objects, - The **addition of statistical features** including zonal statistics (e.g., statistics per vector geometry), grouped statistics (e.g., binning with other variables) and spatial statistics (variogram and kriging) through optional dependencies. - The **addition of filtering and gap-filling features** natively robust to nodata and working similarly for all type of geospatial objects. diff --git a/doc/source/scalability_concept.md b/doc/source/scalability_concept.md index bd0d34a45..07c8fce8f 100644 --- a/doc/source/scalability_concept.md +++ b/doc/source/scalability_concept.md @@ -1,7 +1,19 @@ +--- +file_format: mystnb +jupytext: + formats: md:myst + text_representation: + extension: .md + format_name: myst +kernelspec: + display_name: geoutils-env + language: python + name: geoutils +--- (scalability-concept)= # Concept definitions -This section describes scalability concepts important to grasp to manipulate our objects. +This section describes scalability concepts that are important to grasp to efficiently manipulate geospatial objects. Scalable execution relies on three complementary mechanisms: @@ -9,27 +21,55 @@ Scalable execution relies on three complementary mechanisms: - **Chunked execution:** Operations that process data tile-by-tile to limit memory usage, - **Lazy execution:** Operations whose computation is deferred until explicitly requested. -These mechanisms are often combined (e.g., Dask operations are always chunked **and** lazy) but are conceptually independent. +These mechanisms are often combined (e.g., Dask operations are always **chunked and lazy**) but are conceptually independent. -Finally, one should note that the above concepts only apply to operations that interact with the underlying **data arrays or geometries** of GeoUtils objects. -Naturally, all **metadata operations** (e.g., accessing {attr}`~geoutils.Raster.crs`, {attr}`~geoutils.Raster.bounds`, or {meth}`~geoutils.Raster.info`) have no effect on the array, and do not trigger any loading. +```{note} +The above concepts only apply to operations that interact with the underlying **data arrays or geometries** of objects. +Naturally, all **metadata operations** (e.g., accessing {attr}`~geoutils.Raster.crs`, {attr}`~geoutils.Raster.bounds`, or {meth}`~geoutils.Raster.info`) have no effect on the array, and therefore do not trigger any loading or scalable execution. +``` ## Deferred I/O and implicit loading -**Deferred input/output** refers to operations that modify only **internal I/O metadata**, avoid reading the data entirely and postponing loading. +**Deferred input/output** refers to operations that modify only **internal I/O metadata**, avoiding reading the data entirely and postponing loading. Typical examples include {meth}`~geoutils.Raster.crop`, {meth}`~geoutils.Raster.copy`, and {meth}`~geoutils.Raster.translate`, which behave similarly as Xarray's {meth}`~xarray.DataArray.sel`, {meth}`~xarray.DataArray.copy`, or {meth}`~xarray.DataArray.assign_coords`. -When using the Xarray `rst` accessor, this behavior follows the **native Xarray deferred I/O model**. The {class}`~geoutils.Raster` class implements the -same behavior so that both APIs have consistent semantics. - -This behaviour pairs intrinsically with **implicit loading:** When an object is opened, only metadata is loaded. -Accessing {attr}`~geoutils.Raster.data`, or calling operations that require the array, will **implicitly load the data into memory**. +When using the Xarray {class}`rst ` accessor, this behavior follows the **native Xarray deferred I/O model**. The {class}`~geoutils.Raster` class implements the +same behavior so that both APIs have consistent semantics. An important aspect of **deferred I/O** is that it works with both **in-memory** (NumPy) and **scalable backends** (Dask), allowing to extract parts of large files without any chunked or lazy considerations. +```{code-cell} +import geoutils as gu + +# We open the dataset without Dask backend (same behaviour with Raster class) +filename_rast = gu.examples.get_path("exploradores_aster_dem") +ds = gu.open_raster(filename_rast) + +# We crop the data +ds_cropped = ds.rst.icrop((0, 0, 100, 100)) + +# Neither input nor output dataset are loaded yet +print(f"Input loaded by deferred I/O? {ds.rst.is_loaded}") +print(f"Output loaded by deferred I/O? {ds_cropped.rst.is_loaded}") +``` + +This behaviour pairs intrinsically with **implicit loading:** When an object is opened, only metadata is loaded. +Accessing {attr}`~geoutils.Raster.data`, or calling operations that require the underlying array or geometries will **implicitly load the data into memory**. + +```{code-cell} +# Is the above dataset loaded? +print(f"Loaded before data operation? {ds.rst.is_loaded}") + +# We compute statistics, which loads the array +ds.rst.get_stats() + +# The dataset is now loaded +print(f"Loaded after data operation? {ds.rst.is_loaded}") +``` + ## Chunked execution **Chunked execution** refers to processing raster data **tile-by-tile** instead of loading the full array into memory. @@ -38,21 +78,64 @@ This enables **out-of-core execution**, allowing datasets larger than available In GeoUtils, chunked execution is implemented through two backends: -- **Dask**, used through the Xarray `rst` accessor, +- **Dask**, used through the Xarray {class}`rst ` accessor, - **Multiprocessing**, used through the {class}`~geoutils.Raster` object. Both backends read and process raster chunks sequentially, keeping peak memory usage proportional to the chunk size rather than the full dataset size. Chunked execution therefore allows GeoUtils to scale to large datasets while maintaining a **predictable memory footprint**. For a list of expected memory usage per operation, see the {ref}`scalability-support` page. +```{code-cell} +# Open raster (data is not loaded) +rast = gu.Raster(filename_rast) + +# Create Multiprocessing config, output filepath optional (temporary file by default) +mp_config = gu.multiproc.MultiprocConfig(chunk_size=200) + +# Filter raster with a gaussian in a chunked manner through Multiprocessing +rast_filt = rast.filter(method="gaussian", sigma=4, mp_config=mp_config) + +# The operation happened out-of-memory in chunk-by-chunk +print(f"Temporary raster file created during operation: {rast_filt.name}") +print(f"Is input raster loaded? {rast.is_loaded}") +print(f"Is output raster loaded? {rast_filt.is_loaded}") +``` + ## Lazy execution Lazy execution refers to **deferring computation until results are explicitly requested**. -In GeoUtils, lazy execution is available through the Xarray `rst` accessor with **Dask-backed arrays**. +In GeoUtils, lazy execution is available through the Xarray {class}`rst ` accessor with **Dask-backed arrays**. Operations build a **Dask computation graph** instead of executing immediately. The computation is triggered only when required, for example when calling `compute()` or when writing results to disk. It is particularly useful when **chaining multiple raster operations**, because intermediate results do not need to be materialized or written/read from disk (which costs extra I/O time, often much longer than compute time). -Lazy execution always relies on **chunked execution**, but the reverse is not true: chunked processing can also run eagerly, as in the Multiprocessing backend. \ No newline at end of file +Lazy execution always relies on **chunked execution**, but the reverse is not true: chunked processing can also run eagerly, as in the Multiprocessing backend. + +```{code-cell} +# Open raster lazily with chunks (enables Dask) +ds = gu.open_raster(filename_rast, chunks={"x": 200, "y": 200}) + +print("Input is lazy (Dask arrays):\n") +ds +``` + +```{code-cell} +# Interpolate 30 points from array in chunk-by-chunk +import numpy as np +rng = np.random.default_rng(seed=42) +x = rng.uniform(ds.rst.bounds.left, ds.rst.bounds.right, size=30) +y = rng.uniform(ds.rst.bounds.bottom, ds.rst.bounds.top, size=30) +ds_interp = ds.rst.interp_points((x, y), as_array=True) + +# Result is still lazy +print("Result is still lazy after raster interpolation:\n") +ds_interp +``` + +We can materialize it with `compute()`: + +```{code-cell} +ds_interp.compute() +``` \ No newline at end of file diff --git a/doc/source/scalability_logic.md b/doc/source/scalability_logic.md index 635051415..a66cb9f09 100644 --- a/doc/source/scalability_logic.md +++ b/doc/source/scalability_logic.md @@ -1,4 +1,127 @@ (scalability-logic)= # Implementation strategies -TODO LAST \ No newline at end of file +Implementing **chunked execution** requires developping substantial internal logic often invisible to the user, making it difficult to understand what is happening in the background and how to potentially address a scalability issue. + +Additionally, for certain methods, there is often no single best solution. We therefore propose several **strategies** to further improve performance based on the nature of the input data. + +Below, we detail the logic behind our **chunked execution** implementations, both as an educational resource and to help optimize your code (such as memory usage). + +## Summary + +Several operations are easy to support for **chunked execution** as they directly re-use existing **Dask** methods: +- The {meth}`~geoutils.Raster.filter` function uses {func}`~dask.array.map_overlap` with a `depth` (overlap) half the `size` of the filter, +- The {meth}`~geoutils.Raster.proximity` function uses {func}`~dask.array.map_overlap` with a `max_distance` parameter. + +Other operations are more complex and require specific logic {ref}`specific logic described further below` and summarized as: +- The {meth}`~geoutils.Raster.reproject` function **maps the intersection of projected source grid chunks for each destination chunk** (with potentially different CRS, resolution and bounds), and defines default output chunksizes based on resolution change to avoid unexpected memory blowup, +- The {meth}`~geoutils.Raster.polygonize` function polygonizes implements **three chunk-boundary reconciliation strategies** with different considerations for connected-component labeling and stitching of geometries, +- The {meth}`~geoutils.Vector.rasterize` function performs a geometry subsetting then directly **maps rasterized output blocks** only utilizing these subset geometries. +- The {meth}`~geoutils.Raster.interp_points` function performs a **fast regular-grid mapping of point locations in raster chunks**, expanding raster chunks by a few pixels depending on the resampling method, then performing ordered concatenation of outputs, +- The {meth}`~geoutils.Raster.subsample` function performs an initial chunk-by-chunk sum of valid values to define the requested sample size, then samples values per chunk through **reproducible chunk-invariant seeding** (default) or **faster chunk-dependent seeding**. + +(specific-logic)= +## Logic of implementations + +### Chunked reprojection + +**Reprojection with chunked execution** requires mapping **destination chunks**—defined on a regular grid in the output CRS and potentially with different resolution or bounds—to **projected source chunks**, which become geometrically deformed after reprojection. + +The diagram below illustrates this mapping procedure, with a new CRS and a downsampling of 2: + +```{eval-rst} +.. plot:: code/diagram_chunked_reproject.py + :width: 100% +``` + + +During reprojection, destination chunks are **expanded by 3 pixels** to ensure adequate resampling at chunk boundaries. As a result of this and deformations, the number of source chunks used at once in memory is therefore **always 4-9**. + +Unlike most chunked operations, **reprojection does not preserve identical input and output chunk sizes**. When the resolution changes, maintaining the same chunk size would cause the number of intersecting source chunks to scale approximately with **the square of the downsampling factor** (because reprojection operates in two spatial dimensions). For example, a downsampling factor of 2 would cause four times more source chunks to be accessed per destination chunk, which would quickly increase memory usage during coarse reprojection. + +To prevent this, GeoUtils automatically **scales the output chunk size according to the resolution change**, keeping the number of source chunks involved in each operation bounded. + +Finally, note that **chunked reprojection with GCPs or RCPs is currently not supported**. + +### Chunked polygonization + +**Polygonization with chunked execution** requires identifying raster regions that may extend across **chunk boundaries**. Because chunks are processed independently, connected regions intersecting a boundary must later be **reconciled across neighboring chunks** to produce correct polygons. + +The diagram below illustrates the three strategies implemented in GeoUtils for performing this reconciliation: + +```{eval-rst} +.. plot:: code/diagram_chunked_polygonize.py + :width: 100% +``` + +All strategies begin by processing individual raster chunks independently, then reconstruct continuous polygons that span chunk boundaries. + +Conceptually, the methods differ in how cross-chunk regions are reconstructed: +- **`label_union`** labels values in each chunk, then finds the **union of matching labels across chunk seams** before polygonization, avoiding vector comparisons and requiring only a final dissolve step, +- **`label_stitch`** labels values as in **`label_union`**, then polygonizes each chunk independently and **stitches polygons afterward in vector space**, avoiding the need for a union–find structure, +- **`geometry_stitch`** bypasses labeling entirely by performing **polygonization on halo-expanded chunks** (1-pixel overlap), then stitches polygons similarly as in **`label_stitch`** after clipping. + +Connectivity assumptions influence how many neighboring chunks must be used in memory, which remains small and bounded: + +- **4-connectivity** (4 cardinal directions): typically 2–4 chunks, +- **8-connectivity** (adding 4 diagonals): up to 4–9 chunks. + +Finally, note that **polygon stitching occurs only for polygons touching chunk boundaries**. Polygons fully contained within a chunk are produced directly without additional processing. + +### Chunked rasterization + +**Rasterization with chunked execution** distributes the burn operation across **output raster chunks**, allowing vector datasets to be converted into rasters without materializing the entire output array in memory. + +The diagram below illustrates the execution strategy: + +```{eval-rst} +.. plot:: code/diagram_chunked_rasterize.py + :width: 100% +``` + +For every chunk, GeoUtils first performs a **spatial query on the vector geometries**, selecting only those whose bounding boxes intersect the chunk bounds. These candidate geometries are then rasterized into the chunk-local array using {func}`~dask.array.map_blocks`. If no geometries intersect the chunk, the block function exits early and directly returns an array filled with the background value. + +The memory footprint during chunked rasterization is therefore approximately limited to: + +- **One output chunk array**, and +- **The subset of geometries intersecting that chunk**. + +Unlike chunked reprojection or polygonization, rasterization does not require overlap or cross-chunk reconciliation, since each pixel value is determined independently of the vector intersections within the chunk. + +### Chunked interpolation at points + +**Interpolation at points with chunked execution** evaluates raster values at point coordinates by combining **1D chunking of the input points** with **2D chunking of the raster grid**. + +The diagram below illustrates the workflow: + +```{eval-rst} +.. plot:: code/diagram_chunked_interp_points.py + :width: 100% +``` + +First, the input points are **chunked along their 1D sequence**. For each point chunk, GeoUtils uses a fast regular-grid mapping to **raster chunks containing the corresponding point coordinates**. Only those raster chunks are processed, avoiding loading the full raster into memory. + +Interpolation is then performed independently on each required raster chunk. To ensure correct interpolation near chunk boundaries, the raster chunk is **expanded by an overlap depth equal to the half-interpolation-order rounded up + 1** (for example, 1 pixels for nearest, 2 for linear and 3 for cubic). + +Because points within a point chunk may fall into different raster chunks, interpolation proceeds by **looping over the intersecting raster chunks** for that point chunk. The resulting interpolated values are then **concatenated and reordered** to match the original point order. + +This approach keeps memory usage low, as only the raster chunks needed for the current point chunk are loaded at any given time. + +Finally, note that **using eager (in-memory) point coordinates is typically much faster** when the number of points is moderate. Unless necessary, keep points in memory to avoid additional chunking overhead. + +### Chunked subsampling + +**Chunked subsampling of valid raster values** enables selecting representative pixels from very large rasters without materializing all valid values in memory. + +The diagram below illustrates the workflow: + +```{eval-rst} +.. plot:: code/diagram_chunked_subsample.py + :width: 100% +``` + +Subsampling proceeds in **two stages**. First, GeoUtils performs a lightweight pass that counts the number of **valid pixels in each chunk**. These counts are summed to determine the final subsample size requested by the user, which may be specified either as a **fraction of valid pixels** or as an **absolute number**. This pass is inexpensive because it only requires scanning chunk masks and does not materialize the valid-value array. + +Once the target sample size is known, pixels are selected using one of two strategies: +- With **`topk` (chunk-invariant)** sampling, each valid pixel is assigned a deterministic pseudo-random key derived from the random seed and its global linear index (row + col * number of rows). The **k smallest keys** are selected globally. Because the key depends only on the pixel index and the seed, the resulting subsample is **independent of chunk layout** and reproducible with any chunking or fully in-memory raster. +- With **`sequential` (chunk-dependent)** sampling, pixels are drawn from the **flattened sequence of valid values encountered during chunk traversal**. This strategy is typically **faster** because it avoids computing deterministic keys, but the result **depends on the chunk structure and valid-value ordering**. diff --git a/doc/source/scalability_usage.md b/doc/source/scalability_usage.md index 9f786ce25..3e16c6a66 100644 --- a/doc/source/scalability_usage.md +++ b/doc/source/scalability_usage.md @@ -17,7 +17,7 @@ GeoUtils supports scalable execution for most of its **raster** and (soon) **poi It relies on two execution backends: -- **Dask**, through its `rst` Xarray accessor and `pc` Pandas accessor (**lazy** and **chunked** execution), +- **Dask**, through its {class}`rst ` Xarray accessor and `pc` Pandas accessor (**lazy** and **chunked** execution), - **Multiprocessing**, through its {class}`~geoutils.Raster` and {class}`~geoutils.PointCloud` objects (**chunked** execution only) . Both backends mirror the **exact same object operations and chunked logic**, and yield **identical** results as in-memory operations. @@ -26,13 +26,17 @@ For details on scalability concepts, see the {ref}`scalability-concept` page. As a rule of thumb: -- Use **Dask** to work on **Xarray and GeoPandas objects** through our accessors `rst` and `pc`, and if you want to chain several operations lazily. +- Use **Dask** to work on **Xarray and GeoPandas objects** through our accessors {class}`rst ` and `pc`, and if you want to chain several operations lazily. - Use **Multiprocessing** to work with our {class}`~geoutils.Raster` and {class}`~geoutils.PointCloud` objects, and if you are fine with intermediate writing/reading between steps. - Use standard **in-memory execution** to work efficiently on small rasters, which is possible even if those were loaded from larger rasters (use {class}`~geoutils.Raster.crop`). +```{note} +GeoUtils currently targets **scalable CPU execution**. However, as many of our numerical operations rely on **NumPy**, **SciPy** or **Numba**, those are planned to be linked to their **GPU** counterparts (**CuPy** and **Numba CUDA**). +``` + ## Using Dask through accessors -With Dask, raster operations are both **chunked** and **lazy**. +With **Dask**, raster operations are both **chunked** and **lazy**. This behavior is enabled by opening a raster with the `chunks` argument, which returns an Xarray object backed by Dask arrays. ```{code-cell} python @@ -45,7 +49,7 @@ ds = gu.open_raster(filename_rast, chunks={"x": 200, "y": 200}) ds ``` -GeoUtils, through the `rst` accessor, automatically detects the **Dask** input and switches to a chunked implementation. +GeoUtils, through the {class}`rst ` accessor, automatically detects the **Dask** input and switches to a chunked implementation. ```{code-cell} python # Change output resolution @@ -60,7 +64,7 @@ ds_reproj = ds.rst.reproject( ds_reproj ``` -The resulting raster remains **lazy**. Computation only happens when explicitly requested with `compute()`. +The resulting raster remains **lazy**. Computation only happens when explicitly requested with {meth}`~dask.array.Array.compute()`. For a raster output, one typically wants to write to file lazily to avoid loading it in-memory: @@ -81,7 +85,7 @@ sub_ds = ds_reproj.rst.subsample( sub_ds ``` -The output array is again lazy, and in this case we can use `compute()` to return the in-memory NumPy array: +The output array is again lazy, and in this case we can use {meth}`~dask.array.Array.compute()` to return the in-memory NumPy array: ```{code-cell} python sub_ds.compute() ``` @@ -132,6 +136,12 @@ samp_rast_mp = rast_reproj_mp.subsample( samp_rast_mp ``` +```{code-cell} ipython3 +:tags: [remove-cell] +import os +os.remove(mp_config.outfile) +``` + This backend is convenient when working directly with {class}`~geoutils.Raster` objects and performing **step-by-step processing**. ## Good practices with chunked and lazy operations @@ -139,6 +149,7 @@ This backend is convenient when working directly with {class}`~geoutils.Raster` - If **memory** is the limitating factor for you, use a **single-threaded scheduler** through Dask (```dask.config.set(scheduler='single-threaded')```) or Multiprocessing (default cluster), - If **speed** is the limiting factor for you, use **parallelized processes** through Dask (see [Dask scheduler configuration](https://docs.dask.org/en/stable/scheduler-overview.html#scheduler-overview)) or Multiprocessing (see our Cluster configuration), - Choose chunk sizes large enough to reduce scheduling overhead, but **small enough to fit comfortably in memory**, +- Check that your data files have **on-disk chunksizes** (otherwise loads everything) and use a multiple of it for optimal **in-memory chunking**, - Keep chunk sizes **consistent across operations** to avoid unnecessary rechunking, - Insert **breakpoints** (for example by writing intermediate results to disk) to prevent building overly large Dask graphs. diff --git a/doc/source/summary.md b/doc/source/summary.md deleted file mode 100644 index 921efecda..000000000 --- a/doc/source/summary.md +++ /dev/null @@ -1,216 +0,0 @@ -(method-summary)= - -# Feature overview - -GeoUtils provides a unified API for manipulating **raster**, **vector**, and **point-cloud** data, with **scalable execution** for most raster operations. - -The **tables below** summarize the core operations of GeoUtils, their scalability and backends. - -If you are interested in converting from GDAL/OGR, see our {ref}`cheatsheet-osgeo` page. - -## Summary of methods and scalability - -Methods of GeoUtils are shared across object types and expose a **consistent API** for clarity (similarly as the recent [GDAL CLI overhaul](https://gdal.org/en/stable/programs/index.html)). They also support convenient inputs such as **match-reference arguments** (e.g., matching a grid for reprojection or rasterization, matching bounds for cropping, matching point coordinates for interpolation). See the {ref}`core-match-ref` page for details. - -Nearly all **raster operations** support **scalable execution** through [Dask](https://www.dask.org/) or Multiprocessing, allowing large datasets to be processed chunk-by-chunk without loading the full array into memory. -While the table below provide a scalability summary, details on exact **supported operations** relative to inputs/outputs are available on the {ref}`scalability-support` section. - -Some operations also support multiple computational **backends** (for example SciPy or Numba implementations for numerical routines). - -All methods are tested to ensure they produce **identical results** whether executed in-memory, using chunked processing, or through alternative computational backends. - -## Data operations - -We first describe GeoUtils' core **data operations**, which operate on underlying arrays or geometries and therefore benefit from **scalable execution**. - -**Legend:** **“/”** indicates methods **shared across object types**, while **“⟷”** indicates methods **interfacing between two object types**. - -```{list-table} Common API for data operations -:widths: 3 5 1 2 -:header-rows: 1 -:align: left -:class: tight-table - -* - Method - - Notes - - Scalable - - Backend - -* - Raster / Vector / Point - - - - - - - -* - {meth}`~geoutils.Raster.reproject()` - - Reproject to other CRS. Default tolerance parameters ensure chunk-invariance. - - ✅ - - Rasterio / PyProj - -* - {meth}`~geoutils.Raster.crop()` - - Crop to bounds. For vectors, can return geometries intersecting (untouched) or clipped. - - ✅ - - Rasterio / GeoPandas - -* - {meth}`~geoutils.Raster.translate()` - - Apply a grid shift to object. - - ✅ - - NumPy / GeoPandas - -* - {meth}`~geoutils.Raster.proximity()` - - Estimate proximity distance to target values or geometries. - - ❌ - - SciPy - -* - {meth}`~geoutils.Raster.plot()` - - Visualization helper. - - ❌ - - Matplotlib - -* - Raster / Point - - - - - - - -* - {meth}`~geoutils.Vector.create_mask()` - - Create boolean mask of a vector geometries over raster or point. - - ✅ - - Rasterio / GeoPandas - -* - {meth}`~geoutils.Raster.get_stats()` - - Compute statistics of valid values over a valid mask. - - ❌ - - NumPy / SciPy - -* - {meth}`~geoutils.Raster.subsample()` - - Randomly sample valid values. Chunk-invariant seed ensures reproducibility. - - ✅ - - NumPy - -* - {meth}`~geoutils.Raster.filter()` - - Filter over window. Fast vectorized logic with NaN support. - - ✅ - - SciPy - -* - Raster ⟷ Vector - - - - - - - -* - {meth}`~geoutils.Raster.polygonize()` - - Convert raster regions to vector polygons. Multiple chunked strategies for performance. - - ✅ - - Rasterio / GeoPandas - -* - {meth}`~geoutils.Vector.rasterize()` - - Burn vector geometries onto a raster grid. - - ✅ - - Rasterio - -* - Raster ⟷ Point - - - - - - - -* - {meth}`~geoutils.Raster.interp_points()` - - Interpolate raster at point locations. Fast regular-grid logic with added NaN propagation. - - ✅ - - SciPy - -* - {meth}`~geoutils.Raster.reduce_points()` - - Aggregate raster values around points. - - ❌ - - NumPy - -* - {meth}`~geoutils.PointCloud.grid()` - - Grid irregular points onto a raster grid. Multiple approaches with added NaN propagation. - - ❌ - - SciPy - -* - {meth}`~geoutils.Raster.from_pointcloud_regular()` - - Direct conversion when points lie on a regular grid. - - ❌ - - NumPy - -* - {meth}`~geoutils.Raster.to_pointcloud()` - - Conversion to point cloud. - - ❌ - - NumPy -``` - -## Metadata properties and operations - -In addition to data operations, GeoUtils exposes **metadata** properties and methods consistently across geospatial objects. -These operate only on metadata and therefore **do not load or modify underlying data arrays**. - -```{list-table} Common API from metadata operations -:widths: 3 7 -:header-rows: 1 -:align: left -:class: tight-table - -* - Attribute / Method - - Description - -* - Raster / Vector / Point - - - -* - {attr}`~geoutils.Raster.crs` - - Coordinate reference system (CRS) of object. - -* - {attr}`~geoutils.Raster.bounds` - - Bounding box of object. - -* - {attr}`~geoutils.Raster.footprint` - - Footprint polygon geometry of object. - -* - {attr}`~geoutils.Raster.is_loaded` - - Whether geospatial object is loaded in-memory. - -* - {attr}`~geoutils.Raster.name` - - Filename of object on disk, if it exists. - -* - {meth}`~geoutils.Raster.get_bounds_projected()` - - Bounds projected in other CRS. - -* - {meth}`~geoutils.Raster.get_footprint_projected()` - - Footprint polygon geometry in other CRS. - -* - {meth}`~geoutils.Raster.get_metric_crs()` - - Get metric CRS suitable for this object. - -* - {meth}`~geoutils.Raster.info()` - - Summary of attributes for geospatial object. - -* - Raster / Point - - - -* - {attr}`~geoutils.Raster.data` - - Data array (2D grid for raster, 1D for point cloud). - -* - {attr}`~geoutils.Raster.shape` - - Shape of data array. - -* - {attr}`~geoutils.Raster.is_mask` - - Whether object is a mask. Clarifies ambiguity of raster/point file types often not supporting boolean types. - -* - Raster - - - -* - {attr}`~geoutils.Raster.transform` - - Geotransform to map raster indices to spatial coordinates. - -* - {attr}`~geoutils.Raster.nodata` - - Nodata value used to represent missing data on disk. - -* - {attr}`~geoutils.Raster.area_or_point` - - Pixel interpretation of raster values, either center point or area average. - -* - Point - - - -* - {attr}`~geoutils.PointCloud.point_count` - - Number of points in the point cloud. -``` - - - diff --git a/doc/source/vector_class.md b/doc/source/vector_class.md index e9552300b..1255120e1 100644 --- a/doc/source/vector_class.md +++ b/doc/source/vector_class.md @@ -12,7 +12,7 @@ kernelspec: --- (vector-class)= -# The georeferenced vector ({class}`~geoutils.Vector`) +# The georeferenced vector Below, a summary of the {class}`~geoutils.Vector` object and its methods. From 1c86ffaf1dad882a6473e8e5e6e994f54a6af5fd Mon Sep 17 00:00:00 2001 From: Romain Hugonnet Date: Mon, 9 Mar 2026 22:26:17 -0800 Subject: [PATCH 3/3] Remove commented code + add other changes than doc --- doc/source/feature_overview.md | 229 ---------------------------- geoutils/__init__.py | 2 +- geoutils/interface/vectorization.py | 2 +- geoutils/raster/__init__.py | 1 + 4 files changed, 3 insertions(+), 231 deletions(-) diff --git a/doc/source/feature_overview.md b/doc/source/feature_overview.md index c17f365d1..87857c838 100644 --- a/doc/source/feature_overview.md +++ b/doc/source/feature_overview.md @@ -222,232 +222,3 @@ These operate only on metadata and therefore **do not load or modify underlying - Number of points in the point cloud. ``` -[//]: # ((examples-overview)=) - -[//]: # (## Examples) - -[//]: # () -[//]: # (The following presents a descriptive example show-casing core features of GeoUtils.) - -[//]: # () -[//]: # (```{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) - -[//]: # (```) - -[//]: # () -[//]: # (### Core objects and accessors) - -[//]: # () -[//]: # (GeoUtils operates on **rasters, vectors and point clouds**.) - -[//]: # () -[//]: # (Below, we demonstrate features through **accessors**: {class}`rst ` for **Xarray** rasters, `vct` for **GeoPandas** vectors, `pc` for **GeoPandas** point clouds.) - -[//]: # (The exact same operations are available through the {class}`~geoutils.Raster`, {class}`~geoutils.Vector` and {class}`~geoutils.PointCloud` objects.) - -[//]: # () -[//]: # (**Accessors** operations always return another **Xarray** or **GeoPandas** object, while **GeoUtils object** operations always return a **GeoUtils object**.) - -[//]: # (However, you can pass **any object** as an argument to a function.) - -[//]: # () -[//]: # (First, we open datasets with {meth}`~geoutils.open_raster` and {meth}`~geoutils.open_vector`, passing a `chunks` argument to trigger out-of-memory **Dask** behaviour for the raster.) - -[//]: # () -[//]: # (```{code-cell} ipython3) - -[//]: # (import geoutils as gu) - -[//]: # (import geopandas as gpd) - -[//]: # () -[//]: # (# Example files: infrared band of Landsat and glacier outlines) - -[//]: # (filename_rast = gu.examples.get_path("everest_landsat_b4")) - -[//]: # (filename_vect = gu.examples.get_path("everest_rgi_outlines")) - -[//]: # () -[//]: # (# Open datasets) - -[//]: # (ds = gu.open_raster(filename_rast, chunks={"x": 200, "y": 200})) - -[//]: # (gdf = gpd.read_file(filename_vect)) - -[//]: # () -[//]: # (# Raster and vector objects) - -[//]: # (ds) - -[//]: # (gdf) - -[//]: # (```) - -[//]: # () -[//]: # (GeoUtils accessors expose spatial methods directly on these objects while keeping compatibility with the broader **NumPy, Xarray and GeoPandas ecosystem**.) - -[//]: # () -[//]: # (### Plotting) - -[//]: # () -[//]: # (For plotting, GeoUtils includes lightweight plotting helpers to simplify common tasks such as overlaying rasters and vectors with multiple colorbars, while passing through standard plotting arguments.) - -[//]: # () -[//]: # (```{code-cell} ipython3) - -[//]: # (# Plot raster and vector) - -[//]: # (ds.rst.plot(cbar_title="Distance to glacier outline")) - -[//]: # (gdf.plot(rast_proximity_to_vec, fc="none")) - -[//]: # (```) - -[//]: # () -[//]: # (```{seealso}) - -[//]: # (For publication-quality cartographic figures, see) - -[//]: # (Cartopy (https://scitools.org.uk/cartopy/docs/latest/) or) - -[//]: # (GeoPlot (https://residentmario.github.io/geoplot/index.html).) - -[//]: # (```) - -[//]: # () -[//]: # (### Referencing, transformations and interfacing) - -[//]: # () -[//]: # (Most transformations and interfacing methods can accept another dataset as a **reference to match** during the operation (for example matching bounds, CRS or grid resolution). See the {ref}`core-match-ref` page for details.) - -[//]: # () -[//]: # (Here, we crop the vector to the raster extent.) - -[//]: # () -[//]: # (```{code-cell} ipython3) - -[//]: # (# Print initial bounds) - -[//]: # (print(gdf.vct.bounds)) - -[//]: # () -[//]: # (# Crop vector to raster extent) - -[//]: # (gdf_crop = gdf.vct.crop(ds, clip=True)) - -[//]: # () -[//]: # (print(gdf_crop.bounds)) - -[//]: # (```) - -[//]: # () -[//]: # (As GeoUtils uses **consistent method names across object types** whenever possible, the same {meth}`~geoutils.Raster.bounds` attribute ) - -[//]: # (and {meth}`~geoutils.Raster.crop` method would also work on a **raster** or **point cloud**.) - -[//]: # () -[//]: # (### Numerical operations and NumPy interface) - -[//]: # () -[//]: # (GeoUtils integrates naturally with the **scientific Python stack**.) - -[//]: # () -[//]: # (Objects, whether {class}`~geoutils.Raster` and {class}`~geoutils.PointCloud` or **Xarray/GeoPandas** data structures, can be used directly with **NumPy or SciPy**. ) - -[//]: # (Standard Python operations such as arithmetic, indexing or logical comparisons behave as expected.) - -[//]: # () -[//]: # (```{code-cell} ipython3) - -[//]: # (import numpy as np) - -[//]: # () -[//]: # (# Normalize raster values) - -[//]: # (ds = (ds - np.min(ds)) / (np.max(ds) - np.min(ds))) - -[//]: # (```) - -[//]: # () -[//]: # (### Masking and implicit behaviour) - -[//]: # () -[//]: # (Logical operations automatically create **boolean raster masks**:) - -[//]: # () -[//]: # (```{code-cell} ipython3) - -[//]: # (mask_aoi = np.logical_and(ds > 0.6, rast_proximity_to_vec > 200)) - -[//]: # (```) - -[//]: # () -[//]: # (These masks can be used for indexing, or to trigger implicit behaviour in certain spatial operations such as proximity calculation or polygonization where target pixels are clearly defined.) - -[//]: # () -[//]: # (```{code-cell} ipython3) - -[//]: # (vect_aoi = mask_aoi.rst.polygonize()) - -[//]: # (```) - -[//]: # () -[//]: # (Finally, GeoUtils contains its own {attr}`is_mask` attribute and option on file opening to manipulate masks conveniently under-the-hood (opening, reading, writing), as many geospatial file formats ) - -[//]: # (do not support **boolean types** and thus expect conversion to an integer type.) - -[//]: # () -[//]: # (### Higher-level analysis tools) - -[//]: # () -[//]: # (GeoUtils also provides higher-level spatial analysis tools.) - -[//]: # () -[//]: # (For example, one can compute the distance to vector geometries on a raster grid:) - -[//]: # () -[//]: # (```{code-cell} ipython3) - -[//]: # (# Compute proximity to vector on raster grid) - -[//]: # (rast_prox = vect.vct.proximity(rast)) - -[//]: # (```) - -[//]: # () -[//]: # (```{note}) - -[//]: # (The raster array may still not be loaded in memory.) - -[//]: # (GeoUtils loads data only when required by an operation.) - -[//]: # (```) - -[//]: # () -[//]: # (## Saving results) - -[//]: # () -[//]: # (GeoUtils objects and derived datasets can easily be written to disk.) - -[//]: # () -[//]: # (```{code-cell} ipython3) - -[//]: # (vect_aoi.to_file("myaoi.gpkg")) - -[//]: # (```) - - - diff --git a/geoutils/__init__.py b/geoutils/__init__.py index 7fac605a0..b1e94f70d 100644 --- a/geoutils/__init__.py +++ b/geoutils/__init__.py @@ -24,7 +24,7 @@ from geoutils._config import config # noqa from geoutils.raster import Raster, xr_accessor # noqa isort:skip -from geoutils.raster.xr_accessor import open_raster # noqa isort:skip +from geoutils.raster.xr_accessor import open_raster, RasterAccessor # noqa isort:skip from geoutils.vector import Vector # noqa isort:skip from geoutils.pointcloud import PointCloud # noqa isort:skip diff --git a/geoutils/interface/vectorization.py b/geoutils/interface/vectorization.py index 257d4e8c1..75492a43f 100644 --- a/geoutils/interface/vectorization.py +++ b/geoutils/interface/vectorization.py @@ -1677,7 +1677,7 @@ def _chunked_polygonize_core( The function supports three chunked strategies: - 1) ``label_union`` (exact, label-based) + 1) ``label_union`` (label-based) - Build connected-component labels per block. - Scan seams between neighboring blocks and build a union-find mapping that merges labels that touch across seams (and have identical raster value). diff --git a/geoutils/raster/__init__.py b/geoutils/raster/__init__.py index 328a819b2..4098396c7 100644 --- a/geoutils/raster/__init__.py +++ b/geoutils/raster/__init__.py @@ -16,6 +16,7 @@ # See the License for the specific language governing permissions and # limitations under the License. +from geoutils.raster.base import RasterBase from geoutils.raster.raster import Raster, RasterType, handled_array_funcs # noqa isort:skip from geoutils.raster.array import * # noqa from geoutils.raster.multiraster import * # noqa