Displaying satellite imagery on a web map b0d1b3884d9740a49f4da79979d5318f

**Keywords** :index:`data used; sentinel-2`, :index:`analysis; interactive map`,index:`ipyleaflet`,

Background

Leaflet is the leading open-source JavaScript library for mobile-friendly interactive maps. Functionality it provides is exposed to Python users by ipyleaflet. This library enables interactive maps in the Jupyter notebook/JupyterLab environment.

Description

This notebook demonstrates how to plot an image and dataset footprints on a map.

  1. Load packages

  2. Find a location

  3. Find some datasets to load

  4. Load pixel data in EPSG:3857 projection, same as used by most web maps

  5. Create dataset footprints to display on a map

  6. Create an opacity control to display on the same map

  7. Display image loaded from the datasets on the same map

Getting started

To run this analysis, run all the cells in the notebook, starting with the “Load packages” cell.

Load packages

[1]:

import os
import ipyleaflet
import numpy as np
from ipywidgets import widgets as w
from IPython.display import display
import matplotlib.pyplot as plt
import matplotlib as mpl

import datacube
import odc.ui
from odc.ui import with_ui_cbk

from deafrica_tools.plotting import display_map
from deafrica_tools.plotting import rgb
from deafrica_tools.datahandling import load_ard

/env/lib/python3.8/site-packages/geopandas/_compat.py:106: UserWarning: The Shapely GEOS version (3.8.0-CAPI-1.13.1 ) is incompatible with the GEOS version PyGEOS was compiled with (3.9.1-CAPI-1.14.2). Conversions between both will be slow.
  warnings.warn(
/env/lib/python3.8/site-packages/datacube/storage/masking.py:7: DeprecationWarning: datacube.storage.masking has moved to datacube.utils.masking
  warnings.warn("datacube.storage.masking has moved to datacube.utils.masking",

Connect to the datacube

[2]:

dc = datacube.Datacube(app='Imagery_web_map')

Find a location

The selected latitude and longitude will be displayed as a red box on the map below the next cell. This map can be used to find coordinates of other places, simply scroll and click on any point on the map to display the latitude and longitude of that location.

[3]:

# Set the central latitude and longitude
central_lat = 7.656
central_lon = 0.021

# Set the buffer to load around the central coordinates
buffer = 0.2

# Compute the bounding box for the study area
lats = (central_lat - buffer, central_lat + buffer)
lons = (central_lon - buffer, central_lon + buffer)

display_map(x=lons, y=lats, margin=-0.2)

[3]:

Find datasets

Use the Digital Earth Africa Explorer or dc.list_products() to find avaliable datasets. For more information on using dc.list_products(), see the Products and measurements notebook.

In this example we are using the Sentinel-2A ARD product. We will be visualizing a portion of the swath taken by Sentinel-2A on 12-Jan-2018.

[4]:

# Define products
products = 's2_l2a'

# Specify the parameters to pass to the load query
query = {
    "x": lons,
    "y": lats,
    "time": ('2018-01-12'),
    "measurements":['red', 'green', 'blue'],
    "output_crs": 'EPSG:6933',
    "resolution": (-10, 10),
    "group_by": "solar_day"
}

# Load the data
ds = load_ard(dc, products=products, **query)

Using pixel quality parameters for Sentinel 2
Finding datasets
    s2_l2a
Applying pixel quality/cloud mask
Loading 1 time steps

Create Leaflet Map with dataset footprints

We want to display dataset footprints as well as captured imagery. Therefore we use dss = dc.find_datasets(..) to obtain a list of datacube.Dataset objects overlapping with our query first.

Then we convert list of dataset objects into a GeoJSON of dataset footprints, while also computing bounding box. We will use the bounding box to set initial viewport of the map.

[5]:

dss = dc.find_datasets(product=products, **query)

polygons, bbox = odc.ui.dss_to_geojson(dss, bbox=True)

Create ipyleaflet.Map with full-screen and layer visibility controls. Set initial view to be centered around dataset footprints. We will not be displaying the map just yet.

[6]:

zoom = odc.ui.zoom_from_bbox(bbox)
center = (bbox.bottom + bbox.top) * 0.5, (bbox.right + bbox.left) * 0.5

m = ipyleaflet.Map(
    center=center,
    zoom=round(zoom*1.2),
    scroll_wheel_zoom=True,  # Allow zoom with the mouse scroll wheel
    layout=w.Layout(
        width='600px',   # Set Width of the map to 600 pixels, examples: "100%", "5em", "300px"
        height='600px',  # Set height of the map
    ))

# Add full-screen and layer visibility controls
m.add_control(ipyleaflet.FullScreenControl())
m.add_control(ipyleaflet.LayersControl())

Now we add footprints to the map.

[7]:

m.add_layer(ipyleaflet.GeoJSON(
    data={'type': 'FeatureCollection',
          'features': polygons},
    style={
        'opacity': 0.3,      # Footprint outline opacity
        'fillOpacity': 0     # Do not fill
    },
    hover_style={'color': 'tomato'},  # Style when hovering over footprint
    name="Footprints"                 # Name of the Layer, used by Layer Control widget
))

Create Leaflet image layer

Under the hood mk_image_layer will:

  1. Convert 16-bit rgb xarray to an 8-bit RGBA image

  2. Encode RGBA image as PNG data odc.ui.to_rgba

  3. Render PNG data to “data uri”

  4. Compute image bounds

  5. Construct ipyleaflet.ImageLayer with uri from step 3 and bounds from step 4

JPEG compression can also be used (e.g fmt="jpeg"), useful for larger images to reduce notebook size in bytes (use quality=40 to reduce size further), downside is no opacity support unlike PNG.

Satellite imagery is often 12-bit and higher, but web images are usually 8-bit, hence we need to reduce bit-depth of the input imagery such that there are only 256 levels per color channel. This is where clamp parameter comes in. In this case we use clamp=2000. Input values of 2000 and higher will map to value 255 (largest possible 8-bit unsigned value), 0 will map to 0 and every other value in between will scale linearly.

[8]:

img_layer = odc.ui.mk_image_overlay(
    ds,
    clamp=2000,  # 2000 -- brightest pixel level
    bands=['red','green','blue'],
    fmt='png')   # "jpeg" is another option

# Add image layer to a map we created earlier
m.add_layer(img_layer)

Add opacity control

  • Add Vertical Slider to the map

  • Dragging slider down lowers opacity of the image layer

  • Use of jslink from ipywidgets ensures that this interactive behaviour will work even on a pre-rendered notebook (i.e. on nbviewer)

[9]:

slider = w.FloatSlider(min=0, max=1, value=1,        # Opacity is valid in [0,1] range
                       orientation='vertical',       # Vertical slider is what we want
                       readout=False,                # No need to show exact value
                       layout=w.Layout(width='2em')) # Fine tune display layout: make it thinner

# Connect slider value to opacity property of the Image Layer
w.jslink((slider, 'value'),
         (img_layer, 'opacity') )
m.add_control(ipyleaflet.WidgetControl(widget=slider))

Finally display the map

[10]:

display(m)

Sharing notebooks online

Unlike notebooks with matplotlib figures, saving a notebook after running it is not enough to have interactive maps displayed when sharing rendered notebooks online. You also need to make sure that “Widget State” is saved. In JupyterLab make sure that Save Widget State Automatically setting is enabled. You can find it under Settings menu.

Additional information

License: The code in this notebook is licensed under the Apache License, Version 2.0. Digital Earth Africa data is licensed under the Creative Commons by Attribution 4.0 license.

Contact: If you need assistance, please post a question on the Open Data Cube Slack channel or on the GIS Stack Exchange using the open-data-cube tag (you can view previously asked questions here). If you would like to report an issue with this notebook, you can file one on Github.

Compatible datacube version:

[11]:

print(datacube.__version__)

1.8.5

Last Tested:

[12]:

from datetime import datetime
datetime.today().strftime('%Y-%m-%d')

[12]:

'2021-08-27'