<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 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
Step1:  Download a raster file representing the flow accumulation for South America. This gives an idea of the river network.
Step2:  We transform the data a bit so that rivers appear thicker.
Step3:  The original data is in the WGS 84 projection, but Leaflet uses Web Mercator, so we need to reproject.
Step4:  Let's convert our NumPy array to an image. For that we must specify a colormap (here plt.cm.jet).
Step5:  The image is embedded in the URL as a PNG file, so that it can be sent to the browser.
Step6:  Finally we can overlay our image and if everything went fine it should be exactly over South America.
Step7:  You can play with the opacity slider and check that rivers from our data file match the rivers on OpenStreetMap.
<ASSISTANT_TASK:>
Python Code:
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
Explanation: 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.
End of explanation
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')]
adffile = name + '/' + name + '/w001001.adf'
if not os.path.exists(adffile):
    r = requests.get(url, stream=True)
    with open(filename, 'wb') as f:
        total_length = int(r.headers.get('content-length'))
        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('.')
Explanation: Download a raster file representing the flow accumulation for South America. This gives an idea of the river network.
End of explanation
dataset = rasterio.open(adffile)
acc_orig = dataset.read()[0]
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
circle = np.zeros((2*radius+1, 2*radius+1)).astype('uint8')
y, x = np.ogrid[-radius:radius+1,-radius:radius+1]
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 = acc[::shrink, ::shrink]
Explanation: We transform the data a bit so that rivers appear thicker.
End of explanation
# 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
Explanation: The original data is in the WGS 84 projection, but Leaflet uses Web Mercator, so we need to reproject.
End of explanation
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))
acc_mask = np.where(np.isfinite(acc_web), 255, 0)
mask = PIL.Image.fromarray(np.uint8(acc_mask), mode='L')
im = PIL.Image.new('RGBA', acc_norm.shape[::-1], color=None)
im.paste(acc_im, mask=mask)
Explanation: Let's convert our NumPy array to an image. For that we must specify a colormap (here plt.cm.jet).
End of explanation
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
Explanation: The image is embedded in the URL as a PNG file, so that it can be sent to the browser.
End of explanation
b = dataset.bounds
bounds = [(b.bottom, b.left), (b.top, b.right)]
io = ImageOverlay(url=imgurl, bounds=bounds)
center = [-10, -60]
zoom = 2
m = Map(center=center, zoom=zoom, interpolation='nearest')
m
tile = basemap_to_tiles(basemaps.Esri.WorldStreetMap)
m.add_layer(tile)
Explanation: Finally we can overlay our image and if everything went fine it should be exactly over South America.
End of explanation
m.add_layer(io)
io.interact(opacity=(0.0,1.0,0.01))
Explanation: You can play with the opacity slider and check that rivers from our data file match the rivers on OpenStreetMap.
End of explanation
<END_TASK>