add row2lat() and geodetic_to_mercator()

Author: andrewgiessel

folium/utilities.py

@@ -344,4 +344,50 @@ def png_pack(png_tag, data):
         png_pack(b'IHDR', struct.pack("!2I5B", width, height, 8, 6, 0, 0, 0)),
         png_pack(b'IDAT', zlib.compress(raw_data, 9)),
         png_pack(b'IEND', b'')])
+
+def _row2lat(row):
+    return 180.0/np.pi*(2.0*np.arctan(np.exp(row*np.pi/180.0))-np.pi/2.0)
+    
+def geodetic_to_mercator(geodetic):
+    """This function takes an 2D array in geodetic coordinates (ie: lat x
+    lon unprojected) and converts it to web mercator.  This is needed
+    to correctly overlay an image on a leaflet map, which uses web
+    mercator by default.
+
+    This code works with arrays that match the relative proportions of
+    latitude and longitude of the earth, meaning that they have twice
+    as many points in longitude as latitude (i.e.: (180, 360) for 1
+    degree resolution).
+
+    The code for this is from:
+    http://stackoverflow.com/questions/25058880/convert-to-web-mercator-with-numpy
+
+    Parameters
+    ----------
+    geodetic: numpy image array
+        Latitude x Longitude array, in mono (NxM), rgb (NxMx3) or RGBA (NxMx4)
+
+    Returns
+    -------
+    meractor: projected numpy image array
+
+    """
+    
+    geo = np.repeat(np.atleast_3d(geodetic), 2, axis=0)
+    merc = np.zeros_like(geo)
+    side = geo.shape[0]
+
+    for row in range(side):
+        lat = _row2lat(180 - ((row * 1.0)/side) * 360)
+        g_row = (abs(90 - lat)/180)*side
+        fraction = g_row-np.floor(g_row)
+
+        # Here I optimized the code by using the numpy vector operations 
+        # instead of the for loop:
+
+        high_row = geo[np.floor(g_row), :, :] * (fraction)
+        low_row = geo[np.ceil(g_row), :, :] * (1-fraction)
+        merc[row, :, :] = high_row + low_row
+
+    return np.squeeze(merc)