Technische Artikel

Beschleunigung einer geospatialen anwendung mit MATLAB-Werkzeugen für verteilte Berechnungen

Von Narfi Stefansson, MathWorks, Kelly Luetkemeyer, MathWorks, und Rob Comer, MathWorks


dct_map_w.jpg
Das Titelbild zeigt verschiedene Aspekte dreier übereinander gelegter Karten eines Abschnitts der Golfküste von Florida, die auf dem U.S. National Land Cover Dataset basieren. Die Farbcodierung und die zu Grunde liegenden Berechnungen warden in dem Artikel „Beschleunigung einer geospatialen Anwendung mit MATLABWerkzeugen für verteilte Berechnungen“ vorgestellt. Die Titelillustration wurde komplett in MATLAB erstellt. Zum Vergrößern auf das Bild klicken.

Seit seinem Ursprung als einfaches “Matrix Laboratory” hat MATLAB Wissenschaftlern und Ingenieuren geholfen, große Datensätze effizient zu verarbeiten. Die Distributed Computing Toolbox und die MATLAB Distributed Computing Engine erweitern MATLAB nun um Werkzeuge, mit denen sich rechenintensive Probleme lösen lassen, die die Möglichkeiten eines einzelnen Rechners deutlich übersteigen. Dieser Artikel beschreibt die Erstellung einer Landbedeckungskarte mit MATLAB-Werkzeugen für verteilte Berechnungen.

Die Methode wurde gemeinsam vom MathWorks-Team für geospatiale Berechnungen und Dr. Molly E. Brown von Science Systems and Applications, Inc, entwickelt. Anlass war die Absicht, die Landbedeckungskarten sämtlicher US-Staaten mit ihrer Einteilung in Feuchtgebiete, Wälder, städtische Bereiche etc. so umzuformulieren, dass der Vergleich mit niedriger aufgelösten U.S.- oder globalen Fernaufklärungsdaten möglich wird.

Diese Problemstellung ist ein hervorragendes Beispiel für eine grobkörnig verteilte Anwendung. Dieser Problemtyp ist auch als „embarrassingly parallel application“ bekannt und wird oft mit Monte Carlo-Methoden oder Parameter-Studien angegangen. Die gleiche Anwendung läuft in unabhängigen Instanzen auf mehreren Rechnerknoten (Nodes), die keine Kommunikation und keinen Datenaustausch unter-einander haben. Ausschlaggebend für die Gesamtdauer der Anwendung von ihrem Start bis zum Ende ist die Laufzeit auf den einzelnen Nodes.

Aggregation der Landbedeckung und Gewinnung der Mosaikkarte

Im ersten Schritt wird das National Land Cover Dataset (NLCD) von der U.S. Geological Survey (USGS) heruntergeladen. Diese Datenbank enthält Klassifikationsgitter für sämtliche US-Staaten, die aus Planquadraten (Zellen) mit 30m Kantenlänge bestehen. Definitionsgemäß wird einer Zelle ausschließlich der Landbedeckungstyp zugewiesen, die in ihr dominiert — sie enthält also z.B. entweder 0% Mischwald oder 100% Mischwald (Abb. 1). Insgesamt gibt es etwa 20 verschiedene Landbe-deckungsklassen.

Die hohe Dichte und der Umfang der NLCD-Daten plus die Tatsache, dass jedes Klassifikationsgitter im individuellen Koordinatensystem seiner jeweiligen Karte abgelegt ist, macht den Vergleich mit vielen anderen nationalen oder globalen Datensätzen so gut wie unmöglich. Wir selektieren darum eine bestimmte Klasse oder eine Gruppe verwandter Klassen (in Planqua-draten zusammenhängende Zellen ähnlicher Klassen innerhalb des Gitters eines Einzelstaates) und konstruieren daraus ein Mosaik, das das gesamte Staatsgebiet der USA in einem neuen Koordinatensystem abdeckt und eine gröbere räumliche Auflösung aufweist.

dct_fig1_w.gif
Abb. 1. Farbcodierte NLCD-Karte der Landbedeckung Floridas (in verringerter Auflösung) und zugehöriger Klassifikationsschlüssel. Karte und Schlüssel wurden, ebenso wie sämtliche Kartenansichten in diesem Artikel, komplett in MATLAB erstellt. Zum Vergrößern auf das Bild klicken.

Dieses Verfahren opfert räumliche Auflösung zugunsten einer besseren Datenauflösung: Anstatt nur hinsichtlich zwei Stufen zu quantisieren (0% oder 100%), enthält das neue Mosaik den Prozentanteil jedes einzelnen Landbedeckungstyps als kontinuierlichen Wert für jede Zelle.Da jedes Klassifikationsgitter sein eigenes Koordinatensystem besitzt, müssen wir für jedes Gitter eine eigene Zusammenfassung durchführen und danach die Ergebnisse auf das Koordinatensystem projizieren, das wir für die endgültige Mosaikkarte der USA verwenden (Abb. 2).

dct_fig2_w.jpg
Abb. 2. Mosaik der Dichte erschlossener Gebiete in den US-Staaten; die Farbcodierung entspricht dem Prozentzsatz erschlossener Gebiete in Planquadraten mit 8km Seitenlänge. Zum Vergrößern auf das Bild klicken.

Zusammenfassung einer Landbedeckungsklasse

MATLABs Funktionen und Fähigkeiten zur Indexierung von Arrays ermöglichen uns die gewünschte Zusammenfassung in einem einzigen, effizienten Verarbeitungsschritt. In den Beispielbildern werden zwei Arten von Feuchtgebieten für ein Areal an der Südwestküste Floridas in 900m-Qua-draten zusammengefasst. Alle drei Bilder sind noch im NLCD-Koordinatensystem für Florida dargestellt.

dct_fig3_w.gif
Zweifarbenbild, in dem Feuchtgebiete weiß dargestellt sind; alle anderen Landschaftsarten sind dunkelblau.
dct_fig4_w.gif
Das zusammengefasste Bild zeigt gleitende Abstufungen des Feuchtgebiets-anteils mit erheblich gröberer räumlicher Auflösung.
dct_fig5_w.gif
Die ursprüngliche NLCD-Klassifikation; der Klassifikationsschlüssel ist identisch mit dem aus Abb. 1. Das Vorherrschen heller Blautöne in der Nordhälfte weist auf eine Konzentration waldiger Feuchtgebiete hin.

Implementierung der Methode — Sequenzielle Fassung

Diese Anwendung bietet sich für eine verteilte Berechnung geradezu an, da der Großteil der Verarbeitung Staat für Staat durchgeführt wird und somit sehr gut für eine grobe Parallelisierung geeignet ist. Weil Anwendungen für verteilte Berechnungen häufig erheblich schwieriger zu realisieren sind als sequenzielle Anwendungen, entwickeln wir zuerst den Programmcode für eine sequenzielle Berechnung. Dabei können wir bereits sämtliche Stärken der interaktiven MATLAB-Umgebung nutzen1. Sobald dieser sequenzielle Code stabil läuft, verteilen wir die Berechnungen. Die Funktionen der sequenziellen Fassung folgen dem natürlich erscheinenden Berechnungsablauf:

  1. Einrichten der korrekten Referenzprojektion und Initialisierung des Ausgabegitters.
  2. Abarbeiten aller Staaten mit folgenden, wiederkehrenden Aufgaben:

    1. Importieren des NLCD-Kartengitters für den aktuellen Staat.
    2. Erzeugen einer zusammengefassten Landbedeckungskarte für diesen Staat.
    3. Reprojektion der Ergebnisse aus diesen Zusammenfassungen auf das Koordinatengitter für die USA.
  3. Zusammenstellung der einzelnen reprojizierten Landbedeckungsbilder aller Staaten in einem einzelnen Mosaikbild.

Diese Struktur hat eine Reihe von Vorteilen: Sie reduziert die Komplexität der Aufgabenverteilung, begrenzt die Ein- und Ausgabedaten jeder Funktion und sorgt sowohl bei der sequenziellen als auch bei der verteilten Version für ein Minimum an Redundanz im Programmcode.

Reprojektion aller Einzelstaaten

Nachdem die Zusammenfassung der Landbedeckungsklassen für einen Einzelstaat berechnet worden ist, muss das resultierende Gitter auf die Koordinaten der Karte für die USA reprojiziert (geometrisch transformiert) und in das Mosaik kopiert werden.

Um eine geometrische Transformation auf ein Gitter oder Bild anzuwenden, wird normalerweise der Mittelpunkt jeder Ergebniszelle durch die inverse Transformation auf den Raum der Ausgangsdaten abgebildet und dann der neue Wert interpoliert. In unserem Beispiel wird diese inverse Transformation in zwei Schritten erreicht: Wir ermitteln die Positionen der Zellmittelpunkte aus den US-Kartenkoordinaten in Längen- und Breitenangaben und reprojizieren sie anschließend mit Hilfe der aus der Mapping Toolbox stammenden Funktionen „projinv“ und „projfwd“ in das Koordinatensystem der Einzelstaaten.

[lat, lon] = projinv(regionProj, usX, usY);
[stateX, stateY] = projfwd(stateProj, lat, lon); 

Schließlich skalieren und verschieben wir stateX und stateY in die ihnen im jeweiligen Gitter der Staaten zugewiesene Zeile und Spalte und geben diese Werte an die Funktion „tformarray“ aus der Image Processing Toolbox weiter. Diese interpoliert die Daten passend zur Mosaikkarte der USA.

dct_fig6_w.gif
Der Umriss von Massachusetts illustriert die Wirkung der Reprojektion durch Anwendung kleiner Drehungen, Maßstabsveränderungen und Ursprungsverschiebungen.

Übergang zum Verteilten Modell

