Mit Rasterkoordinaten aus Bild A DGM-Höhen aus anderem Raster auslesen

Oha, da war ich nicht deutlich in meinem Geschreibsel.

Grundsätzlich habe ich mir, ausgehend von einem DGM, einige Positionen im Gelände in Form bestimmter Pixel herausgearbeitet. Sie stellen eine Art kleine Sprünge bzw "Hopser" im Gelände dar. Von diesen bestimmten Pixeln/Geländepositionen möchte ich nun nicht nur wissen, wo sie liegen, sondern auch, wie sich die Höhe in diesen Punkten ändert.

Ich habe hier also einerseits ein waschechtes DGM mit Landeskoordinaten und Höhen in Metern, wie es sich gehört.

Und ich habe durch viel Rumgebastel ein Rasterbild erzeugt, das in der Lage identisch ist mit dem DGM. In diesem Rastebild befinden sich nur noch einige wenige Pixelhaufen. Diese symbolisieren die für mich interessanten Positionen im Gelände. Sie haben im Laufe meiner Rumbastelei den Pixelwert/Grauwert 254 zugewiesen bekommen.
Der Rest der Pixel im Rasterbild hat den Grauwert 0. Sonst sind keine anderen Grauwerte mehr vorhanden. Ich habe also sozusagen alle für mich uninteresanten Geländepunkte eliminiert.

Ich möchte nun wissen, wie sich die Höhe an Position meiner "Wert 254-Pixel" ändert. Ich dachte an etwas in der Art "Lese Höheninformation des "Wert 254-Pixels" im DGM sowie die Höheninformation des nördlich und südlich gelegenen Pixels im DGM; und gibt mir die Höhenänderung zwischen dem "Wert 254-Pixel" und dem nördlich davon liegenden sowie dem südlich davon liegenden Pixel aus.
Mir würde jedoch tatsächlich die Höhenänderung in Nord-Süd ODER OST-West-Richtung genügen, und eigentlich bräuchte ich sogar nur die stärktste der vier Höhenänderungen. Um sozusagen die stärkste Ändertung der Höhe in meinen Punkten zu erhalten, also wie stark der "hopser" im Gelände in meinen Punkten höchstens ist.

Also so eine Art kleines Höhenprofil in zwei Himmelsrichtungen an meinen "Wert 254-Pixeln" im DGM. Wobei das ja eigentlich schon mehr ist, als ich brauche.

In der ArcToolbox gibt es unter "3dAnalystTools - RasterSurface" das Werkzeug "Slope", das, wenn ich es richtig verstehe, die maximale Höhenänderung / Steigung an einer Rasterzelle zu den umliegenden Rasterzellen bestimmt (für ein komplettes Bild).
Ich bräuchte nun so etwas ähnliches, leider befinden sich im meinen Rasterbild nur noch Pixelwerte 0 sowie 254. Die Höheninformation, die ich benötige, befindet sich im ursprünglichen DGM.
(Ansonsten hätte ich versucht, die Slope-Funktion auf mein Rasterbild anzuwenden und irgendwie damit weiterzubasteln)

Wäre es alternativ eine Lösung, die Slope-Funktion auf das ganze DGM anzuwenden und irgendwie nur von den "Wert 254-Pixeln" in meinem Rasterbild die durch "Slope" entstehenden Informationen bezüglich der maximalen Höhenänderung auszulesen?

Ciao
Andi

Vielen Dank für die ausführliche Beschreibung. Ob ich's begriffen habe, merke ich aber erst, wenn ich es durchprobiert habe. ... ;)

Zunächst mal habe ich jetzt unter den "SpatialAnalystTools - Zonal" die ZonalStatistics ausprobiert. Toll, die lassen sich auch aus einem vba-Skript aus aufrufen! Als "Als InputRaster or feature zone Data" das Shape mit den Polygonen (aus werkzeug RasterToPolygon_conversion) angegeben, als "Zone field" den GRIDCODE, das müssten ja die Pixelwerte sein. Das ""InputVAlueRaster" ist mein DGM, der "StatisticsType" ist Range und "Ignore NoData in Calculations" ist gesetzt.

