7. Mapping with Ipyleaflet Rasters

In the last lesson, we worked with vector data, adding it to the map and displaying changing outputs according to some user-set parameters.

In this lesson, we will create a new notebook and a new map, and use this to work with some raster data. We will be loading in some remote sensing data covering Leeds, and will create an (unscientific!) model of air temperature which depends on solar irradiation. This is to demonstrate the process of building a model based app, rather than to create a model that accurately forecasts urban heat.

If you are used to wokring with geospatial data, you have likely encountered geotiffs before.

If not, geotiffs are .tif images with some geospatial information attached to them, defining their extent, location, coordinates reference system, and so on. This geospatial information allows the image to be drawn on the map in the right place and at the right scale. These images encode raster data.

Normal images on a computer, when opened in Python, are stored in a 3 dimensional array with the shape (Height, Width, Channels). The number of channels is almost always 3, for RGB images, or sometimes 1, for greyscale images. RGB images represent their colours as integer values between 0 and 255 for each colour. So a pixel with a value of 255 for red and 255 for blue, but 0 for green will be displayed as purple on the screen, and so on.

Raster data, in the geospatial context, can encode many kinds of information. This includes satellite or aerial imagery, satellite instrument data (e.g. radar), or gridded data like elevation, air temperature, and so on. Because of the wide variety of types of information that can be stored in .tif format, the number of channels can be very large. Hyperspectral images, for example, can have up to 220 bands each representing a small part of the electromagnetic spectrum. The values can also be very large, or very small, so dont easily map onto the 0-255 range of png or jpeg images.

This makes displaying .tif files on the web a bit problematic. If a tif file has more than 3 channels, and/or values which go beyond 255 (or are negative!), then the values need to be converted to the correct format to be displayed.

In addition to this, even if your .tif file has the right values and format, most web browsers cannot display them, so they must be converted to a .png or .jpg file.

#We will use the rasterio library to load and work with the geotiff data
from ipyleaflet import Map, Marker, ImageOverlay
import rasterio as rio


center = (53.8008, -1.5491)

m = Map(center=center, zoom=15)

ndvi_reader = rio.open('Data/leeds_NDVI_aug_highres.tif')

metadata = ndvi_reader.profile
bounds = ndvi_reader.bounds

# ipyleaflet cant work with geotiffs directly, we have to do some covnersions before we can place the image on the map
# First, we need to convert the bounds to the correct format so the map knows where to draw the image

SW = (bounds.bottom, bounds.left)
NE = (bounds.top, bounds.right)
bounds_tuple = (SW, NE)

Next, we need to convert the image to a format that the map can understand. We will use the rasterio function to convert the image to a numpy array. Then, we will convert it to a png using the PIL library. This is because most browsers cannot show TIFF files, so adding a TIFF file will mean we can't see what is going on. This is a bit of a cumbersome process, but is unfortunately necessary in order to display tiff images in a web browser.

import numpy as np
from PIL import Image

array = ndvi_reader.read()
array = np.moveaxis(array, 0, -1)

nan_mask = ~np.isnan(array) * 1 
nan_mask *= 255
nan_mask = nan_mask.astype(np.uint8)
array = np.nan_to_num(array)
# We need to move the axis of the array so that the color channel is the last axis

# We need to scale the values in the array to be between 0 and 255 to be viewable as an image
array_max = np.max(array)
array_min = np.min(array)
array = np.clip((array - array_min) / (array_max - array_min) * 255, 0, 255)
array = array.astype(np.uint8)

Ok - now we have an array in the right format and shape. We now need to add this to the map.

Leaflet is a web mapping tool first and foremost, so it doesn't natively work with arrays and other Python datatypes. Leaflet can take a URL to an image or file and display that, but URLs can also encode data in string format. If you encode an image in the URL, and pass this to Leaflet, then Leaflet can show the image without requiring you to manage files or storage.

# Now we can convert the array to an image
from io import BytesIO
from base64 import b64encode

image = Image.fromarray(np.squeeze(np.stack([array, array, array, nan_mask], axis=-1)), mode="RGBA")

s = BytesIO()

# We must convert the image to a series of bytes
# and encode this image in a url, which we can then
# pass to the Leaflet map 

