Utility classes and functions

Bases: object

Wrapper class for basic raster handling

Attributes:

x: Read-only attribute for the x-ticks of raster grid
y: Read-only attribute for the y-ticks of raster grid
values: Read-only attribute for the raster grid data
path: Read-only attribute for the path to original raster file
tmpfile: Read-only attribute for the path to working raster file
md5: Read-only attribute for the hash of working raster file content
count: Read-only attribute for the number of bands of raster dataset
is_masked: Read-only attribute indicating whether raster has missing data.
shape: Read-only attribute indicating the shape of raster grid
height: Read-only attribute for the number of rows of raster grid
bbox: Read-only attribute for the bounding box of the raster grid
src: Read-only attribute for access to opened dataset handle
width: Read-only attribute for the number of columns of raster grid
dx: Read-only attribute for grid distance in x direction
dy: Read-only attribute for grid distance in y direction
crs: Read-only attribute for raster CRS
nodatavals
transform: Read-only attribute for raster’s transform
dtypes: Read-only attribute for raster’s data type
nodata
xres: Read-only attribute for grid resolution in x direction
yres: Read-only attribute for grid resolution in y direction
resampling_method: Modifiable attribute for stored raster resampling method
chunk_size: Modfiable attribute for stored square raster window size
overlap: Modfiable attribute for stored raster windows overlap amount

Methods

modifying_raster(use_src_meta=True, kwargs)**	Context for modifying raster and saving to disk after.
get_x(window=None)	Get X values from the raster.
get_y(window=None)	Get Y values from the raster.
get_xy(window=None)	Get raster position tuples.
get_values(window=None, band=None, kwargs)**	Get values at all the points in the raster.
get_xyz(window=None, band=None)	Get raster position tuples and values horizontally stacked.
get_multipolygon(zmin=None, zmax=None, window=None, overlap=None, band=1)	Extract multipolygon from raster data.
get_bbox(crs=None, output_type=’polygon’)	Get the raster bounding box.
contourf(…)	Plot a filled contour from the raster data.
tags(i=None)	Get tags set on raster dataset.
read(i, masked=True, kwargs)**	Call underlying dataset read method.
dtype(i)	Return data type from the source raster.
nodataval(i)	Return the no-data value of the source raster.
sample(xy, i)	Call underlying dataset sample method.
close()	Release source raster object.
add_band(values, tags)**	Add a new band of data with specified tags to the raster.
fill_nodata()	Fill no-data points in the raster dataset.
gaussian_filter(kwargs)**	Apply Gaussian filter on the raster data.
average_filter(size, drop_above, drop_below)	Apply average filter on the raster data.
generic_filter(function, kwargs)**	Apply generic filter on the raster data.
mask(shapes, i=None, kwargs)**	Mask raster data by shape.
read_masks(i=None)	Read source raster masks.
warp(dst_crs, nprocs=-1)	Reproject raster data.
resample(scaling_factor, resampling_method=None)	Resample raster data.
save(path)	Save-as raster data to the provided path.
clip(geom)	Clip raster data by provided shape.
adjust(geom=None, inside_min=-np.inf, outside_max=np.inf, cond=None)	Modify raster values based on constraints and shape.
get_contour(level, window)	Calculate contour of specified level on raster data.
get_channels(level=0, width=1000, tolerance=None)	Calculate narrow areas based on input level and width.
iter_windows(chunk_size=None, overlap=None)	Return raster view windows based on chunk size and overlap.
get_window_data(window, masked=True, band=None)	Return raster values based for the input window.
get_window_bounds(window)	Return window bounds.
get_window_transform(window)	Return window transform for the input raster

Notes

This class makes use of property and descriptor type attributes. Some of the properties even point to the value set by the descriptors. It’s important to not that the properties are set lazily. It’s also noteworthy to mention that this class is currently not picklable due temporary file and file-like object attributes.

Raster data manipulator.

Parameters:

pathstr or os.PathLike: Path to a raster image to work on (.tiff or .nc).
crsstr or CRS, default=None: CRS to use and override input raster data with. Note that no transformation takes place.
chunk_sizeint or None, default=None: Square window size to be used for data chunking
overlapint or None , default=None: Overlap size for calculating chunking windows on the raster

Attributes:

bbox: Read-only attribute for the bounding box of the raster grid
chunk_size: Modfiable attribute for stored square raster window size
count: Read-only attribute for the number of bands of raster dataset
crs: Read-only attribute for raster CRS
dtypes: Read-only attribute for raster’s data type
dx: Read-only attribute for grid distance in x direction
dy: Read-only attribute for grid distance in y direction
height: Read-only attribute for the number of rows of raster grid
is_masked: Read-only attribute indicating whether raster has missing data.
md5: Read-only attribute for the hash of working raster file content
nodata
nodatavals
overlap: Modfiable attribute for stored raster windows overlap amount
path: Read-only attribute for the path to original raster file
resampling_method: Modifiable attribute for stored raster resampling method
shape: Read-only attribute indicating the shape of raster grid
src: Read-only attribute for access to opened dataset handle
tmpfile: Read-only attribute for the path to working raster file
transform: Read-only attribute for raster’s transform
values: Read-only attribute for the raster grid data
width: Read-only attribute for the number of columns of raster grid
x: Read-only attribute for the x-ticks of raster grid
xres: Read-only attribute for grid resolution in x direction
y: Read-only attribute for the y-ticks of raster grid
yres: Read-only attribute for grid resolution in y direction

Methods

