User Tools

Site Tools


projekte:misc:perfectpeak

FIXME This page is not fully translated, yet. Please help completing the translation.
(remove this paragraph once the translation is finished)

The perfect Peak

Das Beispiel des perfekten Gipfels wurde nicht nur aus persönlichem Interesse gewählt, sondern auch weil es ein gelungenes Beispiel für die nahezu beliebige Komplexität, der bei der Beantwortung einer scheinbar einfachen räumlichen Frage zu berücksichtigenden Faktoren, bietet. Gleichzeitig ist es keinen wissenschaftlichen Paradigmen zugeordnet, was die inhaltliche und wissenschaftliche Beschäftigung mit der Thematik deutlich vereinfacht, da so auf die Methodik fokussiert werden kann. Gerade dieser pragmatische Ansatz soll die seltene Kompetenz entwickeln helfen sowohl die iterative (und häufig auch rekursive) Entwicklung wissenschaftlicher Fragestellungen zu trainieren als auch die notwendigen Fähigkeiten vermitteln die konzeptionellen Lösungen auch eigenständig umzusetzen zu können.

Und ganz wichtig- es geht immer besser. Trotz der Methodenlastigkeit geht es hier nicht um “gute Programmierung” sondern darum die Fragen die einem so in den Sinn kommen beantworten zu können und nicht umgekehrt die Antwort meinen Möglichkeiten zu fragen anzupassen.

Ungeachtet dessen sind alle eingeladen klarer und eleganter zu fragen und zu antworten!

Einleitung und Hintergrund

