Anisotropic river bathymetry interpolation

RivBathyInterp.py

# This script interpolates river bathymetry in an anisotropic way.
# Following the stream centerline, the interpolation is weighted longitudinally.
# This is accomplished by projecting the sinuous river to a straight line (s, n) coordinates,
# where s is distance downstream and n is distance left/right perpindicular to centerline.
# In (s, n) coordinates, the river is then compressed longitudinally, and then interpolated normally.
# The compression gives interpolation weight to upstream/downstream points instead of going up the bank.
# For best results the centerline should be smoothed, otherwise the inside bend has overlapping points.
# Input Excel file has three sheets, with simple X, Y, Z columns. X, Y only for centerline.

import pandas as pd
import numpy as np
import transform as tfm #requires shapely
from scipy.interpolate import griddata
from scipy.spatial import distance

longitudinalFactor = 10.0
cellSize = 2

# load bathymetry data
xyz = pd.read_excel('DivPoolXYZ.xlsx',sheetname='bathy')
bank = pd.read_excel('DivPoolXYZ.xlsx',sheetname='banks')
cl = pd.read_excel('DivPoolXYZ.xlsx',sheetname='cline')
pts = xyz.append(bank) #combine survey and bank points
del xyz, bank

# transformation to river coordinates
riv = tfm.SN_CoordinateSystem(cl.X.values,cl.Y.values)
es, en = riv.transform_xy_to_sn(pts.X.values,pts.Y.values)

# set up (s, n) grid system
s = np.arange(int(es.min()), int(es.max())+1, cellSize*5)/longitudinalFactor
n = np.arange(int(en.min()), int(en.max())+1, cellSize*5)
S, N = np.meshgrid(s, n)

# where the magic happens
elevs = griddata((es/longitudinalFactor,en),pts.Z,(S,N),method='linear')
xi, yi = riv.transform_sn_to_xy (S.flatten()*longitudinalFactor, N.flatten())

# filter interpolated points near input data points
delPts = np.min(distance.cdist(np.column_stack((xi,yi)),
                               np.column_stack((pts.X.values,pts.Y.values))),
                axis=1) < cellSize*5
xi, yi, elevs = (xi[~delPts], yi[~delPts], elevs.flatten()[~delPts])

# add in original input data points
np.append(xi,pts.X.values)
np.append(yi,pts.Y.values)
np.append(elevs,pts.Z.values)

# set up (x, y) grid system
x = np.arange(int(pts.X.min()), int(pts.X.max())+1, cellSize)
y = np.arange(int(pts.Y.max())+1, int(pts.Y.min()), -cellSize) #reversed
X,Y = np.meshgrid (x, y)

bathy = griddata((xi, yi), elevs, (X, Y), method='linear',rescale=True)
bathy.tofile('divPool.npy')

# make it easy to import grid to ArcMap
print ('To bring array into in ArcMap as Raster:')
print ('import numpy as np')
print ("bathy = np.fromfile('divPool.npy').reshape(%i, %i)" % bathy.shape)
print ("raster = arcpy.NumPyArrayToRaster(bathy,arcpy.Point(%i,%i),%i,%i,np.nan)" % (
    X.min()- cellSize/2., Y.min() - cellSize/2., cellSize, cellSize))

transform.py

import numpy as np
from shapely.geometry import LineString, Point
#taken from https://github.com/dharhas/smear/blob/master/smear/transform.py

class SN_CoordinateSystem:
    """ 
    Define an SN Coordinate System based on a given centerline
    Convert xy dataset into sn Coordinate System based on a given centerline.
    """

    def __init__(self, cx, cy, slope_distance=0.01, interp=None, interp_params=[]):
        """
        Initialize SN coordinate system. Requires arrays containing cartesian x,y values 
        and arrays containing centerline cx,cy cartesian coordinates 
        """

        if interp is not None:
            cx, cy = self.smooth_centerline(cx, cy, interp, interp_params) 

        self.centerline = LineString(zip(cx.tolist(), cy.tolist()))
        self.slope_distance = slope_distance

    def norm(self, x):
        return np.sqrt(np.square(x).sum())

    def smooth_centerline(self, x, y, interp, interp_params):
        #todo add other interp types like bezier curves

        #spline:

        #spline parameters
        s=3.0 # smoothness parameter
        k=2 # spline order
        nest=-1 # estimate of number of knots needed (-1 = maximal)
        xnew = np.array([])
        ynew = np.array([])

        for start, end, spacing in interp_params:
            if spacing==-1:
                xtmp = x[start:end]
                ytmp = y[start:end]
            else:
                npts = int(LineString(zip(x[start:end].tolist(), y[start:end].tolist())).length / spacing) 
                # find the knot points
                tckp,u = splprep([x[start:end],y[start:end]],s=s,k=k,nest=-1)
                xtmp, ytmp = splev(np.linspace(0,1,npts),tckp)

            xnew = np.concatenate((xnew,xtmp))
            ynew = np.concatenate((ynew,ytmp))

        return (xnew, ynew)
                         
    def transform_xy_to_sn(self, x, y):
        s = np.zeros(x.size)
        n = np.zeros(x.size)
        v = np.zeros((x.size,2))
        vn = np.zeros((x.size,2))
        for ii in range(x.size):
            pt = Point((x[ii],y[ii]))
            s[ii] = self.centerline.project(pt)
            pt_s = self.centerline.interpolate(s[ii])
            pt_s1 = self.centerline.interpolate(s[ii] - self.slope_distance)
            pt_s2 = self.centerline.interpolate(s[ii] + self.slope_distance)
            vn[ii,0] = pt.x - pt_s.x
            vn[ii,1] = pt.y - pt_s.y
            v[ii,0] = pt_s2.x - pt_s1.x
            v[ii,1] = pt_s2.y - pt_s1.y
            n[ii] = pt_s.distance(pt)

        n = -np.sign(np.cross(v,vn)) * n
        return (s,n)

    def transform_sn_to_xy(self,s,n):
        x = np.zeros(s.size)
        y = np.zeros(s.size)
        unit_normal = np.zeros(2)
        xy = np.zeros(2)
        for ii in range(s.size):
            pt_s = self.centerline.interpolate(s[ii])
            xy[0] = pt_s.x
            xy[1] = pt_s.y
            pt_s1 = self.centerline.interpolate(s[ii] - self.slope_distance)
            pt_s2 = self.centerline.interpolate(s[ii] + self.slope_distance)
            unit_normal[0] = -(pt_s2.y - pt_s1.y)
            unit_normal[1] = (pt_s2.x - pt_s1.x)
            unit_normal = unit_normal/self.norm(unit_normal)
            x[ii], y[ii] = xy - unit_normal*n[ii]
            #theta = np.arctan((pt_s2.x - pt_s1.x)/(pt_s2.y - pt_s1.y))
            #x[ii] = pt_s.x - n[ii]*np.cos(theta)
            #y[ii] = pt_s.y + n[ii]*np.sin(theta)
            
        return (x,y)