Dabei kommt ein nettes Rasterbild heraus.
Allerdings habe ich das mit "Maske" in den "Umgebungen noch nicht so ganz begriffen.

So wie ich das verstehe, berechnet mir das Tool auf diese Art innerhalb der Polygone den stärksten Höhenunterschied auf Basis der DGM-Daten. Liege ich da richtig?
Die Polygone enthalten jedoch nur die für mich relevanten Pixel, nicht jedoch die daneben liegenden. Muß ich da nicht noch zuvor die Polygone um ein Pixel "vergrößern", damit tatsächlich der Höhenunterschied zu den Nachbarpixeln bestimmt wird?
Liege ich da richtig?

Jetzt versuche ich mal die Methode "Focale statistics"

Ciao
Andi

Wollte zuerst mal im Polygon-Shape die Polygon in alle Richtungen ein Pixel "dicker" machen.
Irgendwie stelle ich mich dusselig an....
Da bin ich doch bei "AnalysisTools-Proximity -> Buffers" richtig, oder?
Als InputFeature nehme ich mein Polygon-Shape, Sied Type ist Full, DissolveType ist NONE.
Als linear unit nehme ich..... grübel.
nehme ich 25 Meter, was der größe meiner Pixel entspricht, dann habe ich nach der Operation ..... fast nix mehr. Statt um ein Pixel vergrößerter Polygone sind fast alle Polygone verschwunden.
Auch bei anderen EInheiten wird das eher eine Erosion anstelle einer Vergrößerung.
Kann es sein, daß nicht meine Polygone (GRIDCODE 254) vergrößert werden, sondern die Bereiche mit "GRIDCODE 0" ?
Na dann nehme ich doch einfach mal einen negativen Buffer-Radiums.
So gehts. Sorry, aber diese Lösung ist mir erst beim eintippen dieser zeilen eingefallen.....

So, jetz habe ich mein Raster genommen, habe ein Polygon-Shape erzeugt.
Dann habe ich mit "DataManagement -Features - multipart to Singlepart" mein Polygon in einzelne Polygone aufgelöst.
Anschließend mit "AnalysTools - Proximity - Buffers" die Polygone um eine Zelle rundherum vergrößert (hoffe ich).
Anschließend mit "SpatialAnalystTools - Zonal - zonalStatistics" und dem Range-Operator dies berechnung druchgeführt.

Das Ergebnis... ist etwas seltsam?
Wo ich in meinem Rasterbild zuvor "pixelwert 0" hatte, also in der weiten Fläche des Bildes, die sozusagen das "Außenbpolygon" mit "GRIDCODE 0" bildete, habe ich nun flächendeckend "Pixel value 155,221786".

Wo zuvor meine Pixelhaufen mit Grauwert 254 waren, also innerhalb meiner Polygone mit GRIDcode 254.....
Auf dem Rand der vergrößerten Polygone, sozusagen auf der Außenkante des Buffers, habe ich jetzt "NoData". Im inneren der Polygone habe ich ebenfalls "NoData". Mit Außnahme solcher Bereiche, die drei Pixel breit sind, so daß ein Pixel ohne Kontakt zum Buffer war, dort habe ich PixelValue 105,088440. Und außer diesen drei Werten (No data, 155m, 105 Meter) gibt es keine weiteren Werte......

Ja, muß tatsächlich erst die Bereiche mit Gridcode 0 loswerden!
Ist mir zu spät gekommen - völlig verpennt. Sorry.
Kann ja auch irgendwie nicht gehen, denn er hat ja die uninteressanten Bereiche ebenfalls als Polygon vorliegen und rechnet daher darauf rum.

Ist einfach zu warm hier....
Werde erst mal diese uninteressanten Bereiche eliminieren.

Bei "Input raster or feature zone data" habe ich mein Polygonshape eingetragen
Bei "Input value Raster" das ursprüngliche DGM.
Bei "ZoneField" ist GRIDCODE gewählt.
"Statistics Type" ist "Range".

Falls ich heute nicht völlig daneben bin.

Danke nochmal für den Hinweis. Ist nicht mein Tag.

Ciao
Andi

Das mit nach Hause gehen habe ich auch gemacht. Die Idee war gut ;)

