I’m a late-in-life learner of code, so please be gentle.

I’m working with code that needs to compute a number of geometric ranges from arbitrary points to surface heights in HAE from a digital elevation model (DEM). This surface height data is provided in WGS84 Lat/Long/MSL elevation, which I can adjust into HAE from a geoid model and store prior to program execution. So far, so good–I can translate the required DEM data for the needed area around the reference point into an array of local-level points (East/North/Up is what I’m using here).
The reference points data I have is all provided in a local-level reference based on a specific ENU origin point in Lat/Long/Height, which I don’t know until runtime. The most significant compute task is an iteration of various distance computations between (E,N,U) reference points in the local-level frame and up to ~100,000 sampled points in the surface reference contour. Not a big problem for numpy.
The problem I’m having is that the DEM data, translated into the ENU frame, is not regular in East/North alignment. The algorithm I’m implementing requires me to return ranges to regularly spaced East/North points, and return computations specific to a raster of East/North surface points
My ENU-space DEM object is a ~6000x7000x3 ndarray with the relevant search area surface locations. Viewed in East/North, it’s warped a bit like a trapezoid that’s wider and slightly bowed in the direction of the equator as a result of the transformation.

The best I’ve been able to do is:

import numpy as np
from scipy.interpolate import NearestNDInterpolator
# Fake DEM data (don't know how to warp it in the way the TED output is, so just adding some rand) This
# could use a rectilinear interpolator, but I believe my data cannot (warped in E as a function of N and 
# vice versa).
e_dim = 7000
n_dim = 6000
e_rng = np.arange(-e_dim/2, e_dim/2, dtype='float32')*72 + np.random.randint(-100, high=101, size=e_dim)
n_rng = np.arange(-n_dim/2, n_dim/2, dtype='float32')*95 + np.random.randint(-100, high=101, size=n_dim)
dem_u = np.random.randint(-100, high=4500, size=(n_dim,e_dim))
row, col = np.meshgrid(n_rng, e_rng, indexing='ij')
ENU = np.dstack((row,col,dem_u))

# looking for a faster approach to lookup Nx2 arrays of E/N points that are on a regular sample grid for # rastered final output

interp = NearestNDInterpolator(
        (ENU[:,:,0].flatten(),ENU[:,:,1].flatten()),
         ENU[:,:,2].flatten())
# times about 12-13s on my hardware--pretty darn slow. 
interval_m = 11338.7 # Just a coarse starting sample dimension
samp_offsets_m = np.array([-150*1852+i*interval_m for i in
                               range(pts_side) ] )
tiles_x, tiles_y = np.meshgrid(samp_offsets_m, samp_offsets_m, indexing='xy' )      
pTL = np.asarray( [[*tiles_x.flatten()],    E pos 
                   [*tiles_y.flatten()]])   N pos

interp(pTL[0], pTL[1])   # pretty fast, about 4ms to update 'Up' on ~2500 samples.

Is there a better way to create a regularly aligned structure to do DEM lookups into in the local-level plane?
Of note, my target environment can’t run GDAL/Rasterio, etc. I can prep the DEM set with those tools, but it isn’t in the runtime environment.

Thanks for any guidance you can provide!
J