7. Library

7.1. Custom Classes

Classes used for FEM

class IslandsLib.Classes.river(rname, x, y, rtype='b', downstream=True)

container for rivers and by extension the boundaries of the domain, I can see the comments of some people…, :))

class IslandsLib.Classes.topo(dname, x, y, z)

container for DEMs.

7.2. Finite Elements Modeling

IslandsLib.IslandLens(islands='Desirade', fname='../data/Contours/Atlantic/Guadeloupe/Desirade.txt', ttype='pq33a1000', fi=0.0001, sub_sampling=20, clockwise=True, lakes=[], plot=True)

Poisson Resolution for an Island (closed contour with null potential)

Parameters:

  • islands (str, optional) – Island name. Defaults to ‘Desirade’.

  • fname (str, optional) – Coastline contour data file. Defaults to “../../data/Desirade.txt”.

  • ttype (str, optional) – Triangle triangulation parameters. Defaults to “pq33a1000”.

  • fi (interpolator or float, optional) – value of the Poisson coefficient (dimensionless). Defaults to 1e-4

  • sub_sampling (int, optional) – sampling of the contour (1 = all points are preserved, 10 = one point every 10 points is conserved)

  • lakes (boolean, None) – to be used as a dic when lake contours are to be included

  • plot (boolean, True) – if True plots the figures

Returns:

  • array (\(u = \phi^2\) solution of poisson equation)

  • Th object (Th domain of the solution)

  • 2D array (X positions on a grid)

  • 2D array (Y positions on a grid)

  • 2D array (Zm = masked water table height on a (X,Y) grid)

  • object (itp regular grid interpolator for Z)

Finite Element Functions

  • Mesh generation using Triangle and pyFreeFem

  • Functions to solve Laplace or poisson equations

  • Graphic functions

IslandsLib.FEM.add_inner_vbn(vertices, segments, rivers)

Add inner boundaries (rivers) to the list of vertices and segments and create river nodes

Parameters:

  • vertices (list) – outer contour vertices

  • segments (list) – outer contour segments

  • rivers (list of rivers) – inner rivers

Returns:

  • lists of real numbers (vertices = list of contour vertices and inner rivers)

  • list of integers (segments = list of contour segments and inner rivers defined by their endpoint nodes)

  • list of integers (rnodes = contour nodes and inner rivers)

IslandsLib.FEM.add_lakes_vbn(vertices, segments, lakes, len_xpt)

add lakes as internal boundaries and create lake nodes

Parameters:

  • vertices (list) – list of vertices

  • segments (list) – list of segments

  • lakes (list of rivers) – island’s lakes

  • len_xpt (int) – original length of outer border (island contour)

Returns:

  • list – vertices

  • list – segments

  • list – lake nodes

  • list – lake center coordinates for holes calculation in triangle

IslandsLib.FEM.add_mask(ax, px, py, outer=False, extent=[])

Mask to remove the very ugly and especially false extrapolated part

Parameters:

  • ax (axis of the graphical representation)

  • px.py (outline for the mask)

  • outer (if true, use extent to draw an outer mask between the polygons)

  • transitioning (formed by extent and outline. This is mainly used to remove the extrapolations when)

  • interpolation. (from a triangular mesh to a square grid by)

  • extent (extension of the grid.)

Return type:

p (polygon) of the mask

Note

Adapted from https://www.matecdev.com/posts/shapely-point-polygon.html

IslandsLib.FEM.boundary_values(label_number, wanted_value, Th)

Defines a set of value at a specific boundary

  • First step is to pick out the corresponding nodes in the mesh

  • Second step is to assign the values

Parameters:

  • label_number (Label of the boundary to modify (integer defined by FreeFem))

  • wanted_values (List of boundary conditions)

  • Th (Th object containing the FreeFem mesh)

Returns:

array

Return type:

list of boundary values corresponding to each nodes (boundary grammian)

Note

Source: Olivier Devauchelle & Céleste Romon

IslandsLib.FEM.build_border_dictionnary(mesh, nodes, bname_list)

obsolete: Creates a full border dictionary from the triangle-generated mesh

IslandsLib.FEM.create_FreeFem_mesh(xv, yv, mesh, borders, adapt=False, plot=False)

Creates the FreeFem mesh from the mesh obtained with triangle

Parameters:

  • xv (arrays) – coordinates of the triangle nodes

  • yv (arrays) – coordinates of the triangle nodes

  • mesh (triangle mesh)

  • borders (dict) – external and internal boundaries

  • adapt (bool) – ​​if true, then FreeFem creates an adaptmesh

  • plot (bool) – if true, plots the mesh