Zu den Zonen: daß in jedem Polygon der selbe Gridcode drinsteht, stimmt eigentlich. Ich dachte, es ginge nur darum, das "Feld" mit dem Kriterium festzulegen und das Kriterium wäre "es steht ein Wert im Feld (in Gridcode) oder halt nicht. Denn ich habe nicht gesehen, daß man iregndwo den Inhhaltswert von Gridcode in dem Benutzerdialog angeben könnte.

Habe im Rasterbild (sorry, im Polygonshape aus dem Rasterbild) insgesamt ZWEI verschiedenen Gridcodes: Wert 254 für die interessanten Rasterzellen-Bereiche/Polygone und GRIDCODE 0 für den Rest. Und wenn ich die Gridcode-Bereiche mit Wert 0 auch noch zu "NoData" eliminiert habe (Dein Post vom 1.Sep.2009 10:07 ), dann habe ich ja nur noch die Polygonbereiche (mit Gridcode wahlweise 254 oder reklassifiziert 1).
Aber offensichtlich lag ich völlig wieder falsch...

Arghl .... Nur stehe ich jetzt vor dem Problem, jedem Einzelpolygon einen "eindeutigen Bezeichner" zu verpassen.
Wie geht denn das schon wieder ?
In meinen Handbüchern habe ich das bis jetzt noch nicht gefunden...

Irgendwie wird die Aufgabe immer komplizierter. Dabei wollte ich doch von der Idee her "einfach nur" Höheninformtionen aus dem DGM für Landeskoordinatenstellen auslesen, die in meinem Rasterfile festgelegt sind....(und die der nördlich und südlich gelegenen Rasterzelle)

Das ist jetzt schon umfangreicher und aufwändiger als die eigentliche Bestimmung der interessanten Rasterzellen, dabei dachte ich, es wäre eine einfache Standardaufgabe für ArcMap.

Ciao
ein immer verwirrterter Andi ...

Danke, da muß ich mal drüber nachgübeln.

Aber mal anders gefragt... nur so als Idee:

Kann man nicht mit einem ArcTool in einem Rasterbild an der Stelle meiner Wert 245-Pixel (oder reklassifiziert Wert1-Pixel, Rest ist dann NoData) die Landeskoordinaten in eine Tabelle auslesen, dann mit einem anderen ArcTool an diesen Landeskoordinaten aus dem DGM die Höhenwerte auslesen sowie die Höhenwerte der Rasterzelle nördlich und südlich davon und diese in eine Tabelle oder ein neues Rasterbild schreiben?
In so einem Rasterbild bzw. Tabelle könnte ich dann per Eigenbau-Filter die Höhenänderung berechnen.
Ist das vielleicht eine einfachere Methode?

(was ich inzwischen gelernt habe, ist, wie ich im Rasterbild mittels durchlaufender Zeilen-/Spaltennummer auf die Pixelwerte zugreifen und diese Manipulieren kann = eine Art Eigenbau-Pixelwert-Filter. Aber keine Ahnung, wie ich damit auf die Landeskoordinaten der Rasterzelle zugreifen könnte. Außerdem muß ich mal nachschauen, ob jede Zelle eine eigene Landeskoordinate hat, oder ob die nicht sich über die Zellen-Fläche hinweg verändern. Immerhin gibt mir das Identify-Werkzsug drei Nachkommastellen aus, dabei habe ich eine horizontale Rasterweite von 20 Metern.)

Ciao
Andi

Das mit dem Speicherbedarf ist eigentlich nicht unbedintgt ein Problem, da es um nur wenige hundert Punkte geht. "Andere Einschränkungen" düften eher Probleme machen.
Leider bekomme ich es ja aber im Rasterformat auch nicht so einfach hin. Wenn es denn irgendwann einmal funktioniert, sind das seehr viele abzuarbeitende Einzelschritte.

Dabei arbeite ich jetzt ja noch an der ArcMap-Oberfläche mit den Werkzeugen. Das ist ja alles auch noch nicht in VBA/ArcObjects umgesetzt.