`add_band`(values, **tags)	Add a new band for values with tags tags to the raster
`adjust`([geom, inside_min, outside_max, cond])	Adjust raster data in-place based on specified shape.
`average_filter`(size[, drop_above, ...])	Apply average(mean) filter on the raster
`clip`(geom)	Clip raster data in-place, outside the specified shape.
`close`()	Delete source object
`contourf`([band, window, axes, vmin, vmax, ...])	Plot filled contour for raster data.
`dtype`(i)	Raster data type
`fill_nodata`()	Fill missing values in the raster in-place
`gaussian_filter`(**kwargs)	Apply Gaussian filter to the raster data in-place
`generic_filter`(function, **kwargs)	Apply Gaussian filter to the raster data in-place
`get_bbox`([crs, output_type])	Calculate the bounding box of the raster.
`get_channels`([level, width, tolerance])	Calculate narrow width polygons based on specified input
`get_contour`(level[, window])	Calculate contour lines for specified data level.
`get_multipolygon`([zmin, zmax, window, ...])	Calculate and return a multipolygon based on the raster data
`get_values`([window, band])	Return the data stored at each point in the raster grid
`get_window_bounds`(window)	Returns west, south, east, north bounds of the window
`get_window_data`(window[, masked, band])	Return the positions and values of raster data for the window
`get_window_transform`(window)	Returns raster's affine transform calculated for the window
`get_x`([window])	Get X positions of the raster grid.
`get_xy`([window])	Get raster positions tuple array
`get_xyz`([window, band])	Return the data stored at each point in the raster grid
`get_y`([window])	Get Y positions of the raster grid.
`iter_windows`([chunk_size, overlap])	Calculates sequence of windows for the raster
`mask`(shapes[, i])	Mask data based on input shapes in-place
`modifying_raster`([use_src_meta])	Context manager for modifying and storing raster data
`nodataval`(i)	Value used for filling no-data points
`read`(i[, masked])	Read the data from raster opened file
`read_masks`([i])	Read existing masks on the raster data
`resample`(scaling_factor[, resampling_method])	Resample raster data in-place based on a scaling factor
`sample`(xy, i)	Get value of the data in specified positions
`save`(path)	Save-as raster dataset to a new location
`tags`([i])	Return a dictionary of dataset or band's tags
`warp`(dst_crs[, nprocs])	Reproject the raster data to specified dst_crs in-place

modifying_raster(use_src_meta: bool = True, **kwargs: Any) → Generator[DatasetReader, None, None]

Context manager for modifying and storing raster data

This is a helper context manager method that handles creating new temporary file and replacing the old one with it when raster data is successfully modified.

Parameters:

use_src_metabool, default=True: Whether or not to copy the metadata of the source raster when creating the new empty raster file
**kwargsdict, optional: Options to be passed as metadata to raster database. These options override values taken from source raster in case use_src_meta is True

Yields:

rasterio.DatasetReader: Handle to the opened dataset on temporary file which will override the old values if no exception occurs during the context

get_x(window: Window | None = None) → ndarray[Any, dtype[float]]

Get X positions of the raster grid.

Parameters:

windowwindows.Window, default=None: The window over which X positions are to be returned

Returns:

np.ndarray: A vector X positions from minimum to the maximum

get_y(window: Window | None = None) → ndarray[Any, dtype[float]]

Get Y positions of the raster grid.

Parameters:

windowwindows.Window, default=None: The window over which Y positions are to be returned

Returns:

np.ndarray: A vector Y positions from minimum to the maximum

get_xy(window: Window | None = None) → ndarray

Get raster positions tuple array

Parameters:

windowwindows.Window, default=None: The window over which positions are to be returned

Returns:

np.ndarray: A n :math:` imes 2` matrix of positions

get_values(window: Window | None = None, band: int | None = None, **kwargs: Any) → ndarray[Any, dtype[float]]

Return the data stored at each point in the raster grid

Parameters:

windowwindows.Window, default=None: The window over which values are to be returned.
bandint, default=None: The band from which the values should be read. If None return data from band 1.
**kwargsdict, optional: Additional arguments to pass to rasterio.DatasetReader.read.

Returns:

np.ndarray: A 2D matrix of values on the raster grid.

get_xyz(window: Window | None = None, band: int | None = None) → ndarray

Return the data stored at each point in the raster grid

Parameters:

windowwindows.Window, default=None: The window over which positions and values are to be returned.
bandint, default=None: The band from which the values should be read. If None return data from band 1.

Returns:

np.ndarray: A n :math:` imes 3` matrix of x, y and values

get_multipolygon(zmin: float | None = None, zmax: float | None = None, window: Window | None = None, overlap: int | None = None, band: int = 1) → MultiPolygon

Calculate and return a multipolygon based on the raster data

Calculates filled contour from raster data between specified limits on the specified band and creates a multipolygon from the filled contour to return.

Parameters:

zminfloat or None, default=None: Lower bound of raster data for filled contour calculation
zmaxfloat or None, default=None: Upper bound of raster data for filled contour calculation
windowwindows.Window or None, default=None: Window over whose data the multipolygon is calculated
overlapint or None, default=None: Overlap used for generating windows if window is not provided
bandint, default=1: Raster band over whose data multipolygon is calculated

Returns:

MultiPolygon: The calculated multipolygon from raster data

get_bbox(crs: str | CRS | None = None, output_type: Literal['polygon', 'bbox'] = 'polygon') → Polygon | Bbox

Calculate the bounding box of the raster.

Parameters:

crsstr or CRS or None, default=None: The CRS in which the bounding box is requested.
output_type{‘polygon’, ‘bbox’}: The label of the return type for bounding box, either a shapely ‘polygon’ or matplotlib ‘bbox’.

Returns:

box or Bbox: The bounding box of the raster.

Raises:

TypeError: If the label of return type is not valid.

contourf(band: int = 1, window: Window | None = None, axes: Axes | None = None, vmin: float | None = None, vmax: float | None = None, cmap: str = 'topobathy', levels: List[float] | None = None, show: bool = False, title: str | None = None, figsize: Tuple[float, float] | None = None, colors: int = 256, cbar_label: str | None = None, norm=None, **kwargs: Any) → Axes

Plot filled contour for raster data.

Parameters:

