使用 Numpy 转换为 Web Mercator答案

【问题标题】：Convert to Web Mercator With Numpy使用 Numpy 转换为 Web Mercator
【发布时间】：2014-07-31 12:29:12
【问题描述】：

我的程序垂直拉伸了一个 Numpy 数组，表示一个 180 x 360 的地图图像，因此它表示一个 Web 墨卡托地图图像。

我写了一个函数（如下），它可以做我想做的事 - 但它非常慢（大约需要五分钟）。有没有更快更简单的方法来做到这一点？也许使用 Numpy interpolate2d 或 MatPlotLib？

def row2lat(row):
  return 180.0/math.pi*(2.0*math.atan(math.exp(row*math.pi/180.0))-math.pi/2.0)

def mercator(geodetic):
    geo = np.repeat(geodetic, 2, axis=0)
    merc = np.zeros_like(geo)
    side = geo[0].size
    for row in range(side):
        lat = row2lat(180 - ((row * 1.0)/side) * 360)
        g_row = (abs(90 - lat)/180)*side
        fraction = g_row-math.floor(g_row)
        for col in range(side):
            high_row = geo[math.floor(g_row)][col] * (fraction)
            low_row = geo[math.ceil(g_row)][col] * (1-fraction)
            merc[row][col] = high_row + low_row
    return merc

【问题讨论】：

没有示例输入很难理解（至少对我而言）。然而，使 numpy 快速的经验法则不是使用 for 循环，而是对所有内容进行矢量化。
如果措辞令人困惑，我们深表歉意。任务很简单；垂直拉动一个 Numpy 数组，使其在顶部和底部延伸得更多。
您只想在数组的顶部和底部附加零吗？或者调整它的大小（就像你调整图片大小一样）？
顶部和底部像素的拉伸大约是它们高度的 3 倍，中心像素越多拉伸越少。就像印在焦糖上的图片。
我添加了显示函数作用的图像。

标签： python numpy gis matplotlib-basemap

【解决方案1】：

尽量避免内部 for 循环并矢量化您的函数。 Numpy 经过高度优化，可以高效地运行这些东西。然后你的函数会读起来像

def mercator_faster(geodetic):
    geo = np.repeat(geodetic, 2, axis=0)
    merc = np.zeros_like(geo)
    side = geo[0].size
    for row in range(side):
        lat = row2lat(180 - ((row * 1.0)/side) * 360)
        g_row = (abs(90 - lat)/180)*side
        fraction = g_row-math.floor(g_row)

        # Here I optimized the code by using the numpy vector operations 
        # instead of the for loop:

        high_row = geo[math.floor(g_row), :] * (fraction)
        low_row = geo[math.ceil(g_row), :] * (1-fraction)
        merc[row, :] = high_row + low_row

    return merc

如果我在我的机器上运行它只需要不到一秒钟：

%timeit mercator_faster(geo)
1 loops, best of 3: 727 ms per loop

它看起来像这样（我不得不重新调整它，因为它对 SO 来说太大了）：

可能外部 for 循环也可能被矢量化，但我想这要困难得多。

【讨论】：

嘿！我喜欢这个想法，但是这个函数的输出与原始函数的输出不同。
好的。我也添加了一张图片:)