image.save(s, 'png')
data = b64encode(s.getvalue())
data = data.decode('ascii')
imgurl = 'data:image/png;base64,' + data

image_layer = ImageOverlay(url=imgurl, bounds=bounds_tuple)
m.add_layer(image_layer)
display(m)

There we have it - we've added raster data to the map!

There are many web apps which basically just allow the user to load and view many different kinds of data, so these tools can be used to create interesting apps by themselves, but we will add more interaction soon.

def add_tiff_to_map(path, map):
    reader = rio.open(path)

    metadata = reader.profile
    bounds = reader.bounds

    SW = (bounds.bottom, bounds.left)
    NE = (bounds.top, bounds.right)
    bounds_tuple = (SW, NE)

    array = reader.read()
    array = np.moveaxis(array, 0, -1)

    nan_mask = ~np.isnan(array) * 1 
    nan_mask *= 255
    nan_mask = nan_mask.astype(np.uint8)
    array = np.nan_to_num(array)


    array_max = np.max(array)
    array_min = np.min(array)
    array = np.clip((array - array_min) / (array_max - array_min) * 255, 0, 255)
    array = array.astype(np.uint8)

    image = Image.fromarray(np.squeeze(np.stack([array, array, array, nan_mask], axis=-1)), mode="RGBA")
    
    #Convert the image to bytes and encode in the url
    s = BytesIO()
    image.save(s, 'png')
    data = b64encode(s.getvalue())
    data = data.decode('ascii')
    imgurl = 'data:image/png;base64,' + data

    # Now we can add the image to the map
    image_layer = ImageOverlay(url=imgurl, bounds=bounds_tuple)
    map.add_layer(image_layer)
    
    # Close the reader
    reader.close()

Now, using the functions we have defined earlier, we will create a simple app, which shows a map with a dropdown menu that lets users select which image to display.

import ipywidgets as widgets


image_path_dictionary = {'NDVI':"Data/leeds_NDVI_aug_highres.tif",
 'NDBI': "Data/leeds_NDBI_aug_highres.tif", 
 'NDWI': "Data/leeds_NDWI_aug_highres.tif"}


image_type_dropdown = widgets.Dropdown(options=image_path_dictionary.keys(), description='Raster Type:')

def update_image(change):
    
    add_tiff_to_map(image_path_dictionary[image_type_dropdown.value], m)

image_type_dropdown.observe(update_image, names='value')
display(m)
display(image_type_dropdown)

Through using the function which adds an image to the map, along with the observer property which can be applied to the dropdown menu, the map can be updated automatically when the input is changed by the user.

Now you have an interactive map where users can choose which data to view and explore areas interesting to them. One trick we are employing here is to just load and unload data, rather than calculating NDVI or similar indices on the fly using remote sensing images. You can apply this logic to all kinds of analysis, if your model has a relatively small number of parameters, and it is possible for you to pre-compute all the options beforehand, it might be better for your app to do this and just let users view the data interactively. This saves a lot of computing resources on mobile and web devices which might not be very powerful. It also means users dont have to wait a long time for servers to fetch and return data, perform analysis, crunch numbers, and so on.

Now, lets add some modelling of data which we can do on the fly, making this app a bit more dynamic. We are going to define a function, which takes the images we are using as input data and produces an image representing air temperature as an output of the model. We are going to use a model based on solar irradiation, where air temperature of a pixel is a linear combination of NDVI (vegetation index), NDBI (built up index), NDWI (water index), and a constant solar irradiation value. We will encode these predictions in an array, and display it as an image in the same way we displayed the NDVI image before.

import numpy as np

# First, lets define a really simple model for the air temperature of point on the land surface. We will model this as a linear combination of 
# solar irradiance, NDVI, NDBI, and NDWI. We will also add some random noise to the model to make it more realistic. 
# This relies on the intuition of the urban heat island effect, where the emissivity and other characteristics of the land surface
# affect temperature. 

def model_airtemp(solar_irradiance, ndvi, ndbi, ndwi, c=25, ndvi_beta=-3, ndbi_beta=4, ndwi_beta=-2):
    airtemp =  ndvi_beta*ndvi + ndbi_beta*ndbi + ndwi_beta*ndwi + np.random.normal(-1, 1) + c + solar_irradiance
    return airtemp