bandint, default=1: Raster band from which data is used.
windowwindows.Window or None, default=None: Raster window from which data is used.
axesAxes or None, default=None: Matplotlib axes to draw contour on>
vminfloat or None, default=None: Minimum value of the filled contour.
vmaxfloat or None, default=None: Maximum value of the filled contour.
cmapstr, default=’topobathy’: Colormap to use for filled contour.
levelslist of float or None, default=None: Prespecified list of contour levels.
showbool, default=False: Whether to show the contour on creation or not.
titlestr or None, default=None: Title used on the axes of the contour
figsizetuple of float or None, default=None: Figure size used for the contour figure
colorsint, default=256: Contour colors associated with levels
cbar_labelstr or None, default=None: Label of the colorbar
normNormalize or None, default=None: Normalizer object
**kwargsdict, optional: Keyword arguments passed to the matplotlib contourf() function

Returns:

Axes: Axes object from matplotlib library that holds onto the contour plot object

tags(i: int | None = None) → Dict[str, str]

Return a dictionary of dataset or band’s tags

Parameters:

iint or None, default=None: The band from which the tags are read.

Returns:

dict: Dictionary of tags

read(i: int, masked: bool = True, **kwargs: Any) → ndarray[Any, dtype[float]]

Read the data from raster opened file

Parameters:

iint: The index of the band to read the data from
maskedbool, default=True: Whether or not to return a masked array
**kwargsdict, optional: Additional keyword arguments passed to rasterio read()

Returns:

ndarray: Array of raster data

dtype(i: int) → ndarray[Any, dtype[float]]

Raster data type

Parameters:

iint: The index of the band to read the data from

Returns:

Any: Data type of raster values

nodataval(i: int) → float

Value used for filling no-data points

Parameters:

iint: The index of the band to read the data from

Returns:

float: The value to be used for points that have missing data

sample(xy: Iterable, i: int) → ndarray[Any, dtype[float]]

Get value of the data in specified positions

Parameters:

xyiterable: Pairs of xy coordinates for which data is retrieved
iint: The index of the band to read the data from

Returns:

ndarray: Array of values for specified input positions

close() → None: Delete source object

add_band(values: ndarray[Any, dtype[float]], **tags: Any) → int

Add a new band for values with tags tags to the raster

Parameters:

valuesarray-like: The values to be added to the raster, it must have the correct shape as the raster.
**tagsdict, optional: The tags to be added for the new band of data.

Returns:

int: ID of the new band added to the raster.

fill_nodata() → None: Fill missing values in the raster in-place

gaussian_filter(**kwargs: Any) → None

Apply Gaussian filter to the raster data in-place

Parameters:

**kwargsdict, optional: Keyword arguments passed to SciPy gaussian_filter function

Returns:

None

average_filter(size: int | ndarray[Any, dtype[int]], drop_above: float | None = None, drop_below: float | None = None, apply_on_bands: List[int] | None = None) → None

Apply average(mean) filter on the raster

Parameters:

size: int, npt.NDArray[int]: size of the footprint
drop_above: float or None: elevation above which the cells are ignored for averaging
drop_below: float or None: elevation below which the cells are ignored for averaging

Returns:

None

generic_filter(function, **kwargs: Any) → None

Apply Gaussian filter to the raster data in-place

Parameters:

function: callable, LowLevelCallable: Function to be used on the footprint array
**kwargsdict, optional: Keyword arguments passed to SciPy generic_filter function

Returns:

None

mask(shapes: Iterable, i: int | None = None, **kwargs: Any) → None

Mask data based on input shapes in-place

Parameters:

shapesiterable: List of GeoJSON like dict or objects that implement Python geo interface protocol (passed to rasterio.mask.mask).
iint or None, default=None: The index of the band to read the data from.
**kwargsdict, optional: Keyword arguments used to create new raster Dataset.

Returns:

None

read_masks(i: int | None = None) → ndarray[Any, dtype[bool]]

Read existing masks on the raster data

Parameters:

iint or None, default=None

Returns:

np.ndarray or view: Raster band mask from the dataset

warp(dst_crs: CRS | str, nprocs: int = -1) → None

Reproject the raster data to specified dst_crs in-place

Parameters:

dst_crsCRS or str: Destination CRS to which raster must be transformed
nprocsint, default=-1: Number of processors to use for the operation

Returns:

None

resample(scaling_factor: float, resampling_method: str | None = None) → None

Resample raster data in-place based on a scaling factor

Parameters:

scaling_factorfloat: The scaling factor to use for resampling data
resampling_methodstr or None, default=None: Name of the resampling method passed to rasterio.DatasetReader.read

Returns:

None

save(path: str | PathLike) → None

Save-as raster dataset to a new location

Parameters:

pathstr or path-like: The path to which raster must be saved.

Returns:

None

clip(geom: Polygon | MultiPolygon) → None

Clip raster data in-place, outside the specified shape.

Parameters:

geomPolygon or MultiPolygon: Shape used to clip the raster data

Returns:

None

adjust(geom: None | Polygon | MultiPolygon = None, inside_min: float = -inf, outside_max: float = inf, cond: Callable[[ndarray[Any, dtype[float]]], ndarray[Any, dtype[bool]]] | None = None) → None

Adjust raster data in-place based on specified shape.

This method can be used to adjust e.g. raster elevation values based on a more accurate land-mass polygon.

Parameters:

geomNone or Polygon or MultiPolygon: Filled shape to determine which points are considered inside or outside (usually land-mass polygon)
inside_minfloat: The minimum value to truncate raster data that falls inside the specified geom shape
outside_maxfloat: The maximum value to truncate raster data that falls outside the specified geom shape

Returns:

None

get_contour(level: float, window: Window | None = None) → LineString | MultiLineString

Calculate contour lines for specified data level.

This method can be used e.g. to calculated coastline based on raster data.

Parameters:

levelfloat: The level for which contour lines must be calculated
windowwindows.Window or None: The raster window for which contour lines must be calculated

Returns:

