In dieser 4-teiligen Serie zeigen wir Ihnen die Vorzüge von ArcGIS Notebooks in ArcGIS Pro. Der letzte Teil zeigt, wie man Satellitenbild-Daten aus dem Copernicus-Programm der EU direkt in ArcGIS hineinholen und nutzen kann.
Diese Blogserie beschreibt ArcGIS Notebooks in ArcGIS Pro und zeigt anhand einiger Beispiele die Verwendungsmöglichkeiten. Wie im ersten Beitrag beschrieben, sind Notebooks in der gesamten Plattform nutzbar und so vielfältig wie ein Schweizer Taschenmesser. Das nachfolgende Beispiel führt Schritt für Schritt durch einen Anwendungsfall aus dem Bereich der Fernerkundung.
Esri stellt einige Fernerkundungsdaten als Image Services bereit – das sind Dienste, die Nutzer für eigene Analysen verwenden können. Ein solcher Dienst umfasst alle Sentinel-2-Daten aus den jeweils letzten 14 Monaten. Da es sich also um ein „rollling archive“ handelt, werden die Bilddaten in den gezeigten Bildschirmausschnitten bei einem neuem Start anders aussehen.
Die Sentinel-2-Satelliten sind Erdbeobachtungssatelliten des Copernicus Progamms der Europäischen Weltraumorganisation. Die Satelitten liefern Bilddaten aus dem sichtbaren und infratotem Spektrum zwischen 443 und 2190nm. Dieser Spektralbereich wird in 13 Kanälen abgebildet , bei einer Abtastbreite von rund 290 km und einer maximalen Auflösung von 10 m. Die Daten sind ideal um etwa Veränderungen in der Vegetation zu erkennen, Erntevorhersagen zu erstellen, Waldbestände zu kartieren oder das Wachstum von Wild- und Nutzpflanzen zu bestimmen. Esri liefert die Daten des letzten Jahres über einen Image Service aus und stellt dabei bereits erste Verarbeitungsvorlagen zur direkten Anwendung zur Verfügung.
Es gibt verschiedene Möglichkeiten wie diese Daten in ArcGIS Pro genutzt werden können. Eine davon führt über die ArcGIS Notebooks:
1. Neues Notebook öffnen
Über den Abschnitt „Geoverarbeitung“ lässt sich ein neues Python-Notebook erstellen.
2. Importieren der benötigten Module
Die ArcGIS API for Python kann genutzt werden, um eine Verbindung zu einem Portal aufzubauen und einen Image Service ansprechen zu können. Das entsprechende Modul heißt arcgis.
Import arcgis
Das sieht im Notebook so aus:
Die Zeile wird über „Run“ oder die Tastenkombination Shift+Enter abgearbeitet.
Hinweis: So lange in der eckigen Klammer ein Sternchen erscheint, ist der Prozess noch nicht abgeschlossen.
Im Verlauf werden weitere Module benötigt, für die Ansprache des Portals (GIS), die Geocodierung (geocode), die Kartengestaltung (ipywidgets) und die Datenausgabe (Pandas):
import pandas as pd
from arcgis.gis import GIS
from arcgis.geocoding import geocode
from ipywidgets import *
import pandas as pd
3. Verbindung zu ArcGIS Online
Wenn ArcGIS Pro bereits eine Verbindung mit ArcGIS Online hat (das lässt sich im rechten oberen Bereit des Software erkennen), kann auf diese Verbindungsdefinition direkt zurückgegriffen werden, indem folgende Zeilen in die erste Zelle des neuen Notebooks eingegeben werden:
gis= GIS(“pro”)
Alternativ kann natürlich auch ganz dezidiert die Verbindung definiert werden (in der Art gis = GIS(“https://pythonapi.playground.esri.com/portal”, “arcgis_python”, “amazing_arcgis_123”).
4. Suchen des Sentinel-2-Dienstes, Abfragen der Eigenschaften
Mit einer einzigen Zeile wird der Sentinel-2-Dienst angesprochen:
s2 = gis.content.search(‘Sentinel_View’, ‘Imagery Layer’, outside_org=True)[0]
s2
Die Eigenschaften dieses Dienstes sind genau beschrieben und können entsprechend ausgegeben werden, zum Beispiel in einer HTML-Variante:
from IPython.display import HTML
HTML(s2.description)
Um den Dienst in eine Karte zu bringen, wird der Dienst als Layer definiert und einer Karte hinzugefügt. Dabei lässt sich der Kartenausschnitt auch definieren, beispielsweise über einen Ortsnamen (hier wird München > Munich verwendet, es funktioniert aber auch mit anderen Varianten):
s2_layer = s2.layers[0]
map_munich = gis.map(‘Munich’)
map_munich.add_layer(s2_layer)
map_munichsentinel
5. Räumlicher Filter
Dieser Dienst umfasst die Aufnahmen der letzten 14 Monate weltweit. Diese Zeilen zeigen, wie auf einen bestimmten Extent gefiltert werden kann. Wechseln wir dazu den Interessensbereich und gehen von München nach Berlin:
aoi = geocode(‘Berlin’, out_sr=s2_layer.properties.spatialReference)[0][‘extent’]
s2_layer.extent = aoi
s2_layer
6. Wolkenfreie Bilder finden
Gezeigt wird übrigens hier immer das sog. „beste“ Bild, also das Bild aus der jüngsten Vergangenheit mit geringster Wolkenbedeckung.
Jetzt gilt es alle vorhandenen Bilder mit einem maximalen Bewölkungsgrad von 2% herauszufiltern:
s2_filtered = s2_layer.filter_by(where = ‘(Category = 1) AND (CloudCover <= 0.02)’, geometry=arcgis.geometry.filters.within(aoi))
fs = s2_filtered.query(out_fields=”acquisitiondate, objectid, name, CloudCover, Category”, return_geometry=True, return_distinct_values=False, order_by_fields=”acquisitiondate”)
Mit Hilfe von Pandas lassen sich die Ergebnisse schön in einer Tabelle zeigen:
Es wurden also 18 Bilder für den Zeitraum der letzten 14 Monate für den Raum Berlin gefunden, die wolkenfrei bzw. nur sehr gering bewölkt waren Hinweis: Dieses Ergebnis wird bei erneuter Durchführung von dem hier Gezeigten abweichen, da der Sentinel-2-Dienst täglich aktualisiert wird.
Für weitere Analysen können bspw. über die Objekt-ID Inhalte genau angesprochen werden. In diesem Fall soll untersucht werden, wie sich die Bodenfeuchte zwischen Frühjahr und Sommer entwickelt hat.
layer1 = s2_layer.filter_by(‘objectid=10516568’, geometry=arcgis.geometry.filters.within(aoi))
layer2 = s2_layer.filter_by(‘objectid=11641732’, geometry=arcgis.geometry.filters.within(aoi))
Die Bodenfeuchte lässt sich über den Normalized Moisture Index abbilden (das ist der Unterschied zwischen den Messungen im Nahen und fernerem Infrarot):
from arcgis.raster.functions import *
Layer1ndmi= apply(layer1, ‘NDMI Colorized’)
Layer2ndmi= apply(layer2, ‘NDMI Colorized’)
Der Unterschied zwischen beiden Layern lässt sich einfach visualisieren:
diff = stretch(composite_band([Layer1ndmi, Layer2ndmi, Layer1ndmi]), stretch_type=’stddev’, num_stddev=3, min=0, max=255, dra=True, astype=’u8′)
7. Ausgabe der Ergebnisse
Die Ergebnisse sollen in 3 Karten nebeneinander präsentiert werden: Karte 1 enthält den Wasserindex der Frühjahrsszene, Karte 2 den Wasserindex der Sommerszene und Karte 3 die Differenzdarstellung:
m1 = gis.map(‘Berlin’, 11)
m1.add_layer(Layer1ndmi)
m2 = gis.map(‘Berlin’, 11)
m2.add_layer(Layer2ndmi)
m3 = gis.map(‘Berlin’, 11)
m3.add_layer(diff)
m1.layout=Layout(flex=’1 1′, padding=’6px’, height=’300px’)
m2.layout=Layout(flex=’1 1′, padding=’6px’, height=’300px’)
m3.layout=Layout(flex=’1 1′, padding=’6px’, height=’300px’)
box = HBox([m1, m2, m3])
box
An diesem Zeitausschnitt ist übrigens deutlich zu sehen, dass Flächen wie Berlin Tempelhof deutlich trockener geworden sind, siehe violett eingefärbter Bereich im Süden des Ausschnitts.
Nun sind wir am Ende unserer 4-teiligen Blogserie.
Falls Sie einen Teil verpasst haben, finden Sie hier 1-3 zum Nachlesen:
Teil 1: ArcGIS Notebooks in der Plattform
Teil 2: ArcGIS Notebooks in ArcGIS Pro
Teil 3: Datenaufbereitung und Analyse mit Notebooks in ArcGIS Pro