Bei unseren sequenziellen Berechnungen werden sämtliche Staaten nacheinander auf die immer gleiche Weise abgearbeitet. Dazu nutzen wir entweder einen einzigen, großen Task oder aber bis zu 48 kleinere Tasks, einer für jeden Staat. Um flexibel zu bleiben, kann man bei der verteilten Version einstellen, in wie viele Tasks die Berechnung unterteilt werden soll. Die Zeitmessungen basieren einheitlich auf der Aufteilung in 48 Tasks.

Nach Unterteilung der Staaten in Gruppen, eine für jeden Task, erzeugen wir mit der Distributed Computing Toolbox den Berech-nungsjob und seine Tasks und schicken diesen an den Computer-Cluster, auf dem die Landbedeckung berechnet werden soll. Sobald ein Worker (Einzelrechner) seine Aufgabe erledigt hat, wird das Ergebnis vom Job-Manager an den MATLAB-Client übergeben, der es in das Mosaikbild einfügt.

Folgende Überlegungen waren für die gewählte Aufteilung ausschlaggebend:

  • Die mittlere Laufzeit für jeden Staat beträgt etwa 90 Sekunden. Da der Task-Overhead im Vergleich zur Laufzeit mini-mal ist, bleibt die Anwendung auch bei 48 Tasks immer noch grobkörnig.
  • Keiner der Rechner ist ein Fileserver. Um die Leistung des Clusters zu optimieren, kopieren wir daher die gesamten 30 GB an Landbedeckungsdaten auf jeden einzelnen Worker. Der Task selber muss dadurch insgesamt nur einige Dutzend Kilobytes an Daten übertragen.
  • findet keine Lastverteilung statt; die Staaten werden alphabetisch auf die Tasks verteilt, nicht nach Größe.

Speicher-Mapping

Auch für Einzelstaaten sind die NLCD-Datensätze noch sehr groß.MATLABs Funktionen zum Speicher-Mapping machen es möglich, ein Datengitter für einen Staat wie ein normales MATLAB-Array zu behandeln, ohne es ständig komplett im Speicher vorzuhalten. Dazu konstruieren wir ein memmapfile-Objekt:


m = memmapfile(filename, ‘Offset’, 0, ‘Format’,...
{‘uint8’, [arrayWidth, arrayHeight], ‘nlcdArray’});


und geben das m.Data.nlcdArrayan unsere Zusammenfassungs-Funktion weiter, als wäre es ein gewöhnliches MATLAB-Array im uint8 Format:

[aggregatedData, mask] = aggregate(m.Data.nlcdArray, ...
nlcdCategory, scaleFactor);

Hardwareplattform und Ergebnisse

Die fertige Anwendung läuft auf einem Cluster aus vier 2,8 GHz Pentium IV-Rechnern, die durch ein Gigabit- Ethernet verbunden sind, je 512 MB RAM aufweisen und unter Windows XP laufen. Die MATLAB Distributed Computing Engine ist auf jedem Rechner installiert und startet jeweils einen Worker. Auf einem der Rechner läuft außerdem der Job-Manager. Der Clientrechner ist ein Laptop mit nur 512 MB RAM.

Die Laufzeit der Landbedeckungs-Berechnungen lässt sich erheblich verbessern, wenn alle vier Rechner eingebunden sind anstatt nur einer: Ein einzelner Task auf nur einem Worker benötigt 4.651 Sekunden. 48 Tasks auf allen vier Workern dagegen brauchen nur 1.280 Sekunden — angesichts des minimalen Aufwands an zusätzlicher Entwicklungszeit ist das eine bemerkenswerte Verbesserung.

dct_fig7_w.jpg
Abb. 3. Darstellung der Interaktion zwischen dem Clientrechner, auf dem die Distributed Computing Toolbox die Jobs und Tasks definiert, und der MATLAB Distributed Engine. Zum Vergrößern auf das Bild klicken.

Glossar zum Verteilten Rechnen

  • Client—, in den die MATLAB-Befehle eingegeben werden
  • Job—Große MATLAB- oder Simulink-Operation, die auf einem Computer-Cluster ausgeführt wird
  • Tasks—Teile eines Jobs: die Arbeitseinheiten, die den Workern zur Ausführung übergeben werden
  • Job manager—Prozess, der die Tasks koordiniert und asynchron auf die Worker verteilt
  • Worker—Nicht-interaktive MATLAB-Instanz, die vom Job-Manager zugewiesene Tasks verarbeitet und deren Ergebnisse an den Job-Manager zurückgibt

Dynamische Lizenzverwaltung

Die MATLAB Distributed Computing Engine führt automatisch alle Algorithmen aus, die aus Toolboxen oder Blocksets stammen, für die der Client Lizenzen besitzt. Für die einzelnen Worker müssen also keine zusätzlichen Lizenzen erworben werden.

Veröffentlicht 2006

Artikel für ähnliche Einsatzgebiete anzeigen