LineString or MultiLineString: The contour lines calculated for the specified level

get_channels(level: float = 0, width: float = 1000, tolerance: float | None = None) → Polygon | MultiPolygon

Calculate narrow width polygons based on specified input

By using buffer functionality this method finds narrow regions of the domain. The level specifies at which data level domain polygon should be calculated and width describes the narrow region cut-off. tolerance is used for simplifying the polygon before buffering it to reduce computational cost.

Parameters:

levelfloat, default=0: Reference level to calculate domain polygon for narrow region calculation.
widthfloat, default=1000: Cut-off used for designating narrow regions.
tolerancefloat or None, default=None: Tolerance used for simplifying domain polygon.

Returns:

Polygon or MultiPolygon: The calculated narrow regions based on raster data

iter_windows(chunk_size: int | None = None, overlap: int | None = None) → Generator[Window, None, None]

Calculates sequence of windows for the raster

This method calculates the sequence of square windows for the raster based on the provided chunk_size and overlap.

Parameters:

chunk_sizeint or None , default=None: Square window size to be used for data chunking
overlapint or None , default=None: Overlap size for calculating chunking windows on the raster

Yields:

windows.Window: Calculated square window on raster based on the window size and windows overlap values.

get_window_data(window: Window, masked: bool = True, band: int | None = None) → Tuple[ndarray[Any, dtype[float]], ndarray[Any, dtype[float]], ndarray[Any, dtype[float]]]

Return the positions and values of raster data for the window

Returns:

tuple of np.ndarray: The triple of x-ticks, y-ticks and data on the raster grid

get_window_bounds(window: Window) → Tuple[int, int, int, int]

Returns west, south, east, north bounds of the window

Returns:

tuple of int: West, south, east, north bounds of the window

get_window_transform(window: Window) → None | Affine

Returns raster’s affine transform calculated for the window

Returns:

Affine: Affine transform matrix for specified window

property x: ndarray[Any, dtype[float]]

Read-only attribute for the x-ticks of raster grid

This is a read-only property that returns the same results as get_x method

property y: ndarray[Any, dtype[float]]

Read-only attribute for the y-ticks of raster grid

This is a read-only property that returns the same results as get_y method

property values: ndarray[Any, dtype[float]]

Read-only attribute for the raster grid data

This is a read-only property that returns the same results as get_values method

property path: Path: Read-only attribute for the path to original raster file

property tmpfile: Path: Read-only attribute for the path to working raster file

property md5: str

Read-only attribute for the hash of working raster file content

This is a read-only property that is recalculated every time from the content of the temporary working raster file.

property count: int: Read-only attribute for the number of bands of raster dataset

property is_masked: bool

Read-only attribute indicating whether raster has missing data.

This is a read-only property that indicates whether or not the raster has missing data points. The value of this property is recalculated every time the property is retrieved.

property shape: Tuple[int, int]: Read-only attribute indicating the shape of raster grid

property height: int: Read-only attribute for the number of rows of raster grid

property bbox: Polygon | Bbox

Read-only attribute for the bounding box of the raster grid

This is a read-only property that returns the same results as get_bbox method

property src: DatasetReader: Read-only attribute for access to opened dataset handle

property width: int: Read-only attribute for the number of columns of raster grid

property dx: float: Read-only attribute for grid distance in x direction

property dy: float: Read-only attribute for grid distance in y direction

property crs: CRS

Read-only attribute for raster CRS

This read-only property returns CRS as a pyproj.CRS type not rasterio.CRS.

property transform: Affine: Read-only attribute for raster’s transform

property dtypes: Any: Read-only attribute for raster’s data type

property xres: float

Read-only attribute for grid resolution in x direction

This read-only property returns the same value as dx

property yres: float

Read-only attribute for grid resolution in y direction

This read-only property returns the same value as -dy

property resampling_method: Resampling: Modifiable attribute for stored raster resampling method

property chunk_size: int: Modfiable attribute for stored square raster window size

property overlap: int: Modfiable attribute for stored raster windows overlap amount

ocsmesh.utils.must_be_euclidean_mesh(func)

ocsmesh.utils.mesh_to_tri(mesh): mesh is a jigsawpy.jigsaw_msh_t() instance.

ocsmesh.utils.cleanup_isolates(mesh)

ocsmesh.utils.cleanup_duplicates(mesh)

Cleanup duplicate nodes and elements

Notes

Elements and nodes are duplicate if they fully overlapping (not partially)

ocsmesh.utils.cleanup_folded_bound_el(mesh): delete all boundary elements whose nodes (all 3) are boundary nodes

ocsmesh.utils.put_edge2(mesh)

ocsmesh.utils.geom_to_multipolygon(mesh)

ocsmesh.utils.get_boundary_segments(mesh)

ocsmesh.utils.get_mesh_polygons(mesh)

ocsmesh.utils.repartition_features(linestring: LineString, max_verts: int)

ocsmesh.utils.transform_linestring(linestring: LineString, target_size: float)

ocsmesh.utils.needs_sieve(mesh, area=None)

ocsmesh.utils.put_id_tags(mesh)

ocsmesh.utils.finalize_mesh(mesh, sieve_area=None)

ocsmesh.utils.remesh_small_elements(opts, geom, mesh, hfun): This function uses all the inputs for a given jigsaw meshing process and based on that finds and fixes tiny elements that might occur during initial meshing by iteratively remeshing

ocsmesh.utils.sieve(mesh, area=None): A mesh can consist of multiple separate subdomins on as single structure. This functions removes subdomains which are equal or smaller than the provided area. Default behaviours is to remove all subdomains except the largest one.

ocsmesh.utils.sort_edges(edges)

ocsmesh.utils.index_ring_collection(mesh)

ocsmesh.utils.outer_ring_collection(mesh)

ocsmesh.utils.inner_ring_collection(mesh)

ocsmesh.utils.get_multipolygon_from_pathplot(ax)