Returns:

  • pyFreeFem object (Th object for FreeFem)

  • List of Th (list of boundaries retrieved from Th)

Note

in Th x and y are stored according to segment number whereas in triangle they are stored according to node number Boundary segments are stored in a dictionnary whose keys are the border labels and vals are the segments.

IslandsLib.FEM.create_grid_from_u(Th, u, dupuit=True, fval=0, dxy=10)

Create gridded values for the water table height Z from the reconstructed values u at each node of the mesh

Parameters:

  • Th (Th object) – mesh

  • u (array) – nodal values of the water table

  • dupuit (bool, optional) – if yes \(u = Z^2\), by default True

  • fval (int, optional) – fill value for gridding, by default 0

  • dxy (int, optional) – grid step in meters, by default 10 OR defined by 1000x1000 matrix

Returns:

  • arrays – X,Y,Z gridded positions and values of the water table

  • arrays – dx, dy step

  • object – itp regular grid interpolator for Z

Note

If dupuit : \(Z = \sqrt{u}\)

IslandsLib.FEM.create_outer_vbn(borders, for_FF=False)

Create vertices, segments, and nodes from a list of borders for triangulation

Parameters:

  • borders (list) – list of river objects that form the outer boundary of the area

  • for_FF (boolean) – if true, add an end node for FreeFem

Returns:

  • list of floats (vertices = list of contour vertices)

  • list of ints (segments = list of contour segments defined by their end nodes)

  • list of ints (nodes = contour nodes)

Note

  • I Remove the last segment that doesn’t correspond to anything and add a closing segment.

  • A flag should be added in case the contour is naturally closed.

  • be careful to close the borders or else triangle won’t succeed in segments creation and everything crashes.

IslandsLib.FEM.create_triangle_mesh(borders, vertices, segments, holes=None, plot=False, ttype='p')

Create the mesh using the triangle library

Parameters:

  • borders (lits of rivers) – the external and internal borders

  • vertices (the boundary vertices)

  • segments (the boundary segments)

  • plot (bool) – if true, represents the resulting mesh

  • ttype (str) – triangulation type, by default p is the simplest triangulation without adaptation and creation of new nodes. This is typically the one I use before a FreeFem adaptmesh.

Returns:

  • array of real numbers (xv, yv, vertex coordinates)

  • mesh (dictionary) (mesh, mesh produced by triangle)

  • list of rivers (borders, recalculated boundaries with reordered nodes)

Note

triangle allows constrained triangulation with internal borders. Personally, I find that the triangle meshes after adaptation are more balanced than those of FreeFem. O.D. says I just don’t know how to use adaptmesh. He’s probably right.

IslandsLib.FEM.grid_plot_u(fig, ax, Th, X, Y, Z)

Graphical representation of \(Z = \sqrt{u}\) on an interpolated grid with streamlines and contour mask

Parameters:

  • fig (figure) – figure on which to plot

  • ax (axis) – axis on which to plot

  • Th (Th object) – mesh and boundaries

  • X (array) – X position matrix

  • Y (array) – Y position matrix

  • Z (array) – Z (water table) position matrix

Returns:

  • figure (fig)

  • axis (ax)

IslandsLib.FEM.is_inside(node, check_list)

checks if a node is inside a list

Parameters:

  • node (int) – the number of a node of the triangulation

  • check_list (list of int) – a list of node numbers

Returns:

bool

Return type:

True/False

IslandsLib.FEM.is_on(xv, yv, seg1, seg2, pval=False, epsilon=1e-05)

Tests if a segment seg2 belongs to another segment seg1

Used to retrieve additional edge nodes created by triangles when asked for an adapted triangulation (equivalent to a FreeFem adaptmesh) under constraints

Parameters:

  • xv (arrays) – real coordinate vectors

  • yv (arrays) – real coordinate vectors

  • seg1 (list) – containing segment

  • seg2 (list) – segment to test

  • pval (bool) – if true, prints the slope coefficient values ​​of the two segments

  • epsilon (real) – margin of error tolerated in comparing the two segments. epsilon must be non-zero, due to numerical approximations, but very small to avoid confusion (double membership, for example). The default value seems to be suitable…

Returns:

bool

Return type:

true or false

IslandsLib.FEM.list_errors()

Small function to get the description of python errors.

IslandsLib.FEM.mask_data(X, Y, Z, p)

Masks the Z array using X,Y and p. this is useful because the interpolation of u, Th on a regular grids always leads to some smear outside the island contour

Parameters:

  • X (2D array) – X positions

  • Y (2D array) – Y positions

  • Z (2D array) – Unmasked water table

  • p (Shapely polygon) – island contour used as a mask

Returns:

