Journey zur Sichtanalyse in ArcGIS Pro

Im Rahmen der PWRG2 für die Umsetzung des 3-30-300 Konzept von Konijnendijk

In dieser Storymap wird die technische Umsetzung erprobt, wie die Sichtanalyse für das 3-30-300 Konzept von Cecil Konijnendijk für Schweizer Städte umgesetzt werden könnte. Das 3-30-300 Konzept kann einen wichtigen Beitrag für die zukünftige Stadtentwicklung leisten, wie es in der bereits durchgeführten PWRG 1 von Pascal Bieri aufgezeigt wurde. Die Untersuchungen werden zum einen in R stattfinden, wo hauptsächlich das LidR-package von Roussel et al. (2022) untersucht wird. Die daraus entstandenen Datensätze werden anschliessend in der Software CloudCompare rasterisiert und anschliessend wird im ArcGIS Pro die eigentliche Sichtanalyse durchgeführt und beurteilt.

 

Datenmaterial

Für die Untersuchungen wurde das Stadtzürcher Quartier "Alt-Wiedikon" ausgewählt. Das Quartier zeichnet sich durch unterschiedliche Gebäude und Baumstrukturen aus. Es ist wichtig, eine möglichst breite Variabilität an unterschiedlicher Bäumen zu untersuchen, so dass in zukünftigen Anwendungen allfällige Spezialfälle bereits im vornherein detektiert werden können und so ein spezieller Augenmerk auf den bereits im vornherein detektierte Spezialfall gelegt werden kann. Das Gebiet Alt-Wiedikon weisst viele unterschiedlich hohe und unterschiedlich dicht beieinanderstehende Bäume und Baumgruppen auf. Daher ist es ideal geeignet für diese Voruntersuchung.Für die Untersuchungen wurden der LiDAR Punktwolken Datensatz der Stadt Zürich verwendet, welcher in den Jahren 2021 bis 2022 erstellt wurde. Die Punktewolke weist eine mittlere Punktdichte von 16 Pkt/m2, eine Lagegenauigkeit von +/- 20 cm und eine Höhengenauigkeit von +/- 10 cm aus. Es handelt sich dabei um eine hochaufgelöste Punktwolke, welche für diesen Anwendungszweck eine ausreichende Auflösung aufweist.  

Bereits vom Kanton Zürich klassifizierte Klassen der LiDAR Punktwolke

Für diese Arbeit werden die Klassen 2, 5 und 6 verwendet. Die niedrige und mittlere Vegetation wird in dieser Arbeit vernachlässigt, da Bäume kleiner als 3 m einen zu vernachlässigenden Einfluss auf die mentale Gesundheit und nur einen marginalen Einfluss auf das Stadtklima haben. 

Die Vektordaten der Gebäude stammen aus dem schweizweit verfügbaren Datensatz SwissBUILDING3D. Darin enthalten sind alle permanent vorhandenen Gebäuden der Schweiz als 3D Modell mit einer Lagegenauigkeit in alle Dimensionen von +/- 30 cm - 50 cm.

Um einen visuellen Vergleich durchführen zu können wurde ebenfalls ein lokaler Ausschnitt des Orthofotomosaik "SWISSIMAGE" verwendet. Für eine bessere räumliche visuelle Überprüfung wurde Google Streetview beigezogen.

Methode

LidR package (Roussel et al.,2022)

Das LidR-package ist für das Bearbeiten von luftgestützten Laser Scans entwickelt worden. Es ermöglicht, LiDAR Daten zu lesen, zu klassifizieren, zu manipulieren und zu exportieren. Da die LiDAR Daten der Stadt Zürich bereits klassifiziert sind, wurde LidR lediglich dazu verwendet um ein canopy height model (chm) und eine Baumsegmentierung durchzuführen. Im folgenden Abschnitt wird nun das Vorgehen erläutert, wie aus den LiDAR Grunddaten der Stadt Zürich einen LiDAR Baumdatensatz enthält, in dem jeder einzelne Baum detektiert und segmentiert ist und eine eigene Baumidentifikationsnummer aufweist.

Als erstes wird der LiDAR Datensatz in die Arbeitsumgebung von R geladen. Dabei können .las sowie .laz Dateiformate gelesen werden. Da es sich bei Punktwolken um grosse Datenmengen handelt, ist das Gebiet in mehrere 225 ha grossen Kacheln eingeteilt, welche anschliessend separat geladen werden müssen. Durch das Laden der Daten geht zumindest bei den Daten der Stadt Zürich die Zuweisung des Koordinatensystems verloren, dies kann jedoch mittels Funktion «sf::st_crs» wieder zugewiesen werden. In den Metadaten der jeweiligen Punktdaten ist notiert, in welchem Referenzsystem die Daten vorliegen. Anschliessend können die Kacheln mit der Funktion «rbind» zusammengefügt werden. 