# Fun question - why do I add the solar_arradiance 
# to the model at the end?

Now, we will create a map as we have done before, but this time we will use the dynamically created array as the data for the image, and add that to the map.

# First lets create a function which is a bit more general than the 
# tiff reading function from earlier, which allows any array
# to be encoded as a URL

def encode_array(array):
    array = np.moveaxis(array, 0, -1)
    nan_mask = ~np.isnan(array) * 1 
    nan_mask *= 255
    nan_mask = nan_mask.astype(np.uint8)
    array_max = np.nanmax(array)
    array_min = np.nanmin(array)


    array = np.nan_to_num(array)


    array = np.clip((array - array_min) / (array_max - array_min) * 255, 0, 255)
    array = array.astype(np.uint8)
    
    image = Image.fromarray(np.squeeze(np.stack([array, array, array, nan_mask], axis=-1)), mode="RGBA")

    
    #Convert the image to bytes and encode in the url
    s = BytesIO()
    image.save(s, 'png')
    data = b64encode(s.getvalue())
    data = data.decode('ascii')
    imgurl = 'data:image/png;base64,' + data
    
    return imgurl

Now, we can use this function to create and add images to maps on the fly, updating dynamically as the user interacts with the application.

#Load the rasters into memory so they can be more easily accessed
ndvi_array = rio.open('Data/leeds_NDVI_aug_highres.tif').read()
ndwi_array = rio.open('Data/leeds_NDWI_aug_highres.tif').read()
ndbi_array = rio.open('Data/leeds_NDBI_aug_highres.tif').read()

solar_irradiance_slider = widgets.FloatSlider(value=14, min=10, max=20, step=0.1, description='Solar Irradiance')
ndvi_slider = widgets.FloatSlider(value=0, min=-10, max=10, step=0.01, description='NDVI Coefficient')
ndbi_slider = widgets.FloatSlider(value=0, min=-10, max=10, step=0.01, description='NDBI Coefficient')
ndwi_slider = widgets.FloatSlider(value=0, min=-10, max=10, step=0.01, description='NDWI Coefficient')

# As before, we will define a callback function which attach
# to the observer of the leaflet map


def updateMap(change):

    new_airtemp = model_airtemp(solar_irradiance_slider.value, ndvi_array, ndbi_array, ndwi_array, 13, ndvi_slider.value, 
                  ndbi_slider.value, ndwi_slider.value)
    
    path = encode_array(new_airtemp)



    imageLayer = ImageOverlay(url=path, bounds=bounds_tuple)
    m.add_layer(imageLayer)
    
# Here we add observers to each of the control widgets we have
# created, so that whenever the user updates a parameter
# of the map, the data displayed on the map is updated

solar_irradiance_slider.observe(updateMap, 'value')
ndvi_slider.observe(updateMap, 'value')
ndbi_slider.observe(updateMap, 'value')
ndwi_slider.observe(updateMap, 'value')


# Now we lay out the widgets, and display them

irradiance_slider_container = widgets.Box([solar_irradiance_slider])
parameters_container = widgets.VBox([ndvi_slider, ndbi_slider, ndwi_slider])

controls_container = widgets.VBox([irradiance_slider_container, parameters_container])
mapDisplay = widgets.Output()

mapLayout = widgets.HBox([m, controls_container])

display(mapLayout)

And that's it! You have created a map which will update and dislay data in real time according to user input. This is the basis of how web apps can work, and by adding more layers of interaction and computation you can build more complex and useful applications. In the next lesson, we'll learn how to take your applications and upload them for others to share. Before that though, there are some exercises for you to complete:

Exercise 1

You might notice that this app slows down over time, have a look at how the code is adding information to the map and re-write the update map function to avoid this problem. Hint: use the m.layers property and the m.remove_layer method.

Exercise 1 Solution

def updateMap(change):

    new_airtemp = model_airtemp(solar_irradiance_slider.value, ndvi_array, ndbi_array, ndwi_array, 13, ndvi_slider.value, 
                  ndbi_slider.value, ndwi_slider.value)
    
    path = encode_array(new_airtemp)

    for l in m.layers[1:]:
        m.rmove_layer(l)
        
    imageLayer = ImageOverlay(url=path, bounds=bounds_tuple)
    m.add_layer(imageLayer)