masked water table

Return type:

2D array

IslandsLib.FEM.order_nodes(xv, yv, mesh, borders, plot=False)

Classifies the contour nodes.

Parameters:

  • xv (arrays) – vectors of the triangle mesh node coordinates

  • yv (arrays) – vectors of the triangle mesh node coordinates

  • mesh (dictionary of vertices, segments, and nodes) – the triangle mesh

  • rivers (borders; list of) – the boundaries

  • plot (bool) – if true, represents the graph of the oriented nodes

Returns:

list or rivers

Return type:

borders whose nodes have been reclassified

Note

This operation is very important because it is a necessary condition for successfully converting a triangle mesh to a FreeFem mesh.

IslandsLib.FEM.plot_psi(Th, psi, ax)

Tricontour plot de u la solution de poisson obtenue par solve_fem

Parameters:

  • Th (objet Th)

  • u (potential solution to be plotted)

Note

The difference with plot_u is that the solution was obtained by solving for the symmetric boundary conditions hence for the flow lines

IslandsLib.FEM.plot_u(Th, u, ax)

Tricontour plot of u the solution of poisson equaiton obtained with solve_fem

Parameters:

  • Th (Th object)

  • u (solution for the potential)

IslandsLib.FEM.prepare_boundaries(borders, lakes=None, rivers=None, for_FF=False)

Create the border contour from borders

Parameters:

  • borders (list of rivers) – external boundaries

  • rivers (list of rivers) – internal rivers

  • lakes (list of rivers) – island’s lakes

  • for_FF (bool) – if true, arranges segments and end nodes to match FreeFem mesh requirements

Returns:

  • list of rivers (complete external and internal borders)

  • list of pairs of real numbers (vertices of border nodes)

  • list of pairs of integer numbers (border segments)

  • list of pairs of intergers (holes coordinates)

Note

this supposes that the borders are oriented in a correct continuous fashion

IslandsLib.FEM.solve_fem(Th, bcs=None, fi=0.0001)

Solve Poisson equation

Parameters:

  • Th (Th object,)

  • bcs (dic,) – boundary conditions dictionary,

  • fi (array, optional,) – right-hand side of the Poisson equation \((\Delta\rho/\rho) (P/K)\),

Returns:

  • array (u, the solution for the potential (here, h^2))

  • object (Th, the Th object that contains all the mesh information)

Note

For now, \(fi = (\Delta\rho/\rho) (P/K)\) is a constant, but nothing prevents non-uniform precipitation distributions

IslandsLib.FEM.update_node_list(borders, rivers, nodes, mesh)

loop to add the boundary nodes created during mesh generation in the node lists

Parameters:

  • borders (list) – list of border classes (outer borders)

  • rivers (list) – list of rivers (inner borders)

  • nodes (list) – list of nodes

  • mesh (dictionnary) – mesh dictionnary (be clearer…)

Returns:

  • list (list of nodes)

  • list (list of new_nodes)

  • list of list (list of border names)

Note

pb d’affectation des noeuds. Il faut empêcher la double affectation sinon on est dans le caca !

7.3. IGN and Topographic Data treatment

IGN Data Processing Functions

IslandsLib.IGN.IGN_to_segment(lg={}, list_key=[], ex=[], inverted=False, exclusion=[])

Create order in the IGN mess (unordered segments) to obtain continuous contours

Warning

OBSOLETE CRASHES REGULARLY

Note

  1. Create a list of successive segments by noting their direction relative to the first segment in the list.

    • Decide on a direction (end=True means we expect the segment closest to the current segment to be oriented like this one)

    • Take the endpoints of a segment

    • Look for the segment closest to the end (if end=True) or the beginning (if end=False)

    • Test whether the end or the beginning of the segment is closest and modify the end function

    • Do the same with the next segment

  2. Take each segment from the list, flip it if necessary, and extract its coordinates to create a single contour.

Parameters:

  • lg (dictionary of segments under shape code:(x,y))

  • list_key (list of segment codes)

  • ex (list of segment endpoints (x1,y1,x2,y2))

  • inverted (False by default, if True returns inverted coordinates)

  • exclusion ([] by default, list of exclusion codes)

Returns:

arrays

Return type:

x,y contour or segment coordinates

IslandsLib.IGN.IGN_to_segment_v2(lg, list_ori, ex_dic, plot_res=False)

Transforms a series of segments into a single contour.

Parameters:

  • (dict) (ex_dic)

  • (list) (list_ori)

  • (dict)

  • (bool (plot_res)

  • optional) (plot true Defaults to False.)

Returns:

arrays

Return type:

x,y

IslandsLib.IGN.clip_ign_shapefile(fname, px, py)

