Displaying satellite imagery on a web map 5f41a5bb9b834bf587c36d8017f43ad4

  • Sign up to the DEA Sandbox to run this notebook interactively from a browser

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

  • Products used: ga_s2bm_ard_3

Background

Leaflet is the leading open-source JavaScript library for mobile-friendly interactive maps. Functionality it provides is exposed to Python users by ipyleaflet and folium. 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 some pixel data in EPSG:3857 projection, same as used by most web maps

  2. Display image loaded from these datasets on a map

  3. Display dataset footprints on a map along with 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

import sys
sys.path.insert(1, '../Tools/')
from dea_tools.maps import folium_map, folium_dual_map
from dea_tools.maps import ipyleaflet_map

Connect to the datacube

[2]:
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 the captured imagery, but later on we will also need the dataset footprints. 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.

[3]:
# Define product and red/green/blue bands in the given product
product = 'ga_s2bm_ard_3'
RGB = ('nbart_red', 'nbart_green', 'nbart_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.

[4]:
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
[4]:
<xarray.Dataset>
Dimensions:      (time: 1, y: 2436, x: 989)
Coordinates:
  * time         (time) datetime64[ns] 2018-01-13T00:57:00.462000
  * 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
    spatial_ref  int32 3857
Data variables:
    nbart_red    (time, y, x) int16 -999 -999 -999 -999 ... -999 -999 -999 -999
    nbart_green  (time, y, x) int16 -999 -999 -999 -999 ... -999 -999 -999 -999
    nbart_blue   (time, y, x) int16 -999 -999 -999 -999 ... -999 -999 -999 -999
Attributes:
    crs:           EPSG:3857
    grid_mapping:  spatial_ref

Place data on a map

We can display this data on a folium map. We will use Datacube OWS styles to render the datasets into an image.

[5]:
# datacube OWS style configuration (see: https://datacube-ows.readthedocs.io/en/latest/styling_howto.html)
rgb_ows_cfg = {
    "components": {
        "red": {"nbart_red": 1.0},
        "green": {"nbart_green": 1.0},
        "blue": {"nbart_blue": 1.0},
    },
    "scale_range": (50, 3000),
}

folium_map(rgb, ows_style_config=rgb_ows_cfg, width=600, height=600)
[5]:
Make this Notebook Trusted to load map: File -> Trust Notebook