From 6186e30d854ec58930b80a33c09412abbf56d745 Mon Sep 17 00:00:00 2001 From: Romain Hugonnet Date: Thu, 13 Jun 2024 17:48:25 -0800 Subject: [PATCH] Incremental commit on doc --- doc/source/api.md | 6 +-- doc/source/config.md | 37 +++++++++++++++++ doc/source/distance_ops.md | 8 ++-- doc/source/feature_overview.md | 10 +++++ doc/source/georeferencing.md | 53 ++++++++++++++++-------- doc/source/geotransformations.md | 43 +++++++++++-------- doc/source/index.md | 5 ++- doc/source/quick_start.md | 10 +++++ doc/source/raster_vector_point.md | 69 ++++++++++++++++++++++++------- geoutils/config.ini | 7 +++- geoutils/raster/raster.py | 10 +++-- geoutils/vector.py | 36 ++++++++++------ tests/test_vector.py | 2 +- 13 files changed, 217 insertions(+), 79 deletions(-) create mode 100644 doc/source/config.md diff --git a/doc/source/api.md b/doc/source/api.md index 13f62d75..60c59b9e 100644 --- a/doc/source/api.md +++ b/doc/source/api.md @@ -84,8 +84,8 @@ documentation. Raster.reproject Raster.polygonize Raster.proximity - Raster.value_at_coords Raster.interp_points + Raster.reduce_points ``` ### Plotting @@ -121,7 +121,7 @@ documentation. Raster.load Raster.save Raster.to_pointcloud - Raster.from_regular_pointcloud + Raster.from_pointcloud_regular Raster.to_rio_dataset Raster.to_xarray ``` @@ -135,7 +135,7 @@ documentation. Raster.xy2ij Raster.ij2xy Raster.coords - Raster.shift + Raster.translate Raster.outside_image ``` diff --git a/doc/source/config.md b/doc/source/config.md new file mode 100644 index 00000000..79b83c38 --- /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 +``` \ No newline at end of file diff --git a/doc/source/distance_ops.md b/doc/source/distance_ops.md index 1330fd19..c7a6912a 100644 --- a/doc/source/distance_ops.md +++ b/doc/source/distance_ops.md @@ -14,7 +14,7 @@ kernelspec: # Distance operations Computing distance between sets of geospatial data or manipulating their shape based on distance is often important -for later analysis. To facilitate this type of operations, GeoUtils, implements distance-specific functionalities +for later analysis. To facilitate this type of operations, GeoUtils implements distance-specific functionalities for both vectors and rasters. ```{code-cell} ipython3 @@ -57,7 +57,7 @@ vect = gu.Vector(gu.examples.get_path("everest_rgi_outlines")) ```{code-cell} ipython3 # Compute proximity to vector outlines proximity = vect.proximity(rast) -proximity.plot(cmap="viridis") +proximity.plot(cmap="viridis", cbar_title="Proximity (m)") ``` ## Buffering without overlap @@ -71,6 +71,6 @@ between shapes, which is sometimes undesirable. Using Voronoi polygons, we provi # Compute buffer without overlap on glacier outlines vect_buff_nolap = vect.buffer_without_overlap(buffer_size=500) # Plot with color to see that the attributes are retained for every feature -vect_buff_nolap.plot(ax="new", column="Area") -vect.plot(ec="k", column="Area", alpha=0.5) +vect.plot(ax="new", ec="k", column="Area", alpha=0.5, add_cbar=False) +vect_buff_nolap.plot(column="Area", cbar_title="Glacier area (km)") ``` diff --git a/doc/source/feature_overview.md b/doc/source/feature_overview.md index 52b655de..09ee87bf 100644 --- a/doc/source/feature_overview.md +++ b/doc/source/feature_overview.md @@ -25,6 +25,16 @@ All pages of this documentation containing code cells can be **run interactively Alternatively, start your own notebook to test GeoUtils at [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/GlacioHack/geoutils/main). ``` +```{code-cell} ipython3 +:tags: [remove-cell] + +# To get a good resolution for displayed figures +from matplotlib import pyplot +pyplot.rcParams['figure.dpi'] = 600 +pyplot.rcParams['savefig.dpi'] = 600 +pyplot.rcParams['font.size'] = 9 +``` + ## The core {class}`~geoutils.Raster` and {class}`~geoutils.Vector` classes In GeoUtils, geospatial handling is object-based and revolves around {class}`~geoutils.Raster` and {class}`~geoutils.Vector`. diff --git a/doc/source/georeferencing.md b/doc/source/georeferencing.md index 57489ecd..9a28d2f0 100644 --- a/doc/source/georeferencing.md +++ b/doc/source/georeferencing.md @@ -50,23 +50,31 @@ vect = gu.Vector(gu.examples.get_path("exploradores_rgi_outlines")) ``` ```{code-cell} ipython3 -# Print raster and vector info +# Print raster info rast.info() +``` + +```{code-cell} ipython3 +# Print vector info vect.info() ``` ### Coordinate reference systems -[Coordinate reference systems (CRSs)](https://en.wikipedia.org/wiki/Spatial_reference_system), sometimes also called a +[Coordinate reference systems (CRSs)](https://en.wikipedia.org/wiki/Spatial_reference_system), sometimes also called spatial reference systems (SRSs), define the 2D projection of the geospatial data. They are stored as a -{class}`pyproj.crs.CRS` object in {attr}`~geoutils.Raster.crs`. More information on the manipulation -of {class}`pyproj.crs.CRS` objects can be found in [PyProj's documentation](https://pyproj4.github.io/pyproj/stable/). +{class}`pyproj.crs.CRS` object in {attr}`~geoutils.Raster.crs`. ```{code-cell} ipython3 -# Show CRS attribute of a raster and vector +# Show CRS attribute of raster print(rast.crs) -print(vect.crs) ``` +```{code-cell} ipython3 +# Show CRS attribute of vector as a WKT +print(vect.crs.to_wkt()) +``` + +More information on the manipulation of {class}`pyproj.crs.CRS` objects can be found in [PyProj's documentation](https://pyproj4.github.io/pyproj/stable/). ```{note} 3D CRSs for elevation data are only emerging, and not consistently defined in the metadata. @@ -74,15 +82,19 @@ The [vertical referencing functionalities of xDEM](https://xdem.readthedocs.io/e can help define a 3D CRS. ``` +(bounds)= ### Bounds Bounds define the spatial extent of geospatial data, composed of the "left", "right", "bottom" and "top" coordinates. The {attr}`~geoutils.Raster.bounds` of a raster or a vector is a {class}`rasterio.coords.BoundingBox` object: ```{code-cell} ipython3 -# Show bounds attribute of a raster and vector -print(rast.bounds) -print(vect.bounds) +# Show bounds attribute of raster +rast.bounds +``` +```{code-cell} ipython3 +# Show bounds attribute of vector +vect.bounds ``` ```{note} @@ -103,15 +115,20 @@ reliable computation of extents between CRSs. ```{code-cell} ipython3 # Print raster footprint rast.get_footprint_projected(rast.crs) +``` +```{code-cell} ipython3 # Plot vector footprint vect.get_footprint_projected(vect.crs).plot() ``` ### Grid (only for rasters) -A raster's grid origin and resolution are defined by its geotransform attribute, {attr}`~geoutils.Raster.transform`, and its shape by the data array shape. -From it are derived the resolution {attr}`~geoutils.Raster.res`, the 2D raster shape -{attr}`~geoutils.Raster.shape` made up of its {attr}`~geoutils.Raster.height` and {attr}`~geoutils.Raster.width`. +A raster's grid origin and resolution are defined by its geotransform attribute, {attr}`~geoutils.Raster.transform`. +Comined with the 2D shape of the data array {attr}`~geoutils.Raster.shape` (and independently of the number of +bands {attr}`~geoutils.Raster.bands`), these two attributes define the georeferenced grid of a raster. + +From it are derived the resolution {attr}`~geoutils.Raster.res`, and {attr}`~geoutils.Raster.height` and +{attr}`~geoutils.Raster.width`, as well as the bounds detailed above in {ref}`bounds`. ```{code-cell} ipython3 # Get raster transform and shape @@ -126,7 +143,8 @@ A largely overlooked aspect of a raster's georeferencing is the pixel interpreta [AREA_OR_POINT metadata](https://gdal.org/user/raster_data_model.html#metadata). Pixels can be interpreted either as **"Area"** (the most common) where **the value represents a sampling over the region of the pixel (and typically refers to the upper-left corner coordinate)**, or as **"Point"** -where **the value relates to a point sample (and typically refers to the center of the pixel)**, used often for DEMs. +where **the value relates to a point sample (and typically refers to the center of the pixel)**, the latter often used +for digital elevation models (DEMs). Pixel interpretation is stored as a string in the {attr}`geoutils.Raster.area_or_point` attribute. @@ -140,11 +158,11 @@ interpretation during analysis**, especially for raster–vector–point interfa or re-gridding, and might also be a problem if defined differently when comparing two rasters. ```{important} -By default, **pixel interpretation can induce a half-pixel shift during raster–point interfacing** +By default, **pixel interpretation induces a half-pixel shift during raster–point interfacing for a "Point" interpretation** (mirroring [GDAL's default ground-control point behaviour](https://trac.osgeo.org/gdal/wiki/rfc33_gtiff_pixelispoint)), -but only **raises a warning for raster–raster operations**. +but only **raises a warning for raster–raster operations** if interpretations differ. -This behaviour can be modified at the package-level by using [GeoUtils' configuration parameters]() +This behaviour can be modified at the package-level by using GeoUtils' {ref}`config` `shift_area_or_point` and `warns_area_or_point`. ``` @@ -172,7 +190,7 @@ The metric system returned can be either "universal" (zone of the Universal Tran ```{code-cell} ipython3 # Get local metric CRS -print(rast.get_metric_crs()) +rast.get_metric_crs() ``` ### Re-set georeferencing metadata @@ -183,7 +201,6 @@ To add? {func}`geoutils.Raster.set_transform` {func}`geoutils.Raster.set_area_or_point` - ### Ccoordinates to indexes (only for rasters) Raster grids are notoriously unintuitive to manipulate on their own due to the Y axis being inverted and stored as first axis. diff --git a/doc/source/geotransformations.md b/doc/source/geotransformations.md index ce1c22e0..5d8a3ea5 100644 --- a/doc/source/geotransformations.md +++ b/doc/source/geotransformations.md @@ -43,7 +43,7 @@ For rasters, it can be useful to use {func}`~geoutils.Raster.reproject` in the s for instance when downsampling to a new resolution {attr}`~geoutils.Raster.res`. ```{tip} -Due to the loss information when re-gridding, it is important to **minimize the number of reprojections during the +Due to the loss of information when re-gridding, it is important to **minimize the number of reprojections during the analysis of rasters** (performing only one, if possible). For the same reason, when comparing vectors and rasters in different CRSs, it is usually **better to reproject the vector with no loss of information, which is the default behaviour of GeoUtils in raster–vector–point interfacing**. @@ -82,12 +82,13 @@ rast_reproj = rast.reproject( f, ax = plt.subplots(1, 2) ax[0].set_title("Before reprojection") -rast.plot(ax=ax[0]) +rast.plot(ax=ax[0], cmap="gray", add_cbar=False) vect.plot(rast, ax=ax[0], ec="k", fc="none") ax[1].set_title("After reprojection") -rast_reproj.plot(ax=ax[1]) +rast_reproj.plot(ax=ax[1], cmap="gray", add_cbar=False) vect_reproj.plot(ax=ax[1], ec="k", fc="none") _ = ax[1].set_yticklabels([]) +plt.tight_layout() ``` ```{note} @@ -101,13 +102,16 @@ We can also simply pass another raster as reference to reproject to match the sa and resolution: ```{code-cell} ipython3 +--- +mystnb: + output_stderr: warn +--- # Reproject vector to CRS of raster by simply passing the raster rast_reproj2 = rast.reproject(rast2) ``` -As shown above, if the rasters have different {ref}`Pixel interpretation`, GeoUtils -raises a warning to ensure this is intended. -This warning can be turned off at the package-level using [GeoUtils' global configuration](). +GeoUtils raises a warning because the rasters have different {ref}`Pixel interpretation`, +to ensure this is intended. This warning can be turned off at the package-level using GeoUtils' {ref}`config`. ```{code-cell} ipython3 :tags: [hide-input] @@ -117,16 +121,17 @@ This warning can be turned off at the package-level using [GeoUtils' global conf f, ax = plt.subplots(1, 3) ax[0].set_title("Raster 1") -rast.plot(ax=ax[0]) +rast.plot(ax=ax[0], cmap="gray", add_cbar=False) vect.plot(rast, ax=ax[0], ec="k", fc="none") ax[1].set_title("Raster 2") -rast2.plot(ax=ax[1], cmap="Reds",) +rast2.plot(ax=ax[1], cmap="Reds", add_cbar=False) vect.plot(rast, ax=ax[1], ec="k", fc="none") ax[2].set_title("Match-ref\nreprojection") -rast_reproj2.plot(ax=ax[2]) +rast_reproj2.plot(ax=ax[2], cmap="gray", add_cbar=False) vect_reproj.plot(ax=ax[2], ec="k", fc="none") _ = ax[1].set_yticklabels([]) _ = ax[2].set_yticklabels([]) +plt.tight_layout() ``` ## Crop or pad @@ -156,23 +161,25 @@ vect_clipped = vect_reproj.crop(rast, clip=True) f, ax = plt.subplots(1, 2) ax[0].set_title("Before clipping") -rast.plot(ax=ax[0]) +rast.plot(ax=ax[0], cmap="gray", add_cbar=False) vect_reproj.plot(ax=ax[0], ec="k", fc="none") ax[1].set_title("After clipping") -rast.plot(ax=ax[1]) +rast.plot(ax=ax[1], cmap="gray", add_cbar=False) vect_clipped.plot(ax=ax[1], ec="k", fc="none") +_ = ax[1].set_yticklabels([]) +plt.tight_layout() ``` ## Translate {func}`geoutils.Raster.translate` or {func}`geoutils.Vector.translate`. -Shifting **modifies the georeferencing of the data by a horizontal offset** without modifying the underlying data, +Translations **modifies the georeferencing of the data by a horizontal offset** without modifying the underlying data, which is especially useful to align the data due to positioning errors. ```{code-cell} ipython3 -# Shift the raster by a certain offset +# Translate the raster by a certain offset rast_shift = rast.translate(xoff=1000, yoff=1000) ``` @@ -183,12 +190,14 @@ rast_shift = rast.translate(xoff=1000, yoff=1000) : code_prompt_hide: "Hide the code for plotting the figure" f, ax = plt.subplots(1, 2) -ax[0].set_title("Before shifting") -rast.plot(ax=ax[0]) +ax[0].set_title("Before translation") +rast.plot(ax=ax[0], cmap="gray", add_cbar=False) vect_clipped.plot(ax=ax[0], ec="k", fc="none") -ax[1].set_title("After shifting") -rast_shift.plot(ax=ax[1]) +ax[1].set_title("After translation") +rast_shift.plot(ax=ax[1], cmap="gray", add_cbar=False) vect_clipped.plot(ax=ax[1], ec="k", fc="none") +_ = ax[1].set_yticklabels([]) +plt.tight_layout() ``` :::{admonition} See also diff --git a/doc/source/index.md b/doc/source/index.md index afc36d91..61733987 100644 --- a/doc/source/index.md +++ b/doc/source/index.md @@ -35,9 +35,9 @@ GeoUtils is built on top of core geospatial packages (Rasterio, GeoPandas) and n (NumPy, Xarray, SciPy) to provide **consistent higher-level functionalities at the interface of raster, vector and point cloud objects** (such as match-reference reprojection, point interpolation or gridding). -It is **tailored to perform quantitative analysis that implicitly understands the intricacies of geospatial data**, +It is **tailored to perform quantitative analysis that implicitly understands the intricacies of geospatial data** (nodata values, projection, pixel interpretation), through **an intuitive object-based API to foster accessibility**, -and **to be computationally scalable** through Dask. +and strives **to be computationally scalable** through Dask. If you are looking to **port your GDAL or QGIS workflow in Python**, GeoUtils is made for you! @@ -127,6 +127,7 @@ analysis_examples/index api cli +config background ``` 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 index dec047c1..3b36a886 100644 --- a/doc/source/raster_vector_point.md +++ b/doc/source/raster_vector_point.md @@ -14,7 +14,8 @@ kernelspec: # 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**. +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] @@ -34,7 +35,7 @@ pyplot.rcParams['font.size'] = 9 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 feature index by default. +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 @@ -95,7 +96,7 @@ the targets are implicitly the valid values of the mask. rasterized_vect.set_mask(rasterized_vect == 0) # Polygonize all non-zero values vect_repolygonized = rasterized_vect.polygonize() -vect_repolygonized.plot(ax="new", column="id", cbar_title="Feature index") +vect_repolygonized.plot(ax="new", column="id", fc="none", cbar_title="Feature index") ``` ## Raster–point operations @@ -152,15 +153,10 @@ Point reduction of a raster is **the estimation of the values at point coordinat 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.value_at_coords` +{func}`geoutils.Raster.reduce_points` ```{code-cell} ipython3 -# 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.value_at_coords(x_coords, y_coords, window=5, reducer_function=np.nanmedian) +vals = rast.reduce_points((x_coords, y_coords), window=5, reducer_function=np.nanmedian) ``` ```{code-cell} ipython3 @@ -172,7 +168,7 @@ vals = rast.value_at_coords(x_coords, y_coords, window=5, reducer_function=np.na 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") +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)) @@ -181,6 +177,7 @@ _ = ax[1].set_yticklabels([]) plt.tight_layout() ``` + ### Raster to points {func}`geoutils.Raster.to_pointcloud` @@ -188,14 +185,54 @@ plt.tight_layout() **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. -### Regular point to raster +```{code-cell} ipython3 +pc = rast.to_pointcloud(subsample=10000) +``` -{func}`geoutils.Raster.from_regular_pointcloud` +```{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. Every point of the point cloud is associated to a -pixel in the raster grid, and the values are set to the raster. The point cloud does not need to contain points for -all grid coordinates of the output raster, as missing pixels are set to nodata values. +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 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 7ce51503..69b4c3cb 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 @@ -4243,7 +4245,9 @@ 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 99731eca..14529e60 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()}", @@ -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,9 @@ 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: