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.
Raster–vector operations#
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).
Show 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"))
# Rasterize the vector features based on their glacier ID number
rasterized_vect = vect.rasterize(rast)
Show 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#
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.
# Create a boolean mask from all vector features
mask = vect.create_mask(rast)
Show 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 Mask
, a georeferenced boolean raster (or optionally, a boolean NumPy array), which
can both be used for indexing or index assignment of a raster.
# Mean of values in the mask
np.mean(rast[mask])
172.66101371277432
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 Mask
,
the targets are implicitly the valid values of the mask.
# Mask 0 values
rasterized_vect.set_mask(rasterized_vect == 0)
# Polygonize all non-zero values
vect_repolygonized = rasterized_vect.polygonize()
Show 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#
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 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.
# 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))
Show 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.
geoutils.Raster.reduce_points()
vals = rast.reduce_points((x_coords, y_coords), window=5, reducer_function=np.nanmedian)
Show 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#
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.
pc = rast.to_pointcloud(subsample=10000)
Show 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#
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 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.
rast_from_pc = gu.Raster.from_pointcloud_regular(pc, transform=rast.transform, shape=rast.shape)
Show 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#
geoutils.PointCloud.grid()
Coming with the next release on point clouds!