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