[PYTHON] Erste Satellitendatenanalyse von Tellus

Einführung

Dieser Artikel wurde als Tag 21 von Cisco Adventskalender 2019 zwischen Cisco Systems GK veröffentlicht. Welche Verbindung besteht überhaupt zwischen Cisco und Satellitendaten? Ich denke, es gibt eine Frage, deshalb möchte ich sie auch einführen.

Was können Sie mit diesem Artikel tun?

――Verstehen, was Satellitendaten sind

• Verstehen, was die Satellitendatenplattform Tellus ist
• Verstehen, wie die Tellus-Entwicklungsumgebung verwendet wird
• Satellitendaten können mit der Tellus-Entwicklungsumgebung analysiert werden
• Verstehen, was NDVI und NDWI sind

Am Ende können Sie dies tun

• Zeigen Sie hochauflösende Satellitenbilder an, indem Sie Datum und Position auswählen Di, 31. Juli 2018 12:42:49 GMT, nördlicher Breitengrad 35,675 - 35,605: östlicher Längengrad 139,745 - 139,835

--Visualisieren Sie den Vegetationsstatus (Grünfläche) anhand von Satellitenbildern

--Quantifizierung des Vegetationsstatus (NDVI)

Was ist die Satellitendatenplattform Tellus?

Von der offiziellen Seite

Tellus ist Japans erste offene und kostenlose Satellitendatenplattform, die darauf abzielt, mithilfe von Regierungssatellitendaten einen neuen Geschäftsmarkt zu schaffen. Es multipliziert mehrere Daten und bietet alle Funktionen zur Förderung der Gründung neuer Unternehmen. https://www.tellusxdp.com/ja/about

Bisher musste die Satellitendatenverwendung Bilddaten verarbeiten, die Hunderttausende Yen pro Blatt in einem speziellen Datenformat kosten. Es gibt verschiedene Probleme in Bezug auf Preis, Speicher und Rechenressourcen, und trotz der Fülle an Satellitendaten war die Nutzung sehr begrenzt. Um diese Probleme zu lösen, startete das Ministerium für Wirtschaft, Handel und Industrie bis 2020 ein dreijähriges Projekt mit dem Titel "Open & Free Government Satellite Data and Data Utilization Promotion Project", das die Aussetzung seines Bestehens darstellen wird. Es wird Tellus sein.

Charakteristisch

1. Alle Daten sind kostenlos verfügbar (mit Obergrenze)
2. Es können nicht nur Satellitendaten, sondern auch Bodendaten (menschlicher Fluss usw.) verwendet werden.
3. Es können nicht nur Daten, sondern auch Analyseressourcen (CPU / GPU / Speicher) verwendet werden
4. Mit Marketplace können Sie Daten, Anwendungen und Algorithmen handeln
5. Erweiterte Schulungsprogramme zur Unterstützung der Nutzung

Relevanz für Cisco

In Bezug auf die Beziehung zwischen Cisco und dem Weltraumgeschäft führte Cisco 2003 (experimenteller Satellit) und 2009 (kommerzieller Satellit) häufig IP-Router in den Weltraum ein. Es gibt vier Arten von Satelliten: meteorologische Satelliten, Beobachtungssatelliten, Positionierungssatelliten und Kommunikationssatelliten. Von diesen hat Cisco in der Vergangenheit an Kommunikationssatelliten gearbeitet. https://www.cisco.com/c/ja_jp/about/newsroom/archives/us/2010/prod-072010b.html

Bei Tellus beteiligen sich verschiedene Unternehmen, Forschungsinstitute und Organisationen in verschiedenen Bereichen wie die xData Alliance, und Cisco nimmt ab dem 21. Februar 2019 ebenfalls teil. Dieses Mal wollen wir mit Beobachtungssatelliten der oben genannten Typen ein neues Geschäft aufbauen. https://www.xdataalliance.com

