Basic profile calculations

Profiles/_modidx.py

@@ -5,4 +5,5 @@
                 'doc_host': 'https://civildot.github.io',
                 'git_url': 'https://github.com/civildot/Profiles',
                 'lib_path': 'Profiles'},
-  'syms': {'Profiles.core': {'Profiles.core.foo': ('core.html#foo', 'Profiles/core.py')}}}
+  'syms': { 'Profiles.profiles': { 'Profiles.profiles.profile': ('profiles.html#profile', 'Profiles/profiles.py'),
+                                   'Profiles.profiles.submodel_gis': ('profiles.html#submodel_gis', 'Profiles/profiles.py')}}}

Profiles/profiles.py

@@ -0,0 +1,196 @@
+# AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/00_Profiles.ipynb.
+
+# %% auto 0
+__all__ = ['submodel_gis', 'profile']
+
+# %% ../nbs/00_Profiles.ipynb 3
+from operator import itemgetter
+
+import matplotlib.tri as mtri
+import numpy as np
+import fiona
+from shapely.geometry import shape
+from shapely.geometry import Point, LineString
+
+# %% ../nbs/00_Profiles.ipynb 4
+def submodel_gis(pline: LineString, gis_file: str, 
+                 layer:str ='Triangulation') -> tuple:
+
+    '''
+    Takes a line geometry and TIN terrain in civildot format
+    and returns reduced TIN to cover just the line.
+
+        Parameters:
+            pline (LineString) : A shapely linestring geometry
+            gis_file (str)     : TIN file location (civildot format)
+            layer (str)        : Layer containing the TIN
+
+        Returns:
+            (tuple) : (matplotlib triangular grid format, x (np-array), 
+                       y (np-array), z (np-array))
+    '''
+
+    with fiona.open(gis_file, layer=layer) as source:
+        hits = list(source.items(bbox=pline.bounds))
+
+    triangles = list()
+    for tri in hits:
+        geom = shape(tri[1].geometry)
+        if geom.intersects(pline):
+            triangles.append((tri[1], geom))
+
+    del hits
+
+    vertices = dict()
+    for tri in triangles:
+        props = dict(tri[0].properties)
+        coords = list(tri[1].exterior.coords)
+        vertices[props['pointA']] = [coords[0]]
+        vertices[props['pointB']] = [coords[1]]
+        vertices[props['pointC']] = [coords[2]]
+
+    refs = list(vertices.keys())
+    refs.sort()
+
+    x = list()
+    y = list()
+    z = list()
+    for i, ref in enumerate(refs):
+        vertices[ref].append(i)
+        x.append(vertices[ref][0][0])
+        y.append(vertices[ref][0][1])
+        z.append(vertices[ref][0][2])
+    x = np.array(x)
+    y = np.array(y)
+    z = np.array(z)
+
+    elements = list()
+    for tri in triangles:
+        props = dict(tri[0].properties)
+        pt1_ref = props['pointA']
+        pt2_ref = props['pointB']
+        pt3_ref = props['pointC']
+        pt1 = vertices[pt1_ref][1]
+        pt2 = vertices[pt2_ref][1]
+        pt3 = vertices[pt3_ref][1]
+        elements.append([pt1, pt2, pt3])
+    elements = np.array(elements)
+
+    triang = mtri.Triangulation(x, y, elements)
+
+    return (triang, x, y, z)
+
+
+# %% ../nbs/00_Profiles.ipynb 8
+def profile(line: LineString, gis_file: str, layer: str='Triangulation', 
+            stake_start: float=0.0, stake_end: float=None,
+            ptol: float=0.001, dtol: float=0.001) -> list:
+
+    '''
+    Takes a line geometry and TIN terrain in civildot format
+    and returns the profile data of the line.
+
+        Parameters:
+            line (LineString)    : A shapely linestring geometry
+            gis_file (str)       : TIN file location (civildot format)
+            layer (str)          : Layer containing the TIN
+            stake_start (float)  : Custom start chainage, default 0.0
+            stake_end (float)    : Custom end chainage, 
+                                   default stake_start + line length
+            ptol (float)         : Nodes closer to the line than ptol will 
+                                   be included in the profile.
+            dtol (float)         : Points closer than dtol on the profile
+                                   line is ommited
+
+        Returns:
+            [list] : each record, 
+                     [no, stake, elevation, x, y, 3d distance]
+    '''
+
+    if stake_end is None:
+        stake_end = stake_start + line.length
+        llength = line.length
+    else:
+        llength = stake_end - stake_start
+
+    coords = list(line.coords)
+    spt = coords[0]
+    ept = coords[-1]
+
+    triang, x, y, z = submodel_gis(line, gis_file)
+
+    # Calculate the profile
+    prof = list()
+    surface = mtri.LinearTriInterpolator(triang, z, trifinder=None)
+    # Start
+    elevation = round(float(surface.__call__(spt[0], spt[1])), 3)
+    prof.append([0, stake_start, elevation, round(spt[0], 3), 
+                 round(spt[1], 3), stake_start])
+    # End
+    elevation = round(float(surface.__call__(ept[0], ept[1])), 3)
+    prof.append([None, round(stake_end, 3), elevation, 
+                 round(ept[0], 3), round(ept[1], 3), None])
+
+    # Edges
+    for edge in triang.edges:
+        pt1 = Point(x[edge[0]], y[edge[0]], z[edge[0]])
+        pt2 = Point(x[edge[1]], y[edge[1]], z[edge[1]])
+        cedge = LineString((pt1, pt2))
+        if cedge.crosses(line):
+            xpt = cedge.intersection(line)
+            fraction = line.project(xpt, normalized=True)
+            stake = stake_start+(fraction*llength)
+            stake = round(stake, 3)
+            elevation = round(xpt.z, 3)
+            prof.append([None, stake, elevation, round(xpt.x, 3), 
+                 round(xpt.y, 3), None])
+
+    # Nodes
+    xy = zip(x, y)
+    for nd in xy:
+        pt = Point(nd[0], nd[1])
+        fraction = line.project(pt, normalized=True)
+        stake = stake_start+(fraction*llength)
+        if abs(stake-stake_start) < dtol or abs(stake-stake_end) < dtol:
+            continue
+        pt2 = line.interpolate(fraction, normalized=True)
+        if pt.distance(pt2) < ptol:
+            elevation = round(float(surface.__call__(pt2.x,pt2.y)), 3)
+            stake = round(stake, 3)
+            prof.append([None, stake, elevation, round(pt2.x, 3), 
+                         round(pt2.y, 3), None])
+
+    prof.sort(key=itemgetter(1))
+
+    prof2 = list()
+    count = 0
+    for i, rec in enumerate(prof):
+        if i == 0:
+            rec[0] = count
+            count += 1
+            prof2.append(rec)
+            continue
+        if rec[1]-prof[i-1][1] < dtol:
+            continue
+        rec[0] = count
+        count += 1
+        prof2.append(rec)
+
+    del prof
+
+    cumulative = 0.0
+    for i, rec in enumerate(prof2):
+        if i == 0:
+            continue
+        x1, y1, z1 = prof2[i-1][3], prof2[i-1][4], prof2[i-1][2]
+        x2, y2, z2 = rec[3], rec[4], rec[2]
+        Δx, Δy, Δz = x2-x1, y2-y1, z2-z1
+        dist = (Δx**2+Δy**2)**0.5
+        dist3d = (dist**2+Δz**2)**0.5
+        cumulative += dist3d
+        rec[5] = round(cumulative, 3)
+
+    for item in prof2:
+        print(item)
+
+    return prof2