Project Name | Stars | Downloads | Repos Using This | Packages Using This | Most Recent Commit | Total Releases | Latest Release | Open Issues | License | Language |
---|---|---|---|---|---|---|---|---|---|---|
Martini | 559 | a year ago | 8 | isc | JavaScript | |||||
A JavaScript library for real-time RTIN terrain mesh generation | ||||||||||
Manifoldplus | 225 | 2 years ago | 6 | C++ | ||||||
ManifoldPlus: A Robust and Scalable Watertight Manifold Surface Generation Method for Triangle Soups | ||||||||||
Pixel2meshplusplus | 198 | 3 years ago | 5 | bsd-3-clause | Python | |||||
Pixel2Mesh++: Multi-View 3D Mesh Generation via Deformation. In ICCV2019. | ||||||||||
Bmeshunity | 192 | a year ago | 5 | mit | C# | |||||
A Unity package to make runtime procedural mesh generation more flexible. | ||||||||||
Occupancy_flow | 174 | a year ago | 3 | mit | Python | |||||
This repository contains the code for the ICCV 2019 paper "Occupancy Flow - 4D Reconstruction by Learning Particle Dynamics" | ||||||||||
Seismicmesh | 101 | 5 months ago | 32 | July 18, 2021 | 13 | gpl-3.0 | Python | |||
2D/3D serial and parallel triangular mesh generation tool for finite element methods. | ||||||||||
Voxel2mesh | 95 | 8 months ago | gpl-3.0 | Python | ||||||
Voxel2Mesh: 3D Mesh Model Generation from Volumetric Data | ||||||||||
Unityhexgrid | 91 | 4 years ago | mit | C# | ||||||
A hex grid generation tool, for use in unity | ||||||||||
Voxelman | 70 | a year ago | 1 | mit | C++ | |||||
A voxel engine for godot. I'ts a c++ engine module. | ||||||||||
Mesh2d | 65 | 3 years ago | 2 | other | MATLAB | |||||
MESH2D is a MATLAB-based Delaunay mesh generator for two-dimensional geometries. |
🌊 🌀
Coastal ocean mesh generation from vector and raster GIS data.
The Python version of the algorithm does not yet have a citation; however, similar algorithms and ideas are shared between both version.
[1] - Roberts, K. J., Pringle, W. J., and Westerink, J. J., 2019.
OceanMesh2D 1.0: MATLAB-based software for two-dimensional unstructured mesh generation in coastal ocean modeling,
Geoscientific Model Development, 12, 1847-1868. https://doi.org/10.5194/gmd-12-1847-2019.
Besides posting issues with the code on Github, you can also ask questions via our Slack channel here.
Otherwise please reach out to Dr. Keith Roberts ([email protected]) with questions or concerns!
Please include version information when posting bug reports.
oceanmesh
uses versioneer.
The version can be inspected through Python
python -c "import oceanmesh; print(oceanmesh.__version__)"
or through
python setup.py version
in the working directory.
To see what's going on with oceanmesh
when running scripts, you can turn on logging (which is by default suppressed) by importing the two modules sys
and logging
and then placing one of the three following logging commands after your imports in your calling script. The amount of information returned is the greatest with logging.DEBUG
and leas with logging.WARNING
.
import logging
import sys
logging.basicConfig(stream=sys.stdout, level=logging.WARNING)
# logging.basicConfig(stream=sys.stdout, level=logging.INFO)
# logging.basicConfig(stream=sys.stdout, level=logging.DEBUG)
The notes below refer to installation on platforms other than MS Windows. For Windows, refer to the following section.
For installation, oceanmesh
needs cmake, CGAL:
sudo apt install cmake libcgal-dev
CGAL and can also be installed with conda
:
conda install -c conda-forge cgal
After that, clone the repo and oceanmesh
can be updated/installed using pip.
pip install -U -e .
On some clusters/HPC in order to install CGAL, you may need to load/install gmp and mpfr. For example, to install:
sudo apt install libmpfr-dev libgmp3-dev
Python under Windows can easily experience DLL hell due to version incompatibilities. Such is the case for the combination of support packages required for OceanMesh. For this reason, we have provided 'install_cgal.bat' to build a CGAL development distribution separately as a prerequisite.
Prerequisites to build CGAL using the provided batch file are as follows:
After successful intallation of a CGAL development package, proceed via one of the two options below to generate a python environment with OceanMesh installed.
If you are using a conda-based Python distribution, then 'install_oceanmesh.bat' should take care of everything, provided no package conflicts arise.
If you have a different Python distribution, or if you do not want to use packages from conda forge, then the following process may be required:
⚠️
WARNING: THIS PROGRAM IS IN ACTIVE DEVELOPMENT. INSTALLATION IS ONLY RECOMMENDED FOR DEVELOPERS AT THIS TIME. WHEN A STABLE API IS REACHED, THE PROGRAM WILL BE AVAILABLE VIA pypi
oceanmesh
can mesh a polygonal region in the vast majority of coordinate reference systems (CRS). Thus, all oceanmesh
scripts begin with declaring the extent and CRS to be used and/or transofmring to a different CRS like this.
import oceanmesh as om
EPSG = 32619 # A Python int, dict, or str containing the CRS information (in this case UTM19N)
bbox = (
-70.29637,
-43.56508,
-69.65537,
43.88338,
) # the extent of the domain (can also be a multi-polygon delimited by rows of np.nan)
extent = om.Region(
extent=bbox, crs=4326
) # set the region (the bbox is given here in EPSG:4326 or WGS84)
extent = extent.transform_to(EPSG) # Now I transform to the desired EPSG (UTM19N)
print(
extent.bbox
) # now the extents are in the desired CRS and can be passed to various functions later on
oceanmesh
uses shoreline vector datasets (i.e., ESRI shapefiles) and digital elevation models (DEMs) to construct mesh size functions and signed distance functions to adapt mesh resolution for complex and irregularly shaped coastal ocean domains.
Shoreline datasets are necessary to build signed distance functions which define the meshing domain. Here we show how to download a world shoreline dataset referred to as GSHHG and read it into to oceanmesh
.
import zipfile
import requests
import oceanmesh as om
# Download and load the GSHHS shoreline
url = "http://www.soest.hawaii.edu/pwessel/gshhg/gshhg-shp-2.3.7.zip"
filename = url.split("/")[-1]
with open(filename, "wb") as f:
r = requests.get(url)
f.write(r.content)
with zipfile.ZipFile("gshhg-shp-2.3.7.zip", "r") as zip_ref:
zip_ref.extractall("gshhg-shp-2.3.7")
fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
EPSG = 4326 # EPSG code for WGS84 which is what you want to mesh in
# Specify and extent to read in and a minimum mesh size in the unit of the projection
extent = om.Region(extent=(-75.000, -70.001, 40.0001, 41.9000), crs=EPSG)
min_edge_length = 0.01 # In the units of the projection!
shoreline = om.Shoreline(
fname, extent.bbox, min_edge_length, crs=EPSG
) # NB: the Shoreline class assumes WGS84:4326 if not specified
shoreline.plot(
xlabel="longitude (WGS84 degrees)",
ylabel="latitude (WGS84 degrees)",
title="shoreline boundaries",
)
# Using our shoreline, we create a signed distance function
# which will be used for meshing later on.
sdf = om.signed_distance_function(shoreline)
DEMs are used to build some mesh size functions (e.g., wavelength, enforcing size bounds, enforce maximum Courant bounds) but are not essential for mesh generation purposes. The DEM used below 'Eastcoast.nc' was created using the Python package elevation with the following command:
eio clip - o EastCoast.nc --bounds -74.4 40.2 -73.4 41.2
This data is a clip from the SRTM30 m elevation dataset.
import oceanmesh as om
fdem = "datasets/EastCoast.nc"
# Digital Elevation Models (DEM) can be read into oceanmesh in
# either the NetCDF format or GeoTiff format provided they are
# in geographic coordinates (WGS84)
# If no extents are passed (i.e., the kwarg bbox), then the entire extent of the
# DEM is read into memory.
# Note: the DEM will be projected to the desired CRS automatically.
EPSG = 4326
dem = om.DEM(fdem, crs=EPSG)
dem.plot(
xlabel="longitude (WGS84 degrees)",
ylabel="latitude (WGS84 degrees)",
title="SRTM 30m",
cbarlabel="elevation (meters)",
vmin=-10, # minimum elevation value in plot
vmax=10, # maximum elevation value in plot
)
The domain is defined by a signed distance function. A signed distance function can be automatically generated from a complex coastal ocean domain as such
import oceanmesh as om
fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
EPSG = 4326 # EPSG:4326 or WGS84
extent = om.Region(extent=(-75.00, -70.001, 40.0001, 41.9000), crs=EPSG)
min_edge_length = 0.01 # minimum mesh size in domain in projection
shoreline = om.Shoreline(fname, extent.bbox, min_edge_length)
# build a signed distance functiona automatically
sdf = om.signed_distance_function(shoreline)
In some situations, it is necessary to flip the definition of what's inside
the meshing domain. This can be accomplished via the invert
kwarg for the signed_distance_function
.
sdf = om.signed_distance_function(shoreline, invert=True)
Setting invert=True
will be mesh the 'land side' of the domain rather than the ocean.
In oceanmesh
mesh resolution can be controlled according to a variety of feature-driven geometric and topo-bathymetric functions. In this section, we briefly explain the major functions and present examples and code. Reasonable values for some of these mesh sizing functions and their affect on the numerical simulation of barotropic tides was investigated in Roberts et. al, 2019
All mesh size functions are defined on regular Cartesian grids. The properties of these grids are abstracted by the Grid class.
A high degree of mesh refinement is often necessary near the shoreline boundary to capture its geometric complexity. If mesh resolution is poorly distributed, critical conveyances may be missed, leading to larger-scale errors in the nearshore circulation. Thus, a mesh size function that is equal to a user-defined minimum mesh size h0 along the shoreline boundary, growing as a linear function of the signed distance d from it, may be appropriate.
import oceanmesh as om
fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
EPSG = 4326 # EPSG:4326 or WGS84
extent = om.Region(extent=(-75.00, -70.001, 40.0001, 41.9000), crs=EPSG)
min_edge_length = 0.01 # minimum mesh size in domain in projection
shoreline = om.Shoreline(fname, extent.bbox, min_edge_length)
edge_length = om.distance_sizing_function(shoreline, rate=0.15)
fig, ax, pc = edge_length.plot(
xlabel="longitude (WGS84 degrees)",
ylabel="latitude (WGS84 degrees)",
title="Distance sizing function",
cbarlabel="mesh size (degrees)",
holding=True,
)
shoreline.plot(ax=ax)
One major drawback of the distance mesh size function is that the minimum mesh size will be placed evenly along straight stretches of shoreline. If the distance mesh size function generates too many vertices (or your application can tolerate it), a feature mesh size function that places resolution according to the geometric width of the shoreline should be employed instead (Conroy et al., 2012;Koko, 2015).
In this function, the feature size (e.g., the width of channels and/or tributaries and the radius of curvature of the shoreline) along the coast is estimated by computing distances to the medial axis of the shoreline geometry. In oceanmesh
, we have implemented an approximate medial axis method closely following Koko, (2015).
import oceanmesh as om
fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
EPSG = 4326 # EPSG:4326 or WGS84
extent = om.Region(extent=(-75.00, -70.001, 40.0001, 41.9000), crs=EPSG)
min_edge_length = 0.01 # minimum mesh size in domain in projection
shoreline = om.Shoreline(fname, extent.bbox, min_edge_length)
sdf = om.signed_distance_function(shoreline)
# Visualize the medial points
edge_length = om.feature_sizing_function(
shoreline, sdf, max_edge_length=0.05, plot=True
)
fig, ax, pc = edge_length.plot(
xlabel="longitude (WGS84 degrees)",
ylabel="latitude (WGS84 degrees)",
title="Feature sizing function",
cbarlabel="mesh size (degrees)",
holding=True,
xlim=[-74.3, -73.8],
ylim=[40.3, 40.8],
)
shoreline.plot(ax=ax)
Some mesh size functions will not produce smooth element size transitions when meshed and this can lead to problems with numerical simulation. All mesh size function can thus be graded such that neighboring mesh sizes are bounded above by a constant. Mesh grading edits coarser regions and preserves the finer mesh resolution zones.
Repeating the above but applying a gradation rate of 15% produces the following:
import oceanmesh as om
fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
EPSG = 4326 # EPSG:4326 or WGS84
extent = om.Region(extent=(-75.00, -70.001, 40.0001, 41.9000), crs=EPSG)
min_edge_length = 0.01 # minimum mesh size in domain in projection
shoreline = om.Shoreline(fname, extent.bbox, min_edge_length)
sdf = om.signed_distance_function(shoreline)
edge_length = om.feature_sizing_function(shoreline, sdf, max_edge_length=0.05)
edge_length = om.enforce_mesh_gradation(edge_length, gradation=0.15)
fig, ax, pc = edge_length.plot(
xlabel="longitude (WGS84 degrees)",
ylabel="latitude (WGS84 degrees)",
title="Feature sizing function with gradation bound",
cbarlabel="mesh size (degrees)",
holding=True,
xlim=[-74.3, -73.8],
ylim=[40.3, 40.8],
)
shoreline.plot(ax=ax)
In shallow water theory, the wave celerity, and hence the wavelength λ, is proportional to the square root of the depth of the water column. This relationship indicates that more mesh resolution at shallower depths is required to resolve waves that are shorter than those in deep water. With this considered, a mesh size function hwl that ensures a certain number of elements are present per wavelength (usually of the M2-dominant semi-diurnal tidal species but the frequency of the dominant wave can be specified via kwargs) can be deduced.
In this snippet, as before we compute the feature size function, but now we also compute the wavelength-to-gridscale sizing function using the SRTM dataset and compute the minimum of all the functions before grading. We discretize the wavelength of the M2 by 30 elements (e.g., wl=30)
import oceanmesh as om
fdem = "datasets/EastCoast.nc"
fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
min_edge_length = 0.01
extent = om.Region(extent=(-74.3, -73.8, 40.3,40.8), crs=4326)
dem = om.DEM(fdem, bbox=extent,crs=4326)
shoreline = om.Shoreline(fname, dem.bbox, min_edge_length)
sdf = om.signed_distance_function(shoreline)
edge_length1 = om.feature_sizing_function(shoreline, sdf, max_edge_length=0.05)
edge_length2 = om.wavelength_sizing_function(
dem, wl=100, period=12.42 * 3600
) # use the M2-tide period (in seconds)
# Compute the minimum of the sizing functions
edge_length = om.compute_minimum([edge_length1, edge_length2])
edge_length = om.enforce_mesh_gradation(edge_length, gradation=0.15)
fig, ax, pc = edge_length.plot(
xlabel="longitude (WGS84 degrees)",
ylabel="latitude (WGS84 degrees)",
title="Feature sizing function + wavelength + gradation bound",
cbarlabel="mesh size (degrees)",
holding=True,
xlim=[-74.3, -73.8],
ylim=[40.3, 40.8],
)
shoreline.plot(ax=ax)
The distance, feature size, and/or wavelength mesh size functions can lead to coarse mesh resolution in deeper waters that under-resolve and smooth over the sharp topographic gradients that characterize the continental shelf break. These slope features can be important for coastal ocean models in order to capture dissipative effects driven by the internal tides, transmissional reflection at the shelf break that controls the astronomical tides, and trapped shelf waves. The bathymetry field contains often excessive details that are not relevant for most flows, thus the bathymetry can be smoothed by a variety of filters (e.g., lowpass, bandpass, and highpass filters) before calculating the mesh sizes.
import oceanmesh as om
fdem = "datasets/EastCoast.nc"
fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
EPSG = 4326 # EPSG:4326 or WGS84
bbox = (-74.4, -73.4, 40.2, 41.2)
extent = om.Region(extent=bbox, crs=EPSG)
dem = om.DEM(fdem, crs=EPSG)
min_edge_length = 0.0025 # minimum mesh size in domain in projection
max_edge_length = 0.10 # maximum mesh size in domain in projection
shoreline = om.Shoreline(fname, extent.bbox, min_edge_length)
sdf = om.signed_distance_function(shoreline)
edge_length1 = om.feature_sizing_function(
shoreline,
sdf,
max_edge_length=max_edge_length,
crs=EPSG,
)
edge_length2 = om.bathymetric_gradient_sizing_function(
dem,
slope_parameter=5.0,
filter_quotient=50,
min_edge_length=min_edge_length,
max_edge_length=max_edge_length,
crs=EPSG,
) # will be reactivated
edge_length3 = om.compute_minimum([edge_length1, edge_length2])
edge_length3 = om.enforce_mesh_gradation(edge_length3, gradation=0.15)
After mesh generation has terminated, a secondary round of mesh improvement strategies is applied that is focused towards improving the geometrically worst-quality triangles that often occur near the boundary of the mesh and can make simulation impossible. Low-quality triangles can occur near the mesh boundary because the geospatial datasets used may contain features that have horizontal length scales smaller than the minimum mesh resolution. To handle this issue, a set of algorithms is applied that iteratively addresses the vertex connectivity problems. The application of the following mesh improvement strategies results in a simplified mesh boundary that conforms to the user-requested minimum element size.
Topological defects in the mesh can be removed by ensuring that it is valid, defined as having the following properties:
the vertices of each triangle are arranged in counterclockwise order;
conformity (a triangle is not allowed to have a vertex of another triangle in its interior); and
traversability (the number of boundary segments is equal to the number of boundary vertices, which guarantees a unique path along the mesh boundary).
Here are some of the relevant codes to address these common problems.
# Address (1) above.
points, cells = fix_mesh(points, cells)
# Addresses (2)-(3) above. Remove degenerate mesh faces and other common problems in the mesh
points, cells = make_mesh_boundaries_traversable(points, cells)
# Remove elements (i.e., "faces") connected to only one channel
# These typically occur in channels at or near the grid scale.
points, cells = delete_faces_connected_to_one_face(points, cells)
# Remove low quality boundary elements less than min_qual
points, cells = delete_boundary_faces(points, cells, min_qual=0.15)
# Apply a Laplacian smoother that preserves the element density
points, cells = laplacian2(points, cells)
Mesh generation is based on the DistMesh algorithm and requires only a signed distance function and a mesh sizing function. These two functions can be defined through the previously elaborated commands above; however, they can also be straightforward functions that take an array of point coordinates and return the signed distance/desired mesh size.
In this example, we demonstrate all of the above to build a mesh around New York, United States with an approximate minimum element size of around 1 km expanding linear with distance from the shoreline to an approximate maximum element size of 5 km.
Here we use the GSHHS shoreline here and the Python package meshio
to write the mesh to a VTK file for visualization in ParaView. Other mesh formats are possible; see meshio
for more details
import meshio
import oceanmesh as om
fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
EPSG = 4326 # EPSG:4326 otherwise known as WGS84
extent = om.Region(extent=(-75.00, -70.001, 40.0001, 41.9000), crs=EPSG)
min_edge_length = 0.01 # minimum mesh size in domain in projection
shore = om.Shoreline(fname, extent.bbox, min_edge_length)
edge_length = om.distance_sizing_function(shore, max_edge_length=0.05)
domain = om.signed_distance_function(shore)
points, cells = om.generate_mesh(domain, edge_length)
# remove degenerate mesh faces and other common problems in the mesh
points, cells = om.make_mesh_boundaries_traversable(points, cells)
points, cells = om.delete_faces_connected_to_one_face(points, cells)
# remove low quality boundary elements less than 15%
points, cells = om.delete_boundary_faces(points, cells, min_qual=0.15)
# apply a Laplacian smoother
points, cells = om.laplacian2(points, cells)
# write the mesh with meshio
meshio.write_points_cells(
"new_york.vtk",
points,
[("triangle", cells)],
file_format="vtk",
)
The major downside of the DistMesh algorithm is that it cannot handle regional domains with fine mesh refinement or variable datasets due to the intense memory requirements. The multiscale mesh generation technique addresses these problems and enables an arbitrary number of refinement zones to be incorporated seamlessly into the domain.
Areas of finer refinement can be incorporated seamlessly by using the generate_multiscale_mesh
function. In this case, the user passes lists of signed distance and edge length functions to the mesh generator but besides this the user API remains the same to the previous mesh generation example. The mesh sizing transitions between nests are handled automatically to produce meshes suitable for FEM and FVM numerical simulations through the parameters prefixed with "blend".
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
import oceanmesh as om
fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
EPSG = 4326 # EPSG:4326 or WGS84
extent1 = om.Region(extent=(-75.00, -70.001, 40.0001, 41.9000), crs=EPSG)
min_edge_length1 = 0.01 # minimum mesh size in domain in projection
bbox2 = np.array(
[
[-73.9481, 40.6028],
[-74.0186, 40.5688],
[-73.9366, 40.5362],
[-73.7269, 40.5626],
[-73.7231, 40.6459],
[-73.8242, 40.6758],
[-73.9481, 40.6028],
],
dtype=float,
)
extent2 = om.Region(extent=bbox2, crs=EPSG)
min_edge_length2 = 4.6e-4 # minimum mesh size in domain in projection
s1 = om.Shoreline(fname, extent1.bbox, min_edge_length1)
sdf1 = om.signed_distance_function(s1)
el1 = om.distance_sizing_function(s1, max_edge_length=0.05)
s2 = om.Shoreline(fname, extent2.bbox, min_edge_length2)
sdf2 = om.signed_distance_function(s2)
el2 = om.distance_sizing_function(s2)
# Control the element size transition
# from coarse to fine with the kwargs prefixed with `blend`
points, cells = om.generate_multiscale_mesh(
[sdf1, sdf2],
[el1, el2],
)
# remove degenerate mesh faces and other common problems in the mesh
points, cells = om.make_mesh_boundaries_traversable(points, cells)
# remove singly connected elements (elements connected to only one other element)
points, cells = om.delete_faces_connected_to_one_face(points, cells)
# remove poor boundary elements with quality < 15%
points, cells = om.delete_boundary_faces(points, cells, min_qual=0.15)
# apply a Laplacian smoother that preservers the mesh size distribution
points, cells = om.laplacian2(points, cells)
# plot it showing the different levels of resolution
triang = tri.Triangulation(points[:, 0], points[:, 1], cells)
gs = gridspec.GridSpec(2, 2)
gs.update(wspace=0.5)
plt.figure()
bbox3 = np.array(
[
[-73.78, 40.60],
[-73.75, 40.60],
[-73.75, 40.64],
[-73.78, 40.64],
[-73.78, 40.60],
],
dtype=float,
)
ax = plt.subplot(gs[0, 0]) #
ax.set_aspect("equal")
ax.triplot(triang, "-", lw=1)
ax.plot(bbox2[:, 0], bbox2[:, 1], "r--")
ax.plot(bbox3[:, 0], bbox3[:, 1], "m--")
ax = plt.subplot(gs[0, 1])
ax.set_aspect("equal")
ax.triplot(triang, "-", lw=1)
ax.plot(bbox2[:, 0], bbox2[:, 1], "r--")
ax.set_xlim(np.amin(bbox2[:, 0]), np.amax(bbox2[:, 0]))
ax.set_ylim(np.amin(bbox2[:, 1]), np.amax(bbox2[:, 1]))
ax.plot(bbox3[:, 0], bbox3[:, 1], "m--")
ax = plt.subplot(gs[1, :])
ax.set_aspect("equal")
ax.triplot(triang, "-", lw=1)
ax.set_xlim(-73.78, -73.75)
ax.set_ylim(40.60, 40.64)
plt.show()
Using oceanmesh is now possible for global meshes. The process is done in two steps:
import os
import numpy as np
import oceanmesh as om
from oceanmesh.region import to_lat_lon
import matplotlib.pyplot as plt
# utilities functions for plotting
def crosses_dateline(lon1, lon2):
return abs(lon1 - lon2) > 180
def filter_triangles(points, cells):
filtered_cells = []
for cell in cells:
p1, p2, p3 = points[cell[0]], points[cell[1]], points[cell[2]]
if not (crosses_dateline(p1[0], p2[0]) or crosses_dateline(p2[0], p3[0]) or crosses_dateline(p3[0], p1[0])):
filtered_cells.append(cell)
return filtered_cells
# Note: global_stereo.shp has been generated using global_tag() function in pyposeidon
# https://github.com/ec-jrc/pyPoseidon/blob/9cfd3bbf5598c810004def83b1f43dc5149addd0/pyposeidon/boundary.py#L452
fname = "tests/global/global_latlon.shp"
fname2 = "tests/global/global_stereo.shp"
EPSG = 4326 # EPSG:4326 or WGS84
bbox = (-180.00, 180.00, -89.00, 90.00)
extent = om.Region(extent=bbox, crs=4326)
min_edge_length = 0.5 # minimum mesh size in domain in meters
max_edge_length = 2 # maximum mesh size in domain in meters
shoreline = om.Shoreline(fname, extent.bbox, min_edge_length)
sdf = om.signed_distance_function(shoreline)
edge_length = om.distance_sizing_function(shoreline, rate=0.11)
# once the size functions have been defined, wed need to mesh inside domain in
# stereographic projections. This is way we use another coastline which is
# already translated in a sterographic projection
shoreline_stereo = om.Shoreline(fname2, extent.bbox, min_edge_length, stereo=True)
domain = om.signed_distance_function(shoreline_stereo)
points, cells = om.generate_mesh(domain, edge_length, stereo=True, max_iter=100)
# remove degenerate mesh faces and other common problems in the mesh
points, cells = om.make_mesh_boundaries_traversable(points, cells)
points, cells = om.delete_faces_connected_to_one_face(points, cells)
# apply a Laplacian smoother
points, cells = om.laplacian2(points, cells, max_iter=100)
lon, lat = to_lat_lon(points[:, 0], points[:, 1])
trin = filter_triangles(np.array([lon,lat]).T, cells)
fig, ax, pc = edge_length.plot(holding = True, plot_colorbar=True, cbarlabel = "Resolution in °", cmap='magma')
ax.triplot(lon, lat, trin, color="w", linewidth=0.25)
plt.tight_layout()
plt.show()
See the tests inside the testing/
folder for more inspiration. Work is ongoing on this package.
To run the oceanmesh
unit tests (and turn off plots), check out this repository and type tox
. tox
can be installed via pip.
This software is published under the GPLv3 license