ocsmesh.utils.signed_polygon_area(vertices)

ocsmesh.utils.vertices_around_vertex(mesh)

ocsmesh.utils.get_surrounding_elem_verts(mesh, in_vert): Find vertices of elements connected to input vertices

ocsmesh.utils.get_lone_element_verts(mesh): Also, there might be some dangling triangles without neighbors, which are also missed by path.contains_point()

ocsmesh.utils.get_verts_in_shape(mesh: jigsaw_msh_t, shape: box | Polygon | MultiPolygon, from_box: bool = False, num_adjacent: int = 0) → Sequence[int]

ocsmesh.utils.select_adjacent(mesh, in_indices, num_layers)

ocsmesh.utils.get_incident_edges(mesh, *args, **kwargs)

ocsmesh.utils.get_cross_edges(mesh, *args, **kwargs)

ocsmesh.utils.clip_mesh_by_shape(mesh: jigsaw_msh_t, shape: box | Polygon | MultiPolygon, use_box_only: bool = False, fit_inside: bool = True, inverse: bool = False, in_place: bool = False, check_cross_edges: bool = False, adjacent_layers: int = 0) → jigsaw_msh_t

ocsmesh.utils.remove_mesh_by_edge(mesh: jigsaw_msh_t, edges: Sequence[Tuple[int, int]], in_place: bool = False) → jigsaw_msh_t

ocsmesh.utils.clip_mesh_by_vertex(mesh: jigsaw_msh_t, vert_in: Sequence[int], can_use_other_verts: bool = False, inverse: bool = False, in_place: bool = False) → jigsaw_msh_t

ocsmesh.utils.get_mesh_edges(mesh, *args, **kwargs)

ocsmesh.utils.calculate_tria_areas(mesh, *args, **kwargs)

ocsmesh.utils.calculate_edge_lengths(mesh, *args, **kwargs)

ocsmesh.utils.elements(mesh, *args, **kwargs)

ocsmesh.utils.faces_around_vertex(mesh, *args, **kwargs)

ocsmesh.utils.get_boundary_edges(mesh): Find internal and external boundaries of mesh

ocsmesh.utils.get_pinched_nodes(mesh): Find nodes through which fluid cannot flow

ocsmesh.utils.has_pinched_nodes(mesh)

ocsmesh.utils.cleanup_pinched_nodes(mesh)

ocsmesh.utils.interpolate(src: jigsaw_msh_t, dst: jigsaw_msh_t, **kwargs)

ocsmesh.utils.interpolate_euclidean_mesh_to_euclidean_mesh(src: jigsaw_msh_t, dst: jigsaw_msh_t, method='linear', fill_value=nan)

ocsmesh.utils.interpolate_euclidean_grid_to_euclidean_mesh(src: jigsaw_msh_t, dst: jigsaw_msh_t, bbox=None, kx=3, ky=3, s=0)

ocsmesh.utils.tricontourf(mesh, ax=None, show=False, figsize=None, extend='both', colorbar=False, **kwargs)

ocsmesh.utils.triplot(mesh, axes=None, show=False, figsize=None, color='k', linewidth=0.07, **kwargs)

ocsmesh.utils.reproject(mesh: jigsaw_msh_t, dst_crs: str | CRS)

ocsmesh.utils.limgrad(mesh, dfdx, imax=100): See https://github.com/dengwirda/mesh2d/blob/master/hjac-util/limgrad.m for original source code.

ocsmesh.utils.msh_t_to_grd(msh: jigsaw_msh_t) → Dict

ocsmesh.utils.grd_to_msh_t(_grd: Dict) → jigsaw_msh_t

ocsmesh.utils.msh_t_to_2dm(msh: jigsaw_msh_t)

ocsmesh.utils.sms2dm_to_msh_t(_sms2dm: Dict) → jigsaw_msh_t

ocsmesh.utils.msh_t_to_utm(mesh, *args, **kwargs)

ocsmesh.utils.estimate_bounds_utm(bounds, crs='EPSG:4326')

ocsmesh.utils.estimate_mesh_utm(mesh, *args, **kwargs)

ocsmesh.utils.get_polygon_channels(polygon, width, simplify=None, join_style=3)

ocsmesh.utils.merge_msh_t(*mesh_list, out_crs='EPSG:4326', drop_by_bbox=True, can_overlap=True, check_cross_edges=False)

ocsmesh.utils.add_pool_args(func)

ocsmesh.utils.drop_extra_vertex_from_line(lstr: LineString) → LineString

ocsmesh.utils.drop_extra_vertex_from_polygon(mpoly: Polygon | MultiPolygon) → MultiPolygon

ocsmesh.utils.remove_holes(poly: Polygon | MultiPolygon) → Polygon | MultiPolygon

Remove holes from the input polygon(s)

Given input Polygon or MultiPolygon, remove all the geometric holes and return a new shape object.

Parameters:

polyPolygon or MultiPolygon: The input shape from which the holes are removed

Returns:

Polygon or MultiPolygon: The resulting (multi)polygon after removing the holes

See also

remove_holes: Remove all the whole from the input shape

Notes

For a Polygon with no holes, this function returns the original object. For MultiPolygon with no holes, the return value is a unary_union of all the underlying `Polygon`s.

If rel_size=1 is specified the result is the same as remove_holes function, except for the additional cost of calculating the areas.

ocsmesh.utils.get_element_size_courant(characteristic_velocity_magnitude: float | ndarray[Any, dtype[float]], timestep: float, target_courant: float = 1) → float | ndarray[Any, dtype[float]]

Calculate the element size based on the specified Courant number.

Calculate the element size based on the specified Courant number and input value for timestep and characteristic velocity

Parameters:

target_courantfloat: The Courant number to be achieved by the calculated element size
characteristic_velocity_magnitudefloat or array of floats: Magnitude of total velocity used for element size calculation (\(\frac{m}{sec}\))
timestepfloat: Timestep size (\(seconds\)) to