Eine Normalisierung kann durchgeführt werden. Dies ist je nach Verwendung sinnvoll, da die Daten damit einfacher zu interpretieren werden. Durch die Normalisierung entstehen absolute Zahlenwerte, welche beispielsweise die Höhe eines Baumes direkt wiedergibt.  Werden die Daten jedoch als 3D-Daten in ArcGIS Scene verwendet, ist eine Normalisierung nicht zielführend, da sonst der 3D-Bezug verloren geht.

Da die Punktdaten bereits klassifiziert sind, kann nun via die Filterfunktion «lidR::filter_poi» die gewünschte Klasse ausgewählt werden. Die für diese Arbeit verwendeten Klassen sind die Klasse 2, 5 und 6 (Bodenpunkte, Vegetation Hoch (> 3 m) und Gebäude). Die drei Klassen ergeben das digitale Oberflächen Modell (DOM), welches später in der Sichtbarkeitsanalyse verwendet wird. Für die Segmentierung der Bäume wird nur mit der Klasse 5 weitergearbeitet.

Canopy Height Model (chm)

Das canopy height model (chm) stellt einen Rasterlayer dar, bei dem jeweils der höcste Punkt der Baumbedeckung den Zellenwert der Rasterzelle ergibt. Es gibt drei Möglichkeiten ein chm zu erhalten, welche sich in ihrer Erstellung jeweils leicht unterscheiden undje nach Datenlage bessere oder schlechtere Ergebnisse liefern. 

Ich habe im Rahmen dieser Arbeit diese Methoden ausprobiert und einander gegenübergestellt. Die p2r-Methode ergibt ein solides Abbild. Da die Baumdichte innerhalb des bebauten Gebietes klein ist, sind die leeren Bereiche nicht erstaunlich. Einzelne Bäume werden bereits gut sichtbar dargestellt. Die TIN-Methode ist ohne ein Max_edge Argument nicht aussagekräftig, da die grossen Leerflächen zu sehr starken Anomalien im Bereich dieser Leerflächen führen. Um dieses Problem zu beheben wurden unterschiedliche Max_edge Längen ausprobiert (1,2,3,4,5). Ist die Länge zu gering, werden zu viele Bereiche im chm nicht wiedergegeben. Sind die Dreieckslängen zu lang (4 und 5), so wird das chm von Anomalien verfälscht und es ergibt eine grosse Unsicherheit.

Ist die Länge jedoch zwischen 2 und 3, so erhält man ähnliche Ergebnisse wie bei der p2r Methode. Wird zum Beispiel das Tin_3 chm mit dem p2r_1 Modell verglichen, so sieht man etwa ähnlich grosse Leerstellen in beiden chm. Beim Tin_3 chm fällt jedoch auf, dass es weniger 1x1 grosse detektierte Flächen hat. Daraus kann geschlussfolgert werden, dass das Tin_3 leicht kleinere Baumflächen marginalisiert, bzw die p2r Methode kleine Baumflächen leicht übervertreten darstellt. 

LidR individuelle Baum-Erkennung (ITD)

Die individuelle Baum-Erkennung wird mithilfe eines lokalen Maximum Filters durchgeführt. Es handelt sich dabei um eine Nachbarschaftsanalyse, welche die umliegenden Höhenwerte aus dem LiDAR-Datensatz miteinander vergleicht und der jeweils Höchste Wert ermittelt wird. 

ttops <- locate_trees(las, lmf(ws = 5))

Im Rahmen dieser Arbeit wurden unterschiedliche Baumkronen-Datensätze generiert. Dabei wurden unterschiedliche fixe Fenstergrössen, sowie unterschiedliche variable Fenstergrössen verwendet. Für eine bessere Visualisierung wurden diese jeweils über ein p2r-chm gelegt. Ebenfalls wurden mithilfe «SWISSIMAGE» und Google Streetview die generierten Baumkronen verglichen. Dadurch konnte effektiv abgeschätzt werden, ob der generierte Baumkronenlayer die Wirklichkeit wahrheitsgetreu abbildet oder nicht. 

Genau für dieses Problem wurde die variable Fenstergrösse entwickelt. Die in den Beispielsunterlagen von LidR (Roussel et al., 2022) erwähnte Funktion für die variable Fenstergrösse sieht folgendermassen aus: 