Auf diesen zweiten Weg bin ich gekommen durch den Thread "mit Hilfe eines Raster aus vorgegebenen Koordinaten die Geländehöhe ermitteln" von Nils Landmeyer etwas weiter unten in der Übersicht. Er versucht wohl für eine Punkte-Tabelle mit xy-Koordinaten die Geländehöhe für die Punkte aus einem Rasterfile zu ermitteln.
Dazu erzeugt er ein Layer aus seiner Tabelle (DataManagementTools-MakeXYEventLayer), wandelt diese zu einem shapefile um (ConversionTools - FeatureClasstoShape) und holt sich die Höhe der Punkte in Form der Zellwerte aus einer Rasterdatei (SpatialAnalystTools - Extraction - Extract Value to Point). Damit will er eine Tabelle mit den Punktkoordinaten und der zugehörigen Höheninformation erhalten. Allerdings bittet er hier im Forum um Hilfe, weil es wohl nicht so ganz funktioniert.

Ciao
Andi

Es muß doch ohne solche enormen Verrenkungen möglich sein, für Punkte mit vorhandenen XY-Landeskoordinaten die korrespondierenden Z-Werte/Höheninformation aus einen anderen Rasterbild/DGM auszulesen.

Mal schrittweise und ohne das mit dem "nördliches und südliches Pixel"-Zeugs
Ich kann ja mein Rasterfile auch z.B. mit "ConversionTools - Raster to Point" in ein Shape /PointFeature umwandeln. das gibt zwar eine "dicke"Datei, aber das wäre zu handeln....

Kann man über eine Verschneidung der beiden Rasterbilder, also dem DGM und dem Rasterbild mit den "interessanten Pixeln" mit Wert 1 (der ganze Rest jetzt mal per Reklassifizierung zu "nicht existent/DoData" gewandelt) nicht ein neues Rasterbild erzeugen?
Also sozusagen, immer wenn im Rasterbild Wert = 1 ist, soll im neuen Rasterbild das entsprechende Pixel aus dem DGM übernommen werden. Ansonsten (Wert = NoData) soll nix übernommen werden?

Oder aus dem DGM genau die Pixel ausschneiden, die im Reklassifizierten Raster wert 1 haben (interessante Werte=1, Rest=0) und daraus ein neues Raster erzeugen?

He, bin gerade über folgendes gestoplert:

SpatialAnalystTools - Extraction - Extract by mask.

Wenn ich hier als InputRaster das DGM angebe und als Maske mein reklassifiziertes Raster (in dem alle uninteressanten Bereiche = NoDate und alle interessanten Pixel = 1 sind)....

Habe ich damit ein neues Bild erzeugt, in dem an den Stellen, die im Raster = 1 sind, die DGM-Höhen zugeordnet sind?

Ciao
Andi

Ups, Entschuldigung Hannes17.
Dachte, ich zerlege meine Probleme in Themenbereiche. Die Grundlage ist zwar insgesamt natürlich mein Rumgebastel als solches, aber es fallen halt jede Menge Teilprobleme an, wie "ArcObject hillshade_3d bekomme ich nicht zum laufen" oder "wie zerlege ich eine Rasterdatei in Einzelgebiete" oder wie hier "Mit Rasterkoordinaten aus Bild A DGM-Höhen aus anderem Raster auslesen".
Aber schon klar, es ergeben sich auf diese Art eine ganze Menge einzelner Threads, insbesondere wenn ein Anfänger wie ich relativ planlos rumbastelt.

Werde versuchen, das in Zukunft etwas einzuschränken und besser zu sammeln.

Meinst Du mit "Convert Feature to z-feature" das Tool
3DAnalystTools - Conversion - FromFeatureClass - feature Class z to ASCII ?

Ciao
Andi

Hi Rena,
eigentlich hätte das ja gar nicht so übermäßig umfangreich werden sollen. Was ich da bastle, stellt eigentlich nur die Aneinanderreihung von ein paar ArcObjects-Aufrufen rund um ein eigenes Filtersystem zur Auffindung meiner "interessanten Punkte"dar. Der Bau dieser Filteroperationen ist eigentlich das, worum es mir geht.

Was ich "drumherum" machen will, ist mir eigentlich relativ klar, ich hänge nur immer wieder daran, wie ich das konkret mit ArcMap mache und wie die ArcObjects tatsächlich funktionieren. Das habe ich mir eigentlich mehr als "einfache, per einfachem Werkzeugaufruf durchführbare Standardaufgabe" vorgestellt.