Returns:

float or array of floats: The calculated element size(s) to achieve the given Courant #

ocsmesh.utils.can_velocity_be_approximated_by_linear_wave_theory(depth: float | ndarray[Any, dtype[float]], wave_amplitude: float = 2) → bool | ndarray[Any, dtype[bool]]

Checks whether the particle velocity can be appoximated.

Based on the input depth, checks whether or not the velocity can be approximated from the linear wave theory

Parameters:

depthfloat or array of float: Depth of the point for which the approximation validation is checked
wave_amplitudefloat, default=2: Free surface elevation (\(meters\)) from the reference (i.e. wave height)

Returns:

bool or array of bool: Whether or not the value at given input depth can be approximated

Notes

Linear wave theory approximation breaks down when \(\nu \sim h\) overland. So this method just returns whether depth is below or over depth of wave_amplitude magnitude.

References

Based on OceanMesh2D approximation method https://doi.org/10.13140/RG.2.2.21840.61446/2.

ocsmesh.utils.estimate_particle_velocity_from_depth(depth: float | ndarray[Any, dtype[float]], wave_amplitude: float = 2) → float | ndarray[Any, dtype[float]]

Approximate particle velocity magnitude based on depth of water

Estimate the value of particle velocity magnitude based on the linear wave theory as \(\left|u\right| = \nu \sqrt{\frac{g}{h}}\) for ocean and \(\left|u\right| = \sqrt{gH}\) for overland region where \(\nu \sim h\) so instead of linear wave theory we take \(H \approx \nu\)

Parameters:

depthfloat or array of floats: The depth of still water (e.g. from DEM)
wave_amplitudefloat: Wave amplitude for approximation as defined by linear wave theory

Returns:

float or array of floats: Estimated water particle velocity

References

Based on OceanMesh2D approximation method https://doi.org/10.13140/RG.2.2.21840.61446/2.

ocsmesh.utils.approximate_courant_number_for_depth(depth: float | ndarray[Any, dtype[float]], timestep: float, element_size: float | ndarray[Any, dtype[float]], wave_amplitude: float = 2) → float | ndarray[Any, dtype[float]]

Approximate the Courant number for given depths

Approximate the value of Courant number for the input depth, timestep and element size. The velocity is approximated based on the input depth.

Parameters:

depthfloat or array of floats
timestepfloat: Timestep size (\(seconds\)) to
element_sizefloat or array of floats: Element size(s) to use for Courant number calculation. Must be scalar otherwise match the dimension of depth
wave_amplitudefloat, default=2: Free surface elevation (\(meters\)) from the reference (i.e. wave height)

Returns:

float or array of floats: The approximated Courant number for each input depth

ocsmesh.utils.create_rectangle_mesh(nx, ny, holes, x_extent=None, y_extent=None, quads=None)

Note:

x = x-index y = y-index

holes or quads count starting at 1 from bottom corner square

node-index(node-id)

25(26) 29(30)

5 ————*

| | | |

4 ————*

| | | |

3 ————*

| | | |

2 ————*

| | # | |

1 ————*

| | | |

0 ————*
0(1) 4(5)

0 1 2 3 4

ocsmesh.utils.raster_from_numpy(filename, data, mgrid, crs=<Geographic 2D CRS: EPSG:4326> Name: WGS 84 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: World. - bounds: (-180.0, -90.0, 180.0, 90.0) Datum: World Geodetic System 1984 ensemble - Ellipsoid: WGS 84 - Prime Meridian: Greenwich) → None

ocsmesh.utils.msht_from_numpy(coordinates, *, edges=None, triangles=None, quadrilaterals=None, values=None, crs=<Geographic 2D CRS: EPSG:4326> Name: WGS 84 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: World. - bounds: (-180.0, -90.0, 180.0, 90.0) Datum: World Geodetic System 1984 ensemble - Ellipsoid: WGS 84 - Prime Meridian: Greenwich)

ocsmesh.utils.shape_to_msh_t(shape: Polygon | MultiPolygon) → jigsaw_msh_t

Calculate vertex-edge representation of polygon shape

Calculate jigsawpy vertex-edge representation of the input shapely shape.

Parameters:

shapePolygon or MultiPolygon: Input polygon for which the vertex-edge representation is to be calculated

Returns:

jigsaw_msh_t: Vertex-edge representation of the input shape

Raises:

NotImplementedError

ocsmesh.utils.shape_to_msh_t_2(shape: Polygon | MultiPolygon | GeoDataFrame | GeoSeries) → jigsaw_msh_t

Triangulate the input shape, with additional points provided

Use triangle package to triangulate the input shape. If provided, use the input additional points as well. The input opts controls how triangle treats the inputs. See the documentation for triangle more information

Parameters:

shapePolygon or MultiPolygon or GeoDataFrame or GeoSeries: Input polygon to triangulate
aux_ptsnumpy array or Point or MultiPoint or GeoDataFrame or GeoSeries: Extra points to be used in the triangulation
optsstring, default=’p’: Options for controlling triangle package
type_tint, default=2: Options for controlling the workflow of the triangulation (see Notes)

Returns:

jigsaw_msh_t: Generated triangulation

Raises:

ValueError: If input shape is invalid or the point input cannot be used to generate point shape

Notes

type_t was added to encompass polygons that fail to be triangulated the default type_t = 2 is equal to the standard triangulate_polygon prior to the implementation of the type_t variable. type_t = 1 uses shape_to_msh_t instead of shape_to_msh_t_2 type_t = 3 places the point for the hole at the corned of the polygon

it should never be used in case there are multiple holes

Triangulate the input shape smoothly

Parameters:

shapePolygon or MultiPolygon or GeoDataFrame or GeoSeries: Input polygon to triangulate
min_int_anginteger, default=30: Minimal internal angle for triangulation

Returns:

jigsaw_msh_t: Generated triangulation

ocsmesh.utils.filter_el_by_area(mesh, l_limit=0, u_limit=1e-07)

