Models and Maps with Ipyleaflet

Contents

Models and Maps with Ipyleaflet#

In this next section, we’re going to use all of what you have learned so far to create a heatwave forecasting map which lets users explore future events.

Please note, the model for air temperature in this section is basically totally made up. It makes some intuitive sense but the parameters are basically random, and probably will not give any reasonable estimate, it is purely to demonstrate how to use a model to create data which will go onto a map.

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 ad the end?

import ipywidgets as widgets
from ipyleaflet import Map, ImageOverlay
import rasterio as rio

def add_tiff_to_map(path, map, filename="temp.png"):
    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")
    image.save(filename)


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

import os
center = (53.8008, -1.5491)
m = Map(center=center, zoom=12)

# First, we will load the spectral indices data from the raster files. 

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()

# All three of these arrays have the same shape, so we can use the shape and extent of one of them to define the extent of air temperature data.

image_extent = rio.open('Data/leeds_NDVI_aug_highres.tif').bounds
SW = (image_extent.bottom, image_extent.left)
NE = (image_extent.top, image_extent.right)
bounds_tuple = (SW, NE)

# Now, lets create the airtemperature array which holds the data
airtemp_array = model_airtemp(3, ndvi_array, ndbi_array, ndwi_array)

from PIL import Image

array = np.moveaxis(airtemp_array, 0, -1)
nan_mask = ~np.isnan(array) * 1 
nan_mask *= 255
array_min = np.nanmin(array)

nan_mask = nan_mask.astype(np.uint8)
array = np.nan_to_num(array)


array_max = np.max(array)
array = np.clip((array - array_min) / (array_max - array_min) * 255, 0, 255)
array = array.astype(np.uint8)
try:
    os.remove("airtemp.png")
except:
    pass

image = Image.fromarray(np.squeeze(np.stack([array, array, array, nan_mask], axis=-1)), mode="RGBA")
image.save('airtemp.png')
##NB
## If you pass an image (even if it is different) to the ImageOverlay, it will not update the image if the URL is the same. 
## So, we need to delete the old image if it exists, and then save the new image with a new name if we want to keep updating the map. 

# Now we can add the image to the map
image_layer = ImageOverlay(url="airtemp.png", bounds=bounds_tuple)
m.add_layer(image_layer)
display(m)

Now that we have added the data to the map, its time to put it all together!

Using what you have learned in the previous notebooks, your job is to:

Create a slider to alter solar_irradiance between 10 and 23
Use an observer to calculate a new air temperature raster each time the slider is changed
Clear the map of its old layers, and add the new raster once it is calculated
Create an additional panel for tuning the coefficients of the model: ndvi_beta, ndwi_beta, and ndbi_beta between -5 and 5
Use containers to compose these widgets together so the map is on one side, and controls are on the other

To get started, I will provide a function which automatically converts the air temperature array to a png file

Hint - the full solution is available in the file ‘voila_demo.ipynb’#

def array_to_png(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")
    image.save(filename)
    return filename

## Bonus Exercise! 
# The image added to the map is currently gray-scale, and doesn't contain any information about the scale of the values. 
# Your task is to update the widgets you have created above to generate colour images using a colour map of your choice. 
# We have to manually convert the single-band array to a 3-band colour image.
# We can do that using the following functions
import matplotlib
import matplotlib.cm as cm

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 

def greyscale_to_colour(greyscale_array, colour_map='magma'):
    array = np.nan_to_num(greyscale_array) # Convert NaNs to 0
    array_max = np.max(array)
    array_min = np.min(array)
    # Matplotlib wants the data in range [0, 1], so we need to normalise the data
    normaliser = matplotlib.colors.Normalize(vmin=array_min, vmax=array_max, clip=True)
    scalar_mapper = cm.ScalarMappable(norm=normaliser, cmap=colour_map)
    colour_array = scalar_mapper.to_rgba(np.squeeze(greyscale_array))
    colour_array = image_to_8bit(colour_array)
    return colour_array 

#The rest is up to you

##Bonus Exercise Solution
def array_to_png(array, filename='temp.png', colour_map='None'):
    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)

    #Remove and overwrite the file if it already exists
    try:
        os.remove(filename)
    except:
        pass

    if colour_map != 'None':
        array = greyscale_to_colour(array, colour_map=colour_map)
        array [:,:,-1] = np.squeeze(nan_mask)
        image = Image.fromarray(np.squeeze(array), mode="RGBA")
    else:
        image = Image.fromarray(np.squeeze(np.stack([array, array, array, nan_mask], axis=-1)), mode="RGBA")
   
    image.save(filename)
    return filename

import os
center = (53.8008, -1.5491)
m = Map(center=center, zoom=12)

# First, we will load the spectral indices data from the raster files. 

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()

# All three of these arrays have the same shape, so we can use the shape and extent of one of them to define the extent of air temperature data.

image_extent = rio.open('Data/leeds_NDVI_aug_highres.tif').bounds
SW = (image_extent.bottom, image_extent.left)
NE = (image_extent.top, image_extent.right)
bounds_tuple = (SW, NE)

# Now, lets create the airtemperature array which holds the data
airtemp_array = model_airtemp(3, ndvi_array, ndbi_array, ndwi_array)

from PIL import Image

array = np.moveaxis(airtemp_array, 0, -1)
filename = 'new_airtemp.png'
filename = array_to_png(array, filename, colour_map='magma')


##NB
## If you pass an image (even if it is different) to the ImageOverlay, it will not update the image if the URL is the same. 
## So, we need to delete the old image if it exists, and then save the new image with a new name if we want to keep updating the map. 

# Now we can add the image to the map
image_layer = ImageOverlay(url=filename, bounds=bounds_tuple)
m.add_layer(image_layer)
display(m)