Die Programmiererei in VBA ist eigentlich gar kein Problem, mehr als eine Handvoll grundlegender Befehle (Schleifen, Variablendeklaration usw.) braucht's irgendwie gar nicht.
Problematischer ist da schon eher, wenn die ESRI-Beispiele für den Aufruf der ArcObjects Fehler enthalten oder gar nicht nicht für vba beschrieben sind, so daß man dann doch rumbasteln muß. Wenn man sich dann noch dusselig anstellt und murks eintipt....
Und natürlich muß ich Anfänger die passenden Werkzeuge auch erst mal finden.

Öhm... ein ModelBuilder-Modell kann man als vba-Skript exportieren? Dann brauche ich ein neues ArcMap Handbuch. Hatte eines aus der Bücherei, da stand nichts von Export-Möglichkeiten.


Hm, wenn ich nun ein Rasterbild erzeugt habe, das meine interessanten Punkte samt zugrhöriger DGM-Höhen enthält, dann muß ich mich jetzt darum kümmern, wie ich an das "nördlich und südlich gelegene Pixel" rankomme.

Als Ausgangs-Status habe ich jetzt also (so hoffe ich) mein DGM sowie ein Rasterbild mit den interssanten Pixeln, die mit den zugehörigen DGM-Höhen versehen sind.
Mal sehen, ob ich etwas finde, das mir beispielsweise den Höhenunterschied zwischen (den Pixeln in) meinem neuen Rasterbild und den "Umgebungspixeln" im DGM bestimmt.

Bin natürlich für Tips dankbar.

Ciao
Andi

Hi,
was haltet ihr von folgendem Lösungsweg:

1) das Ursprungs-DGM mit "ArcToolbox - SpartialAnalystTools - Neighborhood - FocalStatistics" unter Einsatz des Range-Operators bearbeiten. (=> Datei focal_aus_dgm.img)
Damit müßte doch ein neues Rasterbild entstehen, welches in den Rasterzellen den jeweils größten Zellwertunterschied (= Höhenunterschied beim DGM ) zu den Nachbarzellen in der Maskengeometrie stehen hat. Wenn man statt der 3x3- eine 1x3-Umgebung verwendet, sollte in den neuen Rasterzellen dann der größte Zellwertunterschied (= Höhenunterschied) zwischen der betrachteten Zelle und dem "oberhalb" bzw. "unterhalb" liegenden Pixels stehen. Oder?

2) Das Rasterbild mit meinen "interessanten Pixeln" mittels
"ArcToolbox - SpatialAnalystTools - Reclass - Reclassify"
reklassifizieren, so daß meine "interessanten Pixel" den Wert 0 bekommen und alles andere zu NoData wird.

3) Nun mit "ArcToolbox - SpatialAnalystTools - Extraction - Extract by mask" diese beiden Zwischenergebnisse verschneiden: Mittels dem reklassifiziretn Rasterbild aus focal_aus_dgm.img ein neues Raster ausschneiden.

Dann sollte doch ein Rasterfile herauskommen, welches an den "interessanten Stellen" den größeren der beiden Höhendifferenzen zwischen betrachtetem Pixel/"interessanter Stelle" sowie dem oberhalb bzw. unterhalb liegenden Pixel stehen hat. Und alles andere ist NoData?

Wäre das eine Lösung oder ist das falsch gedacht?

Und noch zwei zusätzliche Fragen.
- ist die 1x3-geometrie tatsächelich die "Nord-Süd-Richtung" und 3x1 die "ost-west-Richtung" (meine Ursprungsbilder sind sozusagen genordet).
- wenn ich statt dessen eine 1x5-Geometrie nehme, was wird dann an die neue Rasterzelle geschrieben? Ich nehme an, der höchste vorkommende Wertunterschied zwischen dem "zentralen Pixel" und den anderen 4 Pixeln?

Wäre klasse, wenn jemand diesen Gedankengang nachprüfen könnte.
Denn so langsam sehe ich vor ArcTools keine Bäume mehr .... oder so ähnlich.

Ciao
Andi