stompy.grid — Grid reading, writing and manipulating

Modules for generating, converting, plotting, and editing unstructured grids as used in hydrodynamic models (e.g. UnTRIM, SUNTANS, DFlow-FM).

For generating grids, the only stable and robust code for this is in tom.py (triangular orthogonal mesher) and paver.py. tom.py is a command line interface to paver.py. A simple example of calling tom is in stompy/tests/test_tom.sh, and invoking tom.py -h will show the other options available.

For most other grid-related tasks, the best module to use is unstructured_grid.py, as it supports non-triangular meshes (such as mixed triangles/quads), and is actively developed. Grid generation methods built on unstructured_grid.py are in front.py, but these are not stable and generally should not be used.

Note that for any significant amount of grid modification, the CGAL python bindings are essential. A backup pure python implementation is included, but will be orders of magnitude slower and likely less robust numerically. Watch for error messages near the start of using tom.py to see whether there are issues loading CGAL.

Submodules

stompy.grid.depth_connectivity module

depth_connectivity.py —

Assign bathymetry to a grid based on cell-to-cell connectivity derived from a high-resolution DEM.

Adapted from Holleman and Stacey, JPO, 2014.

Primary entry point: edge_depths=edge_connection_depth(g,dem,edge_mask=None,centers=’lowest’)

see end of file

stompy.grid.depth_connectivity.cell_mean_depth(g, dem)

Calculate “true” mean depth for each cell, at the resolution of the DEM. This does not split pixels, though.

stompy.grid.depth_connectivity.edge_connection_depth(g, dem, edge_mask=None, centers='circumcenter')

Return an array g.Nedges() where the selected edges have a depth value corresponding to the minimum elevation at which adjacent cells are hydraulically connected, evaluated on the dem. g: instance of UnstructuredGrid dem: field.SimpleGrid instance, usually GdalGrid edge_mask: bitmask for which edges to calculate, defaults to bounds of dem.

centers controls the reference point for each cell.

‘circumcenter’: use cell circumcenter ‘centroid’: use cell centroid ‘lowest’: use lowest point within the cell.

stompy.grid.depth_connectivity.greedy_edge_mean_to_node(g, orig_node_depth=None, edge_depth=None, n_iter=100)

Return node depths such that the mean of the node depths on each edge approximate the provided edge_mean_depth. The approach is iterative, starting with the largest errors, visiting each edge a max of once.

Still in development, has not been tested.

stompy.grid.depth_connectivity.greedy_edgemin_to_node(g, orig_node_depth, edge_min_depth)

A simple approach to moving edge depths to nodes, when the hydro model (i.e. DFM) will use the minimum of the nodes to set the edge. It sounds roundabout because it is, but there is not a supported way to assign edge depth directly.

For each edge, want to enforce a minimum depth in two sense: 1. one of its nodes is at the minimum depth 2. neither of the nodes are below the minimum depth and.. 3. the average of the two nodes is close to the average DEM depth

of the edge

Not yet sure of how to get all of those. This method focuses on the first point, but in some situations that is still problematic. The 3rd point is not attempted at all, but in DFM would only be relevant for nonlinear edge depths which are possibly not even supported for 3D runs.

stompy.grid.depth_connectivity.min_connection_elevation(ijs, min_depth, max_depth, F)

stompy.grid.depth_connectivity.min_graph_elevation_for_edge(g, dem, j, starts='lowest')

g: unstructured_grid j: edge index dem: a Field subclass which supports extract_tile().

starts:

‘circumcenter’ connections are between voronoi centers ‘centroid’ connections are between cell centroids ‘lowest’ connections are between lowest point in cell

returns: the minimum edge elevation at which the cells adjacent

to j are hydraulically connected. nan if j is not adjacent to two cells (i.e. boundary).

stompy.grid.depth_connectivity.points_inside_poly(points, ijs)

stompy.grid.depth_connectivity.points_to_mask(hull_ijs, nx, ny)

stompy.grid.depth_connectivity.poly_mean_elevation(dem, pnts)

stompy.grid.exact_delaunay module

stompy.grid.front module

stompy.grid.merge_ugrid_subgrids module

stompy.grid.orthogonalize module

stompy.grid.orthomaker module

stompy.grid.paver module

stompy.grid.tom module

stompy.grid.trigrid module

Class for representing and manipulating triangular grids. This has largely been superceded by unstructured_grid.py, but remains here for compatibility with some SUNTANS code, as well as the foundation for tom.py, which is still better than the unstructured_grid.py-based grid generation.

exception stompy.grid.trigrid.NoSuchCellError

Bases: stompy.grid.trigrid.TriGridError

exception stompy.grid.trigrid.NoSuchEdgeError

Bases: stompy.grid.trigrid.TriGridError

class stompy.grid.trigrid.TriGrid(sms_fname=None, tri_basename=None, suntans_path=None, processor=None, suntans_reader=None, tec_file=None, gmsh_basename=None, edges=None, points=None, cells=None, readonly=False)

Bases: object

Ncells()

Nedges()

Npoints()

add_cell(c)

