Displaying satellite imagery on a web map

  • Compatibility: Notebook currently compatible with both the NCI and DEA Sandbox environments

  • Products used: s2b_ard_granule

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. Find some datasets to load

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

  3. Display dataset footprints on a map

  4. Display image loaded from these 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 datacube
import odc.ui
from odc.ui import with_ui_cbk

Configure datacube for efficient S3 access

This step is optional. Hopefully in the future it will happen automatically. Cell below configures internal libraries used by Open Data Cube (GDAL and rasterio) in a way that allows more efficient data loading in the cloud environment (i.e. reading Cloud Optimized GeoTIFFs from S3).

[2]:
from datacube.utils.rio import set_default_rio_config

# Only run if executing in the cloud, will successfully do nothing on NCI
if 'AWS_ACCESS_KEY_ID' in os.environ:
    set_default_rio_config(aws={'region_name': 'auto'},
                           cloud_defaults=True)

Connect to the datacube

[3]:
dc = datacube.Datacube(app='Imagery_on_web_map')

Find datasets

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

We want to display dataset footprints as well as captured imagery. Rather than calling dc.load(..) directly with the time and spatial bounds we first use

dss = dc.find_datasets(..)

to obtain a list of datacube.Dataset objects overlapping with our query first.

[4]:
# Define product and red/green/blue bands in the given product
product = 's2b_ard_granule'
RGB = ('nbar_red', 'nbar_green', 'nbar_blue')

# Region and time of interest
query = dict(
    lat=(-30, -36),
    lon=(137, 139),
    time='2018-01-13',
)

dss = dc.find_datasets(product=product, **query)
print(f"Found {len(dss)} datasets")
Found 7 datasets

Load red/green/blue bands

Since we already have a list of datasets (dss) we do not need to repeat the query, instead we supply datasets to dc.load(.., datasets=dss, ..) along with other parameters used for loading data. Note that since we do not supply lat/lon bounds we will get all the imagery referenced by the datasets found earlier and the result will not be clipped to a lat/lon box in the query above.

We will load imagery at 200 m per pixel resolution (1/20 of the native) in the Pseudo-Mercator (EPSG:3857) projection, same as used by most webmaps.

[5]:
rgb = dc.load(
    product=product,             # dc.load always needs product supplied, this needs to be fixed in `dc.load` code
    datasets=dss,                # Datasets we found earlier
    measurements=RGB,            # Only load red,green,blue bands
    group_by='solar_day',        # Fuse all datasets captured on the same day into one raster plane
    output_crs='EPSG:3857',      # Default projection used by Leaflet and most webmaps
    resolution=(-200, 200),      # 200m pixels (1/20 of the native)
    resampling='bilinear',       # Use bilinear resampling when scaling down
    progress_cbk=with_ui_cbk())  # Display load progress
rgb
[5]:
<xarray.Dataset>
Dimensions:     (time: 1, x: 989, y: 2436)
Coordinates:
  * time        (time) datetime64[ns] 2018-01-13T00:56:59.027000
  * y           (y) float64 -3.478e+06 -3.478e+06 ... -3.965e+06 -3.965e+06
  * x           (x) float64 1.514e+07 1.514e+07 ... 1.534e+07 1.534e+07
Data variables:
    nbar_red    (time, y, x) int16 -999 -999 -999 -999 ... -999 -999 -999 -999
    nbar_green  (time, y, x) int16 -999 -999 -999 -999 ... -999 -999 -999 -999
    nbar_blue   (time, y, x) int16 -999 -999 -999 -999 ... -999 -999 -999 -999
Attributes:
    crs:      EPSG:3857

Create Leaflet Map with dataset footprints

First 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.

[6]:
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.

[7]:
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=zoom,
    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.

[8]:
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=3000. Input values of 3000 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.

[9]:
img_layer = odc.ui.mk_image_overlay(
    rgb,
    clamp=3000,  # 3000 -- brightest pixel level
    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)

[10]:
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

[11]:
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 Australia 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.

Last modified: December 2019

Compatible datacube version:

[12]:
print(datacube.__version__)
1.7+142.g7f8581cf.dirty

Tags

Browse all available tags on the DEA User Guide’s Tags Index

Tags: NCI compatible, sandbox compatible, sentinel 2, widgets, ipyleaflet, interactive