Kürzlich ein Service-Design, das eine neue urbane Erfahrung namens [Tokyo Moonshot Challenge] bietet (https://www.cisco.com/c/m/ja_jp/training-events/events-webinars/2019-moonshotchallenge.html) Wir veranstalten auch einen Wettbewerb zur Ideenfindung / -realisierung. Wenn Sie interessiert sind, nehmen Sie bitte an uns teil.

Verwenden wir Tellus OS

Verwenden wir Tellus tatsächlich, um Satellitendaten zu nutzen. Das Verfahren zur Verwendung ist wie folgt.

1. Benutzerregistrierung bei Tellus
2. Melden Sie sich bei Tellus OS an
3. Zugriff auf die Tellus-Entwicklungsumgebung
4. Wählen Sie die zu verwendende Sprache aus --R und Python können verwendet werden, aber dieses Mal werden wir Python verwenden
5. Schreiben Sie Code in Jupyter Notebook
6. Ausführung und Bestätigung der Ergebnisse

Registrieren wir uns als Benutzer bei Tellus

Greifen Sie auf Tellus zu, geben Sie die erforderlichen Informationen ein und registrieren Sie sich als Benutzer. Möglicherweise gibt es einige Orte, an denen Sie eintreten können, aber lassen Sie sich nicht entmutigen.

Verwenden wir die Tellus-Entwicklungsumgebung

Lassen Sie uns nach Abschluss der Benutzerregistrierung Tellus verwenden. Es gibt zwei Arten von Nutzungsumgebungen für Tellus. Melden Sie sich zunächst über die Schaltfläche oben rechts oder unterhalb der Tellus-Webseite bei Tellus OS an.

1. Tellus OS Wir haben eine Umgebung vorbereitet, in der Sie über die GUI intuitiv arbeiten können. Die Daten werden als Voreinstellungen vorbereitet, sodass Sie diese grob für die Analyse verwenden können. Nachdem Sie das Analyseziel hier gegebenenfalls eingegrenzt haben, führen Sie gegebenenfalls eine detaillierte Analyse in der folgenden Entwicklungsumgebung durch.

2. Tellus-Entwicklungsumgebung

APIs, die von Python und R verwendet werden können, sind verfügbar. Es werden nicht nur APIs, sondern auch Umgebungen wie Computing, Storage und GPU bereitgestellt. Auf die Entwicklungsumgebung kann vom Tellus-Betriebssystem aus zugegriffen werden. Sie können auf die Entwicklungsumgebung zugreifen, indem Sie unter Tellus OS auf die Schaltfläche unten klicken.

Für den Zugriff auf die Entwicklungsumgebung ist eine Anwendung erforderlich.

Anwendungsmethode für die Entwicklungsumgebung

1. Gehen Sie zu Meine Seite

2. Wählen Sie die gewünschte Entwicklungsumgebung in der Entwicklungsumgebungsanwendung auf Meine Seite aus und beantragen Sie sie. ――Es kann einige Zeit dauern, da die Anwendungen derzeit konzentriert sind.

API TOKEN Registrierung

Wenn die Entwicklungsumgebung verfügbar wird, registrieren Sie API TOKEN. Sie können sich über Meine Seite registrieren. Geben Sie dem Schlüssel einen beliebigen Namen.

Starten Sie das Jupyter-Notizbuch

Melden Sie sich bei Tellus OS an und rufen Sie die Entwicklungsumgebung auf.

Schreiben wir den Code

Gehen Sie nach dem Anmelden in das Arbeitsverzeichnis. Lassen Sie uns ein neues Python3-Notizbuch erstellen.

Ich werde das Jupyter-Notizbuch hier nicht erklären, aber es ist eine vertraute Umgebung für diejenigen, die es verwendet haben. Die diesmal verwendeten Module sind wie folgt. Importieren Sie diese zunächst.

#Importieren Sie die erforderlichen Module und Bibliotheken
import numpy as np 
from skimage import io, color, img_as_ubyte, filters, transform
import requests
from io import BytesIO
from PIL import Image
import matplotlib.pyplot as plt
%matplotlib inline
import math

Hier möchten wir den diesmal verwendeten Satelliten vorstellen.

Geben Sie das Zugriffstoken an.

#API-Zugriffstoken-API_In TOKEN aufbewahren
API_TOKEN = 'YOUR TOKEN HERE'

Dieses Mal möchte ich ein Bild von der Gegend um die Bucht von Tokio machen, also geben Sie den Breiten- und Längengrad wie folgt ein.

#Mindestspielraum
min_lat = 35.675
#Minimaler Längengrad
min_lon = 139.745
#Maximaler Breitengrad
max_lat = 35.605
#Maximaler Längengrad
max_lon = 139.835

Geben Sie die URL, den Header und die Parameter für die Erfassung von Satellitendaten mithilfe der API an.

#Geben Sie die Server-URL, den Header und die Parameter des zu erfassenden Satellitenbilds an
# ASNARO-URL von 1
urls_ASNARO1 = 'https://gisapi.tellusxdp.com/api/v1/asnaro1/scene'
#Header
headers = {'Authorization': 'Bearer ' + API_TOKEN}
#Geben Sie Parameter an (Maximal- und Minimalwerte für Längen- und Breitengrad).
parameter = { 'min_lat': min_lat, 'min_lon': min_lon, 'max_lat': max_lat, 'max_lon': max_lon}

Die Server-URL und -Parameter finden Sie in der Tellus-API-Referenz. API-Referenz: https://www.tellusxdp.com/ja/dev/api

Es gibt zwei Haupttypen von Satellitendaten-APIs.

1. Welche Art von Szene hat jeder Satellit?
2. Bilddaten für jede Szene, Bilddaten für jede Wellenlänge, SAR-Daten

Lassen Sie uns zunächst überprüfen, welche Art von Szene ASNARO-1 hat.

#Durchsuchen Sie die Szene, die den im Anforderungsmodul angegebenen Bereich beobachtet hat, und rufen Sie die Liste im JSON-Format ab
r_ASNARO1 = requests.get(urls_ASNARO1, params=parameter, headers=headers).json()

#Lassen Sie uns den Inhalt anzeigen
r_ASNARO1[0]
{'acquisitionDate': 'Wed, 31 Oct 2018 12:31:16 GMT',
 'clat': 35.6602740014593,
 'clon': 139.701407725281,
 'cloudCover': 20.0,
 'entityId': '20190115064426419_AS1',
 'max_lat': 35.7191765745438,
 'max_lon': 139.786374532312,
 'min_lat': 35.6013160074949,
 'min_lon': 139.61632039109,
 'path': 0,
 'productId': '20190115064426419_AS1_D01_L1B',
 'row': 0,
 'thumbs_url': 'https://tile.tellusxdp.com/thums/ASNARO-1/20190115064426419_AS1.PNG'}

Lassen Sie uns überprüfen, wie viele Aufnahmedaten sich in diesem Breiten- und Längengradbereich befinden

len(r_ASNARO1)
28

Ich fand, dass es 28 Blätter gibt. Lassen Sie uns den Inhalt anzeigen.

# ASNARO-1 Suchergebnisse anzeigen
print("ASNARO-Liste von 1[Aufnahmedatum und -zeit, entityId, productId, cloudCover]")
for res in r_ASNARO1:
    print(", ".join([res['acquisitionDate'],
        str(res['entityId']), str(res['productId']), str(res['cloudCover'])]))

ASNARO-Liste von 1[Aufnahmedatum und -zeit, entityId, productId, cloudCover]
Wed, 31 Oct 2018 12:31:16 GMT, 20190115064426419_AS1, 20190115064426419_AS1_D01_L1B, 20.0
Sat, 20 Oct 2018 12:52:04 GMT, 20181224064142828_AS1, 20181224064142828_AS1_D01_L1B, 19.0
Tue, 31 Jul 2018 12:42:49 GMT, 20181224064003777_AS1, 20181224064003777_AS1_D01_L1B, 0.0
Sat, 30 Jun 2018 12:17:46 GMT, 20181224062516884_AS1, 20181224062516884_AS1_D01_L1B, 13.0
Fri, 08 Jun 2018 12:25:21 GMT, 20190115062640786_AS1, 20190115062640786_AS1_D01_L1B, 20.0
Sat, 02 Jun 2018 12:36:46 GMT, 20181224062257563_AS1, 20181224062257563_AS1_D01_L1B, 18.0
Sat, 21 Apr 2018 12:38:02 GMT, 20181225033937173_AS1, 20181225033937173_AS1_D01_L1B, 0.0
Sun, 25 Mar 2018 12:27:00 GMT, 20181224062059947_AS1, 20181224062059947_AS1_D01_L1B, 0.0
Wed, 03 Jan 2018 15:30:59 GMT, 20181224061906786_AS1, 20181224061906786_AS1_D01_L1B, 0.0
Tue, 02 Jan 2018 02:16:36 GMT, 20181224061858253_AS1, 20181224061858253_AS1_D01_L1B, 0.0
Fri, 29 Dec 2017 15:24:58 GMT, 20181224061846634_AS1, 20181224061846634_AS1_D01_L1B, 1.0
Sun, 05 Nov 2017 13:55:39 GMT, 20181224061816204_AS1, 20181224061816204_AS1_D01_L1B, 0.0
Sun, 09 Jul 2017 15:40:45 GMT, 20181224061651885_AS1, 20181224061651885_AS1_D01_L1B, 0.0
Mon, 22 May 2017 14:06:36 GMT, 20181224061523296_AS1, 20181224061523296_AS1_D01_L1B, 1.0
Sat, 06 May 2017 15:29:27 GMT, 20181224061506752_AS1, 20181224061506752_AS1_D01_L1B, 8.0
Tue, 04 Apr 2017 15:27:01 GMT, 20181224061413592_AS1, 20181224061413592_AS1_D01_L1B, 1.0
Fri, 24 Mar 2017 13:57:55 GMT, 20181224061303679_AS1, 20181224061303679_AS1_D01_L1B, 0.0
Sat, 18 Mar 2017 15:39:14 GMT, 20181224061254113_AS1, 20181224061254113_AS1_D01_L1B, 2.0
Mon, 24 Oct 2016 14:10:19 GMT, 20181224060959386_AS1, 20181224060959386_AS1_D01_L1B, 0.0
Fri, 17 Jun 2016 15:39:32 GMT, 20181224045149110_AS1, 20181224045149110_AS1_D01_L1B, 29.0
Thu, 05 May 2016 14:08:24 GMT, 20181223152759390_AS1, 20181223152759390_AS1_D01_L1B, 0.0
Fri, 18 Mar 2016 15:31:58 GMT, 20181223152702426_AS1, 20181223152702426_AS1_D01_L1B, 0.0
Mon, 30 Nov 2015 14:31:57 GMT, 20190117124820618_AS1, 20190117124820618_AS1_D01_L1B, 33.0
Tue, 03 Nov 2015 14:32:11 GMT, 20181225123811442_AS1, 20181225123811442_AS1_D01_L1B, 0.0
Sun, 25 Oct 2015 03:35:55 GMT, 20181223152305607_AS1, 20181223152305607_AS1_D01_L1B, 0.0
Thu, 06 Aug 2015 03:12:35 GMT, 20181223060118912_AS1, 20181223060118912_AS1_D01_L1B, 13.0
Wed, 05 Aug 2015 14:11:45 GMT, 20181223060109423_AS1, 20181223060109423_AS1_D01_L1B, 0.0
Sun, 28 Jun 2015 12:41:27 GMT, 20181225122707116_AS1, 20181225122707116_AS1_D01_L1B, 10.0

Hier wird "Wolkendecke" als Wolkenbedeckungsrate bezeichnet und ist ein numerischer Wert dafür, wie viele Wolken im Bild enthalten sind. Diesmal handelt es sich nicht um ein Radar mit zusammengesetzter Apertur, das nicht von SAR-Wolken betroffen ist, sondern um ein normales optisches Satellitenbild. Daher wählen wir die neuesten Daten aus, die keine Wolkendecke aufweisen.

Tue, 31 Jul 2018 12:42:49 GMT, 20181224064003777_AS1, 20181224064003777_AS1_D01_L1B, 0.0

Als Nächstes erhalten wir die Bilddaten für jede Szene mithilfe der zweiten API. ASNARO-1 verwendet productId, um das Bild abzurufen. Speichern wir hier auch die entityId.

#Wählen Sie die Szene aus, die Sie erhalten möchten
#Ausgewählter ASNARO-1 Speichern Sie die Produkt- und Entitäts-ID des Satellitenbilds
ASNARO1_entityId = "20181224064003777_AS1"
ASNARO1_productId = "20181224064003777_AS1_D01_L1B"

Hier erkläre ich die Bildkachel und die Vergrößerung. Es wird ausgedrückt als (z = Zoomstufe, x = Kachel_x, y = Kachel_y), aber in Tellus ist es möglich, 256 x 256 Kachelbilder für jede Zoomstufe zwischen Z = 0-18 mit API zu erhalten. .. Hierbei wird die Web Mercator-Projektion verwendet, bei der es sich um einen Ausschnitt der Weltkarte handelt, die von der Mercator-Projektion im Bereich von 180 ° Ost bis 180 ° West und 85,0511 ° Nord bis 85,0511 ° Süd projiziert wird.

Wenn Z = 1

Wenn Z = 1

Wenn Z = 2

Wenn Z = 7

Wenn Z = 12

Übrigens ist es möglich, den Bereich, den Sie dieses Mal anzeigen möchten, mit Kacheln von Tellus OS wie diesem zu finden. In diesem Fall wählen Sie jedoch "Datenauswahl" → "Kachelbild" wie unten gezeigt.

Wenn es jedoch jedes Mal schwierig ist, so etwas zu tun, werde ich überlegen, wie eine Methode implementiert werden kann, um Breiten- und Längengrade in Kachelkoordinaten in Python umzuwandeln. Ich werde den detaillierten Algorithmus nicht verwenden, da er lang sein wird, aber der folgende Code kann für die Konvertierung verwendet werden.

#Koordinaten der Kachel einschließlich des angegebenen Breiten- und Längengrads (x, y,z) berechnet wird.
#Eine Funktion, die die entsprechenden Kachelkoordinaten mit Breiten- / Längengrad und Zoomstufe z als Eingaben berechnet.
#Die Zoomstufe z, die Sie in der API-Referenz auf der Seite Tellus Developers angeben können, finden Sie.
# https://www.tellusxdp.com/ja/dev/api
def LatLon2Tile(lat,lon, z):
    L = 85.05112878
    tileX = int( (2**(z+7)) * ((lon/180) + 1) )
    tileY = int( (2**(z+7)/math.pi) * (-1 * np.arctanh(math.sin(lat * math.pi/180)) + np.arctanh(math.sin(L * math.pi/180))) )
    return int(tileX/256), int(tileY/256)

#Zentralkoordinaten dieses ausgewiesenen Bereichs(Längengrad Breitengrad)
center_lat = (min_lat + max_lat) / 2
center_lon = (min_lon + max_lon) / 2
#Stellen Sie diesmal die Zoomstufe z auf 12 ein.
z = 12
#Holen Sie sich Kachelkoordinaten
x, y = LatLon2Tile(center_lat,center_lon, z)

Lassen Sie uns mit dieser Funktion tatsächlich ein Satellitenbild aufnehmen. Dieses Mal werden wir Bilder für jedes RGB- und NIR-Band (nahes Infrarot) aufnehmen.

# ASNARO-Holen Sie sich 1 Satellitenbild
#Geben Sie die Server-URL an
url_ASNARO1 =  'https://gisapi.tellusxdp.com/blend/asnaro1'
#Bestimmen Sie die Bandzuordnung
query_ASNARO1 = "r=4&g=3&b=2" 
res_ASNARO1 = requests.get("{}/{}/{}/{}/{}.png?{}".format(
                url_ASNARO1,
                ASNARO1_productId,
                z, x, y,
                query_ASNARO1),
                headers={"Authorization": "Bearer " + API_TOKEN})
RGB_ASNARO1 = np.asarray(Image.open((BytesIO(res_ASNARO1.content))))
#Ändern Sie den Datentyp in eine Ganzzahl
RGB_ASNARO1 = RGB_ASNARO1.astype(np.uint8)

Importieren als R = Band4, G = Band3, B = Band2 durch Bandzuweisung. Das aufgenommene Bild wird als Numpy-Array in 4 Dimensionen mit der Bezeichnung RGBA (A = Alpha Chennel) aufgenommen.

# 7. ASNARO-Holen Sie sich 1 Einzelbandbild(Rot, Blau, Grün, nahes Infrarot)
#Konvertieren Sie in ein zweidimensionales Einzelband-Array aus Rot, Blau und Grün
R_ASNARO1 = RGB_ASNARO1[:,:,0]
G_ASNARO1 = RGB_ASNARO1[:,:,1]
B_ASNARO1 = RGB_ASNARO1[:,:,2]

# ASNARO-Ändern Sie die Bandzuordnung von 1 (rot: Band 4, grün: Band 4, blau: Band 4)
query_ASNARO1 = "r=6&g=6&b=6" 
# ASNARO-1 Holen Sie sich ein Satellitenbild
res_ASNARO1 = requests.get("{}/{}/{}/{}/{}.png?{}".format(
                url_ASNARO1, 
                ASNARO1_productId, 
                z, x, y, 
                query_ASNARO1),
                headers=headers)
#NumPy zweidimensionales Array von Einzelbandbildern im nahen Infrarot
NIR_ASNARO1 = np.asarray(Image.open((BytesIO(res_ASNARO1.content))))
NIR_ASNARO1 = NIR_ASNARO1[:,:,0].astype(np.uint8)

Als nächstes erfassen Sie NIR = Band6 im nahen Infrarot.

Das aufgenommene Bild wird in Echtfarbenkomposition angezeigt, einer Kompositionsmethode in der Nähe des menschlichen Auges. Unten sehen Sie ein Bild der Echtfarbenzusammensetzung.

# 9.Echtfarbenzusammensetzung
#R: Rotes Band, G: Grünes Band, B: Blaues Band
# ASNARO-1
true_ASNARO1 = np.dstack((R_ASNARO1, G_ASNARO1, B_ASNARO1))
# ASNARO-1 Anzeige des Satellitenbildes
plt.figure(figsize=(20,20)) #Anzeigegröße in Zoll 20 x 20 Zoll angegeben
plt.imshow(true_ASNARO1)
plt.show()

Da es mit 256 x 256 Pixel erfasst wird, wird es erheblich unscharf, wenn es in diesem Ausmaß gedehnt wird. Es scheint jedoch, dass einige Informationen wie der Zuflussstatus von Ablagerungen aus dem Hama Rikyu-Garten und Flüssen erfasst werden können.

Ist es mit Google Earth nutzlos?

Bisher habe ich einen ziemlich problematischen Code geschrieben, aber ich denke, es ist möglich, ihn bequemer über GCP in Google Earth zu verwenden, ohne dies zu tun. Ich denke, einige Leute denken das. Ich tat. Es ist wirklich das, was Sie anstreben. Es ist also vollkommen in Ordnung, wenn Sie Ihre Ziele mit Google Earth oder einer anderen Datenquelle erreichen können.

Satellitendaten weisen jedoch eine Reihe von Funktionen auf, die andere Daten nicht haben. Obwohl es sich um eine unabhängige Untersuchung handelt, habe ich die Merkmale der einzelnen unten zusammengefasst.

Wie oben erwähnt, sind die Vorteile von Satellitendaten die Fülle vergangener Daten, Sensoren verschiedener Wellenlängen und die aktive Erfassung durch ein Radar mit synthetischer Apertur, das als SAR bezeichnet wird. Durch die Verwendung dieser wird es möglich, Informationen zu erhalten, die mit sichtbarem Licht nicht erhalten werden konnten.

Darüber hinaus liefert Tellus auch Daten, die aus anderen Quellen als Satelliten stammen, wie z. B. Temperatur- und künstliche Flussdaten, was ebenfalls vorteilhaft ist.

Was sind NDVI und NDWI?

Was sind also die Informationen, die mit sichtbarem Licht nicht erhalten werden konnten? https://sorabatake.jp/5192/

Da der Beobachtungssatellit mit optischen Sensoren verschiedener Wellenlängen ausgestattet ist, unterscheiden sich die Wellenlängen, die erfasst werden können, für jeden Satelliten. Wie oben erwähnt, sind die Wellenlängen der vom Objekt absorbierten und reflektierten elektromagnetischen Wellen unterschiedlich, so dass durch Kombinieren dieser verschiedene Informationen erhalten werden können. Typische Kombinationen sind wie folgt bekannt.

https://sorabatake.jp/5192/

Unter diesen ist NDVI, ein Index dafür, wie viel Pflanzen überwachsen sind, ein häufig verwendeter Index. Es gibt einen Index namens NDWI, der angibt, wie viel Wasser vorhanden ist. Dieses Mal wird das Pixel jedes Bildes durch uint8 $ (2 ^ 8 = 256) $ angezeigt, also durch einen numerischen Wert von 0 bis 255. NDVI und NDWI werden normalisiert zwischen (-1 → 1) angezeigt.

Normalisierter Vegetationsindex NDVI (Kombination von Rot (R) und Nahinfrarot (NIR))

NDVI = \frac{R - NIR}{R + NIR}

Normalisierter Wasserindex NDWI (Kombination aus rotem (R) und kurzwelligem Infrarot (SWIR))

NDWI = \frac{R - SWIR}{R + SWIR}

Als nächstes wird ein Zusammensetzungsverfahren durchgeführt, das als natürliche Farbzusammensetzung bezeichnet wird. Dies ist eine Farbzusammensetzungsmethode, bei der R = NIR (nahes Infrarot) anstelle von R = R verwendet wird. Sie können die Pflanze grün anzeigen. Das spezifische Bild ist wie folgt.

# 10.Natürliche Farbzusammensetzung
#R: Rotes Band, G: Nahes Infrarotband, B: Grünes Band
# ASNARO-1
natural_ASNARO1 = np.dstack((R_ASNARO1, NIR_ASNARO1, G_ASNARO1))

# ASNARO-1 Anzeige des Satellitenbildes
plt.figure(figsize=(20,20)) #Anzeigegröße in Zoll 20 x 20 Zoll angegeben
plt.imshow(natural_ASNARO1)
plt.show()

Der Vegetationsteil wird grün dargestellt, und der Vegetationsbereich kann deutlicher dargestellt werden als bei Betrachtung mit Echtfarbenzusammensetzung.

Als nächstes visualisieren wir den oben beschriebenen NDVI. Geben Sie R und NIR in die Formel ein.

# 1. NDVI
#Ändern Sie den Datentyp von uint8 in float32
NIR_ASNARO1 = NIR_ASNARO1.astype(np.float32)
R_ASNARO1 = R_ASNARO1.astype(np.float32)

#Berechnen Sie den NDVI
NDVI_ASNARO1 = (NIR_ASNARO1 - R_ASNARO1) / (NIR_ASNARO1 + R_ASNARO1)

plt.figure(figsize=(20,20)) #Anzeigegröße in Zoll 6 Zoll x 6 Zoll angegeben
plt.imshow(NDVI_ASNARO1, cmap='Reds')
plt.colorbar()
plt.show()

Ich konnte den NDVI berechnen und den Vegetationsindex quantifizieren. Leider verfügt ASNARO-1 nicht über einen SWIR-Sensor (Short Wave Infrared), sodass der NDWI nicht berechnet werden kann. Selbst in einem solchen Fall ist es jedoch möglich, den NDWI unter Verwendung der Daten von Landsat 8 zu berechnen, das mit SWIR ausgestattet ist. Im Fall von Landsat8 wird jedoch nur die Vergrößerung von Z = 12 unterstützt, sodass die Rauheit dieses Bildes die Grenze darstellt. Auf diese Weise kann Tellus je nach Verwendungszweck Daten von verschiedenen Satelliten wechseln und verwenden.

Aber ist das Bild doch nicht rau?

Ich denke, die Auflösung ist niedriger als das am Anfang des Artikels gezeigte Bild, aber das bedeutet nicht, dass der Webbrowser nicht geladen werden konnte. Jede Zoomstufe erfasst ein Bild mit einer Auflösung von 256 x 256. Wenn Sie also beispielsweise ein Bild mit höherer Auflösung aufnehmen möchten, müssen Sie die Zoomstufe erhöhen. Die Zoomstufe variiert von Satellit zu Satellit, aber selbst innerhalb desselben Satelliten variiert sie von Sensor zu Sensor. In ASNARO-1 beträgt die maximale Zoomstufe 18 für einen optischen Sensor, der RGB kombiniert, aber für einen optischen Sensor für jedes Band wie dieses Mal beträgt die maximale Zoomstufe bis zu 16. Die technischen Details für jeden Satelliten werden auf den folgenden Seiten ausführlich beschrieben. Bitte beziehen Sie sich bei Bedarf auf diese. https://www.tellusxdp.com/ja/dev/data

Wenn ich jetzt versuche, ein hochauflösendes Bild (Z = 16) zu erhalten, wird der Bereich eingeschränkt, da er 256 x 256 Pixel beträgt. Wenn Sie jedoch versuchen, den Bereich zu erweitern (Z = 12), verringert sich die Auflösung aufgrund von 256 x 256 Pixel. In einem solchen Fall ist es möglich, das Bild mit einer Auflösung von Z = 16 im Bereich von Z = 12 zu erfassen, indem die Bilder von Z = 16 verkettet und angezeigt werden.

def get_combined_image(get_image_func, z, topleft_x, topleft_y, size_x=1, size_y=1, option={}, params={}):
    """
Holen Sie sich das kombinierte Kachelbild
    Parameters
    ----------
    get_image_func : function 
Kachelbildaufnahmeverfahren
Argument ist z, x, y, option, params
    z : int 
Kachel-Zoomfaktor
    topleft_x : int 
X-Koordinate der oberen linken Kachel
    topleft_y : int 
Y-Koordinate der oberen linken Kachel
    size_x : int 
Anzahl der Kacheln, die in Längsrichtung verbunden werden können
    size_y : int 
Anzahl der Kacheln, die in Breitengradrichtung verbunden werden können
    option : dict 
API-Pfadoptionen(z,x,y ausgeschlossen)
    params : dict 
Parameter abfragen
    Returns
    -------
    combined_image: ndarray
Kombiniertes Kachelbild
    """
    rows = []
    blank = np.zeros((256, 256, 4), dtype=np.uint8)
    
    for y in range(size_y):
        row = []
        for x in range(size_x):
            try:
                img = get_image_func(z, topleft_x + x, topleft_y + y, option)
            except Exception as e:
                img = blank
                
            row.append(img)
        rows.append(np.hstack(row))
    return  np.vstack(rows)

Hier stellen wir einen Algorithmus vor, der zuerst eine leere Numpy-Matrix erstellt und dann Bilder nacheinander einfügt. Wenn Z um 1 erhöht wird, wird das Bild vervierfacht und vervierfacht. Wenn Sie also von Z = 12 zu Z = 16 wechseln, werden $ 4 ^ {4} = 256 $ Kacheln in einer Reihe angeordnet. Dies ist unten dargestellt.

Der obige Code ist die tatsächliche Implementierung davon.

def get_asnaro1_image(z, x, y, option={}, params={}, query=''):
    """
    ASNARO-Holen Sie sich 1 Kachelbild
    Parameters
    ----------
    z : int 
Kachel-Zoomfaktor
    x : int 
X-Koordinate der Kachel
    y : int 
Y-Koordinate der Kachel
    option : dict 
API-Pfadoptionen(productId)
    params : dict 
Parameter abfragen
    Returns
    -------
    img: ndarray
Fliesenbild
    """
    url = 'https://gisapi.tellusxdp.com/blend/asnaro1/{}/{}/{}/{}.png?{}'.format(option['productId'], z, x, y, query)
    headers = {
        'Authorization': 'Bearer' + API_TOKEN
    }
    r = requests.get(url, headers=headers, params=params)
    if not r.status_code == requests.codes.ok:
        r.raise_for_status()
    return io.imread(BytesIO(r.content))

Dieses Mal habe ich den Algorithmus geschrieben, um das 256x256-Bild zu erfassen, das das Originalbild ist, das wie oben zusammengefügt werden soll. Im Folgenden werden Bilder als R = Band4, G = Band3, B = Band2 für die Echtfarbenzusammensetzung aufgenommen.

# ASNARO-Holen Sie sich 1 Satellitenbild
def rgb_asnaro1(z, x, y, option):
    #Geben Sie die Server-URL an
    url_ASNARO1 =  'https://gisapi.tellusxdp.com/blend/asnaro1'
    #Bestimmen Sie die Bandzuordnung
    query_ASNARO1 = "r=4&g=3&b=2" 
    res_ASNARO1 = requests.get("{}/{}/{}/{}/{}.png?{}".format(
                url_ASNARO1,
                option['productId'],
                z, x, y,
                query_ASNARO1),
                headers={"Authorization": "Bearer " + API_TOKEN})
    RGB_ASNARO1 = np.asarray(Image.open((BytesIO(res_ASNARO1.content))))
    #Ändern Sie den Datentyp in eine Ganzzahl
    #RGB_ASNARO1 = RGB_ASNARO1.astype(np.uint8)
    return RGB_ASNARO1.astype(np.uint8)

Im Folgenden werden Bilder als R = Band6, G = Band6, B = Band6 für die natürliche Farbzusammensetzung aufgenommen.

# ASNARO-Holen Sie sich 1 Satellitenbild
def nir_asnaro1(z, x, y, option):
    #Geben Sie die Server-URL an
    url_ASNARO1 =  'https://gisapi.tellusxdp.com/blend/asnaro1'
    #Bestimmen Sie die Bandzuordnung
    query_ASNARO1 = "r=6&g=6&b=6" 
    res_ASNARO1 = requests.get("{}/{}/{}/{}/{}.png?{}".format(
                url_ASNARO1,
                option['productId'],
                z, x, y,
                query_ASNARO1),
                headers={"Authorization": "Bearer " + API_TOKEN})
    RGB_ASNARO1 = np.asarray(Image.open((BytesIO(res_ASNARO1.content))))
    #Ändern Sie den Datentyp in eine Ganzzahl
    #RGB_ASNARO1 = RGB_ASNARO1.astype(np.uint8)
    return RGB_ASNARO1.astype(np.uint8)

Unten werden die Koordinaten von Top Left X und Top Left Y aus den Mittelkoordinaten berechnet. +1 ist die vorgenommene Korrektur, da die Koordinaten zu viel zurückgeben, wenn 1/2 von den Mittelkoordinaten abgezogen wird.

z = 16
size_x = 16
size_y = 16
x, y = LatLon2Tile(center_lat,center_lon, z)
topleft_x = int(x - size_x / 2 + 1)
topleft_y = int(y - size_y / 2 + 1)
option = {
    'productId': ASNARO1_productId
}
query_ASNARO1 = "r=4&g=3&b=2" 

Geben Sie die obigen Parameter für die Funktion ein, mit der die zusammengesetzten Echtfarbenbilder tatsächlich zusammengefügt werden.

combined = get_combined_image(rgb_asnaro1, z, topleft_x, topleft_y, size_x, size_y, option, params_ASNARO1)

Wir werden auch die obigen Parameter zu der Funktion hinzufügen, die die zusammengesetzten Bilder in natürlichen Farben tatsächlich verbindet.

combined_nir = get_combined_image(nir_asnaro1, z, topleft_x, topleft_y, size_x, size_y, option, query_ASNARO1)

Bilder werden in einem Array für jedes Band gespeichert.

# ASNARO-Holen Sie sich 1 Einzelbandbild(Rot, Blau, Grün, nahes Infrarot)
#Konvertieren Sie in ein zweidimensionales Einzelband-Array aus Rot, Blau und Grün
COMR_ASNARO1 = combined[:,:,0].astype(np.uint8)
COMG_ASNARO1 = combined[:,:,1].astype(np.uint8)
COMB_ASNARO1 = combined[:,:,2].astype(np.uint8)

#NumPy zweidimensionales Array von Einzelbandbildern im nahen Infrarot
COMNIR_ASNARO1 = combined_nir[:,:,0].astype(np.uint8)

Zuallererst Echtfarbenzusammensetzung.

#Echtfarbenzusammensetzung
#R: Rotes Band, G: Grünes Band, B: Blaues Band
# ASNARO-1
true_COMASNARO1 = np.dstack((COMR_ASNARO1, COMG_ASNARO1, COMB_ASNARO1))
# ASNARO-1 Anzeige des Satellitenbildes
plt.figure(figsize=(20,20)) #Anzeigegröße in Zoll 20 x 20 Zoll angegeben
plt.imshow(true_COMASNARO1)
plt.show()

Sie konnten ein ziemlich hochauflösendes Bild erhalten. In Bezug auf Pixel ist es ein Bild von 256 x 16 = 4096 Pixel.

Als nächstes kommt die natürliche Farbzusammensetzung.

#Natürliche Farbzusammensetzung
#R: Rotes Band, G: Nahes Infrarotband, B: Grünes Band
# ASNARO-1
natural_COMASNARO1 = np.dstack((COMR_ASNARO1, COMNIR_ASNARO1, COMG_ASNARO1))

# ASNARO-1 Anzeige des Satellitenbildes
plt.figure(figsize=(20,20)) #Anzeigegröße in Zoll 20 x 20 Zoll angegeben
plt.imshow(natural_COMASNARO1)
plt.show()

Ich konnte hier ein ziemlich hochauflösendes Bild bekommen.

# NDVI
#Ändern Sie den Datentyp von uint8 in float32
COMNIR_ASNARO1 = COMNIR_ASNARO1.astype(np.float32)
COMR_ASNARO1 = COMR_ASNARO1.astype(np.float32)

#Berechnen Sie den NDVI
COMNDVI_ASNARO1 = (COMNIR_ASNARO1 - COMR_ASNARO1) / (COMNIR_ASNARO1 + COMR_ASNARO1)

plt.figure(figsize=(20,20)) #Anzeigegröße in Zoll 6 Zoll x 6 Zoll angegeben
plt.imshow(COMNDVI_ASNARO1, cmap='Reds')
plt.colorbar()
plt.show()

Ein ziemlich hochauflösendes Bild wird auch für NVDI erzeugt.

Bonus

Dieser Artikel beschreibt hauptsächlich die Verwendung von Satellitendaten. Wenn es jedoch um Bilder geht, möchten Sie die Bilderkennung nicht ausprobieren. Ich werde es in diesem Artikel nicht erwähnen, aber natürlich ist eine Bilderkennung auch auf Tellus möglich. Die Ergebnisse der Erkennung eines festgemachten Schiffes mithilfe des maschinellen Lernmodells [^ 3] sind wie folgt.

Wir haben diejenigen mit einer Sicherheit von 0,8 oder höher extrahiert, konnten aber insgesamt 326 verankerte Schiffe finden. Wir haben auch Schiffe entdeckt, die im Fluss in Richtung Land festgemacht haben. Da die Schiffe, die navigieren, vom Erfassungsziel ausgeschlossen sind, ist ersichtlich, dass die winkenden Schiffe nicht erkannt werden.

Reference

• Sorahata
• Umgang mit Neuigkeiten über die Verwendung von Tellus und das Raumfahrtgeschäft
• Tellus
• Klicken Sie hier, um auf das Tellus-Betriebssystem und die Entwicklungsumgebung zuzugreifen

Haftungsausschluss

Die auf dieser Website geäußerten Meinungen und die entsprechenden Kommentare sind die persönlichen Meinungen des Posters und nicht die Meinungen von Cisco. Der Inhalt dieser Website dient nur zu Informationszwecken und darf nicht von Cisco oder einer anderen Partei unterstützt oder zum Ausdruck gebracht werden. Durch das Posten auf dieser Website ist jeder Benutzer allein für den Inhalt aller veröffentlichten, verknüpften oder anderweitig hochgeladenen Informationen verantwortlich und schließt Cisco von jeglicher Haftung für die Nutzung dieser Website aus. Genau.

[^ 1]: Basierend auf zukünftigen Satellitendaten an Bord [^ 2]: Es gibt häufige Methoden wie die Satellitenkonstellation