add_edge(nodeA, nodeB, marker=0, cleft=-2, cright=-2, coerce_boundary=None)

returns the number of the edge for cells that are marked -2, this will check to see if a new cell can be made on that side with other unmeshed edges

add_node(P)

angles()

returns [Nc,3] array of internal angles, in radians

animate_scalar(scalar_frames, post_proc=None)

areas()

returns signed area, CCW is positive

boundary_angle(pnt_i)

returns the interior angle in radians, formed by the boundary at the given point

bounds()

build_edge_index()

build_index()

carve_thalweg(depths, threshold, start, mode, max_count=None)

Ensures that there is a path of cells from the given start edge to deep water with all cells of at least threshold depth.

start: edge index

depths and threshold should all be as soundings - i.e. positive

mode is ‘cells’ - cell-centered depths

or ‘edges’ - edge-centered depths

max_count: max number of cells/edges to deepen along the path (starting

at start).

Modifies depths in place.

cell2edges(cell_i)

cell_centers()

cell_divergence_of_edge_flux(edge_flux)

edge_flux is assumed to be depth integrated, but not horizontally integrated - so something like watts per meter

cell_edge_map()

cell2edges for the whole grid return an integer valued [Nc,3] array, where [i,k] is the edge index opposite point self.cells[i,k]

N.B. this is not kept up to date when modifying the grid.

cell_neighbors(cell_id, adjacent_only=0)

return array of cell_ids for neighbors of this cell. here neighbors are defined by sharing a vertex, not just sharing an edge, unless adjacent_only is specified. (in which case it only returns cells sharing an edge)

cells2edge(nc1, nc2)

cells_on_line(xxyy)

Return cells intersecting the given line segment cells are found based on having vertices which straddle the line, and cell centers which are within the segment’s extent

checkpoint()

closest_cell(p, full=0, inside=False)

full=0: return None if the closest point is not in a cell on this subdomain full=1: exhaustively search all cells, even if the nearest point is not on this subdomain

inside: require that the returned cell contains p, otherwise return None

DANGER: this method is not robust! in particular, a nearby but disconnected cell could have a vertex close to the query point, and we’d never see the right cell.

closest_edge(p)

closest_point(p, count=1, boundary=0)

Returns the count closest nodes to p boundary=1: only choose nodes on the boundary.

commit()

created_cell(c)

created_edge(e)

created_node(i)

default_clip = None

delete_cell(c, replace_with=-2, rollback=1)

replace_with: the value to set on edges that used to reference this cell. -2 => leave an internal hole -1 => create an ‘island’

delete_edge(e, rollback=1)

for now, just make it into a degenerate edge specify rollback=0 to skip recording the undo information

delete_node(i, remove_edges=1)

delete_node_and_merge(n)

For a degree 2 node, remove it and make one edge out its two edges. this used to be in paver, but I don’t think there is any reason it can’t be here in trigrid.

delete_unused_nodes()

any nodes which aren’t in any cells or edges will be removed.

deleted_cell(c)

deleted_edge(e)

deleted_node(i)

edge_centers()

edge_index = None

edges_to_polygons(edgemask)

use the edges (possibly masked by given edgemask) to create a shapely.geometry.Polygon() for each top-level polygon, ordered by decreasing area

edges_to_rings(edgemask=None, ccw=1)

using only the edges for which edgemask is true, construct rings. if edgemask is not given, use all of the current edges

if ccw is 1, only non-intersecting ccw rings will be return if ccw is 0, only non-intersecting cw rings will be return

edges_to_rings_and_holes(edgemask)

using only the edges for which edgemask is true, construct polygons with holes. if edgemask is not given, use all of the current edges

This calls edges_to_rings to get both ccw rings and cw rings, and then determines which rings are inside which.

returns a list [ [ outer_ring_nodes, [inner_ring1,inner_ring2,…]], … ]

faces(i)

file_path(conf_name)

find_cell(nodes)

return the cell (if any) that is made up of the given nodes depends on pnt2cells

find_edge(nodes)

from_data(points, edges, cells)

ghost_cells()

Return a bool array, with ghost cells marked True Ghost cells are determined as any cell with an edge that has marker 6

index = None

init_listeners()

interp_cell_to_edge(F)

given a field [Nc,…], linearly interpolate to edges and return [Ne,…] field.

interp_cell_to_node(F)

listen(event, cb)

listener_count = 0

make_edges_from_cells()

merge_edges(e1, e2)

