# From NumPy to Leaflet¶

This notebook shows how to display some raster geographic data in IPyLeaflet. The data is a NumPy array, which means that you have all the power of the Python scientific stack at your disposal to process it.

The following libraries are needed:

• requests
• tqdm
• rasterio
• numpy
• scipy
• pillow
• matplotlib
• ipyleaflet

The recommended way is to try to conda install them first, and if they are not found then pip install.

In [ ]:
import requests
import os
from tqdm import tqdm
import zipfile
import rasterio
from affine import Affine
import numpy as np
import scipy.ndimage
from rasterio.warp import reproject, RESAMPLING
import PIL
import matplotlib.pyplot as plt
from base64 import b64encode
try:
from StringIO import StringIO
py3 = False
except ImportError:
from io import StringIO, BytesIO
py3 = True
from ipyleaflet import Map, ImageOverlay, basemap_to_tiles, basemaps


Download a raster file representing the flow accumulation for South America. This gives an idea of the river network.

In [ ]:
url = 'http://earlywarning.usgs.gov/hydrodata/sa_30s_zip_grid/sa_acc_30s_grid.zip'
filename = os.path.basename(url)
name = filename[:filename.find('_grid')]

r = requests.get(url, stream=True)
with open(filename, 'wb') as f:
for chunk in tqdm(r.iter_content(chunk_size=1024), total=(total_length/1024) + 1):
if chunk:
f.write(chunk)
f.flush()
zip = zipfile.ZipFile(filename)
zip.extractall('.')


We transform the data a bit so that rivers appear thicker.

In [ ]:
dataset = rasterio.open(adffile)
acc = np.where(acc_orig<0, 0, acc_orig)

shrink = 1 # if you are out of RAM try increasing this number (should be a power of 2)
radius = 5 # you can play with this number to change the width of the rivers
index = x**2 + y**2 <= radius**2
circle[index] = 1
acc = np.sqrt(acc)
acc = scipy.ndimage.maximum_filter(acc, footprint=circle)
acc[acc_orig<0] = np.nan
acc = np.array(acc[::shrink, ::shrink])


The original data is in the WGS 84 projection, but Leaflet uses Web Mercator, so we need to reproject.

In [ ]:
# At this point if GDAL complains about not being able to open EPSG support file gcs.csv, try in the terminal:
# export GDAL_DATA=gdal-config --datadir

with rasterio.Env():
rows, cols = acc.shape
src_transform = list(dataset.affine)
src_transform[0] *= shrink
src_transform[4] *= shrink
src_transform = Affine(*src_transform[:6])
src_crs = {'init': 'EPSG:4326'}
source = acc

dst_crs = {'init': 'EPSG:3857'}
dst_transform, width, height = rasterio.warp.calculate_default_transform(src_crs, dst_crs, cols, rows, *dataset.bounds)
dst_shape = height, width

destination = np.zeros(dst_shape)

reproject(
source,
destination,
src_transform=src_transform,
src_crs=src_crs,
dst_transform=dst_transform,
dst_crs=dst_crs,
resampling=RESAMPLING.nearest)

acc_web = destination


Let's convert our NumPy array to an image. For that we must specify a colormap (here plt.cm.jet).

In [ ]:
acc_norm = acc_web - np.nanmin(acc_web)
acc_norm = acc_norm / np.nanmax(acc_norm)
acc_norm = np.where(np.isfinite(acc_web), acc_norm, 0)
acc_im = PIL.Image.fromarray(np.uint8(plt.cm.jet(acc_norm)*255))
im = PIL.Image.new('RGBA', acc_norm.shape[::-1], color=None)


The image is embedded in the URL as a PNG file, so that it can be sent to the browser.

In [ ]:
if py3:
f = BytesIO()
else:
f = StringIO()
im.save(f, 'png')
data = b64encode(f.getvalue())
if py3:
data = data.decode('ascii')
imgurl = 'data:image/png;base64,' + data


Finally we can overlay our image and if everything went fine it should be exactly over South America.

In [ ]:
b = dataset.bounds
bounds = [(b.bottom, b.left), (b.top, b.right)]
io = ImageOverlay(url=imgurl, bounds=bounds)

In [ ]:
center = [-10, -60]
zoom = 2
m = Map(center=center, zoom=zoom)
m

In [ ]:
tile = basemap_to_tiles(basemaps.Esri.WorldStreetMap)


You can play with the opacity slider and check that rivers from our data file match the rivers on OpenStreetMap.

In [ ]:
m.add_layer(io)
io.interact(opacity=(0.0,1.0,0.01))