f <- function(x) {
  y <- 2.6 * (-(exp(-0.08*(x-2)) - 1)) + 3
  y[x < 2] <- 3
  y[x > 20] <- 5
  return(y)
}

Ich habe daraufhin die Schwellenwerte sukzessiv angepasst, um ein besseres Resultat zu erhalten. Zwischen den Schwellenwerten versuchte ich es zuerst ebenfalls mit einer komplexen Exponentialfunktion mit angepassten Schwellenwerten. Das resultierende Ergebnis war ebenfalls nicht zufriedenstellend. Damit ich mit weniger Aufwand die Schwellenwerte anpassen konnte, habe ich anstelle der komplexen Exponentialfunktion eine lineare Funktion zwischen den Schwellenwerten verwendet.

variable Fenstergrösse ttop_mw_58_2095
f <- function(x) {
  y <- 0.095*x +7.55
  y[x < 5] <- 8
  y[x > 20] <- 9.5
  return(y)
}

So konnte diese einfacher angepasst werden. Der Fokus lag dabei darauf, möglichst die reale Baumanzahl zu detektieren, diese jedoch nicht zu überschätzen. Das beste Ergebnis lieferten dabei die Schwellenwerte 20 m Baumhöhe mit 9.5 m Fenstergrösse. Mit diesen Schwellenwerte wurden sieben von acht Bäumen richtig detektiert. Im Vergleich werden mit einer fixen Fenstergrösse von 9 m nur sechs von acht Bäumen detektiert, hingegen wird bei 8.5 m bereits der erste Baum doppelt detektiert wird und somit findet eine Überschätzung statt. Auch wenn nicht eine vollständig korrekte Detektierung aller Baumkronen stattfindet, wurde die Variante MW58_2095 als die vielversprechendste Annäherung an die Wirklichkeit betitelt und für die weiteren Berechnungen verwendet.

links fixe Fenstergrösse mit ws = 9, rechts variable Fenstergrösse (ttop_mw_58_2095), Hintergrund: chm p2r_1

Neben der eben erläuterter Methode, welche direkt mit Punktdaten arbeitet, gibt es noch die Möglichkeit die Baumkronen mithilfe eines chm zu detektieren. Dies ist eine rasterbasierte Anwendung, welche keine LiDAR Punktwolkendaten als Basis benötigt. Durch die reduzierte Datenstruktur ist dies ein einfacheres und somit schnelleres Verfahren. 

links ttop nur mit chm_p2r und ws = 9, rechts ttop mit Punktwolkendaten und ws = 9, Hintergrund: chm_p2r_1

Individuelle Baum Segmentation (ITS)

In einem weiteren Abschnitt werden nun aus den Punktwolkendaten die Bäume jeweils in einzelne Bäume mit einer eigenen Identifikationsnummer segmentiert. Dieser Output ergibt schlussendlich die Grundlage für die noch folgende Sichtanalyse in ArcGis.

ArcGIS

Nachdem die Punktdaten in R vorbereitet wurden, findet die eigentliche Analyse in ArcGIS Pro statt. Dazu wurde in der Version 2.9.5 gearbeitet. 

Da ich bis anhin noch nie mit LiDAR Daten in ArcGIS gearbeitet habe, musste ich mich zuerst mit der neuen Datenstruktur auseinandersetzen. Die aus R generierten .las Daten können ohne Konvertierung geladen werden. Der Workspace wurde mit einem Orthofoto ergänzt, sowie wurde der SwissBUILDING3D Vektordatensatz geladen. Dieser enthält alle Gebäude der Schweiz als detaillierte 3D-Struktur. Das Dateiformat ist Multipatch. Um die Daten in einer 3D-Ansicht darzustellen wird eine lokale Szene erstellt.

lokale Szene mit ITS-Baumdatensatz, SwissBUILDING3D und Orthofoto

Die grossen Datenmengen des Multipatch-Datensatzes sowie der LiDAR Daten machten sich schnell durch lange Ladezeiten bemerkbar. Deswegen wurde das bereits kleine Sample nochmals verkleinert, bzw. auf einen einzelnen Baum mit seinen umliegenden Gebäude reduziert. 

Der erste Gedanke die Sichtanalyse durchzuführen war eine line-of-sight-Analyse. Dazu konnten die Punkte gewählt werden, wohin der Betrachter «blickt». So könnte eine präzise lokale Analyse durchgeführt werden, da der Baum detailliert als Punktwolke vorhanden ist und die Gebäude präzise dargestellt sind.

line-of-sight-Analyse mit reduziertem Datensatz