Der perfekte Gipfel kann - folgt man Christian Rauch im Panorama Artikel 2/2012, S. 112 [(http://www.alpenverein.de/dav-services/panorama-magazin/dominanz-prominenz-eigenstaendigkeit-eines-berges_aid_11186.html)], unter Vernachlässigung ästhetischer oder kultureller Einflüsse, mit dem Eigenständigkeitswert berechnet werden. Unsere rein praktisch ausgerichtete Problemstellung lautet daher:

Berechne für jeden Gipfel den Eigenständigkeitswert.

Zunächst gilt es den Artikel von Rauch zu analysieren und der teilweise prosaischen Betrachtungsweise die Begriffe und Parameter zu entnehmen, die für die Berechnung der Eigenständigkeit notwendig sind. Folgt man wie der Autor Rauch der Definition von Bernhard [( Bernhard, W., Eigenständigkeit von Gipfeln, http://www.thehighrisepages.de/bergtouren/na_orogr.htm#BK , Letzter Zugriff: 1.5.2012)] erscheint dies zunächst einfach. So schreibt Bernhard:

“Um Gipfel auf ihre geografische, genauer: orografische Bedeutung hin zu betrachten, nimmt der Autor eine Berechnung der Eigenständigkeit (E) nach Höhe, Dominanz und relativer Prominenz vor. Diese drei Kennzahlen dürften die Bedeutung eines Berges am ehesten erfassen. Die Prominenz gibt an, wie tief ein Gipfel herausgearbeitet ist, die Dominanz, wie weit dieser von einem höheren entfernt ist, und die Höhe, wie weit er selbst andere zu überragen vermag. Diese drei Größen sind messbar und damit objektiv und einheitlich, einfach zu erheben, und nahtlos anwendbar für die gesamte Erde, vom Achttausender bis zum Küstenhügel”.

Aus wissenschaftshistorischer und geographischer Sicht greifen die Betrachtungen von Bernhard gewiss zu kurz. Es gibt eine lange fachübergreifende Diskussion zur kognitiv-ästhetischen, orographischen, morphologischen und morphometrischen Bestimmung von Geländeformen im Relief. Dies gilt um so mehr für die Bestimmung der Begriffe Berge, Scharten, Gipfel etc.. Eine für unsere Zwecke guten Überblick bietet Rasemann (2004:8-10) [(Rasemann, S. (2004), Geomorphometrische Struktur eines mesoskaligen alpinen Geosystems, Asgard-Verlag, St. Augustin, Bonner Geographische Abhandlungen, Heft Nr. 111, URL: http://hss.ulb.uni-bonn.de/2003/0211/0211.htm)]. Für eine konzeptionelle Analyse zur Umsetzung der Bernhard'schen Berechnung müssen diese kontroversen Definitionen (zunächst) nicht näher in Betracht gezogen werden. Wir halten uns an Bernhard's Postulat, dass die Berechnung objektiv und universell anwendbar sei. Jedoch wird sich zeigen dass nicht nur die Berücksichtigung der ästhetischen und kognitiven Aspekte der von Rauch geführten Diskussion sondern auch die objektive Klassifikation von Gipfeln, Scharten etc. eine deutlich weitergehende Auseinandersetzung mit diesem Thema erforderlich machen.

Letztlich wird in der Spannweite von dem Ansatz “Gipfelkoordinaten und Namen werden als als gegeben betrachtet” bis hin zu Konzept “Gipfel werden ausschließlich aus der Geländeform der verfügbaren Daten bestimmt” eine Definition notwendig werden. Diese Festlegung wird im wesentlichen von ihrer Zweckmäßigkeit zur Beantwortung der Fragestellung abhängen.

Konzeption

Die Fragestellung ist im Detail einiges komplexer als der erste Eindruck es vermuten lässt. Daher ist es sinnvoll zunächst nur die scheinbar einfachste, weil objektiv berechenbare Vorgehensweise zur Bestimmung der Eigenständigkeit eines Gipfels konzeptionell umzusetzen.

Definitionen und Annahmen

Analysiert man den Text bzw. die Formel zur Berechnung der Eigenständigkeit (E) verwendet der Autor folgende Größen:

h = Höhe                        (absolute Höhe des betrachteten Gipfels)
d = Dominanz                    (horizontale Distanz zum nächsten, höheren Geländepunkt)
p = Prominenz/Schartendifferenz (Höhendifferenz zur tiefsten Scharte im Verbindungsgrat zum nächsthöheren Gipfel)

Diese Begriffe zureichend zu diskutieren kann in diesem Zusammenhang nicht in erschöpfender Weise geschehen. Dennoch müssen einige Überlegungen getätigt werden:

  1. Auf welcher Datengrundlage (räumliche Auflösung, Datenqualität) sollen die obigen Größen berechnet werden?
  2. Wie sollen Gipfel identifiziert bzw. definiert werden?
  3. Wie kann auf dieser Grundlage die Schartendifferenz definiert werden?

Es werden folgende Annahmen getätigt:

  • Zur räumlichen, quasi-kontinuierlichen Bestimmung von Höhenwerten sind rastermodell-basierte Geländemodelle geeignet. Dank der mittlerweile guten Verfügbarkeit in den unterschiedlichsten räumlichen Auflösungen sind rasterbasierte digitale Geländemodells (DGM) für eine automatisierte Bearbeitung vorteilhaft. Generell kann die Berechnung der Eigenständigkeit mit allen verfügbaren Höhenmodelldaten durchgeführt werden. Die Berechnung erfolgt mit frei verfügbaren DGM-Daten. Es bieten sich vor allem die ASTER GDEM[(Aster: http://www.jspacesystems.or.jp/ersdac/GDEM/E/)] oder SRTM [(SRTM: http://srtm.csi.cgiar.org/)] Datensätze an.
  • Gipfel sind:
    • die Maximalwerte in einer definierten räumlichen Umgebung und/oder
    • die Zentroidpunkte einer morphometrisch abgeleiteten Formstruktur und/oder
    • Namenenskodierte Orte
  • Die maximale Schartendifferenz ist die Differenzhöhe zum niedrigsten Punkt der höchstmöglichen Verbindung vom betrachteten Gipfel zum nächsten höheren Gipfel.

Mögliche Algorithmen

Datenvorbereitung

  • DGM Daten einlesen (Formate und Umwandlung sind abhängig von der verwendeten Quelle)
  • Ergebnistabelle mit den folgenden Attributen anlegen
Gipfelkennung Hochwert Rechtswert Dominanzwert Prominenzwert Eigenständigkeitswert
ID

Gipfel identifizieren

Zur Identifikation von Gipfel stehen sehr unterschiedliche Möglichkeiten zur Verfügung. Prinzipiell lassen sie sich in zwei Kategorien unterscheiden. Sie können aus (1) den Daten selbst berechnet werden oder man kann sie (2) aus externen Quellen verfügbar machen. Daraus leiten sich die drei nachfolgenden Ideen ab.

  1. Gipfel nach der morphometrischen Kriterien bestimmen
    • Zentrale Koordinate der Gipfelposition (Zentroid der Fläche oder Zentroidkoordinate des Pixel) der morphometrischen Gipfelbestimmung extrahieren GipfelID zuweisen
    • Gipfel als lokales Höhenmaximun bestimmen (Moore-Nachbarschaft), Zentroidkoordinate des Pixel extrahieren GipfelID zuweisen
  2. Gipfelnamen aus Tabellen und Karten extrahieren
    • Namen von Gipfeln geokodieren bzw. Namen und Koordinaten von Gipfeln verwenden und räumlich möglichst am korrespondierender Position im DGM reverses Geokodieren verorten.

GipfelID, Koordinaten und Höhenwert → Tabelle schreiben

Dominanz berechnen

Die folgenden Schritte müssen für jeden Berg einzeln durchgeführt werden.

  • Für jede GipfelID ein Raster der euklidischen Distanz zu den größeren Höhenwerten berechnen
  • Minimalen Entfernungswert extrahieren und der Ergebnistabelle hinzufügen

Dominanz → Tabelle schreiben

Prominenz berechnen

Zur Berechnung der Prominenz stehen ebenfalls unterschiedliche Möglichkeiten zur Verfügung. Hier werden zwei Ideen vorgestellt. (1) Mit Hilfe eines statischen “Flutmodells” werden jeweils 2 Gipfel miteinander verbunden und (2) es wird eine Identifikation der die Gipfel verbindenden Grate als morphologischer Struktur durchgeführt. Im einzelnen leiten sich die beiden nachfolgenden Ideen ab:

  1. Flutmodell
    • Schrittweise Flächen zulassen, die niedrigere Höhenwerte haben als der Gipfel, bis die Schnittfläche Pixel nur zwei Gipfel enthält
    • Den Mindestwert dieses Rasters als Schartenhöhe extrahieren = Prominenz → Tabelle schreiben
  2. Gratmodell
    • von jedem Gipfel alle Grate zu höheren Gipfeln identifizieren
    • Jeweils den niedrigsten Höhenwert in eine Liste schreiben
    • maximaler Höhenwert dieser Liste = größte Schartenhöhe = Prominenz → Tabelle schreiben

Prominenz → Tabelle schreiben

Eigenständigkeit berechnen

Die Eigenständigkeit berechnet sich nach Bernhard [(#2)] wie folgt:

Für d <100.000m:

E = -((log2  (h / 8848) + log2 (d / 100000) + log2(p / h)) / 3) 

Für d > 100.000m:

E = -((log2 (h / 8848) + log2(p / h)) / 3)

Eigenständigkeit → Tabelle schreiben

Umsetzung der Konzepte mit R, SAGA, GRASS

Wie sich herausstellen wird gibt es gute Gründe ein (willkürliches) Mix an (GIS-)Software zu nutzen. Die Implementierung ist recht umfangreich dokumentiert und kommentiert. Übergeordnet geht es jedoch um die folgenden Kompetenzen:

  1. Erwerb struktureller Kompetenzen zur technische Umsetzung
    • Es soll exemplarisch gezeigt werden wie unterschiedlichste Bibliotheken verwendet werden können
    • Es soll gezeigt werden dass konzeptionelle Vorstellungen immer überprüfbar sein müssen und natürlich auch evaluiert gehören
    • Es soll ausführlich kommentiert und erläutert werden (inkl. wissenschaftlicher und technischer Verweise)
  2. Ganz wichtig ist die (zwar nur minimale) aber sehr wichtige fachliche Diskussion und Definition der Begriffe
    • Geländeformen, Gebirge
    • Regeln zur ästhetischen Klassifikation

Grundlage der gezeigten Skripte sind Arbeiten der MEGaGIS-Kursteilnehmer (Danke an Hanna). Für Einsteiger gibt es technische Starthilfe unterNützliche Links zu GRASS & Python und Warum RSAGA?.

Sehr hilfreich ist das WinR-GIS2Go Open Source GIS Paket für Windows.

Implementierung

Für die nachfolgenden Arbeitsschritte müssen die freien Softwarepakete R, GRASS und SAGA installiert sein. Das Script ist unter Ubuntu 11.10 mit der R Version 2.15.0 (64-bit), GRASS 6.4.2 und SAGA Version: 2.0.8 getestet. Alle aufgerufenen Bibliotheken sind aktualisiert.

Vorbereitung der Projekt-Arbeitsumgebung

Zunächst müssen für R die GRASS und SAGA Schnittstellenpakete und einige R-packages installiert werden. Auf der R-Konsole z.B.

install.packages(c("RSAGA", "spgrass6", "rgdal", "maptools", "raster", "sp"))

Dann geht es mit dem Import der packages und dem Festlegen der Arbeitspfade etc. weiter:

##Laden der grundlegenden libraries für die Verarbeitung räumlicher Daten
#-------.
library(sp)         # R spatial package
library(rgdal)      # R GDAL interface
library(maptools)   # R maptools package
library(spgrass6)  	# R GRASS interface
library(RSAGA)      # R SAGA interface
library(raster)     # R raster package
 
 
## Arbeitsumgebung setzten
#-------.
 
# Arbeitsverzeichnis R
#-------.
setwd("/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA")
Rhome="/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA"
 
 
# Pfade und SETUP für GRASS 
#-------.
 
initGRASS(gisBase="/usr/lib/grass64", home=tempdir(), 
          gisDbase="/home/creu/workspace/grassprojects", 
          location="GIPFEL", mapset="SRTM",
          override=TRUE) 
 
# INFOS zum MAPSET
gmeta6()
str(gmeta6())
 
# Pfade und SETUP für SAGA 
#-------.
 
rsaga.env(workspace="/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA")
myenv=rsaga.env(workspace="/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA",path="/usr/local/bin",modules="/usr/local/lib/saga")
SHome="/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA"

Daten

Die Pfade des persönlichen Arbeitsverzeichnisses sind natürlich wie immer anzupassen. Zunächst gilt es sich die Daten zu besorgen und für eine Berechnung verfügbar zu machen. Für das aktuelle Beispiel wurden SRTM Geländemodelldaten verwendet. Das Skript zeigt mehrere Möglichkeiten auf diese einzulesen. Ebenso wird beschrieben wie die Gipfelliste (teil automatisiert geladen und eingelesen werden kann.

#-------.
#   Daten Import/Export
#-------.
#  Es werden einige Möglichkeiten aufgezeigt verschiedene Daten für die weitere Verarbeitung verfügbar zumachen.
# 1) Datenimport aus XYZ Daten des Landesvermessungsamtes
# 2) Datenimport aus SRTM GEOTIFF Binärdaten 
#    jeweils in SAGA, GRASS, R 
#-------.
# Anmerkungen:
# Für alle weiteren Bedürfnisse muss leider auf die Modulbeschreibungen etc. verwiesen werden. 
# GEO-Datenhandling ist dank GDAL/OGR VIEL leichter geworden gehört aber immer noch mit zum Lästigsten überhaupt  
#-------.
 
 
 
## Import nach GRASS von XYZ ASCII Daten mit SPACE " " getrennt und Punkt "." als Dezimaltrennzeichen Achtung x=RW y=HW Z=h GAUSSKRüger DHDN 3 (so kommts vom LVA)
#execGRASS("r.in.xyz", flags=c("overwrite"), input="/home/creu/workspace//grassprojects/GIPFEL/SRTM/DGM10_3436996_5583209.xyz", output="DGM_1",  fs=" ", y=as.integer(2),x=as.integer(1), z=as.integer(3),pth=as.integer(1), zscale=1.0,percent=as.integer(100))
 
## Import nach SAGA von XYZ ASCII Daten mit SPACE " " getrennt und Punkt "." als Dezimaltrennzeichen Achtung x=RW y=HW Z=h GAUSSKRüger DHDN 3 (so kommts vom LVA)
#rsaga.geoprocessor("io_grid",6,env=myenv,list(FILENAME="DGM10_3436996_5583209.xyz", GRID="DG", CELLSIZE="10",SEPARATOR="space"))
 
## Import nach R von XYZ ASCII Daten mit SPACE " " getrennt und Punkt "." als Dezimaltrennzeichen Achtung kein header!
#dor <- read.table("/home/creu/workspace/grassprojects/GIPFEL/SRTM/DGM10_3436996_5583209.xyz", header=F)
 
## Import nach GRASS von GDAL GeoTIFF Daten Achtung dise müssen zuvor projiziert sein
execGRASS("r.in.gdal", flags=c("overwrite"), input="/home/creu/workspace/grassprojects/GIPFEL/SRTM/srtm_39_03/srtm_39_03utm.tif", output="DGM")
 
 
## Export von GRASS als GEOTIFF 
## (Workaround ist z.B. für SAGA statt direktem Auschreiben notwendig wegen möglicher NAN Fehler in der GDAL Bib)
## Wird hier gemacht6 weils unten je nach Algorithmus mit SAGA Formaten weitergeht 
execGRASS("r.out.gdal", flags=c("c"),input="DGM", type="Float32", output="/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/DGM.tif")
 
## Import SAGA ACHTUNG Bennennung muss GRIDS und FILES heissen NICHT GRID/FILENAME
rsaga.geoprocessor("io_gdal",0,env=myenv,list(GRIDS="DGM.sgrd" , FILES="DGM.tif"))
 
if (gipfelliste==T){
## Tourenwelt Bergliste enthält den wohl aktuellsten Datensatz zu Gipfeln 
## Der wird im folgenden geladen eingelesen reprojiziert und auf das Untersuchunggebiet ausgeschitten
## benötigt Internet
## TODO "shapes_grid",10 benötigt noch die von Hand eingetragfenen Parameter des Grids 
## Die Projektion ist auf UTM festverdrahtet muss also für jede U-Gebiet angepasst werden  
## Die Bergliste muss noch geparst werden da so nur die Koordinaten + Namen (gpx import) oder alles nur keine namen Cvs import
## verfügbar sind (warum auch immer) ALTERNATIVE reverse Geocoding
## Natürlich könnte man auch ein bisschen tricksen. 
## 1) konvertieren als CVS dann austausch der HTM Zeichen von hand oder mit kommandozeile
## 2) Import des ganzen als gpx alles wegschmeisen bis auf koordinaten und Namen
## 3) Merge von beiden
# download der datei
system("wget http://www.tourenwelt.info/commons/download/bergliste-komplett.kmz.php")
# umbenenen da gunzip nur .gz verarbeitet
system("mv /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.kmz.php /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.gz")
# entzippen da es sich um kmz handelt pgsbabel aber nur kml kann
system("gunzip -f /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.gz")
#anhängen der  kml endung
system("mv -f /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.kml")
# export als gpx datei
system("gpsbabel -i kml -f /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.kml -o gpx -F /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.gpx")
#export als csv datei
system("gpsbabel -i kml -f /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.kml -o unicsv -F /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.csv")
}

Gipfelerkennung

Gipfel können als Koordinaten einer morphometrisch analysierten räumlichen Umgebung betrachtet die eine gewisse Form aufwiesen betrachtet werden. Dies ermöglicht eine begründbare Ableitung von Gipfeln aufgrund von räumlichen Formmerkmalen. Hinsichtlich einer morphologischen und räumlichen Diskussion des Begriffs Gipfel und letztlich auch hinsichtlich der kognitiv-ästhetischen Betrachtungsweise des die Fragestellung initierenden Artikels erscheint dies als erfolgversprechender Ansatz. Wood (1996) diskutiert [( Wood, J.D. (1996) The geomorphological characterisation of digital elevation models PhD Thesis, University of Leicester, UK, http://www.soi.city.ac.uk/~jwo/phd)]) das Potenzial, dass in einer multiskaligen, morphometrischen Analyse von Geländeformen liegt. Allerdings sind vor allem aufgrund der räumlichen Skalenübergange des betrachteten Formenschatzes erhebliche Unschärfen zu erwarten (vgl. Fisher et al. 2004 [(Fisher, P., Wood, J. and Cheng, T. (2004): Where is Helvellyn? Fuzziness of Multiscale Landscape Morphometry, Transactions of the Institute of British Geographers, 29(1), pp.106-128.)]. Die von Fisher et. al. 2004 diskutierten Unschärfen finden sich in der zuvor angesprochen geographisch geprägten kontroversen Entwicklung der Begrifflichkeit und Benennung solcher Formen in den unterschiedlichen orographischen Strukturen wieder (vgl. Rasemann (2003)[(#1)]. Berücksichtigt man diese Problematik erscheint es zumindest aus pragmatischer informationstechnischer Sicht als sinnvoll und berechtigt nur nach lokalen Maxima auf Grundlage der direkten Moore-Nachbarschaft Gipfel zu bestimmen.

Letztlich wird allerdings bereits hier deutlich, dass es je nach Zielsetzung, wesentlich sinnvoller sein könnte, nicht aus der Geländeform auf die Gipfel zu schliessen, sondern verfügbare Gipfelkoordinaten zb. “Harry's Gipfelliste” [( Breitkreutz, H.: Gipfelliste http://www.tourenwelt.info/bergliste/bergliste.php. URL:http://www.tourenwelt.info/commons/download/bergliste-komplett.kmz.php. (Zugriff: 14.06.2012))] (natürlich mit namentlicher Zuweisung) im Relief zu verorten und das DGM nur zur Berechnung der Parameter heranzuziehen.

#-------. # Gipfelerkennung #-------. # Variante A # 1) Berechnung morhophometrischer Gipfel-Gridzellen *) # 2) Gipfelgridzellen werden als Punktkkordinaten (Zentroid) mit ihrem Höhenwert extrahiert und # 3) als R-Tabelle verfügbar gemacht. # # Varante B # 1) Filterung der Daten (am einfachsten Mean oder Gauss) # 2) Bestimmung aller lokalen Maxima # 3) MAximum-Gridzellen werden als Punktkkordinaten mit ihrem Höhenwert extrahiert und # 4) als R-Tabelle verfügbar gemacht. # # Variante C # 1) Gipfelliste (z.B. http://www.tourenwelt.info) importieren und bereinigen (siehe Datenimport) # # # *) Wood, J.D. (1996) The geomorphological characterisation of digital elevation models, http://www.soi.city.ac.uk/~jwo/phd #-------. if (wood==T) { # Gipfelerkennung nach Wood (1996) ## Der Parameter s_tol bestimmt die Krümmungstoleranz der Gipfelfäche, c_tol der Krümmungstoleranz von Ebenen ## size gibt die zur Berechnung verwendete Kernel-Größe (Moore-Nachbarschaft) in Pixel an execGRASS("r.param.scale", flags=c("overwrite"), input="DGM", output="gipfel", s_tol=1.0, c_tol=0.00001, size=as.integer(3), param="feature") # Export aus GRASS des Gipfelerkennung als GEOTIFF execGRASS("r.out.gdal", input="gipfel", output="/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/peak_form.tif") ## Import GEOTIFF Gipfelerkennung nach SAGA ACHTUNG Bennennung muss GRIDS und FILES heissen NICHT GRID/FILENAME rsaga.geoprocessor("io_gdal",0,env=myenv,list(GRIDS="peak_form.sgrd" , FILES="peak_form.tif")) ## Wood bestimmt 6 Form-Klassen ## 1=planar,2=pit,3=channel,4=pass,5=ridge,6=peak ## Da uns nur 6=peak interessiert müssen wir das grid ## reklassifizieren. GRASS benötigt hierzu eine Textdatei ## mit Regeln falls die Textdatei "reclass" bereit existiert ## wird sie mit dem Systembefehl "rm reclass" gelöscht system("rm reclass") ## Erzeugen (linux shell) einer Lookuptabelle zur Reklassifizierung der morphometrischen Berechnung ## peak werden zu 1 alle anderen Werte zu 0 system("echo '6 = 1 * = NULL end'
>  reclass" )
 
## reklassifizierung unter Verwendung der erzeugten Datei reclass
execGRASS("r.reclass",flags=c("overwrite"), input="gipfel", output="gipfel_bool", rules="reclass")
 
## Da Gipfel nach Wood anders als lokale Maxima Gebiete ähnlicher Form sind
## können sie z.T. viele einzelene (zusammenhängende) Pixel umfassen
## Ein einfacher Weg aus diesen einen Punkt zu erzeugen ist die Umwandlung der 
## reklassifizierten Raster-Daten in Gipfel-Vektor-Polygone
execGRASS("r.to.vect", flags=c("overwrite"), input="gipfel_bool", output="gipfel_vpoly_bool", feature="area")
 
## Von diesen Poolygonen wird der Zentroid (Flächenschwerpunkt) berechnet und ausschreiben der centroid punkte der Gipfel-Polygone
## um den Gipfel als Punktkoordinate zu erhalten
execGRASS("v.type", flags=c("overwrite"),input="gipfel_vpoly_bool", output="gipfel_vpoint_bool", type="centroid,point")
 
## Zur Speicherung der Höhenwerte aus dem Geländemodell
## muss der erzeugten Vektordatei erst eine neue Attributspalte 
## in der Attributtabelle hinzugefügt werden 
execGRASS("v.db.addcol", map="gipfel_vpoint_bool", columns="heights double precision")
 
## Jetzt noch die Höhenwerte an den Centroid-Koordinaten 
## aus dem DGM Raster extrahieren und in die neue Spalte schreiben
execGRASS("v.what.rast", vector="gipfel_vpoint_bool", raster="DGM", column="heights")
gipfel=read.csv("gipfel.csv", header = FALSE, sep = "", dec=".")
# Löschen der Variable 3 (ID der Vektordatei)
gipfel$V3 <- NULL
 
} # Ende Wood -----------------------------------------------------------------
 
 
# Beginn MaxPix -----------------------------------------------------------------
if (max==T){ #maxPix
  # Gipfelerkennung "lokale Maxima" METHOD=0 ist der MEAN Filter
 
  ## Auschreiben der Reclass Tabelle für SAGA. 
  ## Achtung die DateiDARF NICHT beim Einlesen nach SAGA die Endung csv haben!
  ## AUSSERDEM MUSS sie 6 Ziffern nach dem Dezimalpunkt aufweisen (nsmall = 6)
#  hoehe =minGipfelhoehe
#  rec <- structure(list(MIN = c(0.0), MAX = c(hoehe), CODE = c(0.0)), .Names = c("MIN","MAX", "CODE"), row.names = c(NA,1L), class = "data.frame")
#  write.table(format(rec, nsmall = 6), "rec.txt",row.names = FALSE, quote = F, sep ="\t")
 
  ## Reklassifizieren des GRIDS
#  rsaga.geoprocessor("grid_tools",15,list(INPUT="DGM.sgrd",RESULT="DGMrec",METHOD=2,RETAB="rec.txt",ROPERATOR=1 ))
 
  rsaga.geoprocessor("grid_filter",0,list(INPUT="DGM.sgrd",RESULT="fDGM",MODE=0,METHOD=0,RADIUS=filterumgebung))
 
 
  # Shapes_Grid extrahiert lokale Maxima und Minima Positionen und ihre Höhenwerte (peaks und pits) in shapefiles
  rsaga.geoprocessor("shapes_grid",9,list(GRID="fDGM.sgrd", MAXIMA="localmax",MINIMA="localmin"))       
 
  # Das localmax shapefile enthält in der Spalte Z die Höhenwerte wird hier in eine CSV Textdatei exportiert
  rsaga.geoprocessor("io_shapes",2,list(SHAPES="localmax.shp",FIELD="Z",SEPARATE=0,FILENAME="gipfel.csv"))
 
  ## die csv datei der Gipfel (gleichgültig welche)  wird in die R-Tabelle gipfel geschrieben
  gipfel=read.csv("gipfel.csv", header = FALSE, sep = "")
  # es werden Spaltennamen zuweisen
  colnames(gipfel)=c("X","Y","Z")
 
} # Ende max--------------------------------------------------------------------
 
 
 
# Beginn Gipfelliste -----------------------------------------------------------------
if (gipfelliste==T){ 
# ermitteln des Extend damit ein polygon zum clippen erzeugt werden kann
rsaga.geoprocessor("grid_tools",15,list(INPUT="fDGM.sgrd",RESULT="extent", METHOD=0, OLD=9999.,NEW=100,SOPERATOR=1))
# erzeugen des Polygons
rsaga.geoprocessor("shapes_grid",10,list(SHAPES="umfang.shp", CELLS="0", PARAMETERS_GRID_SYSTEM_NX="1277",  PARAMETERS_GRID_SYSTEM_NY="1120", PARAMETERS_GRID_SYSTEM_X="577537.5",PARAMETERS_GRID_SYSTEM_Y="5211037.5",PARAMETERS_GRID_SYSTEM_D="75"))
 
#### einmal der durchgang 1) Einlesen 2) projizieren 3) clippen für die (von hand korrigierte) CSV liste
rsaga.geoprocessor("io_shapes",3,list(SHAPES="bergliste5.shp", X_FIELD=2,Y_FIELD=1,FILENAME="bergliste-komplett_16.05_2012.csv"))
rsaga.geoprocessor("pj_proj4",5,list(SOURCE_PROJ="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0",TARGET_PROJ="+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0",SOURCE="bergliste5.shp",TARGET="bergliste5_utm.shp"))
rsaga.geoprocessor("shapes_points",8,list(POINTS="bergliste5_utm.shp" ,POLYGONS="umfang.shp" ,CLIPS="bergliste5_cut.shp", METHOD="0"))
 
#### einmal der durchgang 1) Einlesen 2) projizieren 3) clippen für die gpx liste
rsaga.geoprocessor("io_shapes",13,list(SHAPES="bergliste.shp",FILE="bergliste-komplett.gpx"))
rsaga.geoprocessor("pj_proj4",5,list(SOURCE_PROJ="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0",TARGET_PROJ="+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0",SOURCE="bergliste.shp",TARGET="bergliste_utm.shp"))
rsaga.geoprocessor("shapes_points",8,list(POINTS="bergliste_utm.shp" ,POLYGONS="umfang.shp" ,CLIPS="bergliste_cut.shp", METHOD="0"))
 
### Todo Merge der beiden listen
 
# in eine CSV Textdatei exportiert
rsaga.geoprocessor("io_shapes",2,list(SHAPES="bergliste5_cut.shp",  ALL=T, SEPARATE=1,FILENAME="gipfel.txt"))
 
## die csv datei der Gipfel (gleichgültig welche)  wird in die R-Tabelle gipfel geschrieben
gipfel=read.csv("gipfel.csv", header = FALSE, sep = ",")
# es werden Spaltennamen zuweisen
colnames(gipfel)=c("X","Y","Z")
 
} # Ende Gipfelliste

<imgcaption abb1| Beispiel Gipfelerkennung Lahnberge im Raum Marburg> [ ] </imgcaption>

Berechnung der Dominanz

Hier wird für jeden Gipfel ein “proximity Grid” erstellt, was für jeden Pixel die Entfernung zum Gipfel enthält. Zusätzlich wird ein binäres Raster erstellt, welches eine 1 für höhere Pixel und eine 0 für niedrigere Pixel als der Gipfel enthält. Diese Grids werden multipliziert. Der niedrigste Wert davon ist dann die Dominanz. Diese wird als neue Spalte in die Tabelle “gipfel” geschrieben.

<imgcaption abb2| Berechnung der Dominanz>

[] </imgcaption>

#-------.
#	Dominanz
#-------.
# 1) Für jede Gipfel-Koordionate wird mit "proximity" ein Distanz-Grid erstellt
#    (Resultat ist die Entfernung jedes Pixels zum Gipfel. 
# 2) Zusätzlich15wird ein Maske erstellt (1 height>gipfel, 0 height <= Gipfel
# 3) Dominanz=minimum(proximity*maske) 
# 4) Dominaz wird als neue Variable in die Tabelle "gipfel" eingetragen 
#-------.
 
for (i in 1: nrow(gipfel))                                   				  # Für jeden Gipfel
{
 
  # Zeile vom Gipfel i als Text Tabelle speichern
  write.table(gipfel[i,],"tmp1.txt",row.names=FALSE)                     
 
  # Tabelle in SAGA importieren
  rsaga.geoprocessor("io_table",1,env=myenv, list(TABLE="peak",HEADLINE=TRUE,SEPARATOR="space", FILENAME="tmp1.txt"))
 
  # Tabellenimport in Shapedatei konvertieren
  rsaga.geoprocessor("shapes_points",0,env=myenv,list(POINTS="peak",TABLE="peak", X="X",Y="Y")) 
 
  # Erzeugte Shapedatei in eine Rasterdatei konvertieren ACHTUNG die Datei peak.sgrd muss bereits existieren!
  rsaga.geoprocessor("grid_gridding",0,list(INPUT="peak.shp", FIELD="X",LINE_TYPE=2,GRID_TYPE=3,TARGET=1,GRID_GRID="peak.sgrd"))    
 
 
  ## Höhe vom Gipfel aufrunden (simpler workaround damit Pixelwert sicher höher ist als seine Nachbarpixel) 
  ## und in Variable hoehe schreiben
  hoehe=ceiling(gipfel$Z[i])       	
 
  ## SAGA Single Reclass (Method = 0): Alle die <Z =0 (SOPERATOR=3 heißt <)
  #rsaga.geoprocessor("grid_tools",15,list(INPUT="fDGM.sgrd",RESULT="high_px", METHOD=2,SOPERATOR=1))
  ## SAGA Single Reclass (Method = 0): Alle Pixel die  >= Z = 1 (SOPERATOR=3 heißt >=)
  #rsaga.geoprocessor("grid_tools",15,list(INPUT="low_pix.sgrd", RESULT="high_pix",METHOD=0, OLD=hoehe,NEW=1,SOPERATOR=3))
 
  ## Auschreiben der Reclass Tabelle für SAGA. 
  ## Achtung die DateiDARF NICHT beim Einlesen nach SAGA die Endung csv haben!
  ## AUSSERDEM MUSS sie 6 Ziffern nach dem Dezimalpunkt aufweisen (nsmall = 6)
  rec <- structure(list(MIN = c(0.0, hoehe), MAX = c(hoehe, 9999.0), CODE = c(0.0,1.0)), .Names = c("MIN","MAX", "CODE"), row.names = c(NA,2L), class = "data.frame")
  write.table(format(rec, nsmall = 6), "rec.txt",row.names = FALSE, quote = F, sep ="\t")
 
  ## Reklassifizieren des GRIDS
  rsaga.geoprocessor("grid_tools",15,list(INPUT="fDGM.sgrd",RESULT="localmax",METHOD=2,RETAB="rec.txt",ROPERATOR=1 ))
 
  ## SAGA Proximity grid erstellen
  rsaga.geoprocessor("grid_tools",26,list(FEATURES="peak.sgrd",DISTANCE="proximity"))
 
  ## Entfernunggrid (0/1) mit Höhengrid  multiplizieren
  rsaga.geoprocessor("grid_calculus",1,list(GRIDS=list("proximity.sgrd;high_pix.sgrd"), RESULT="dominanz", FORMULA="a*b"))
 
  ## Ergebnis in R verfügbar machen
  dominanz_grd=raster("dominanz.sgrd", band=1,native=TRUE)  
  tmp=getValues(dominanz_grd)
  if (any(tmp>0)){                                        #geringste Distanz zu höheren =Dominanz
    gipfel$dominanz[i]=min(tmp[tmp>0])                      #wird in gipfel-tabelle geschrieben
 
 
  } else {
    gipfel$dominanz[i]=max(tmp[tmp>0])  	                   #Fall für höchsten Gipfel   
  }
 
 
}

Berechnung der Prominenz

Teil 1 Graterkennung und Schartenhöhenmatrix

Hier soll für jeden Gipfel die Schartenhöhe zu jedem anderen Gipfel ermittelt werden. zunächst wird festgestellt, ob der Gipfel mit dem “Ziel-Gipfel” über einen Grat verbunden ist. Wenn das der Fall ist, wird ein Höhenprofil aus dieser Verbindung erstellt.Die minimale Höhe in diesem Profil ist dann die Schartenhöhe. Wenn die Gipfel nicht verbunden sind,wird ein Höhenprofil aus der direkten Verbindung erstellt und der min Wert davon ist die Schartenhöhe.

<imgcaption abb3| Berechnung der Schartenhöhe> [ ] </imgcaption>

#-------.
#	 SCHARTENHÖHE bzw. Prominenz   Variante B                                                                                                       
#-------.
#  1) Für jeden Gipfel wird die Schartenhöhe zu jedem höheren Gipfel ermittelt.                                           
#  2) Falls Gipfel1 und Gipfeln über einen Grat verbunden sind wird das Grat-Höhenprofil ermittelt                        
#  3) Die minimale Höhe des Profils ist die Schartenhöhe.                                                                 
#  4) Das Maximum aller ermittelten Schartenhöhen ist die gesuchte maximal Schartenhöhe                                   
#  5) Falls Gipfel und Zielgipfel NICHT über einen Grat miteinander verbunden sind wird ein Höhenprofil                   
#     aus der Sichtlinie zwischen den Gipfelpunkten erstellt dann weiter mit3) und 4)                                     
#-------.
 
##  Damit mit einer Flow-Analyse das Routing entlang der Grate Sattfinden kann
##  wird das gefilterte DGM invertiert
rsaga.geoprocessor("grid_calculus",1,list(GRIDS="fDGM.sgrd",RESULT="ifDGM",FORMULA="-a"))
 
# Hydrologische Korrektur - Sink Removal Output: pitlessifDGM
rsaga.geoprocessor("ta_preprocessor",2,list(DEM="ifDGM.sgrd",DEM_PREPROC="pitlessifDGM", METHOD=1))
 
# Flow Analysis. Output: direction
rsaga.geoprocessor("ta_channels",5,list(DEM="pitlessifDGM.sgrd",DIRECTION="direction",CONNECTION="tmp", SEGMENTS="segments",BASINS="basins",THRESHOLD=0))
 
# Erzeugen einer Matrix für die Schartenhöhen
schartenhoehe=matrix(ncol=nrow(gipfel),nrow=nrow(gipfel)) 
 
# Füllen der Matrix mit NA
schartenhoehe=schartenhoehe==NA                            
 
 
for (i in 1:nrow(schartenhoehe)) 				                   # für jeden Gipfel
{
 
  write.table(gipfel[i,],"tmp1.txt",row.names=F)    		   # Gipfelkoordinate die geroutet werden soll
  rsaga.geoprocessor("io_table",1, list(TABLE="DerGipfel",HEADLINE=TRUE,SEPARATOR="space",FILENAME="tmp1.txt"))
  rsaga.geoprocessor("shapes_points",0,list(POINTS="ZielGipfel",TABLE="DerGipfel",X="X",Y="Y"))
  if (i==1){schartenhoehe[i,1]=NA} 
  else     {if (i==2){
                      k=1
                      write.table((gipfel[k,]),"tmp2.txt",row.names=F)
                      rsaga.geoprocessor("io_table",1, list(TABLE="AndereGipfel",HEADLINE=TRUE,SEPARATOR="space",FILENAME="tmp2.txt"))
                      rsaga.geoprocessor("shapes_points",0,list(POINTS="Gipfelshp",TABLE="AndereGipfel",X="X",Y="Y"))
                      rsaga.geoprocessor("shapes_tools",2,list(OUT="merged",MAIN="ZielGipfel.shp",LAYERS="Gipfelshp.shp"))
                      rsaga.geoprocessor("grid_gridding",0,list(INPUT="merged.shp",FIELD="X",LINE_TYPE=2,GRID_TYPE=3,TARGET=1,GRID_GRID="tmp.sgrd"))
                      rsaga.geoprocessor("ta_channels",0,list(ELEVATION="pitlessifDGM.sgrd",SINKROUTE="direction.sgrd",CHNLNTWRK="tmp",CHNLROUTE="tmp",SHAPES="line",INIT_GRID="tmp.sgrd",INIT_METHOD=2,INIT_VALUE=0,DIV_CELLS=5,MINLEN=10))    
                      rsaga.geoprocessor("shapes_polygons",6,list(SHAPES="line.shp",POINTS="networkpoints"))
                      rsaga.geoprocessor("io_shapes",2,list(SHAPES="networkpoints.shp",FIELD="Z",SEPARATE=0,FILENAME="Networkpoints.txt"))
                      netpoints=read.table("Networkpoints.txt",header=F,sep="\t")   	
 

An dieser Stelle wurde jetzt versucht die Gipfel i und k über einen Grat zu verbinden. Falls das gelingt wird exakt ein Linienvektor zwischen den Gipfeln erzeugt. In der Tabelle “netpoints” ist die ID dieser Linie enthalten. Es folglich nur eine ID geben, falls dies zutrifft, dann geht es wie folgt weiter:

               if (substr(netpoints[1,3],1,2)==substr(netpoints[nrow(netpoints),3],1,2)) 
                         { # es darf nur eine durchgehende linie sein!!
                           rsaga.geoprocessor("ta_profiles",4,list(DEM="DGM.sgrd",LINES="line.shp",PROFILE="profile")) 
                           rsaga.geoprocessor("io_shapes",2,list(SHAPES="profile.shp",FIELD="Z",SEPARATE=0,FILENAME="Elevprofile.txt"))
                           elevprofile= read.table("Elevprofile.txt",header=F,sep="\t")   # lade die nach R
                           schartenhoehe[i,k]=min (elevprofile[,3]) 	
                         }					

Falls kein Grat zwischen den beiden Gipfeln i und k vorhanden ist, wird schlicht die Luftlinie verwendet:

                      else
                          {    
                           for (k in 1:ncol(schartenhoehe[,1:(i-1)]))						  # für jeden anderen Berg (nicht doppelt berechnen!)
                               {
                                write.table((gipfel[k,]),"tmp2.txt",row.names=F)
                                rsaga.geoprocessor("io_table",1, list(TABLE="AndereGipfel",HEADLINE=TRUE,SEPARATOR="space",FILENAME="tmp2.txt"))
                                rsaga.geoprocessor("shapes_points",0,list(POINTS="Gipfelshp",TABLE="AndereGipfel",X="X",Y="Y")) 				 
                                rsaga.geoprocessor("shapes_tools",2,list(OUT="merged",MAIN="ZielGipfel.shp",LAYERS="Gipfelshp.shp"))  		     
                                rsaga.geoprocessor("grid_gridding",0,list(INPUT="merged.shp",FIELD="X",LINE_TYPE=2,GRID_TYPE=3,TARGET=1,GRID_GRID="tmp.sgrd"))  
                                rsaga.geoprocessor("ta_channels",0,list(ELEVATION="pitlessifDGM.sgrd",SINKROUTE="direction.sgrd",CHNLNTWRK="tmp",CHNLROUTE="tmp",SHAPES="line",INIT_GRID="tmp.sgrd",INIT_METHOD=2,INIT_VALUE=0,DIV_CELLS=5,MINLEN=10))      
                                rsaga.geoprocessor("shapes_polygons",6,list(SHAPES="line.shp",POINTS="networkpoints")) 
                                rsaga.geoprocessor("io_shapes",2,list(SHAPES="networkpoints.shp",FIELD="Z",SEPARATE=0,FILENAME="Networkpoints.txt"))
                                netpoints=read.table("Networkpoints.txt",header=F,sep="\t")  
                                if (substr(netpoints[1,3],1,2)==substr(netpoints[nrow(netpoints),3],1,2))        
                                   {        
                                    rsaga.geoprocessor("ta_profiles",4,list(DEM="DGM.sgrd",LINES="line.shp",PROFILE="profile"))
                                    rsaga.geoprocessor("io_shapes",2,list(SHAPES="profile.shp",FIELD="Z",SEPARATE=0,FILENAME="Elevprofile.txt"))
                                    elevprofile= read.table("Elevprofile.txt",header=F,sep="\t")   											
                                    schartenhoehe[i,k]=min (elevprofile[,3]) 	
                                    cat("Gipfel",i,"von",nrow(gipfel),"wird bearbeitet(",100*(k/ncol(schartenhoehe[,1:(i-1)])),"%). Methode: Grat")
                                   }					
                                 else
                                    {
                                     rsaga.geoprocessor("shapes_lines",1,list(LINES="line",POINTS="merged.shp",ORDER="Y")) 					         
                                     rsaga.geoprocessor("ta_profiles",4,list(DEM="DGM.sgrd",LINES="line.shp",PROFILE="profile")) 		         
                                     rsaga.geoprocessor("io_shapes",2,list(SHAPES="profile.shp",FIELD="Z",SEPARATE=0,FILENAME="Elevprofile.txt"))     
                                    elevprofile= read.table("Elevprofile.txt",header=F,sep="\t")   												     
                                    schartenhoehe[i,k]=min (elevprofile[,3]) 
                                    cat("Gipfel",i,"von",nrow(gipfel),"wird bearbeitet(",100*(k/ncol(schartenhoehe[,1:(i-1)])),"%). Methode: Sichtlinie")
                                    }
                              } 
                          }
       }
  cat("Gipfel",i,"von",nrow(gipfel),"fertig!")
  }

Jetzt gibt es für jeden Berg zu jedem anderen Berg eine Schartenhöhe. Die Matrix wird jetzt gespiegelt, damit die Abfrage nacher einfacher ist:

tmp=t(schartenhoehe)
schartenhoehe[is.na(schartenhoehe)]=tmp[is.na(schartenhoehe)]
write.table(schartenhoehe,"MatrixSchartenhoehe.txt",row.names=F)                               # damit man es nicht immer neu rechnen muss ;)
Teil 2 - Abfrage Prominenz

Hier wird für jeden Gipfel abgefragt, welche Schartenhöhe zu einem höheren Gipfel die höchste ist. Diese wird dann von der höhe des Gipfel subtrahiert und entspricht damit der Prominenz des Gipfels. Die Prominenz wird als weitere Spalte der Tabelle “gipfel” angehängt. Der Fall für den höchsten Berg ist nicht optimal: Da es keinen höheren gibt kann auch keine Prominenz zugewiesen werden. Um den bedeutsamen Berg nicht auszuschließen wird daher hier die minimale Schartenhöhe, die überhaupt in der Matrix vorkommt verwendet.

for (i in 1:nrow(gipfel))
{
  hoehere=which(gipfel$Z>gipfel$Z[i])
  ifelse ((length(hoehere)==0),                                           #...wenn es keine höhere Gipfel gibt
         (gipfel$prominenz[i]=gipfel$Z[i]-min(schartenhoehe[i,],na.rm=T)),#nehme tiefste Verbindung überhaupt
         (gipfel$prominenz[i]=gipfel$Z[i]-max(schartenhoehe[i,hoehere]))) #sonst:maximale Schartenhöhe 
}
gipfel$prominenz[gipfel$prominenz<0]=0.000001  #für den Fall das die Prominenz negativ ist weil Scharte größer
                                               #als Berg...Mit 0 kann man nachher nicht weiterrechnen, deshalb
                                               #wird eine Zahl fast 0 gewählt

Berechnung des Eigenständigkeitswertes

Die Spalte “E” wird zu der Tabelle “gipfel” hinzugefügt und der Eigenständigkeitswert entsprechend der Formeln berechnet.

gipfel$E[gipfel$dominanz<100000]=-1*((logb((gipfel$Z/8848),base=2)+
                                  logb((gipfel$dominanz/100000),base=2)+
                                  logb((gipfel$p/gipfel$Z),base=2))/3)
gipfel$E[gipfel$dominanz>=100000]=-1*((logb((gipfel$Z/8848),base=2)+
                                  logb((gipfel$p/gipfel$Z),base=2))/3)

<imgcaption abb3| Ergebnis der Eigenständigkeitsberechnung> [ ] </imgcaption>

Vollständiges R GRASS/SAGA Skript

Das Skript hat sich von den oben beschriebenen Teilskripten einigermassen fortentwickelt. Sollte bei Gelegenheit angepasst werden. Allerdings sind die Grundkonzepte gleich geblieben. Die Daten und die Verzeichnisstruktur gibts hier. Achtung alle Pfade etc. sollten überprüft werden. Fehler könnten an Schreibrechten liegen. bzw. das SAGA GRid-Dateien fehlen (Saga benötigt diese beim rastern von Vektordateien) Download der Daten fürs Skript

# GipfelWert
# Version 0.9 
# Datum 16.06.2012
# Zweck:      R Skript nutzt SAGA und GRASS zur
#             Berechnung des Eigenständigkeitswert eines Gipfels nach Leonhard mit: 
#             h = Höhe                        (absolute Höhe des betrachteten Gipfels)
#             d = Dominanz                    (horizontale Distanz zum nächsten, höheren Geländepunkt)
#             p = Prominenz/Schartendifferenz (Höhendifferenz zur tiefsten Scharte im Verbindungsgrat 
#                                              zum nächsthöheren Gipfel)
# Quellen:    Rauch. C. (2012): Der perfekte Gipfel.  Panorama, 2/2012, S. 112 
#             http://www.alpenverein.de/dav-services/panorama-magazin/dominanz-prominenz-eigenstaendigkeit-eines-berges_aid_11186.html
#             (Zugriff: 20.05.2012)
#             Leonhard, W. (2012): Eigenständigkeit von Gipfeln. - 
#             Bergtouren. - http://www.thehighrisepages.de/bergtouren/na_orogr.htm
#             (Zugriff: 28.04.2012)
#             Wood, J.D. (1996) The geomorphological characterisation of digital elevation models, 
#             http://www.soi.city.ac.uk/~jwo/phd
#             (Zugriff: 20.05.2012)
#             Breitkreutz, H.: Gipfelliste http://www.tourenwelt.info/bergliste/bergliste.php. 
#             URL:http://www.tourenwelt.info/commons/download/bergliste-komplett.kmz.php 
#             (Zugriff: 14.06.2012)
#             SAGA Modulbeschreibungen:  http://sourceforge.net/apps/trac/saga-gis/wiki/Module%20Descriptions  
#             GRASS Modulbeschreibungen: http://wgbis.ces.iisc.ernet.in/grass/grass64/manuals/html64_user/index.html
#
# Autor(en):  Hanna Meyer (Basisskript R & SAGA (Grate))
#             (c) Chris Reudenbach ( Erweiterung, Überarbeitung, Zusammenführung und Kommentare)
#             Kontakt: giswerk@gis-ma.org
# 
# Anmerkung:  Es wird nahezu ausschließlich mit den beiden generic wrappern rsaga.geoprocessor 
#             und execGRASS gearbeitet. Die packages RSAGA und spgrass6 bieten wesentlich effizientere 
#             z.T. Schnittstellen!
#             Vorteil der generischen Aufrufe ist, dass die Kommandozeilenaufrufe
#             der Pakete fast 1:1 genutzt werden können => steile Lernkurve
#             Achtung wer Rechner mit wenig Hauptspeicher hat sollte mit kleinen Datenmengen operieren. 
#             R ist da ganz schön lahm  :)
#-------.
# TODO:       Ne ganze Menge ;) Viel technische Optimierung siehe Kommentare
#             Vor allem aber konzeptionelle Verbesserungen
#             - Implementierung einer skaligen Gipfelsuche nach Wood (Schleife einbauen ;)
#             - Implementierung einer geeigneten Filtersuche             
#             - Optimierung der Gratsuche
#             - Automatische Korrektur der Gipfelliste 
#             - Implementierung von Regeln und Abfragen auf die "Art" des Gipfels
#                 - Landnutzung, Touren, Exposition Hütten etc. Hier gibts einen 
#                 - riesigen Datenfundus
#-------.
# Copyright: by-nc-sa 3.0 http://creativecommons.org/licenses/by-nc-sa/3.0/
#-------.
 
 
##Laden der grundlegenden libraries für die Verarbeitung räumlicher Daten
#-------.
library(sp)         # R spatial package
library(rgdal)      # R GDAL interface
library(maptools)   # R maptools package
library(spgrass6)  	# R GRASS interface
library(RSAGA)      # R SAGA interface
library(raster)     # R raster package
 
 
## Arbeitsumgebung setzten
#-------.
 
# Arbeitsverzeichnis R
#-------.
setwd("/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA")
Rhome="/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA"
 
 
# Pfade und SETUP für GRASS 
#-------.
 
initGRASS(gisBase="/usr/lib/grass64", home=tempdir(), 
          gisDbase="/home/creu/workspace/grassprojects", 
          location="GIPFEL", mapset="SRTM",
          override=TRUE) 
 
# INFOS zum MAPSET
gmeta6()
str(gmeta6())
 
# Pfade und SETUP für SAGA 
#-------.
 
rsaga.env(workspace="/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA")
myenv=rsaga.env(workspace="/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA",path="/usr/local/bin",modules="/usr/local/lib/saga")
SHome="/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA"
 
# Globale Variablen 
#-------.
filterumgebung=45
#minGipfelhoehe=900
 
 
## Schalter für Gipfelvarianten 
## max = lokale Maxima 
## wood = morphometrische Berechnung nach Wood
## gipfelliste = Gipfelinfos von http://www.tourenwelt.info/
max = T
wood = F
gipfelliste = F
 
## Schalter für Prominez 
## flut = Flutungskonzept  (noch nicht implementiert)
## grat = Graterkennung
flut = T
if (flut==T){grat = F} else {grat==T}
 
 
#-------.
#   Daten Import/Export
#-------.
#  Es werden einige Möglichkeiten aufgezeigt verschiedene Daten für die weitere Verarbeitung verfügbar zumachen.
# 1) Datenimport aus XYZ Daten des Landesvermessungsamtes
# 2) Datenimport aus SRTM GEOTIFF Binärdaten 
#    jeweils in SAGA, GRASS, R 
#-------.
# Anmerkungen:
# Für alle weiteren Bedürfnisse muss leider auf die Modulbeschreibungen etc. verwiesen werden. 
# GEO-Datenhandling ist dank GDAL/OGR VIEL leichter geworden gehört aber immer noch mit zum Lästigsten überhaupt  
#-------.
 
 
 
## Import nach GRASS von XYZ ASCII Daten mit SPACE " " getrennt und Punkt "." als Dezimaltrennzeichen Achtung x=RW y=HW Z=h GAUSSKRüger DHDN 3 (so kommts vom LVA)
#execGRASS("r.in.xyz", flags=c("overwrite"), input="/home/creu/workspace//grassprojects/GIPFEL/SRTM/DGM10_3436996_5583209.xyz", output="DGM_1",  fs=" ", y=as.integer(2),x=as.integer(1), z=as.integer(3),pth=as.integer(1), zscale=1.0,percent=as.integer(100))
 
## Import nach SAGA von XYZ ASCII Daten mit SPACE " " getrennt und Punkt "." als Dezimaltrennzeichen Achtung x=RW y=HW Z=h GAUSSKRüger DHDN 3 (so kommts vom LVA)
#rsaga.geoprocessor("io_grid",6,env=myenv,list(FILENAME="DGM10_3436996_5583209.xyz", GRID="DG", CELLSIZE="10",SEPARATOR="space"))
 
## Import nach R von XYZ ASCII Daten mit SPACE " " getrennt und Punkt "." als Dezimaltrennzeichen Achtung kein header!
#dor <- read.table("/home/creu/workspace/grassprojects/GIPFEL/SRTM/DGM10_3436996_5583209.xyz", header=F)
 
## Import nach GRASS von GDAL GeoTIFF Daten Achtung dise müssen zuvor projiziert sein
execGRASS("r.in.gdal", flags=c("overwrite"), input="/home/creu/workspace/grassprojects/GIPFEL/SRTM/srtm_39_03/srtm_39_03utm.tif", output="DGM")
 
 
## Export von GRASS als GEOTIFF 
## (Workaround ist z.B. für SAGA statt direktem Auschreiben notwendig wegen möglicher NAN Fehler in der GDAL Bib)
## Wird hier gemacht6 weils unten je nach Algorithmus mit SAGA Formaten weitergeht 
execGRASS("r.out.gdal", flags=c("c"),input="DGM", type="Float32", output="/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/DGM.tif")
 
## Import SAGA ACHTUNG Bennennung muss GRIDS und FILES heissen NICHT GRID/FILENAME
rsaga.geoprocessor("io_gdal",0,env=myenv,list(GRIDS="DGM.sgrd" , FILES="DGM.tif"))
 
if (gipfelliste==T){
## Tourenwelt Bergliste enthält den wohl aktuellsten Datensatz zu Gipfeln 
## Der wird im folgenden geladen eingelesen reprojiziert und auf das Untersuchunggebiet ausgeschitten
## benötigt Internet
## TODO "shapes_grid",10 benötigt noch die von Hand eingetragfenen Parameter des Grids 
## Die Projektion ist auf UTM festverdrahtet muss also für jede U-Gebiet angepasst werden  
## Die Bergliste muss noch geparst werden da so nur die Koordinaten + Namen (gpx import) oder alles nur keine namen Cvs import
## verfügbar sind (warum auch immer) ALTERNATIVE reverse Geocoding
## Natürlich könnte man auch ein bisschen tricksen. 
## 1) konvertieren als CVS dann austausch der HTM Zeichen von hand oder mit kommandozeile
## 2) Import des ganzen als gpx alles wegschmeisen bis auf koordinaten und Namen
## 3) Merge von beiden
# download der datei
system("wget http://www.tourenwelt.info/commons/download/bergliste-komplett.kmz.php")
# umbenenen da gunzip nur .gz verarbeitet
system("mv /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.kmz.php /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.gz")
# entzippen da es sich um kmz handelt pgsbabel aber nur kml kann
system("gunzip -f /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.gz")
#anhängen der  kml endung
system("mv -f /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.kml")
# export als gpx datei
system("gpsbabel -i kml -f /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.kml -o gpx -F /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.gpx")
#export als csv datei
system("gpsbabel -i kml -f /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.kml -o unicsv -F /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/bergliste-komplett.csv")
}
 
 
#-------.
#   Gipfelerkennung
#-------.
# Variante A
# 1) Berechnung morhophometrischer Gipfel-Gridzellen *)
# 2) Gipfelgridzellen werden als Punktkkordinaten (Zentroid) mit ihrem Höhenwert extrahiert und 
# 3) als R-Tabelle verfügbar gemacht. 
#
# Varante B
# 1) Filterung der Daten (am einfachsten Mean oder Gauss)
# 2) Bestimmung aller lokalen Maxima
# 3) MAximum-Gridzellen werden als Punktkkordinaten mit ihrem Höhenwert extrahiert und
# 4) als R-Tabelle verfügbar gemacht.
#
# Variante C
# 1) Gipfelliste (z.B. http://www.tourenwelt.info) importieren und bereinigen (siehe Datenimport)
#   
#
# *) Wood, J.D. (1996) The geomorphological characterisation of digital elevation models, http://www.soi.city.ac.uk/~jwo/phd
#-------.
 
if (wood==T) { # Gipfelerkennung nach Wood (1996)
## Der Parameter s_tol bestimmt die Krümmungstoleranz der Gipfelfäche, c_tol der Krümmungstoleranz von Ebenen 
## size gibt die zur Berechnung verwendete Kernel-Größe (Moore-Nachbarschaft) in Pixel an
execGRASS("r.param.scale", flags=c("overwrite"), input="DGM", output="gipfel", s_tol=1.0, c_tol=0.00001, size=as.integer(3), param="feature")
 
# Export aus GRASS des Gipfelerkennung als GEOTIFF
execGRASS("r.out.gdal", input="gipfel", output="/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/peak_form.tif")
 
## Import GEOTIFF  Gipfelerkennung nach SAGA ACHTUNG Bennennung muss GRIDS und FILES heissen NICHT GRID/FILENAME
rsaga.geoprocessor("io_gdal",0,env=myenv,list(GRIDS="peak_form.sgrd" , FILES="peak_form.tif"))
 
 
## Wood bestimmt 6 Form-Klassen
## 1=planar,2=pit,3=channel,4=pass,5=ridge,6=peak
## Da uns nur 6=peak interessiert müssen wir das grid
## reklassifizieren. GRASS benötigt hierzu eine Textdatei
## mit Regeln falls die Textdatei "reclass" bereit existiert 
## wird sie mit dem Systembefehl "rm reclass" gelöscht
system("rm reclass")
 
## Erzeugen (linux shell) einer Lookuptabelle zur Reklassifizierung der morphometrischen Berechnung
## peak werden zu 1 alle anderen Werte zu 0
system("echo '6 = 1  
* = NULL 
                 end'   >>  reclass" )
 
## reklassifizierung unter Verwendung der erzeugten Datei reclass
execGRASS("r.reclass",flags=c("overwrite"), input="gipfel", output="gipfel_bool", rules="reclass")
 
## Da Gipfel nach Wood anders als lokale Maxima Gebiete ähnlicher Form sind
## können sie z.T. viele einzelene (zusammenhängende) Pixel umfassen
## Ein einfacher Weg aus diesen einen Punkt zu erzeugen ist die Umwandlung der 
## reklassifizierten Raster-Daten in Gipfel-Vektor-Polygone
execGRASS("r.to.vect", flags=c("overwrite"), input="gipfel_bool", output="gipfel_vpoly_bool", feature="area")
 
## Von diesen Poolygonen wird der Zentroid (Flächenschwerpunkt) berechnet und ausschreiben der centroid punkte der Gipfel-Polygone
## um den Gipfel als Punktkoordinate zu erhalten
execGRASS("v.type", flags=c("overwrite"),input="gipfel_vpoly_bool", output="gipfel_vpoint_bool", type="centroid,point")
 
## Zur Speicherung der Höhenwerte aus dem Geländemodell
## muss der erzeugten Vektordatei erst eine neue Attributspalte 
## in der Attributtabelle hinzugefügt werden 
execGRASS("v.db.addcol", map="gipfel_vpoint_bool", columns="heights double precision")
 
## Jetzt noch die Höhenwerte an den Centroid-Koordinaten 
## aus dem DGM Raster extrahieren und in die neue Spalte schreiben
execGRASS("v.what.rast", vector="gipfel_vpoint_bool", raster="DGM", column="heights")
gipfel=read.csv("gipfel.csv", header = FALSE, sep = "", dec=".")
# Löschen der Variable 3 (ID der Vektordatei)
gipfel$V3 <- NULL
 
} # Ende Wood -----------------------------------------------------------------
 
 
# Beginn MaxPix -----------------------------------------------------------------
if (max==T){ #maxPix
  # Gipfelerkennung "lokale Maxima" METHOD=0 ist der MEAN Filter
 
  ## Auschreiben der Reclass Tabelle für SAGA. 
  ## Achtung die DateiDARF NICHT beim Einlesen nach SAGA die Endung csv haben!
  ## AUSSERDEM MUSS sie 6 Ziffern nach dem Dezimalpunkt aufweisen (nsmall = 6)
#  hoehe =minGipfelhoehe
#  rec <- structure(list(MIN = c(0.0), MAX = c(hoehe), CODE = c(0.0)), .Names = c("MIN","MAX", "CODE"), row.names = c(NA,1L), class = "data.frame")
#  write.table(format(rec, nsmall = 6), "rec.txt",row.names = FALSE, quote = F, sep ="\t")
 
  ## Reklassifizieren des GRIDS
#  rsaga.geoprocessor("grid_tools",15,list(INPUT="DGM.sgrd",RESULT="DGMrec",METHOD=2,RETAB="rec.txt",ROPERATOR=1 ))
 
  rsaga.geoprocessor("grid_filter",0,list(INPUT="DGM.sgrd",RESULT="fDGM",MODE=0,METHOD=0,RADIUS=filterumgebung))
 
 
  # Shapes_Grid extrahiert lokale Maxima und Minima Positionen und ihre Höhenwerte (peaks und pits) in shapefiles
  rsaga.geoprocessor("shapes_grid",9,list(GRID="fDGM.sgrd", MAXIMA="localmax",MINIMA="localmin"))       
 
  # Das localmax shapefile enthält in der Spalte Z die Höhenwerte wird hier in eine CSV Textdatei exportiert
  rsaga.geoprocessor("io_shapes",2,list(SHAPES="localmax.shp",FIELD="Z",SEPARATE=0,FILENAME="gipfel.csv"))
 
  ## die csv datei der Gipfel (gleichgültig welche)  wird in die R-Tabelle gipfel geschrieben
  gipfel=read.csv("gipfel.csv", header = FALSE, sep = "")
  # es werden Spaltennamen zuweisen
  colnames(gipfel)=c("X","Y","Z")
 
} # Ende max--------------------------------------------------------------------
 
 
 
# Beginn Gipfelliste -----------------------------------------------------------------
if (gipfelliste==T){ 
# ermitteln des Extend damit ein polygon zum clippen erzeugt werden kann
rsaga.geoprocessor("grid_tools",15,list(INPUT="fDGM.sgrd",RESULT="extent", METHOD=0, OLD=9999.,NEW=100,SOPERATOR=1))
# erzeugen des Polygons
rsaga.geoprocessor("shapes_grid",10,list(SHAPES="umfang.shp", CELLS="0", PARAMETERS_GRID_SYSTEM_NX="1277",  PARAMETERS_GRID_SYSTEM_NY="1120", PARAMETERS_GRID_SYSTEM_X="577537.5",PARAMETERS_GRID_SYSTEM_Y="5211037.5",PARAMETERS_GRID_SYSTEM_D="75"))
 
#### einmal der durchgang 1) Einlesen 2) projizieren 3) clippen für die (von hand korrigierte) CSV liste
rsaga.geoprocessor("io_shapes",3,list(SHAPES="bergliste5.shp", X_FIELD=2,Y_FIELD=1,FILENAME="bergliste-komplett_16.05_2012.csv"))
rsaga.geoprocessor("pj_proj4",5,list(SOURCE_PROJ="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0",TARGET_PROJ="+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0",SOURCE="bergliste5.shp",TARGET="bergliste5_utm.shp"))
rsaga.geoprocessor("shapes_points",8,list(POINTS="bergliste5_utm.shp" ,POLYGONS="umfang.shp" ,CLIPS="bergliste5_cut.shp", METHOD="0"))
 
#### einmal der durchgang 1) Einlesen 2) projizieren 3) clippen für die gpx liste
rsaga.geoprocessor("io_shapes",13,list(SHAPES="bergliste.shp",FILE="bergliste-komplett.gpx"))
rsaga.geoprocessor("pj_proj4",5,list(SOURCE_PROJ="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0",TARGET_PROJ="+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0",SOURCE="bergliste.shp",TARGET="bergliste_utm.shp"))
rsaga.geoprocessor("shapes_points",8,list(POINTS="bergliste_utm.shp" ,POLYGONS="umfang.shp" ,CLIPS="bergliste_cut.shp", METHOD="0"))
 
### Todo Merge der beiden listen
 
# in eine CSV Textdatei exportiert
rsaga.geoprocessor("io_shapes",2,list(SHAPES="bergliste5_cut.shp",  ALL=T, SEPARATE=1,FILENAME="gipfel.txt"))
 
## die csv datei der Gipfel (gleichgültig welche)  wird in die R-Tabelle gipfel geschrieben
gipfel=read.csv("gipfel.csv", header = FALSE, sep = ",")
# es werden Spaltennamen zuweisen
colnames(gipfel)=c("X","Y","Z")
 
} # Ende Gipfelliste
 
 
# Weiter mit Dominanz -----------------------------------------------------------------
## Nun weiter mit der Verarbeitung
 
 
 
 
#-------.
#	Dominanz
#-------.
# 1) Für jede Gipfel-Koordionate wird mit "proximity" ein Distanz-Grid erstellt
#    (Resultat ist die Entfernung jedes Pixels zum Gipfel. 
# 2) Zusätzlich wird ein Maske erstellt (1 height>gipfel, 0 height <= Gipfel
# 3) Dominanz=minimum(proximity*maske) 
# 4) Dominaz wird als neue Variable in die Tabelle "gipfel" eingetragen 
#-------.
 
for (i in 1: nrow(gipfel))                                   				  # Für jeden Gipfel
{
 
  # Zeile vom Gipfel i als Text Tabelle speichern
  write.table(gipfel[i,],"tmp1.txt",row.names=FALSE)                     
 
  # Tabelle in SAGA importieren
  rsaga.geoprocessor("io_table",1,env=myenv, list(TABLE="peak",HEADLINE=TRUE,SEPARATOR="space", FILENAME="tmp1.txt"))
 
  # Tabellenimport in Shapedatei konvertieren
  rsaga.geoprocessor("shapes_points",0,env=myenv,list(POINTS="peak",TABLE="peak", X="X",Y="Y")) 
 
  # Erzeugte Shapedatei in eine Rasterdatei konvertieren ACHTUNG die Datei peak.sgrd muss bereits existieren!
  rsaga.geoprocessor("grid_gridding",0,list(INPUT="peak.shp", FIELD="X",LINE_TYPE=2,GRID_TYPE=3,TARGET=1,GRID_GRID="peak.sgrd"))    
 
 
  ## Höhe vom Gipfel aufrunden (simpler workaround damit Pixelwert sicher höher ist als seine Nachbarpixel) 
  ## und in Variable hoehe schreiben
  hoehe=ceiling(gipfel$Z[i])       	
 
  ## SAGA Single Reclass (Method = 0): Alle die <Z =0 (SOPERATOR=3 heißt <)
  #rsaga.geoprocessor("grid_tools",15,list(INPUT="fDGM.sgrd",RESULT="high_px", METHOD=2,SOPERATOR=1))
  ## SAGA Single Reclass (Method = 0): Alle Pixel die  >= Z = 1 (SOPERATOR=3 heißt >=)
  #rsaga.geoprocessor("grid_tools",15,list(INPUT="low_pix.sgrd", RESULT="high_pix",METHOD=0, OLD=hoehe,NEW=1,SOPERATOR=3))
 
  ## Auschreiben der Reclass Tabelle für SAGA. 
  ## Achtung die DateiDARF NICHT beim Einlesen nach SAGA die Endung csv haben!
  ## AUSSERDEM MUSS sie 6 Ziffern nach dem Dezimalpunkt aufweisen (nsmall = 6)
  rec <- structure(list(MIN = c(0.0, hoehe), MAX = c(hoehe, 9999.0), CODE = c(0.0,1.0)), .Names = c("MIN","MAX", "CODE"), row.names = c(NA,2L), class = "data.frame")
  write.table(format(rec, nsmall = 6), "rec.txt",row.names = FALSE, quote = F, sep ="\t")
 
  ## Reklassifizieren des GRIDS
  rsaga.geoprocessor("grid_tools",15,list(INPUT="fDGM.sgrd",RESULT="localmax",METHOD=2,RETAB="rec.txt",ROPERATOR=1 ))
 
  ## SAGA Proximity grid erstellen
  rsaga.geoprocessor("grid_tools",26,list(FEATURES="peak.sgrd",DISTANCE="proximity"))
 
  ## Entfernunggrid (0/1) mit Höhengrid  multiplizieren
  rsaga.geoprocessor("grid_calculus",1,list(GRIDS=list("proximity.sgrd;high_pix.sgrd"), RESULT="dominanz", FORMULA="a*b"))
 
  ## Ergebnis in R verfügbar machen
  dominanz_grd=raster("dominanz.sgrd", band=1,native=TRUE)  
  tmp=getValues(dominanz_grd)
  if (any(tmp>0)){                                        #geringste Distanz zu höheren =Dominanz
    gipfel$dominanz[i]=min(tmp[tmp>0])                      #wird in gipfel-tabelle geschrieben
 
 
  } else {
    gipfel$dominanz[i]=max(tmp[tmp>0])  	                   #Fall für höchsten Gipfel   
  }
 
 
}
 
 
#-------.
#	 SCHARTENHÖHE bzw. Prominenz   Variante B                                                                                                       
#-------.
#  1) Für jeden Gipfel wird die Schartenhöhe zu jedem höheren Gipfel ermittelt.                                           
#  2) Falls Gipfel1 und Gipfeln über einen Grat verbunden sind wird das Grat-Höhenprofil ermittelt                        
#  3) Die minimale Höhe des Profils ist die Schartenhöhe.                                                                 
#  4) Das Maximum aller ermittelten Schartenhöhen ist die gesuchte maximal Schartenhöhe                                   
#  5) Falls Gipfel und Zielgipfel NICHT über einen Grat miteinander verbunden sind wird ein Höhenprofil                   
#     aus der Sichtlinie zwischen den Gipfelpunkten erstellt dann weiter mit3) und 4)                                     
#-------.
 
##  Damit mit einer Flow-Analyse das Routing entlang der Grate Sattfinden kann
##  wird das gefilterte DGM invertiert
rsaga.geoprocessor("grid_calculus",1,list(GRIDS="fDGM.sgrd",RESULT="ifDGM",FORMULA="-a"))
 
# Hydrologische Korrektur - Sink Removal Output: pitlessifDGM
rsaga.geoprocessor("ta_preprocessor",2,list(DEM="ifDGM.sgrd",DEM_PREPROC="pitlessifDGM", METHOD=1))
 
# Flow Analysis. Output: direction
rsaga.geoprocessor("ta_channels",5,list(DEM="pitlessifDGM.sgrd",DIRECTION="direction",CONNECTION="tmp", SEGMENTS="segments",BASINS="basins",THRESHOLD=0))
 
# Erzeugen einer Matrix für die Schartenhöhen
schartenhoehe=matrix(ncol=nrow(gipfel),nrow=nrow(gipfel)) 
 
# Füllen der Matrix mit NA
schartenhoehe=schartenhoehe==NA                            
 
 
for (i in 1:nrow(schartenhoehe)) 				                   # für jeden Gipfel
{
 
  write.table(gipfel[i,],"tmp1.txt",row.names=F)    		   # Gipfelkoordinate die geroutet werden soll
  rsaga.geoprocessor("io_table",1, list(TABLE="DerGipfel",HEADLINE=TRUE,SEPARATOR="space",FILENAME="tmp1.txt"))
  rsaga.geoprocessor("shapes_points",0,list(POINTS="ZielGipfel",TABLE="DerGipfel",X="X",Y="Y"))
  if (i==1){schartenhoehe[i,1]=NA} 
  else     {if (i==2){
                      k=1
                      write.table((gipfel[k,]),"tmp2.txt",row.names=F)
                      rsaga.geoprocessor("io_table",1, list(TABLE="AndereGipfel",HEADLINE=TRUE,SEPARATOR="space",FILENAME="tmp2.txt"))
                      rsaga.geoprocessor("shapes_points",0,list(POINTS="Gipfelshp",TABLE="AndereGipfel",X="X",Y="Y"))
                      rsaga.geoprocessor("shapes_tools",2,list(OUT="merged",MAIN="ZielGipfel.shp",LAYERS="Gipfelshp.shp"))
                      rsaga.geoprocessor("grid_gridding",0,list(INPUT="merged.shp",FIELD="X",LINE_TYPE=2,GRID_TYPE=3,TARGET=1,GRID_GRID="tmp.sgrd"))
                      rsaga.geoprocessor("ta_channels",0,list(ELEVATION="pitlessifDGM.sgrd",SINKROUTE="direction.sgrd",CHNLNTWRK="tmp",CHNLROUTE="tmp",SHAPES="line",INIT_GRID="tmp.sgrd",INIT_METHOD=2,INIT_VALUE=0,DIV_CELLS=5,MINLEN=10))    
                      rsaga.geoprocessor("shapes_polygons",6,list(SHAPES="line.shp",POINTS="networkpoints"))
                      rsaga.geoprocessor("io_shapes",2,list(SHAPES="networkpoints.shp",FIELD="Z",SEPARATE=0,FILENAME="Networkpoints.txt"))
                      netpoints=read.table("Networkpoints.txt",header=F,sep="\t")   	
                      if (substr(netpoints[1,3],1,2)==substr(netpoints[nrow(netpoints),3],1,2)) 
                         { # es darf nur eine durchgehende linie sein!!
                           rsaga.geoprocessor("ta_profiles",4,list(DEM="DGM.sgrd",LINES="line.shp",PROFILE="profile")) 
                           rsaga.geoprocessor("io_shapes",2,list(SHAPES="profile.shp",FIELD="Z",SEPARATE=0,FILENAME="Elevprofile.txt"))
                           elevprofile= read.table("Elevprofile.txt",header=F,sep="\t")   # lade die nach R
                           schartenhoehe[i,k]=min (elevprofile[,3]) 	
                         }					
                          else  # wenn berg nicht durch grat verbunden
                              {   # linie dazwischen
                                rsaga.geoprocessor("shapes_lines",1,list(LINES="line",POINTS="merged.shp",ORDER="Y"))
                                rsaga.geoprocessor("ta_profiles",4,list(DEM="DGM.sgrd",LINES="line.shp",PROFILE="profile"))
                                rsaga.geoprocessor("io_shapes",2,list(SHAPES="profile.shp",FIELD="Z",SEPARATE=0,FILENAME="Elevprofile.txt"))
                                elevprofile= read.table("Elevprofile.txt",header=F,sep="\t")		   		  # lade die nach R
                                schartenhoehe[i,k]=min (elevprofile[,3]) 					          # min Schartenhöhen zu allen anderen Gipfeln
                              }
                     }
                      else
                          {    
                           for (k in 1:ncol(schartenhoehe[,1:(i-1)]))						  # für jeden anderen Berg (nicht doppelt berechnen!)
                               {
                                write.table((gipfel[k,]),"tmp2.txt",row.names=F)
                                rsaga.geoprocessor("io_table",1, list(TABLE="AndereGipfel",HEADLINE=TRUE,SEPARATOR="space",FILENAME="tmp2.txt"))
                                rsaga.geoprocessor("shapes_points",0,list(POINTS="Gipfelshp",TABLE="AndereGipfel",X="X",Y="Y")) 				 
                                rsaga.geoprocessor("shapes_tools",2,list(OUT="merged",MAIN="ZielGipfel.shp",LAYERS="Gipfelshp.shp"))  		     
                                rsaga.geoprocessor("grid_gridding",0,list(INPUT="merged.shp",FIELD="X",LINE_TYPE=2,GRID_TYPE=3,TARGET=1,GRID_GRID="tmp.sgrd"))  
                                rsaga.geoprocessor("ta_channels",0,list(ELEVATION="pitlessifDGM.sgrd",SINKROUTE="direction.sgrd",CHNLNTWRK="tmp",CHNLROUTE="tmp",SHAPES="line",INIT_GRID="tmp.sgrd",INIT_METHOD=2,INIT_VALUE=0,DIV_CELLS=5,MINLEN=10))      
                                rsaga.geoprocessor("shapes_polygons",6,list(SHAPES="line.shp",POINTS="networkpoints")) 
                                rsaga.geoprocessor("io_shapes",2,list(SHAPES="networkpoints.shp",FIELD="Z",SEPARATE=0,FILENAME="Networkpoints.txt"))
                                netpoints=read.table("Networkpoints.txt",header=F,sep="\t")  
                                if (substr(netpoints[1,3],1,2)==substr(netpoints[nrow(netpoints),3],1,2))        
                                   {        
                                    rsaga.geoprocessor("ta_profiles",4,list(DEM="DGM.sgrd",LINES="line.shp",PROFILE="profile"))
                                    rsaga.geoprocessor("io_shapes",2,list(SHAPES="profile.shp",FIELD="Z",SEPARATE=0,FILENAME="Elevprofile.txt"))
                                    elevprofile= read.table("Elevprofile.txt",header=F,sep="\t")   											
                                    schartenhoehe[i,k]=min (elevprofile[,3]) 	
                                    cat("Gipfel",i,"von",nrow(gipfel),"wird bearbeitet(",100*(k/ncol(schartenhoehe[,1:(i-1)])),"%). Methode: Grat")
                                   }					
                                 else
                                    {
                                     rsaga.geoprocessor("shapes_lines",1,list(LINES="line",POINTS="merged.shp",ORDER="Y")) 					         
                                     rsaga.geoprocessor("ta_profiles",4,list(DEM="DGM.sgrd",LINES="line.shp",PROFILE="profile")) 		         
                                     rsaga.geoprocessor("io_shapes",2,list(SHAPES="profile.shp",FIELD="Z",SEPARATE=0,FILENAME="Elevprofile.txt"))     
                                    elevprofile= read.table("Elevprofile.txt",header=F,sep="\t")   												     
                                    schartenhoehe[i,k]=min (elevprofile[,3]) 
                                    cat("Gipfel",i,"von",nrow(gipfel),"wird bearbeitet(",100*(k/ncol(schartenhoehe[,1:(i-1)])),"%). Methode: Sichtlinie")
                                    }
                              } 
                          }
       }
  cat("Gipfel",i,"von",nrow(gipfel),"fertig!")
  }
 
tmp=t(schartenhoehe)
schartenhoehe[is.na(schartenhoehe)]=tmp[is.na(schartenhoehe)]
write.table(schartenhoehe,"MatrixSchartenhoehe.txt",row.names=F)                               # damit man es nicht immer neu rechnen muss ;)
 
#-------.
# Prominenz
#-------.
# Abfrage auf jeden Gipfel zu welchem der höheren Gipfel die Schartenhöhe am höchsten ist.
# Diese wird dann von der höhe des Gipfel subtrahiert und entspricht damit der Prominenz des Gipfels.
#-------.
for (i in 1:nrow(gipfel))
{
  hoehere=which(gipfel$Z>gipfel$Z[i])
  ifelse ((length(hoehere)==0), (gipfel$prominenz[i]=gipfel$Z[i]-min(schartenhoehe[i,],na.rm=T)),   
          (gipfel$prominenz[i]=gipfel$Z[i]-max(schartenhoehe[i,hoehere])))                   #noch nicht optimal: dem höchsten Berg kann man hier eigentlich keine prominenz zuweisen
}
gipfel$prominenz[gipfel$prominenz<0]=0.000001                                                #für den Fall das die Prominenz negativ ist weil scharte größer als Berg...??
 
 
#-------.
#	Eigenständigkeitswert
#-------.
gipfel$E[gipfel$dominanz<100000]=-1*((logb((gipfel$Z/8848),base=2)+logb((gipfel$dominanz/100000),base=2)+logb((gipfel$p/gipfel$Z),base=2))/3)
gipfel$E[gipfel$dominanz>=100000]=-1*((logb((gipfel$Z/8848),base=2)+logb((gipfel$p/gipfel$Z),base=2))/3)
write.table(gipfel,"gipfel.csv",row.names=F,sep=";")
 
#-------.
#								 Finde "perfekten" Gipfel
#-------.
Ziel=gipfel
Ziel=Ziel[order(Ziel$E,decreasing=F),]
write.table(Ziel,"ziel.csv",row.names=F,sep=";")
 
rsaga.geoprocessor("io_gdal",1,list(
  GRIDS="DGM.sgrd",
  FILE="DGM.rst",
  FORMAT=14,
  TYPE=0))                                                                     # export als rst (warum eigentlich?)
demS=raster("DGM.rst")                                                                # lade in R
cols=c("red","orange","yellow","green","blue")
pdf("/home/creu/workspace/peak/V1/Top5.pdf")
plot(demS,ylab="Y",xlab="X",main="Die 5 'besten' Gipfel",col= terrain.colors(20),zlim=c(0,500),ylim=c(5620000,5637500)) # plot Top 5
points(Ziel[1:5,]$X,Ziel[1:5,]$Y,pch=16,col=cols,cex=2)
par(xpd=T)
text(15.2,50,"Elevation (m)",cex=0.8)
legend("topright",pch=16,col=cols,legend=c("1","2","3","4","5"),bty="n",horiz=T)
dev.off()
 
#-------.  
#								 Dies&Das
#-------.
 
 
# --------- Exkurs Import Export Plotten ----------------------
## export der erzeugten Vektordatei 
 
## Und weils so schön ist hier der Export des  Punkt-Shapefiles g.shp als ASCII XYZ File (Option 2) in die Datei gipfel.csv
#rsaga.geoprocessor("io_shapes",2,env=myenv,list(SHAPES="g.shp",FIELD="heights",SEPARATE=0, FILENAME="gipfel.csv"))
 
## falls die auszuschreibende SHapedatei g.* bereit existiert wird sie mit dem Systembefehl gelöscht
##system("rm /home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/g.*")
 
## aus GRASS ins SHAPEformat
##execGRASS("v.out.ogr",  input="gipfel_vpoint_bool", type="point", dsn="/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/g.shp")
## aus GRASS ins ASCII XYZ Format
###execGRASS("v.out.ascii", input="gipfel_vpoint_bool", output="/home/creu/workspace/grassprojects/GIPFEL/SRTM/RSAGA/gipfel.csv", columns="heights", format="point",fs=" ")
 
# Import SHAPE Punkte in R
##xx <- readShapePoints("/home/creu/workspace/peak/V1/g.shp")
## dimensions(xx)
##summary(xx)
 
## plotten der Punkte
##plot(xx,col="red", add=TRUE)
 
##schreibt man die Ausgabebefehle (plot, image etc)
## zwischen die Anweisungen png() dev.off() wird alles in die Grafikdatei geleitet
##png("gipfel.png")
##plot(xx,col="red", add=TRUE)
##dev.off()
# --------- Exkurs Import Export Plotten ----------------------
 
 
###Plot Zwischenschritte
rsaga.geoprocessor("io_gdal",1,list(
  GRIDS="DGM.sgrd",
  FILE="DGM.rst",
  FORMAT=14,
  TYPE=0))#import demsmooth nach rst
dem=raster("DGM.rst")                                                                       # lade in R
pdf("/home/creu/workspace/peak/V1/DEMS.pdf",width=10,height=5)
par(mfrow=c(1,2),mar=c(4.2,4.2,1,5.2))
plot(dem,ylab="Y",xlab="X",main="Raw DEM",col= terrain.colors(20),zlim=c(0,500))                # plot Top 10
plot(demS,ylab="Y",xlab="X",main="DEM und Lage der Gipfel",col= terrain.colors(20),zlim=c(0,500))
points(gipfel$X,gipfel$Y,pch=16,cex=1)
par(xpd=T)
text(15.2,50,"Elevation (m)",cex=0.8)
dev.off()
 
###Hilfe 
#rsaga.get.modules(); z.B. rsaga.get.usage("io_grid",6)
#-------.
#
#-------.
rsaga.geoprocessor("io_shapes",3,env=myenv,list(SHAPES="GipfelShape",X_FIELD=1,Y_FIELD=2,FILENAME="gipfel.xyz"))

Hanna MeyerChris Reudenbach

projekte/misc/perfectpeak.txt · Last modified: 2018/12/23 19:46 (external edit)