Convert to Web Mercator With Numpy
-
21-12-2019 - |
Pregunta
My program vertically stretches a Numpy array, representing a 180 by 360 map image, so it represents a Web Mercator map image.
I wrote a function (below) that does what I want - but it is crazy slow (takes like five minutes). Is there a much faster and easier way to do this? Maybe using Numpy interpolate2d
or 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
Solución
Try to avoid the inner for loop and vectorize your functions. Numpy is highly optimized to run those things efficient. Your function would then read like
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
If I run it on my machine it takes less then a second:
%timeit mercator_faster(geo)
1 loops, best of 3: 727 ms per loop
And it looks like this (I had to rescale it, because it was too big for SO):
Possibly the outer for loop might be vectorized as well, but I guess this is much harder.
Licenciado bajo: CC-BY-SA con atribución
No afiliado a StackOverflow