Keep only the information from the fname file located inside the px.py coordinate polygon.

Parameters:

  • fname (shapefile)

  • px.py (coordinates of the vertices of the clipping polygon.)

Returns:

GeoDataFrame

Return type:

clipped GeoDataFrame

Note

The coordinates must be in the same system, of course.

IslandsLib.IGN.create_river_from_linestrings(gpd, ilist)

Creates two x,y matrices from a given list of linestrings in the correct order.

Parameters:

  • gpd (geopandas dataframe containing the geometric data)

  • ilist (list of indices (iloc) of the topo objects to be retrieved from gpd)

Returns:

float arrays

Return type:

xr, yr coordinate vectors

Note

IGN shapefiles are not ordered; river or coastline segments must be ordered first to reconstruct a continuous contour.

IslandsLib.IGN.get_alti_lims(fname)

Retrieves information about an alti DB file

Parameters:

fname (file name)

IslandsLib.IGN.get_cross_section(dem, npoints=100)

Extracts cross section from dem

Parameters:

  • dem (class) – dem class defined in myFEMlib to store gridded DEM data

  • npoints (int, optional) – number of points of the cross section, by default 100

IslandsLib.IGN.get_elevation_from_dem(dem, borders, minx=0, miny=0)

Retrieves the elevation of each river in the border from the DEM.

Parameters:

  • borders (list of river-class objects)

  • dem (DEM class obtained from bdalti)

Returns:

list of classes

Return type:

borders with elevation

IslandsLib.IGN.get_nearest_point(dem, p)

Gets the point in dem that is nearest to p Using this function implies that the dem is in cartesion (metric) units not in spherical (lat, lon, z).

Parameters:

  • dem (class) – dem class as defined in myFEMlib

  • p (list) – point position

Returns:

list

Return type:

[x,y,z] coordinates of nearest point

IslandsLib.IGN.get_riverlist_inside(fname, px, py)

Returns the list of rivers located within a polygon. Also returns the length and type of their geometry.

Parameters:

  • fname (name of the IGN shapefile)

  • px.py (coordinates of the polygon used as a mask)

Returns:

  • list (list of river names.)

  • float (length of the path (note: this only makes sense if the coordinates are metric).)

  • str (geometry type (e.g., Linestring).)

IslandsLib.IGN.get_rivers_inside(fname, px, py, rnames, rlist)

Extracting rivers from a shape file, retaining only those within a given contour

Parameters:

  • fname (shapefile name)

  • ps (contour within which rivers are retained)

  • py (contour within which rivers are retained)

  • rnames (list of river names (TOPONYMS in the IGN sense))

  • rlist (ordered list (upstream -> downstream) of node lists allowing you to choose the paths to retain and create a single contour)

Returns:

list of classes

Return type:

list of river objects of the extracted rivers

IslandsLib.IGN.linestring_to_list(geom='')

Transforms a Linestring geometry into a list of x,y positions

Parameters:

  • (str (geom)

  • optional) (geometry to be transformed. Defaults to "".)

Returns:

list

Return type:

x,y coordinates of points making the original geometry

Note

use shapely if > 2 available

IslandsLib.IGN.output_to_file(fname, entete, x, y)

Outputs x,y lists to text file. Used for Islands coastlines

Parameters:

  • (str) (entete)

  • (str)

  • array) (y (list or)

  • array)

IslandsLib.IGN.plot_rivers_inside(fname, px, py, rlist=[], border=True, length=0, index=False, ax=None)

Graphical representation of rivers within a polygon

Parameters:

  • fname (shapefile name)

  • px.py (coordinates of the boundary polygon)

  • rlist (list of place names. If not empty, only rivers whose place names are in the list are represented.)

  • border (if true, the boundary polygon is represented.)

  • length (if the value is > 0, filters rivers whose length is less than the value passed as an argument.)

  • index (if true, then the start of the river is marked. Note: this only works for simple geometries (not multi-something).)

  • ax (if ax, then the figure is provided, otherwise one is created.)

IslandsLib.IGN.polygon_to_list(geom='')

Transforms a Polygon type geometry into a list of x,y positions

Parameters:

  • (str (geom)

  • optional) (geometry to be transformed. Defaults to "".)

Returns:

list

Return type:

x,y coordinates of points making the original geometry

Note

use shapely if > 2 available

IslandsLib.IGN.read_alti(fname)

Opens a 25m alti database file, retrieves the altitude and position data, and returns square X, Y, and Z matrices.

Parameters:

fname (Full name of the alti file)

Returns:

(1000, 1000) nunmpy arrays

Return type:

X, Y, Z position and altitude matrices