Doch für eine flächendeckende Analyse ist dies aufgrund der fehlenden Automatismen nicht zielführend. Dabei wurde schnell klar, dass der Detailgrad zum Erreichen der Beantwortung der Frage zu detailliert und somit zu aufwändig ist. Dementsprechend wurde der LiDAR Datensatz rasterisiert, um so eine flächendeckende Analyse zu ermöglichen. Da ArcGIS sehr spezifische Anforderungen an das .las Format hat, scheiterten interne Konvertierungen und es wurde auf die Software CloudCompare ausgewichen. Damit konnten aus den LiDAR Datensätzen Rasterdatensätze erstellt werden. Da ein Rasterdatensatz jeweils nur eine Valuedarstellen kann, wurde zuerst versucht, die Identifikation und die Höheninformation per Codierung zu bewerkstelligen. Dazu wurden bei den ersten 6 Ziffern für die Baumidentifikation und ab der siebten Ziffer die Höheninformation dargestellt. Dies wurde mittels einer Re-klassifizierung umgesetzt und sieht so aus: 2455412325 Damit wurde der Baum mit der Identifikation 412325 mit einer Höhe von 24.55m im Raster dargestellt. Bei solch grossen Zahlen muss beachtet werden, dass der Zahlenraum von 2^31 nicht überschritten werden kann (Integer/32bit).

Ich habe jedoch schlussendlich trotzdem aus dem gleichen LiDAR Datensatz zwei Rasterdatensätze erstellt. Der Eine enthält die Segmentationsinformation der Bäume, resp. die Baum-Identifikationsnummer und der Andere enthält die Höheninformationen der Bäume. 

Sichtanalyse

ArcGIS bietet mehrere Sichtanalysentools für Rasteranwendungen. Diese sind Viewshed, Observerpoints und Visibility. Aufgrund der scheinbar höchsten Flexibilität habe ich zuerst mit dem Tool Visibility versucht meine Fragestellung umzusetzen.

Anschliessend wurde die Analyse für mehrere OFFSETA Werte durchgeführt. Dazu wurde 1 m für das Parterre verwendet, 4 m für den ersten Stock, 7 m für den zweiten Stock.

Im späteren Verlauf der Analysen stosste ich auf das Problem, dass mit dem Parameter «observers» lediglich nur 16 Observer bearbeitet werden können. Somit ist diese Methode für eine flächenhafte Anwendung nicht geeignet. Das gleiche gilt auch für das Tool «Observer Points».

Um das Problem lösen zu können, wurde die Analysemethode umgekehrt. Nun wird untersucht, ob der Baum das Stockwerk sehen kann. Dazu wurde der Workflow «WF_Visibility_tree» entwickelt.

Workflow «WF_Visibility_tree»

Dafür wird zuerst aus den beiden Baumrastern ein Punkfeature generiert, welches jeweils ein Baum mit seiner Maximalhöhe darstellt und als Observer-Feature im Visibility-Tool verwendet wird. Auch wird nun den Parameter «frequencie» wieder verwendet, so dass mehr als 16 Punkte untersucht werden können. Anschliessend wird das Attribut «Gridcode» der Visibility Analyse binär re-klassifizierten ( >= 3 Beobachter = 1, <= 2 Beobachter = 0 ) und anschliessend ein spatial join mit dem Gebäudefussabdruck und dem re-klassifiziertem Visibility-Feature durchgeführt, wobei das Outputfield «Gridcode» mittels Maximum oder Minimum Regel verbunden wird. Als Resultat erhalten wir ein Gebäudefussabdruckfeature mit dem Attribut «Gridcode», welches beschreibt ob das Gebäude die Bedinung von 3 Bäumen erfüllt oder nicht. Da sich im Testausschnitt in Alt-Wiedikon jedes Gebäude die Bedingung erfüllt, wurde noch ein weiterer Test im Quartier Langstrasse um den Gebäudeblock Militärstrasse/Kanonengasse/Lagerstrasse/Tellstrasse durchgeführt, da hier die Gebäude höher sind und es so eine geringere Sichtbarkeit zu erwarten ist.

Fazit

Bereits vom Kanton Zürich klassifizierte Klassen der LiDAR Punktwolke

links fixe Fenstergrösse mit ws = 9, rechts variable Fenstergrösse (ttop_mw_58_2095), Hintergrund: chm p2r_1

links ttop nur mit chm_p2r und ws = 9, rechts ttop mit Punktwolkendaten und ws = 9, Hintergrund: chm_p2r_1

lokale Szene mit ITS-Baumdatensatz, SwissBUILDING3D und Orthofoto

line-of-sight-Analyse mit reduziertem Datensatz

Workflow «WF_Visibility_tree»