Displaying satellite imagery on a web map

• Products used: s2_l2a

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

Contexte

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

import os
import ipyleaflet
import numpy as np
import geopandas as gpd
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 datacube.utils.geometry import Geometry

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

Connect to the datacube

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

Find a location

To define the area of interest, there are two methods available:

1. By specifying the latitude, longitude, and buffer. This method requires you to input the central latitude, central longitude, and the buffer value in square degrees around the center point you want to analyze. For example, lat = 10.338, lon = -1.055, and buffer = 0.1 will select an area with a radius of 0.1 square degrees around the point with coordinates (10.338, -1.055).

2. By uploading a polygon as a GeoJSON or Esri Shapefile. If you choose this option, you will need to upload the geojson or ESRI shapefile into the Sandbox using Upload Files button in the top left corner of the Jupyter Notebook interface. ESRI shapefiles must be uploaded with all the related files (.cpg, .dbf, .shp, .shx). Once uploaded, you can use the shapefile or geojson to define the area of interest. Remember to update the code to call the file you have uploaded.

To use one of these methods, you can uncomment the relevant line of code and comment out the other one. To comment out a line, add the "#" symbol before the code you want to comment out. By default, the first option which defines the location using latitude, longitude, and buffer is being used.

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.

# Set the area of interest
# Method 1: Specify the latitude, longitude, and buffer
aoi = define_area(lat=7.656, lon=0.021, buffer=0.2)

# Method 2: Use a polygon as a GeoJSON or Esri Shapefile.
# aoi = define_area(vector_path='aoi.shp')

#Create a geopolygon and geodataframe of the area of interest
geopolygon = Geometry(aoi["features"][0]["geometry"], crs="epsg:4326")
geopolygon_gdf = gpd.GeoDataFrame(geometry=[geopolygon], crs=geopolygon.crs)

# Get the latitude and longitude range of the geopolygon
lat_range = (geopolygon_gdf.total_bounds[1], geopolygon_gdf.total_bounds[3])
lon_range = (geopolygon_gdf.total_bounds[0], geopolygon_gdf.total_bounds[2])

display_map(x=lon_range, y=lat_range, margin=-0.2)

Make this Notebook Trusted to load map: File -> Trust Notebook

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.

# Define products
products = 's2_l2a'

# Specify the parameters to pass to the load query
query = {
    "x": lon_range,
    "y": lat_range,
    "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.

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.

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.

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.

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)

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

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:

print(datacube.__version__)

1.8.15

Last Tested:

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

'2023-08-14'