2

我有带有 SRTM 高程数据的 numpy 数组，我需要计算路径配置文件 - 带有从一个点（发射器）到地图上每个点的高程数据的向量。这是我目前的解决方案：

def get_elev_vector(self, coord1, coord2, step=90):
    if coord1 == coord2:
        return 0
    lat1, long1 = coord1
    lat2, long2 = coord2
    dist = self.get_dist(coord1, coord2)
    npt = dist / step
    vect = np.zeros(npt + 3, dtype='float32')
    vect[0] = npt
    vect[1] = dist
    vect[2] = self.get_elevation_direct(coord1)
    vect[npt] = self.get_elevation_direct(coord2)
    g = pyproj.Geod(ellps='WGS84')
    coordvect = g.npts(long1, lat1, long2, lat2, npt)
    i = 3
    for lon, lat in coordvect:
        vect[i] = self.get_elevation((lat, lon))
        i += 1
    return vect


#(tx, ty) transmission coordinates
while lat > latmax:
         long = longmin
         while long<longmax:
               vect =  self.carta.get_elev_vector((tx, ty), (lat, lon))
               arr[xarr, yarr] = 100
               yarr +=1
               lon += dlon
         xarr +=1
         lat +=dx

它可以工作，但对于 20x20 公里的区域大约需要 3 个小时。有没有办法让它更快地工作？

添加

    def get_elevation(self, coord):
    """TODO"""
    lat, lon = self.xyfromlatlon(coord)
    mlat = int(lat)
    mlon = int(lon)
    vicin = self.tile['numpy_array'][mlat: mlat + 2, mlon: mlon + 2]
    temp_bool = np.equal(vicin, -32768)
    np.putmask(vicin, temp_bool, 0)
    elev = billin_interpol(vicin[0][0], vicin[0][1], vicin[1][0],
                                        vicin[1][1], lat - mlat, lon - mlon)
    return int(elev + 0.5)


def billin_interpol(tl, tr, bl, br, x, y):
"""Bilinear interpolation.
   Return:
    float, interpolated value
"""
b1 = tl
b2 = bl - tl
b3 = tr - tl
b4 = tl - bl - tr + br
return b1 + b2 * x + b3 * y + b4 * x * y

完整的代码和数据

4

0 回答 0