returns the id of the new edge, which for now will always be one of e1 and e2 (and the other will have been deleted

move_node(i, new_pnt)

op_stack = None

op_stack_serial = 10

plot(voronoi=False, line_collection_args={}, all_cells=True, edge_values=None, edge_mask=None, vmin=None, vmax=None, ax=None, clip=None)

vmin: if nan, don’t set an array at all for the edges

clip=[xmin,xmax,ymin,ymax]: additionally mask edges which are not within the given rectangle

edge_values: defaults to the edge marker.

plot_bad_bcs()

plot_edge_marks(edge_mask=None, clip=None)

label edges with c[nc1]-j[j],mark-c[nc2], rotated so it reads in the correct orientation for nc1, nc2 edge_mask should be a boolean array of size Nedges() clip can be a list like matplotlib axis() - [xmin,xmax,ymin,ymax]

plot_nodes(ids=None)

plot_scalar(scalar, pdata=None, clip=None, ax=None, norm=None, cmap=None)

Plot the scalar assuming it sits at the center of the cells (i.e. use the voronoi centers) scalar should be a 1d array, with length the same as the number of cells

to mask out values, set scalar to nan

pnt2cells(pnt_i)

pnt2edges(pnt_i)

pop_op()

push_op(meth, *data, **kwdata)

read_gmsh(gmsh_basename)

reads output from gmsh - gmsh_basename.{nod,ele}

read_sms(fname)

read_suntans(use_cache=1)

read_tecplot(fname)

read_triangle(tri_basename)

reads output from triangle, tri_basename.{ele,poly,node}

recursive_delete(c, renumber=1)

refresh_metadata()

Call this when the cells, edges and nodes may be out of sync with indices and the like.

renumber()

removes duplicate cells and nodes that are not referenced by any cell, as well as cells that have been deleted (==-1)

revert(cp)

scalar_contour(scalar, V=10, smooth=True)

Generate a collection of edges showing the contours of a cell-centered scalar.

V: either an int giving the number of contours which will be evenly spaced over the range of the scalar, or a sequence giving the exact contour values.

smooth: control whether one pass of 3-point smoothing is applied.

returns a LineCollection

select_edges_by_polygon(poly)

set_edge_markers(pnt1, pnt2, marker)

Find the nodes closest to each of the two points, Search for the shortest path between them on the boundary. Set all of those edges’ markers to marker

set_edge_neighbors_from_cells()

shortest_path(n1, n2, boundary_only=0, max_cost=inf)

dijkstra on the edge graph from n1 to n2 boundary_only: limit search to edges on the boundary (have a -1 for cell2)

smooth_scalar(cell_value)

simple method for smoothing a scalar field. note that this is not conservative of anything! and the degree of smoothing is per cell, not per area, so results may be misleading.

it does take care not to corrupt valid values with nans during the smoothing

split_edge(nodeA, nodeB, nodeC)

take the existing edge AC and insert node B in the middle of it

nodeA: index to node on one end of the existing edge nodeB: (i) index to new new node in middle of edge,

(ii) tuple (coords, dict for add_node options) may be extended to allow arbitrary options for point

nodeC: index to node on other end of existing edge

trim_to_left(path)

Given a path, trim all cells to the left of it.

unadd_cell(old_length)

unadd_edge(old_length)

unadd_node(old_length)

undelete_cell(c, nodes, edge_updates)

undelete_edge(e, e_data)

undelete_node(i, p)

unlisten(cb_id)

unmerge_edges(e1, e2, e1data, e2data)

unmodify_edge(e, old_data)

unmove_node(i, orig_val)

updated_cell(c)

updated_edge(e)

updated_node(i)

valid_edges()

returns an array of indices for valid edges - i.e. not deleted

vcenters()

verbose = 0

verify_bc(do_plot=True)

check to make sure that all grid boundaries have a BC set

write_Triangle(basename, boundary_nodes=None)

duplicate some of the output of the Triangle program - particularly the .node and .ele files

note that node and cell indices are taken as 1-based.

if boundary_nodes is supplied, it should be an integer valued array of length Npoints, and give the boundary marker for each node (usually 0 for internal, nonzero for boundary). this can be used to specify a subset of the boundary nodes for a BC in SWAN.

if not specified, boundary markers will be 0 for internal, 1 for external nodes.

write_cells_shp(shpname, cell_mask=None, overwrite=False, fields=None)

fields: a structure array of fields to write out - see wkb2shp

write_contours_shp(shpname, cell_depths, V, overwrite=False)

like write_shp, but collects edges for each depth in V.

write_mat(fn, order='ccw')

write_obj(fname)

Output to alias wavefront

• scales points to fall within [0,10]

write_shp(shpname, only_boundaries=1, edge_mask=None, overwrite=0)

Write some portion of the grid to a shapefile. If only_boundaries is specified, write out only the edges that have non-zero marker

For starters, this writes every edge as a separate feature, but at some point it may make polygons out of the edges.

write_sms(fname)

write_suntans(pathname)

create cells.dat, edges.dat and points.dat from the TriGrid instance, all in the directory specified by pathname

write_tulip(fname)

Write a basic representation of the grid to a tulip compatible file

exception stompy.grid.trigrid.TriGridError

Bases: exceptions.Exception

stompy.grid.trigrid.circumcenter(p1, p2, p3)

stompy.grid.trigrid.dist(a, b)

stompy.grid.trigrid.ensure_ccw(points)

stompy.grid.trigrid.ensure_cw(points)

stompy.grid.trigrid.is_ccw(points)

stompy.grid.trigrid.outermost_rings(poly_list)

given a list of Polygons, return indices for those that are not inside any other polygon

stompy.grid.trigrid.rot(angle, pnts)

stompy.grid.trigrid.signed_area(points)

stompy.grid.ugrid module

stompy.grid.unstructured_grid module

Module contents