Displaying satellite imagery on a web map 9ea3729975aa4f95b05cd489ae7c6918

  • Produits utilisés : s2_l2a

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

Aperçu

Leaflet est la bibliothèque JavaScript open source leader pour les cartes interactives adaptées aux mobiles. Les fonctionnalités qu’elle fournit sont exposées aux utilisateurs de Python par ipyleaflet. Cette bibliothèque permet de créer des cartes interactives dans l’environnement Jupyter notebook/JupyterLab.

Description

Ce bloc-notes montre comment tracer une image et des empreintes de jeu de données sur une carte.

  1. Charger des paquets

  2. Trouver un emplacement

  3. Trouver des ensembles de données à charger

  4. Charger les données de pixels dans la projection « EPSG:3857 », identique à celle utilisée par la plupart des cartes Web

  5. Créer des empreintes de jeu de données à afficher sur une carte

  6. Créer un contrôle d’opacité à afficher sur la même carte

  7. Afficher l’image chargée à partir des jeux de données sur la même carte

Commencer

Pour exécuter cette analyse, exécutez toutes les cellules du bloc-notes, en commençant par la cellule « Charger les packages ».

Charger des paquets

[1]:

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

Se connecter au datacube

[2]:

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

Trouver un emplacement

Pour définir la zone d’intérêt, deux méthodes sont disponibles :

  1. En spécifiant la latitude, la longitude et la zone tampon. Cette méthode nécessite que vous saisissiez la latitude centrale, la longitude centrale et la valeur de la zone tampon en degrés carrés autour du point central que vous souhaitez analyser. Par exemple, « lat = 10,338 », « lon = -1,055 » et « buffer = 0,1 » sélectionneront une zone avec un rayon de 0,1 degré carré autour du point avec les coordonnées (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 8f77c70a8cd04c7b824bafa16cc0e83a 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.

Pour utiliser l’une de ces méthodes, vous pouvez décommenter la ligne de code concernée et commenter l’autre. Pour commenter une ligne, ajoutez le symbole "#" avant le code que vous souhaitez commenter. Par défaut, la première option qui définit l’emplacement à l’aide de la latitude, de la longitude et du tampon est utilisée.

La latitude et la longitude sélectionnées seront affichées sous forme de cadre rouge sur la carte sous la cellule suivante. Cette carte peut être utilisée pour trouver les coordonnées d’autres lieux. Faites simplement défiler la carte et cliquez sur n’importe quel point pour afficher la latitude et la longitude de cet endroit.

[3]:

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

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

Trouver des ensembles de données

Utilisez l’explorateur Digital Earth Africa <https://explorer.digitalearth.africa/products>`__ ou dc.list_products() pour trouver les jeux de données disponibles. Pour plus d’informations sur l’utilisation de dc.list_products(), consultez le carnet de produits et mesures <../Beginners_guide/02_Products_and_measurements.ipynb>`__.

Dans cet exemple, nous utilisons le produit Sentinel-2A ARD. Nous allons visualiser une partie de la bande capturée par Sentinel-2A le 12 janvier 2018.

[4]:

# 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

Créer une carte de prospectus avec des empreintes de jeu de données

Nous souhaitons afficher les empreintes des jeux de données ainsi que les images capturées. Par conséquent, nous utilisons dss = dc.find_datasets(..) pour obtenir une liste d’objets datacube.Dataset se chevauchant avec notre requête en premier.

Nous convertissons ensuite la liste des objets du jeu de données en un GeoJSON d’empreintes de jeu de données, tout en calculant également le cadre de délimitation. Nous utiliserons le cadre de délimitation pour définir la fenêtre d’affichage initiale de la carte.

[5]:

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

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

Créez « ipyleaflet.Map » avec des contrôles de visibilité plein écran et par couche. Définissez la vue initiale pour qu’elle soit centrée sur les empreintes du jeu de données. Nous n’afficherons pas encore la carte.

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

Nous ajoutons maintenant des empreintes de pas à la carte.

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

Créer un calque d’image de dépliant

Sous le capot, « mk_image_layer » va :

  1. Convertir un xarray « rgb » 16 bits en une image RGBA 8 bits

  2. Encoder l’image RGBA en tant que données PNG odc.ui.to_rgba

  3. Rendre les données PNG en « data uri »

  4. Calculer les limites de l’image

  5. Construisez « ipyleaflet.ImageLayer » avec l’URI de l’étape 3 et les limites de l’étape 4

La compression JPEG peut également être utilisée (par exemple fmt="jpeg"), utile pour les images plus grandes afin de réduire la taille du bloc-notes en octets (utilisez quality=40 pour réduire davantage la taille), l’inconvénient est l’absence de prise en charge de l’opacité contrairement au PNG.

Les images satellites sont souvent de 12 bits et plus, mais les images Web sont généralement de 8 bits, nous devons donc réduire la profondeur de bits des images d’entrée de sorte qu’il n’y ait que 256 niveaux par canal de couleur. C’est là qu’intervient le paramètre « clamp ». Dans ce cas, nous utilisons « clamp=2000 ». Les valeurs d’entrée de « 2000 » et plus seront mappées à la valeur « 255 » (la plus grande valeur non signée possible sur 8 bits), « 0 » sera mappé à « 0 » et toutes les autres valeurs intermédiaires seront mises à l’échelle de manière linéaire.

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

Ajouter un contrôle d’opacité

  • Ajouter un curseur vertical à la carte

  • Faire glisser le curseur vers le bas diminue l’opacité du calque d’image

  • L’utilisation de « jslink » depuis « ipywidgets » garantit que ce comportement interactif fonctionnera même sur un notebook pré-rendu (c’est-à-dire sur 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))

Enfin afficher la carte

[10]:

display(m)

Partage de carnets en ligne

Contrairement aux notebooks avec des figures « matplotlib », l’enregistrement d’un notebook après son exécution ne suffit pas à afficher des cartes interactives lors du partage de notebooks rendus en ligne. Vous devez également vous assurer que « Widget State » est enregistré. Dans JupyterLab, assurez-vous que le paramètre « Save Widget State Automatically » est activé. Vous pouvez le trouver dans le menu « Settings ».

Informations Complémentaires

Licence : Le code de ce carnet est sous licence Apache, version 2.0 <https://www.apache.org/licenses/LICENSE-2.0>. Les données de Digital Earth Africa sont sous licence Creative Commons par attribution 4.0 <https://creativecommons.org/licenses/by/4.0/>.

Contact : Si vous avez besoin d’aide, veuillez poster une question sur le canal Slack Open Data Cube <http://slack.opendatacube.org/>`__ ou sur le GIS Stack Exchange en utilisant la balise open-data-cube (vous pouvez consulter les questions posées précédemment ici). Si vous souhaitez signaler un problème avec ce bloc-notes, vous pouvez en déposer un sur Github.

Version de Datacube compatible :

[11]:

print(datacube.__version__)

1.8.15

Dernier test :

[12]:

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

[12]:

'2023-08-14'