Get the elements that have area within limits

Parameters:

meshocsmesh.mesh.mesh.EuclideanMesh2D: Input mesh to be filtered
l_limitfloat: lower area limit
u_limitfloat: lower area limit

Returns:

dict: Mapping between selected element IDs and associated node Ids

Notes

Default l_limit and u_limit area optimum for finding small area elements

ocsmesh.utils.create_patch_mesh(mesh_w_problem, sel_el_dict, mesh_for_patch, buffer_size=0.01, crs=<Geographic 2D CRS: EPSG:4326> Name: WGS 84 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: World. - bounds: (-180.0, -90.0, 180.0, 90.0) Datum: World Geodetic System 1984 ensemble - Ellipsoid: WGS 84 - Prime Meridian: Greenwich)

Extracts “patches” from a mesh (mesh_for_patch) based on a buffer of size “buffer_size” around elements (“sel_el_dict”) from another mesh (“mesh_w_problem”)

Parameters:

mesh_w_problemocsmesh.mesh.mesh.EuclideanMesh2D: mesh with erroneous elements
sel_el_dictdict: dict with elements for which a buffer will be created
mesh_for_patchocsmesh.mesh.mesh.EuclideanMesh2D: mesh that will be used for creating the patches

Returns:

jigsaw_msh_t: Patch mesh

Notes

The hfun.2dm can be used as mesh_for_patch Default buffer_size is optimum for finding small area elements

ocsmesh.utils.clip_mesh_by_mesh(mesh_to_be_clipped: ~jigsawpy.msh_t.jigsaw_msh_t, mesh_clipper: ~jigsawpy.msh_t.jigsaw_msh_t, inverse: bool = True, fit_inside: bool = False, check_cross_edges: bool = True, adjacent_layers: int = 2, crs=<Geographic 2D CRS: EPSG:4326> Name: WGS 84 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: World. - bounds: (-180.0, -90.0, 180.0, 90.0) Datum: World Geodetic System 1984 ensemble - Ellipsoid: WGS 84 - Prime Meridian: Greenwich, buffer_size=None) → jigsaw_msh_t

Clip a mesh (“mesh_to_be_clipped”) based on the mesh extent of another mesh (“mesh_clipper”)

Parameters:

mesh_to_be_clippedjigsawpy.msh_t.jigsaw_msh_t: mesh to be clipped
mesh_clipperjigsawpy.msh_t.jigsaw_msh_t: mesh that will be used to clip

Returns:

jigsaw_msh_t: clipped mesh

ocsmesh.utils.create_mesh_from_mesh_diff(domain: ~jigsawpy.msh_t.jigsaw_msh_t | ~shapely.geometry.polygon.Polygon | ~shapely.geometry.multipolygon.MultiPolygon | ~geopandas.geodataframe.GeoDataFrame | ~geopandas.geoseries.GeoSeries | list, mesh_1: ~jigsawpy.msh_t.jigsaw_msh_t, mesh_2: ~jigsawpy.msh_t.jigsaw_msh_t, crs=<Geographic 2D CRS: EPSG:4326> Name: WGS 84 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: World. - bounds: (-180.0, -90.0, 180.0, 90.0) Datum: World Geodetic System 1984 ensemble - Ellipsoid: WGS 84 - Prime Meridian: Greenwich, min_int_ang=None, buffer_domain=0.001, hfun_mesh=None) → jigsaw_msh_t

Create a triangulation for the area correspondent to the difference between “mesh_1” and “mesh_2” for the extent defined by “mesh_for_domain”

Parameters:

domainjigsawpy.msh_t or Polygon or: MultiPolygon or GeoDataFrame or GeoSeries extent of entire domain (mesh_1+mesh_2+gap)
mesh_1jigsawpy.msh_t.jigsaw_msh_t: mesh_1
mesh_2jigsawpy.msh_t.jigsaw_msh_t: mesh_2
min_int_ang: Minimal internal angle for triangulation default = None, i.e. direct connection between vertices recommended = 30

Returns:

jigsaw_msh_t: mesh for the area between mesh_1 and mesh_2 within mesh_for_domain

ocsmesh.utils.merge_neighboring_meshes(*all_msht)

Combine meshes whose boundaries match Adapted from:

https://github.com/SorooshMani-NOAA/river-in-mesh/blob/main/river_in_mesh/mesh

Parameters:

all_mshtjigsawpy.msh_t.jigsaw_msh_t: mesh to be merged, sections of boundaries of these mesh must overlap (i.e. they must share boundary nodes)

Returns:

jigsaw_msh_t: merged mesh

ocsmesh.utils.fix_small_el(mesh_w_problem: jigsaw_msh_t, mesh_for_patch: jigsaw_msh_t, l_limit: float = 0.0, u_limit: float = 1e-07, buffer_size: float = 0.01, adjacent_layers: int = 2) → jigsaw_msh_t

Fix spurious small elements (<u_limit=1e-7)

Parameters:

mesh_w_problemocsmesh.mesh.mesh.EuclideanMesh2D: mesh that has small areas
mesh_for_patchocsmesh.mesh.mesh.EuclideanMesh2D: mesh that will be used for fixing the mesh_w_problem

Returns:

jigsaw_msh_t: fixed mesh with no small area elements

Notes

Other optimal attributes were determined based on testing and they can be changed as needed:

l_limit=0,u_limit=1e-7, buffer_size=0.001, adjacent_layers=2

ocsmesh.utils.merge_overlapping_meshes(all_msht: list, adjacent_layers: int = 0, buffer_size: float = 0.005, buffer_domain: float = 0.01, min_int_ang: int = 30, hfun_mesh=None, crs=<Geographic 2D CRS: EPSG:4326> Name: WGS 84 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: World. - bounds: (-180.0, -90.0, 180.0, 90.0) Datum: World Geodetic System 1984 ensemble - Ellipsoid: WGS 84 - Prime Meridian: Greenwich) → jigsaw_msh_t

