xmsstamper

Introduction

xmsstamper is a library used to generate "feature stamps". Feature stamps refer to the insertion of man-made structures into a natural topography or bathymetry set. In common terms, this means adding an embankment (such as a levy) or dredging a channel. A stamped feature usually a linear object or centerline. More information on feature stamping can be found here https://xmswiki.com/wiki/SMS:Feature_Stamping.

The basic workflow is to first define the land surface or bathymetry as a TIN (triangular surface). Then define a feature stamp consisting of a center line, cross sections along the center line, and end caps. The feature stamp is used to generate a new TIN that is cut off by the land surface or bathemetry.

These objects are demonstrated below. This example shows how to create a proposed embankment (stamp) for a bridge reconstruction.

In [1]:
import param
import numpy as np
import geoviews as gv
import holoviews as hv
import pandas as pd
import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import panel as pn

from shapely.geometry import Polygon
from holoviews.streams import Params

import hvplot.pandas
import hvplot.xarray

import xmsmesh
import xmsstamper
import xmsinterp
import xmsgrid

import ugrid_tools

hv.extension('bokeh')

Create a land surface TIN

In [2]:
tin = xmsinterp.triangulate.Tin(pts=ugrid_tools.pts)
tin.triangulate()

x = [x[0] for x in tin.pts]
y = [y[1] for y in tin.pts]
z = [z[2] for z in tin.pts]

nodes = np.column_stack([x, y, z])

srs = ccrs.UTM(zone=18)
nodes = gv.Points(nodes, vdims=['z'], crs=srs)
In [3]:
# display the elevation points
basemap = gv.tile_sources.ESRI
%opts Polygons Points TriMesh [width=700, height=500 colorbar=True tools=['hover']]
nodes = nodes.opts(
    hv.opts.Points(cmap='viridis', color='z')
)

basemap * nodes
Out[3]:

Create a StamperIo class

This class is used to define the feature stamp and set various options. This class is then passed to the stamp method and an output TIN is created.

In this example we create a "fill" stamp because we are creating a feature that will be placed on top of the land surface. If we wanted to create a new or modified channel we would use the "cut" option.

In [4]:
# type of stamp
io_stamping_type = 'fill'

Next we will define the coordinates of the center line of the proposed abutment.

In [5]:
# center line coordinates
io_centerline = [
    (305004.09408, 4386232.78456, 101.1936),
    (305014.04857444204, 4386233.860927977, 101.1936),
    (305022.8669304807, 4386234.814446314, 101.1936),
    (305030.6788192409, 4386235.659136725, 101.1936),
    (305037.59911207034, 4386236.407420116, 101.1936),
    (305043.7295696494, 4386237.070299156, 101.1936),
    (305049.1603384373, 4386237.657521958, 101.1936),
    (305053.971276156, 4386238.177723205, 101.1936),
    (305058.2331261235, 4386238.638551848, 101.1936),
    (305062.0085574754, 4386239.0467848135, 101.1936),
    (305065.3530866969, 4386239.408425137, 101.1936),
    (305068.31589409104, 4386239.72878975, 101.1936),
    (305070.940546764, 4386240.012590553, 101.1936),
    (305073.26563940616, 4386240.263999957, 101.1936),
    (305075.32536176883, 4386240.486715487, 101.1936),
    (305077.15000137076, 4386240.684011613, 101.1936),
    (305078.7663890244, 4386240.858789724, 101.1936),
    (305080.19829318626, 4386241.013619894, 101.1936),
    (305081.4667695012, 4386241.150778675, 101.1936),
    (305082.590470571, 4386241.272283233, 101.1936),
    (305083.58592000004, 4386241.37992, 101.1936),
]

Now we will define the cross sections along the center line. This is example is simple in that we will have the same cross section along the entire center line.

In [6]:
left = right = [(0, 101.1936), (3.81, 101.1936), (4.1148, 100.889)]
index_left_shoulder = index_right_shoulder = 1
left_max = right_max = 10.668
cs = xmsstamper.stamper.StampCrossSection(
    left=left,
    right=right,
    index_left_shoulder=index_left_shoulder,
    index_right_shoulder=index_right_shoulder,
    left_max=left_max,
    right_max=right_max,
)

cs_list = [cs] * len(io_centerline)

Next we have options for defining the end caps of the abutments. There are 3 different types of end caps: wing wall, sloped abutment, and guidebank. We will create a sloped abutment for the "last_end_cap" (the end cap at the end of the center line). The first_end_cap is at the beginning of the center line.

In [7]:
angle = 0.0
sloped_abutment = xmsstamper.stamper.SlopedAbutment(
    max_x=7.62,
    slope=[(0, 101.1936), (1.0, 100.889)]
)
last_end_cap = xmsstamper.stamper.StamperEndCap(
    angle=angle,
    sloped_abutment=sloped_abutment
)

Finally we will create a StamperIo class that holds the options for the feature stamp and is used to create the output TIN. The first TIN that we generate will have no defined land surface or bathemetry. This will create a TIN that conforms to the cross sections and end cap definitions.

In [8]:
stamp_io = xmsstamper.stamper.StamperIo(
    stamping_type=io_stamping_type,
    centerline=io_centerline,
    cs=cs_list,
    last_end_cap=last_end_cap,
    #bathymetry=tin,
)
xmsstamper.stamper.stamper.stamp(stamp_io)

out_tin = stamp_io.get_out_tin()
In [9]:
# draw the output TIN
tris = [x for x in zip(*[iter(out_tin.tris)] * 3)]

x, y, z = zip(*out_tin.pts)
stamp_nodes = np.column_stack([x, y, z])
stamp_nodes = gv.Points(stamp_nodes, vdims='z', crs=srs)

trimesh = gv.TriMesh((tris, stamp_nodes), crs=srs)
basemap * stamp_nodes * trimesh.opts(
   hv.opts.TriMesh(filled=True, colorbar=True, cmap='viridis', edge_color='z', edge_alpha=0.4, node_size=3)
)
Out[9]:

Now we will specify a land surface that will modify the output TIN. The TIN will be cut off by the elevations specified by the land surface TIN.

In [10]:
stamp_io.bathymetry = tin
xmsstamper.stamper.stamper.stamp(stamp_io)
out_tin = stamp_io.get_out_tin()
In [11]:
# draw the output TIN
tris = [x for x in zip(*[iter(out_tin.tris)] * 3)]

x, y, z = zip(*out_tin.pts)
stamp_nodes = np.column_stack([x, y, z])
stamp_nodes = gv.Points(stamp_nodes, vdims='z', crs=srs)

trimesh = gv.TriMesh((tris, stamp_nodes), crs=srs)
nodes * basemap * stamp_nodes * trimesh.opts(
   hv.opts.TriMesh(filled=True, colorbar=True, cmap='viridis', edge_color='z', edge_alpha=0.4, node_size=3)
)
Out[11]: