Ordering lake shoreline data

[scjames]

When generating a flexible-mesh model, most grid-generation software (e.g., RGFGRID) require that lake shoreline data be ordered so that the boundary is swept continuously in the clockwise direction (even if this yields segments with convex curves). A quick Google search was not productive; it was difficult to find a tool/utility for what would seem a common request. MATPLOT has a "boundary" function that might work (but a MATLAB license is required) and there is Alpha Shapes that might also work (the source code was not readily available). I would like to see you develop a utility code (especially if it were in Python) that facilitates conversion of a raw, arbitrarily sorted, two-column ASCII (or *.shp) file into ordered, clockwise, shoreline data. An important feature to include would be the ability to handle islands within the lake.

I have attached an example dataset that includes convex boundary segments and islands, so this would be a challenging example problem and if you can get this working, it should work for most other datasets, too.

WatereeShoreline.zip

[ssolson]

@scjames thanks for your interest in MHKiT and submitting an issue. As discussed offline I and others at Integral have run into this issue when creating land boundary files for setting up simulations and created a python script that has done some of this that we could adapt to solve this common problem. I believe MHKit is a good fit for this type of function and we have some work related to this coming in Q1. I will assign this task to myself to get implemented in Q1 as soon as I knock out some of the previous open issues!

[ssolson]

Below is an arbitrary alpha shape function and plotting script. As can be seen, shoreline detail is currently being cut off. Increasing the detail leads to many smaller individual polygons. The concavity varies significantly throughout the shoreline and makes a single alpha parameter non-ideal. A couple of ideas I would want to try:

1. partitioning the data into smaller sections and considering individual alphas for each
2. considering a mean or expected difference between points and trying to only increase this distance as needed to get to the next point. In this case, most of the data seems consistently spaced. The southeast area is an example problem area of non evenly spaced data. Also, there are areas where bridges go over the lake and could cut off the shoreline data following this approach

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from descartes import PolygonPatch
from shapely.geometry import Point, MultiLineString, MultiPoint
from scipy.spatial import Delaunay
from shapely.ops import  polygonize, cascaded_union

def add_edge(edges, i, j):
    """
    Add a line between the i-th and j-th points,
    if not in the list already
    """        
    if (i, j) not in edges or (j, i) not in edges:    
        edges.add( (i, j) )
    return edges

def alpha_shape(points, alpha):
    """
    Compute the alpha shape (concave hull) of a set
    of points.

    Parameters
    ----------
    points: array
        Iterable container of points.
    alpha: float
        alpha value to influence the length between 2 points at the border. 

    Returns
    -------
    cascaded_union(triangles): object

    """
    assert len(points) > 3, "Need at least four points"

    tri = Delaunay(points)
    edges = set()

    # loop over triangles:
    # ia, ib, ic = indices of corner points of the triangle
    i=0
    for ia, ib, ic in tri.simplices:
      pa = points[ia]
      pb = points[ib]
      pc = points[ic]
      # Lengths of sides of triangle
      a = np.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
      b = np.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
      c = np.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
      # Semiperimeter of triangle
      s = (a + b + c)/2.0
      # Area of triangle by Heron's formula
      area = np.sqrt(s*(s-a)*(s-b)*(s-c))
      circum_r = a*b*c/(4.0*area)
      # Radius filter.      
      if circum_r < 1.0/alpha:
          edges = add_edge(edges, ia, ib)
          edges = add_edge(edges, ib, ic)
          edges = add_edge(edges, ic, ia)

    edge_points = []
    pts = pd.Series(points)
    for idx in edges:
        pt0=pts[idx[0]]
        pt1=pts[idx[1]]
        edge_points.append([pt0, pt1])

    m = MultiLineString(edge_points)
    triangles = cascaded_union(list(polygonize(m)))

    return triangles

def plot_polygon(polygon):
    fig = plt.figure(figsize=(10,10))
    ax = fig.add_subplot(111)
    margin = .3
    x_min, y_min, x_max, y_max = polygon.bounds
    ax.set_xlim([x_min-margin, x_max+margin])
    ax.set_ylim([y_min-margin, y_max+margin])
    patch = PolygonPatch(polygon, fc='#999999',
                         ec='#000000', fill=True,
                         zorder=-1)
    ax.add_patch(patch)
    return fig

df = pd.read_csv('WatereeShoreline.txt', sep='\t')

points_2d = [(x,y) for x,y in zip(df.X, df.Y)]
concave_hull = alpha_shape(points_2d, 0.0005)

plot_polygon(concave_hull)
plt.plot(df.X, df.Y, '.')
plt.show()

[ssolson]

Closing due to inactivity.