Combine meshes whose boundaries match

Parameters:

*all_mshtjigsaw_msh_t: list of meshes to be merged, the later on the list the higher the priority
adjacent_layersint: # of elements of the low priority mesh that will be deleted passed the overlap
buffer_sizefloat: size of the buffer that will be added around the high priority mesh when carving out the low priority mesh (sometimes this necessary if the mesh is too narrow)
buffer_domainfloat: size of the buffer that will be added around the entire mesh domain. This is necessary for preventing slivers at the mesh boundary intersections
min_int_angint: Minimal internal angle for triangulation default = 30, i.e. triangles will have internal angles of at least 30 deg For no smoothness min_int_ang = None

Returns:

jigsaw_msh_t: final merged mesh

ocsmesh.utils.calc_el_angles(msht)

Adapted from: https://github.com/SorooshMani-NOAA/river-in-mesh/tree/main/river_in_mesh

Calculates the internal angle of each node for element (triangular or quadrangular)

Parameters:

mshtjigsawpy.msh_t.jigsaw_msh_t

Returns:

np.array: internal angles of each element
Notes

ocsmesh.utils.order_mesh(msht, crs=<Geographic 2D CRS: EPSG:4326> Name: WGS 84 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: World. - bounds: (-180.0, -90.0, 180.0, 90.0) Datum: World Geodetic System 1984 ensemble - Ellipsoid: WGS 84 - Prime Meridian: Greenwich) → jigsaw_msh_t

Order mesh nodes counterclockwise (triangles and quads) based on the coordinates

Parameters:

mshtjigsawpy.msh_t.jigsaw_msh_t: mesh with nodes out of order

Returns:

jigsawpy.msh_t.jigsaw_msh_t

np.array: mesh whose nodes within each element are oriented counterclockwise

ocsmesh.utils.quads_from_tri(msht) → jigsaw_msh_t

Partially adapted from:

https://stackoverflow.com/questions/69605766/find-position-of-duplicate-elements-in-list

Combines all triangles that share vertices that are not right angles. right angles is defined as the internal angle closest to 90deg

Parameters:

mshtjigsawpy.msh_t.jigsaw_msh_t: triangular mesh

Returns:

jigsawpy.msh_t.jigsaw_msh_t: quadrangular + triangular mesh

ocsmesh.utils.clip_elements_by_index(msht: jigsaw_msh_t, tria=None, quad=None, inverse: bool = False) → jigsaw_msh_t

Adapted from:

https://github.com/SorooshMani-NOAA/river-in-mesh/tree/main/river_in_mesh/utils

Parameters:

mshtjigsawpy.msh_t.jigsaw_msh_t: mesh to beclipped
tria or quad: array with the element ids to be removed
inverse = default:False

Returns:

jigsaw_msh_t: mesh without skewed elements

ocsmesh.utils.cleanup_skewed_el(mesh: jigsaw_msh_t, lw_bound_tri: float = 1.0, up_bound_tri: float = 175.0, lw_bound_quad: float = 10.0, up_bound_quad: float = 179.0) → jigsaw_msh_t

Removes elements based on their internal angles

Parameters:

mshtjigsawpy.msh_t.jigsaw_msh_t
lw_bound_tridefault=1
up_bound_tridefault=175
lw_bound_quaddefault=10
up_bound_quaddefault=179

Returns:

np.array: internal angles of each element

ocsmesh.utils.cleanup_concave_quads(mesh: jigsaw_msh_t) → jigsaw_msh_t

Removes concave quads that might have been wrongly created during the quadrangulation process

Parameters:

mshtjigsawpy.msh_t.jigsaw_msh_t

Returns:

np.array: internal angles of each element

ocsmesh.utils.quadrangulate_rivermapper_arcs(arcs_shp, _buffer_1: float = 0.001, _buffer_2: float = 0.0001, _batch_size: float = 1000, crs=<Geographic 2D CRS: EPSG:4326> Name: WGS 84 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: World. - bounds: (-180.0, -90.0, 180.0, 90.0) Datum: World Geodetic System 1984 ensemble - Ellipsoid: WGS 84 - Prime Meridian: Greenwich) → jigsaw_msh_t

Creates quads using the RiverMapper diagnostic outputs removes quads with repetitive coords removes all elements that intersect merges all individual quads for each stream into a single mesh removes intersecting quads at junctions

Parameters:

arcs_shp.shp file with the RiverMapper diagnostic arcs
_buffer_1default(0.001), controls the the removal of: overlaping quads within a given stream
_buffer_2default(0.0001), controls the the removal of: overlaping quads at the river junctions
_batch_sizedefault(1000), batch size for merging stream quads

Returns:

jigsaw_msh_t: quad mesh

Notes

Future versions of this function should be parallelized for efficiency

ocsmesh.utils.batched(iterable, n): This function is part of itertools for python 3.12+ This function was added to utils.py to ensure OCSMesh can run on older with python<3.12

ocsmesh.utils.delaunay_within(gdf)

Creates the initial delaunay triangules for a gpd composed of polygons (only). Selects those delaunay triangules that fall within domain.

Parameters:

gdfgpd of polygons

Returns:

gdfgpd of triangulated polygons

ocsmesh.utils.triangulate_shp(gdf)

Fills out the gaps left by the delaunay_within

Parameters:

gdfgpd of polygons

Returns:

gdfgpd of triangulated polygons

ocsmesh.utils.shptri_to_msht(triangulated_shp)

Converts a triangulated shapefile to msh_t

Parameters:

triangulated_shptriangulated gpd

Returns:

jigsaw_msh_t

ocsmesh.utils.triangulate_rivermapper_poly(rm_poly)

Creates triangulated mesh using the RiverMapperoutputs

Parameters:

rm_poly.shp (gpd) file with the RiverMapper outputs

Returns:

jigsaw_msh_t: River Mesh