将 2d DEM 数据重新采样到常规局部级数组中

Question

我需要根据数字高程模型 (DEM) 计算 HAE 中从任意点到表面高度的多个几何范围。该表面高度数据以 WGS84 纬度/经度/MSL 高程提供，我可以将其从大地水准面模型调整为 HAE 并在程序执行之前存储。我可以将参考点周围所需区域所需的 DEM 数据转换为本地点数组（东/北/上是我在这里使用的）。

参考点数据是在基于纬度/经度/高度中的特定 ENU 原点的本地参考中提供的，直到运行时我才知道。最重要的计算任务是局部级框架中的 (E,N,U) 参考点与表面参考轮廓中多达约 100,000 个采样点之间的各种距离计算的迭代。对于 NumPy 来说不是一个大问题。

转换为 ENU 框架的 DEM 数据在东/北对齐方面不规则。我正在实现的算法要求我将范围返回到规则间隔的东/北点，并返回特定于东/北表面点栅格的计算。我的 ENU 空间 DEM 对象是一个 ~6000x7000x3 ndarray，具有相关的搜索区域表面位置。从东/北方向看，它有点像梯形，由于变换而变宽并且在赤道方向上稍微弯曲。

我能做的最好的事情是：

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.

是否有更好的方法来创建规则对齐的结构以在局部平面中进行 DEM 查找？我的目标环境无法运行 GDAL/Rasterio 等。我可以使用这些工具准备 DEM 集，但它不在运行时环境中。

Answer 1

如果您的程序对于您正在执行的任务来说是正确的，并且我们希望保留它，则几乎不需要改进。

正如您所注意到的，构建插值器非常耗时：

%timeit -n 1 -r 30 interp = NearestNDInterpolator((ENU[:,:,0].flatten(), ENU[:,:,1].flatten()), ENU[:,:,2].flatten())
# 11.8 s ± 488 ms per loop (mean ± std. dev. of 30 runs, 1 loop each)

并且在我的机器上花费的时间几乎相同，它们应该具有可比性。

首先，您可以使用

ravel

代替

flatten

，这样可以避免复制数组并节省一些内存。

%timeit -n 1 -r 30 interp = NearestNDInterpolator((ENU[:,:,0].ravel(), ENU[:,:,1].ravel()), ENU[:,:,2].ravel())
# 11.7 s ± 630 ms per loop (mean ± std. dev. of 30 runs, 1 loop each)

但是它并没有带来任何时间上的改善，这只是关于记忆。

现在我们必须考虑这样一个事实：这一步花费的大部分时间都在构建底层的 KD-Tree。幸运的是，

NearestNDInterpolator

让您有机会设置如何构建这个

cKDTree

。

%timeit -n 1 -r 30 interp = NearestNDInterpolator((ENU[:,:,0].ravel(), ENU[:,:,1].ravel()), ENU[:,:,2].ravel(), tree_options={"compact_nodes": False, "balanced_tree": False, "leafsize": 64})
# 3.67 s ± 254 ms per loop (mean ± std. dev. of 30 runs, 1 loop each)

制作较粗的树大大减少了构建时间，但可能会增加预测时间。

如果我们比较的话，预测时间：

%timeit -n 1 -r 30 interp(pTL[0], pTL[1])
# Default tree: 4.36 s ± 149 ms per loop (mean ± std. dev. of 30 runs, 1 loop each)
# Tweaked tree: 5.47 s ± 139 ms per loop (mean ± std. dev. of 30 runs, 1 loop each)

较粗的树需要大约 1 秒来预测，但需要 8 秒来构建。这不是 10 或 100 因子优化，而是保持您的程序尽可能通过调整树配置将执行时间减少约 2 倍。