使用Numpy转换为Web Mercator

Question

我的程序需要把一个Numpy数组进行垂直拉伸，这个数组代表的是一个180乘360的地图图像，最终要变成一个Web Mercator地图图像。

我写了一个函数（见下文），它能完成我想要的效果，但速度非常慢，差不多要花五分钟。有没有更快更简单的方法可以做到这一点呢？比如使用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

Answer 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

结果是这样的（我需要调整一下大小，因为它在StackOverflow上太大了）：

外部的for循环也可能可以进行向量化，不过我觉得这会更难一些。