Example - San Diego Meshing

This example illustrates how run a more complicated meshing example for a real-world case. In this case, you use a polygon and other features defining the extents of San Diego Bay to generate a mesh that can be used for running a hydraulic model of the bay. An overview of the model is shown in the image below.

A python-based class has been created for this example that contains two PolyInput objects that can be used to generate a 2D mesh.

In [1]:
import cartopy.crs as ccrs
import geoviews as gv
import pandas as pd
from shapely.geometry import Polygon

import meshing_tools
import xmsmesh as mesh

gv.extension("bokeh")

Plot initial geometry

In [2]:
san_diego_example = meshing_tools.SanDiegoExample()
poly1 = san_diego_example.poly1
poly2 = san_diego_example.poly2

basemap = gv.tile_sources.StamenTerrainRetina().options(width=900, height=900)

polygons = gv.Polygons(
    data=[
        meshing_tools.polys_as_shapely(poly2), 
        meshing_tools.polys_as_shapely(poly1)
    ], 
    crs=ccrs.UTM(zone=11),
).options(fill_color=None,)

points = gv.Points(
    data=[poly1.outside_poly, poly2.outside_poly, poly2.inside_polys[-1]],
    crs=ccrs.UTM(zone=11),
).options(color='black', size=4)

path = gv.Path(
    data=[pd.DataFrame(data=poly2.inside_polys[-1], columns=['x', 'y', 'z'])], 
    crs=ccrs.UTM(zone=11), 
).options(line_color='black')

basemap * polygons * points * path
Out[2]:

The MultiPolyMesherIo class holds all the options for generating meshes from polygon data

In [3]:
mesh_io = mesh.meshing.MultiPolyMesherIo((poly1, poly2))

Generate the mesh

In [4]:
# Generate Mesh
succeded, errors = mesh.meshing.mesh_utils.generate_mesh(mesh_io=mesh_io)
if succeded:
    print("Meshing was successful")
else:
    print("Meshing errors found:")
    for err in errors:
        print("\t{}".format(err))

Meshing was successful

Plot Mesh

The mesh is converted to a pandas dataframe for plotting.

In [5]:
verts,cells = meshing_tools.xmsmesh_to_dataframe(mesh_io.points, mesh_io.cells)
vert_points = gv.project(gv.Points(verts, vdims=['z'], crs=ccrs.UTM(zone=11)))
trimesh = gv.TriMesh((cells, vert_points))

basemap * polygons * points * path * trimesh.edgepaths.options(line_width=1)
Out[5]: