adding in edge stuff and SDF (#9)

* adding in edge stuff and SDF

Author: keith roberts

oceanmesh/__init__.py

@@ -1,6 +1,9 @@
 from oceanmesh.geodata import DEM, Geodata, Shoreline
-from oceanmesh.edgefx import Grid
+from oceanmesh.grid import Grid
 from oceanmesh.edgefx import distance_sizing_function
+from oceanmesh.signed_distance_function import signed_distance_function
+from oceanmesh.edges import get_poly_edges, draw_edges
+from .inpoly import inpoly
 
 __all__ = [
     "SizeFunction",
@@ -9,4 +12,8 @@
     "DEM",
     "Shoreline",
     "distance_sizing_function",
+    "signed_distance_function",
+    "get_poly_edges",
+    "draw_edges",
+    "inpoly",
 ]

oceanmesh/edgefx.py

@@ -1,7 +1,7 @@
 import numpy
 import skfmm
 
-from .Grid import Grid
+from .grid import Grid
 
 __all__ = ["distance_sizing_function"]
 

oceanmesh/edges.py

@@ -0,0 +1,61 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib import collections as mc
+
+
+nan = np.nan
+
+
+def get_poly_edges(poly):
+    """Given a winded polygon represented as a set of ascending line segments
+    with separated features indicated by nans, this function calculates
+    the edges of the polygon such that each edge indexes the start and end
+    coordinates of each line segment of the polygon.
+
+    Parameters
+    ----------
+    poly: array-like, float
+        A 2D array of point coordinates with features sepearated by NaNs
+
+    Returns
+    -------
+    edges: array-like, int
+        A 2D array of integers containing indexes into the `poly` array.
+
+    """
+    ix = np.argwhere(np.isnan(poly[:, 0]))
+    ix = np.insert(ix, 0, 0)
+
+    edges = []
+    for s in range(len(ix) - 1):
+        col1 = np.arange(ix[s], ix[s + 1] - 1)
+        col2 = np.arange(ix[s] + 1, ix[s + 1])
+        tmp = np.vstack((col1, col2)).T
+        tmp = np.append(tmp, [[ix[s + 1] - 1, ix[s] + 1]], axis=0)
+        edges.append(tmp)
+    return np.concatenate(edges, axis=0)
+
+
+def draw_edges(poly, edges):
+    """Visualizes the polygon as a bunch of line segments
+
+    Parameters
+    ----------
+    poly: array-like, float
+        A 2D array of point coordinates with features sepearated by NaNs.
+    edges: array-like, int
+        A 2D array of integers indexing into the `poly` array.
+
+    Returns
+    -------
+    None
+
+    """
+    lines = []
+    for edge in edges:
+        lines.append([poly[edge[0]], poly[edge[1]]])
+    lc = mc.LineCollection(lines, linewidths=2)
+    fig, ax = plt.subplots()
+    ax.add_collection(lc)
+    ax.autoscale()
+    plt.show()

oceanmesh/Grid.py oceanmesh/grid.pyoceanmesh/Grid.py renamed to oceanmesh/grid.py

@@ -7,8 +7,8 @@
 
 class Grid:
     """Abstracts a structured grid along with
-    operations (e.g., min, project, etc.), stores data
-    `values` defined at each grid point.
+    primitive operations (e.g., min, project, etc.) and
+    stores data `values` defined at each grid point.
 
     Parameters
     ----------
@@ -200,7 +200,7 @@ def project(self, grid2):
         new_values[new_values == FILL] = grid2.values[new_values == FILL]
         return Grid(bbox=grid2.bbox, grid_spacing=grid2.grid_spacing, values=new_values)
 
-    def plot(self, hold=False):
+    def plot(self, hold=False, vmin=0.0, vmax=0.1):
         """Visualize the values in :obj:`Grid`
 
         Parameters
@@ -222,7 +222,7 @@ def plot(self, hold=False):
         x, y = self.create_grid()
 
         fig, ax = plt.subplots()
-        ax.pcolor(x, y, self.values, vmin=0.0, vmax=0.1, shading="auto")
+        ax.pcolor(x, y, self.values, vmin=vmin, vmax=vmax, shading="auto")
         ax.axis("equal")
         if hold is False:
             plt.show()

oceanmesh/inpoly.py

@@ -0,0 +1,205 @@
+import numpy as np
+from numba import jit
+
+
+def inpoly(vert, node, edge=None, ftol=4.9485e-14):
+    """A port of Darren Engwirda's `inpoly` routine into Python.
+
+       Returns the "inside/outside" status for a set of vertices VERT and
+       a polygon {NODE,EDGE}  embedded in a two-dimensional plane. General
+       non-convex and multiply-connected polygonal regions can be handled.
+
+
+    Parameters
+    ----------
+    vert: array-like
+        VERT is an N-by-2 array of XY coordinates to query.
+
+    node: array-like
+        NODE is an M-by-2 array of polygon vertices.
+
+    edge: array-like, optional
+        EDGE is a P-by-2 array of edge indexing. Each
+        row in EDGE represents an edge of the polygon, such that
+        node(edge(KK,1),:) and node(edge(KK,2),:) are the coord-
+        inates of the endpoints of the KK-TH edge. If the argum-
+        ent EDGE is omitted it assumed that the vertices in NODE
+        are connected in ascending order.
+
+    ftol: float, otpional
+        FTOL is a floating-point tolerance for boundary comparisons.
+        By default, FTOL = EPS ^ 0.85.
+
+
+    Returns:
+    --------
+    STAT: array-like
+        STAT is an associated N-by-1 logical array,
+        with STAT(II) = TRUE if VERT(II,:) is an interior point.
+
+    BNDS: array-like
+        N-by-1 logical array BNDS, with BNDS(II) = TRUE if VERT(II,:)
+        lies "on" a boundary segment,
+
+
+    Notes
+    -----
+
+    ALL CREDIT GOES TO THE ORIGINAL AUTHOR: Darren Engwirda
+
+    Darren Engwirda : 2017 --
+    Email           : de2363@columbia.edu
+    Last updated    : 31/03/2020
+
+
+    Author of port:
+
+    Keith J. Roberts (2020)
+    Email: krober@usp.br
+    Date: 2020-13-09
+
+    """
+
+    nnod = len(node)
+    nvrt = len(vert)
+
+    # create edge if not supplied
+    if edge is None:
+        edge = np.vstack((np.arange(0, nnod - 1), np.arange(1, nnod))).T
+        edge = np.concatenate((edge, [[nnod - 1, 0]]))
+
+    STAT = np.zeros((nvrt))
+    BNDS = np.zeros((nvrt))
+
+    # prune points using bbox
+    mask = (
+        (vert[:, 0] >= np.nanmin(node[:, 0]))
+        & (vert[:, 0] <= np.nanmax(node[:, 0]))
+        & (vert[:, 1] >= np.nanmin(node[:, 1]))
+        & (vert[:, 1] <= np.nanmax(node[:, 1]))
+    )
+
+    vert = vert[mask]
+
+    # flip to ensure y-axis is the `long` axis
+    vmin = np.amin(vert, axis=0)
+    vmax = np.amax(vert, axis=0)
+    ddxy = vmax - vmin
+
+    lbar = np.sum(ddxy) / 2.0
+
+    if ddxy[0] > ddxy[1]:
+        vert = vert[:, [1, 0]]
+        node = node[:, [1, 0]]
+
+    # sort points via y-value
+    swap = node[edge[:, 1], 1] < node[edge[:, 0], 1]
+
+    tmp = edge[swap]
+    edge[swap, :] = tmp[:, [1, 0]]
+
+    ivec = np.argsort(vert[:, 1])
+    vert = vert[ivec]
+
+    t_stat, t_bnds = _inpoly(vert, node, edge, ftol, lbar)
+
+    stat = np.ones(len(vert))
+    bnds = np.ones(len(vert))
+    stat[ivec] = t_stat
+    bnds[ivec] = t_bnds
+
+    STAT[mask] = stat
+    BNDS[mask] = bnds
+
+    return STAT, BNDS
+
+
+@jit(nopython=True)
+def _inpoly(vert, node, edge, ftol, lbar):
+
+    feps = ftol * lbar ** 2
+    veps = ftol * lbar ** 1
+
+    _nvrt = len(vert)
+    _nedge = len(edge)
+
+    _stat = np.zeros((_nvrt))
+    _bnds = np.zeros((_nvrt))
+
+    # loop over polygon edges
+    for epos in range(_nedge):
+
+        inod = edge[epos, 0]
+        jnod = edge[epos, 1]
+
+        # calc. edge bounding-box
+        yone = node[inod, 1]
+        ytwo = node[jnod, 1]
+        xone = node[inod, 0]
+        xtwo = node[jnod, 0]
+
+        xmin = np.minimum(xone, xtwo)
+        xmax = np.maximum(xone, xtwo)
+
+        xmax += veps
+
+        ymin = yone - veps
+        ymax = ytwo + veps
+
+        ydel = ytwo - yone
+        xdel = xtwo - xone
+
+        # find top VERT[:,1]<YONE
+        ilow = 0  # for python
+        iupp = _nvrt
+
+        while ilow < iupp - 1:  # binary search
+            imid = int(ilow + np.floor((iupp - ilow) / 2))
+            if vert[imid, 1] < ymin:
+                ilow = imid
+            else:
+                iupp = imid
+
+        if vert[ilow, 1] >= ymin:
+            ilow -= 1
+
+        # calc. edge-intersection
+        for jpos in range(ilow + 1, _nvrt):
+
+            if _bnds[jpos]:
+                continue
+
+            xpos = vert[jpos, 0]
+            ypos = vert[jpos, 1]
+
+            if ypos <= ymax:
+                if xpos >= xmin:
+                    if xpos <= xmax:
+                        # compute crossing number
+                        mul1 = ydel * (xpos - xone)
+                        mul2 = xdel * (ypos - yone)
+                        if feps >= np.abs(mul2 - mul1):
+                            # BNDS -- approx. on edge
+                            _bnds[jpos] = 1
+                            _stat[jpos] = 1
+                        elif (ypos == yone) & (xpos == xone):
+                            # BNDS -- match about ONE
+                            _bnds[jpos] = 1
+                            _stat[jpos] = 1
+                        elif (ypos == ytwo) & (xpos == xtwo):
+                            # BNDS -- match about TWO
+                            _bnds[jpos] = 1
+                            _stat[jpos] = 1
+                        elif mul1 < mul2:
+                            if (ypos >= yone) & (ypos < ytwo):
+                                # advance crossing number
+                                _stat[jpos] = 1 - _stat[jpos]
+
+                else:
+                    if (ypos >= yone) & (ypos < ytwo):
+                        # advance crossing number
+                        _stat[jpos] = 1 - _stat[jpos]
+            else:
+                # done -- due to the sort
+                break
+    return _stat, _bnds

oceanmesh/signed_distance_function.py

@@ -0,0 +1,44 @@
+import numpy
+import skfmm
+
+from .grid import Grid
+from .inpoly import inpoly
+from . import edges
+
+__all__ = ["signed_distance_function"]
+
+
+def signed_distance_function(shoreline):
+    """Takes a `shoreline` object containing segments representing islands and mainland boundaries
+    and calculates a signed distance function with it (assuming the polygons are all closed)
+
+    Parameters
+    ----------
+
+    Returns
+    -------
+
+    """
+
+    print("Building a signed distance function...")
+    grid = Grid(bbox=shoreline.bbox, grid_spacing=shoreline.h0)
+    phi = numpy.ones(shape=(grid.nx, grid.ny))
+    lon, lat = grid.create_grid()
+    poly = numpy.vstack((shoreline.inner, shoreline.mainland))
+    indices = grid.find_indices(poly, lon, lat)
+    phi[indices] = -1.0
+    # call Fast Marching Method
+    print("Calculating the distance...")
+    dis = skfmm.distance(phi, grid.grid_spacing)
+    # now sign it
+    e = edges.get_poly_edges(poly)
+    print("Signing the distance...")
+    lon = numpy.reshape(lon, -1)
+    lat = numpy.reshape(lat, -1)
+    qry = numpy.vstack((lon, lat)).T
+
+    inside, _ = inpoly(qry, poly, e)
+
+    grid.values = dis * inside.reshape((dis.shape))
+    grid.build_interpolant()
+    return grid

requirements.txt

@@ -5,3 +5,4 @@ netcdf4
 pyshp
 Pillow==7.1.0
 scikit-fmm==2019.1.30
+numba

setup.py

@@ -28,5 +28,6 @@
         "matplotlib",
         "Pillow",
         "scikit-fmm",
+        "numba",
     ],
 )