multiple linearly projected inverse distance weighting

mlpIDW.py

import numpy as np
import pandas as pd

#~~~~USER INPUT PARAMETERS~~~~
maxNearPoints = 20
anisotropyFactor = 15
searchTransv = 200 #meters
searchLongitud = searchTransv*anisotropyFactor #meters
gridSize = 3 #meters
invDistPower = 1
angularPower = 0.5
arcpy.env.overwriteOutput = True
pyAngle = r"def pyAngle(x1, y1, x2, y2):\n  return math.degrees(math.atan2(x1-x2,y1-y2))"

def angleSimilarity(deg1,deg2):
        return np.maximum(np.cos(np.deg2rad(deg1-deg2)),0)

def calcSdist(nearTable, lnGeoms):
    arcpy.management.AddField(nearTable,'S_DIST', 'DOUBLE')
    arcpy.management.CalculateField(nearTable, "NEAR_ANGLE", 
        "pyAngle(!FROM_X!,!FROM_Y!,!NEAR_X!,!NEAR_Y!)", "PYTHON", pyAngle)
    updateCur = arcpy.da.UpdateCursor(nearTable,['NEAR_FID','NEAR_X','NEAR_Y','S_DIST'])
    for r in updateCur:
        alongLnPct = lnGeoms[r[0]].measureOnLine(arcpy.Point(r[1],r[2]), True)
        if alongLnPct > 0 and alongLnPct < 1:
            r[3] = alongLnPct * lnGeoms[r[0]].length
            updateCur.updateRow(r)
        else:
            updateCur.deleteRow()

def calcLPIDW(pointDF,gridPt):
    query = 'NEAR_FID == %i and S_DIST >= %f and S_DIST <= %f and NEAR_DIST >= %f and NEAR_DIST <= %f' % (
            gridPt[0], gridPt[2] - searchLongitud, gridPt[2] + searchLongitud,
            gridPt[3] - searchTransv, gridPt[3] + searchTransv)
    intPoints = pointDF.query(query)
    if len(intPoints) > 0:
        intPoints['AngleSim'] = angleSimilarity(intPoints.NEAR_ANGLE,gridPt[1])**angularPower
        intPoints = intPoints.query('AngleSim > 0')
        intPoints['SNdist'] = np.fmax(1, ((intPoints.S_DIST - gridPt[2])**2 + 
                            ((intPoints.NEAR_DIST - gridPt[3]) * anisotropyFactor)**2)**invDistPower)
        intPoints['weights'] = intPoints.AngleSim / intPoints.SNdist
        intPoints = intPoints.nlargest(maxNearPoints,'weights').to_records()
        if intPoints.weights.sum() > 0:
            gridPt[4] = (intPoints.weights * intPoints.elevation).sum()/intPoints.weights.sum()
    return gridPt[4]

#~~~Main function for interpolating a surface, can be limited in extent using optional 'xHull' polygon
def mlpidw(inLines, inPoints, outSurface, xHull='calc'):
    ptTabl = arcpy.analysis.GenerateNearTable(inPoints, inLines,'ptSnap',location='LOCATION',angle='ANGLE',
                                    closest='ALL',closest_count=4, method='PLANAR')
    sr = arcpy.Describe(inPoints).spatialReference
    arcpy.env.outputCoordinateSystem = sr
    oidField = arcpy.Describe(inLines).OIDFIeldName
    lnGeoms = {ln[0]:ln[1].projectAs(sr) for ln in arcpy.da.SearchCursor(inLines,[oidField,'SHAPE@'])}
    
    #add s-distance (distance along line) field
    print ('Preparing s-n distance index table...')
    arcpy.management.JoinField(ptTabl, 'IN_FID', inPoints, "OBJECTID", "elevation")
    calcSdist(ptTabl,lnGeoms)
    
    #set up output raster
    ptXY = arcpy.da.FeatureClassToNumPyArray(inPoints,['SHAPE@X','SHAPE@Y'])
    llc = arcpy.Point(int(ptXY['SHAPE@X'].min())-gridSize/2, int(ptXY['SHAPE@Y'].min())-gridSize/2)
    exes = np.arange(llc.X+gridSize/2,int(ptXY['SHAPE@X'].max())+gridSize/2,gridSize)
    whys = np.arange(llc.Y+gridSize/2,int(ptXY['SHAPE@Y'].max())+gridSize/2,gridSize)[::-1]
    outGrid = np.empty((len(whys),len(exes)))
    outGrid *= np.nan
    X, Y = np.meshgrid(exes,whys)
    
    #set up calculation points in raster from convex hull
    print ('Preparing convex hull of input points...')
    if xHull =='calc':
        xHull = arcpy.Multipoint(arcpy.Array([arcpy.Point(a[0],a[1]) for a in ptXY]), sr).convexHull()
    else:
        xHull = arcpy.CopyFeatures_management(xHull,arcpy.Geometry())[0]
    npHull = np.vectorize(xHull.contains)
    npPoint = np.vectorize(arcpy.Point)
    ptGrid = npPoint(X,Y)
    hullGrid = npHull(ptGrid)
    print ('Preparing grid points...')
    gridPts = [arcpy.PointGeometry(pt, sr) for pt in ptGrid[hullGrid]]
    gridTabl = arcpy.analysis.GenerateNearTable(gridPts, inLines,'gridSnap',location='LOCATION',angle='ANGLE',
                                    closest='ALL',method='PLANAR',closest_count=2)
    calcSdist(gridTabl,lnGeoms)
    arcpy.management.AddField(gridTabl,'mlpidw','DOUBLE')
    allPoints = arcpy.da.TableToNumPyArray(ptTabl,['NEAR_FID','NEAR_ANGLE','S_DIST','NEAR_DIST','elevation'])
    pointDF = pd.DataFrame(allPoints)
    gridPoints = arcpy.da.UpdateCursor(gridTabl,['NEAR_FID','NEAR_ANGLE','S_DIST','NEAR_DIST','mlpidw'])
    print ('Calculating grid points...')
    for pt in gridPoints:
        pt[4] = calcLPIDW(pointDF,pt)
        gridPoints.updateRow(pt)
    gridDF = pd.DataFrame(arcpy.da.TableToNumPyArray(gridTabl,['IN_FID','NEAR_DIST','mlpidw'],null_value=np.nan))
    vals = []
    print ('Inverse Distance Weighting grid points...')
    for p in range(len(gridPts)):
        results = gridDF[gridDF.IN_FID==(p+1)]
        if results.mlpidw.count() == 2:
            elevs = results.mlpidw.tolist()
            dists = results.NEAR_DIST.tolist()
            idw = (elevs[0]/dists[0] + elevs[1]/dists[1])/(1/dists[0] + 1/dists[1])
            vals.append(idw)
        elif results.mlpidw.count() == 1:
            vals.append(results[~(results.mlpidw == np.nan)].mlpidw.tolist()[0])
        else:
            vals.append(np.nan)
    outGrid[hullGrid] = vals
    outSurf = arcpy.NumPyArrayToRaster(outGrid,llc,gridSize,gridSize,np.nan)
    outSurf.save(outSurface)