8. Advanced mapping with Ipyleaflet

In the previous lessons, we have covered how to interact with and display both raster and vector data, generating new data on the fly and using user-set parameters to determine what will be displayed.

We are now going to use both raster and vector data together, creating interaction between different kinds of data. We will use the vector data that we worked with earlier to calcluate zonal statistics of the heatwave forecast raster, and then allow users to interact with this data in table format, showing some charts along the way.

This might sound much more complicated than what we have done so far, but we're just going to layer the tools and methods we have already used together in order to do this. The goal will be to create a map which allows users to explore heatwave hazard by census area. We will adjust the average temperature of an area by its population, to see the total heat wave hazard in that area.

First, lets re-define the raster functions from the last lesson here:

from io import BytesIO
from base64 import b64encode
import numpy as np
import rasterio as rio
import matplotlib
import matplotlib.cm as cm
import ipywidgets as widgets
from ipyleaflet import Map, Marker, ImageOverlay
import rasterio as rio
from PIL import Image

center = (53.8008, -1.5491)
m = Map(center=center, zoom=15)

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()
raster_reader = rio.open('Data/leeds_NDBI_aug_highres.tif')
raster_transform = raster_reader.transform
bounds = raster_reader.bounds
SW = (bounds.bottom, bounds.left)
NE = (bounds.top, bounds.right)
bounds_tuple = (SW, NE)




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


def encode_array(array, filename='temp.png'):
    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


def image_to_8bit(array):
    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)
    return array 

Ok. In the previous lessons we created two different systems:

1. A system for loading and displaying vector data, and varying this by parameter
2. A system for doing the same with raster data
3. A system for generating new raster data, and displaying the results

So, if we want to show the total heat exposure for each census area, we need to also:

1. Take the average value of temperatue for each vector feature
2. Multiply this by the population of each feature
3. Update the display and show a table

Lets begin by visualising and interacting with some tabular data - the population and other attributes from the census vector.

Ipywidgets doesn't have a built in way to easily display dataframes, but can make use of the flexible and general purpose ipython display system using the output widget, in much the same way we used the output widget to display matplotlib charts in the earlier lessons.

import geopandas as gpd
vector = gpd.read_file("Data/leeds_census_final.geojson", driver='GeoJSON')

frameDis = widgets.Output()

def display_chart(vector):
    with frameDis:
        display(vector, clear=True)  

Ok, this is great - we can also update the display and show new data using the same widgets.

import random
random_data_button = widgets.Button(description= 'Randomise!')
def randomise(b):
    vector['Random'] = vector['Population'] * random.random()
    display(display_chart(vector))


random_data_button.on_click(randomise)

display(random_data_button)
display(frameDis)

If we want to summarise a raster dataset by geometries, we need calculate what are called 'Zonal Statistics', we can use the rasterstats package to do this. In this project we have been working on creating a heatwave forecasting app, so lets start to examine the differences by census area. We're going to do some really overly simplistic risk modelling, by adjusting the average temperature of an area by its population. We can do this by taking the average temperature value for an area, multiplying it by the population, and then dividing by the maximum value for this new column to normalise the data.

from rasterstats import zonal_stats
def hazard_score(vector, raster, transform):
    stats = zonal_stats(vector, np.squeeze(raster), affine=transform)
    mean_list = []
    for s in stats:
        mean_list.append(s['mean'])
        
    vector['Mean Temp'] = mean_list
    vector['Hazard Score'] = vector['Mean Temp']*vector['Population']
    # Now we'll normalise the hazard score 
    vector['Hazard Score'] /= vector['Hazard Score'].max()
    return vector

    

Lets test this out with a randomly generated air temperature raster.

airtemp_array = model_airtemp(10, ndvi_array, ndbi_array, ndwi_array, 13)

hazard_score(vector, airtemp_array, raster_transform)

Because we have stored (as in the earlier lesson) the vector as a geodataframe, we can simply add a new column containing the calculated mean to the dataframe, making it easy to view later.

Now - lets combine this with our previous map, and use an observer to update the data of the hazard score with the new air temperature data that we display to the map.

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, but now we add
# the frame display widget as well 

frameDis = widgets.Output()


def update_map(change):
    #Remove old layers
    try:
        for layer in m.layers[1:]:
            m.remove_layer(layer)
    except:
        pass    

    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)

    ## Now we add the code to update the dataframe
    # If we declare a variable in a function, we cant
    # access a global variable with the same name
    # so we are assigning the new vector to _vector
    _vector = hazard_score(vector, new_airtemp, raster_transform)
    display_chart(_vector)

    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(update_map, 'value')
ndvi_slider.observe(update_map, 'value')
ndbi_slider.observe(update_map, 'value')
ndwi_slider.observe(update_map, '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])
#We will add the frame display beneath the map
mapChartLayout = widgets.VBox([mapLayout, frameDis])
display(mapChartLayout)

And that's it! We have added automatic interaction between the two different systems that we have built over the entire course, and displayed it in a multi-model format, displaying both a raster map and a table.

The next step is up to you! Alongside the table, there are a few things for you to add to make the display of data more interesting.

Exercise 1

Using the skills you learned in the earlier lessons of this course, please replace the table with a histogram chart which shows the distribution of average temperature for each feature in the vector dataset, and update it when parameters of the function are changed.

Exercise 1 Solution

Update the display_chart function to the following:

import matplotlib.pyplot as plt

def display_chart(_v):
    with frameDis:
        frameDis.clear_output()
        _h = plt.hist(_v['Hazard Score'])
        plt.show(_h)