Project xyz Data im Batch Modus

Hallo Forenmitglieder,

ich hab hier einen ganzen Haufen an xyz-Daten im *.txt Format.
Die möchte ich jetzt erstmal als Event-Layer laden, dann als Shapefile exportieren und schließlich im Rasterformat abspeichern und dann ganz am Schluss zu einem zusammenhängenen DGM mosaikieren.

Aber: Geht das irgendwie im Batch Modus, sonst muss ich Weihnachten noch dran rumbasteln....

Hab mal versucht die Funktion in der ArcToolbox zu finden, war aber bislang erfolglos. Python kann ich nicht, also hab ich auf den Modelbuilder geschaut, aber da muss ich ja auch erstmal wissen, wo man die Funktion finden kann.

Vielen Dank für Hilfe
Beste Grüße
Jochen
Hallo Jochen,

ich hatte vor Jahren das gleiche Problem mit über 400 xyz-Dateien im txt-Format. Eine einzelne Datei (2km*2km im 1m-Raster) als Event-Layer einzuladen und daraus ein Grid zu erzeugen, hat etwa 5 Minuten gedauert. Aufgrund der vielen Dateien sah ich damals nur die Lösung, mittels Python-Skript die über 400 Einzel-Grids zu erzeugen. Diese Einzel-Grids zu einem zusammenzufügen, das ging dann schnell per Hand im ArcGIS.

Mit dem Python-Skript habe ich immer einzelne Pakete rechnen lassen (tagsüber wenn ein Rechner frei war – nachts ging nicht). Die Gesamtrechenzeit bei mir lag so bei 35 Stunden. Vielleicht schaffst Du es doch, Python dafür zu nutzen. Bei mir setzten sich die Dateinamen aus dem verkürzten linken Rechtswert und dem verkürzten unteren Hochwert zusammen, also z.B. „3624.txt".

Mein Python-Skript lautete wie folgt:

# Create the Geoprocessor object
import arcgisscripting, sys, string, os
gp = arcgisscripting.create()

# overwrite the output if it already exists
gp.OverWriteOutput = 1

try:
for i in range (36, 42, 2): # Rechtswert, durchlaeuft die Zahlen ...
for j in range (20, 22, 2): # Hochwert, durchlaeuft die Zahlen ...

# Set local variables
xmin = i
ymin = j

Kachel= str(xmin) + str(ymin)
Textfile = "D:/DGM_2012/Lieferungen/update_2011_06_30/TXT/" + str(Kachel) +".txt"

Ausdehnung = "45" + str(xmin) + "000 57" + str(ymin) + "000 45" + str(xmin+2) + "000 57" + str(ymin+2) + "000"
print Ausdehnung

# Make the XY event ... input is a table view, result is a temporary point feature layer
gp.MakeXYEventLayer(Textfile, "x", "y", "Punktwolke")

# vielleicht erforderlich, vielleicht auch nicht
gp.Extent = Ausdehnung

# FeatureToRaster_conversion ... input is a point feature class, output without extension is a grid
gp.FeatureToRaster_conversion("Punktwolke", "z", "D:/DGM_2012/Lieferungen/update_2011_06_30/Test/"+ Kachel, "1")

print "done: ", Kachel

except:
# Print error message if an error occurs
print gp.GetMessages()


Vielleicht gibt es ja hier jemanden, der das Skript noch optimieren kann. Meine Python-Kenntnisse sind immer noch rudimentär.
Gigi
Hallo Gigi,

vielen Dank für deine ausführliche Antwort. Ich hatte schon gar nicht mehr damit gerechnet und mittlerweile auch eine Lösung gefunden, auch wenn die etwas umständlich ist.

Ich beutze dafür den Model Builder und lasse ihn mit der Funktion "Iterate" über den Ordner mit den Original *.xyz Daten laufen. Dabei werden folgende Funktionen ausgeführt:

ASCII 3d to Feature (zu finden unter: 3D Analyst Tools --> Conversion --> From File --> ASCII 3d to feature)

dann

Add xy-coordinates (zu finden unter: Data-Management Tools --> Features --> Add xy coordinates)

dann haben die Daten x,y,z Spalten

dann

feature to raster (zu finden unter: conversion tools --> to raster --> feature to raster).

Ich habe das mal gestestet und das tut's soweit. Allerdings habe ich die Testdaten nicht zu einem Mosaik zusammengefügt, es sollte aber stimmen, die haben ja alle das gleiche Koordinatensystem zugeordnet bekommen.

Das dauert ziemlich lange schon... 35 Stunden hätte ich auch gern.

Beste Grüße
Jochen
Hallo Jochen,
beim Import der XYZ-Dateien mit dem Tool 3D-Analyst Tools / Konvertierung / von Datei / 3D ASCII zu Feature Class ist es möglich, als Eingabe einen ORDNER anzugeben. Dann muss man die Dateiendung eintragen, das verwendete Kommazeichen und ggf. noch die durchschnittliche Punktdichte.
Es werden dann ALLE Dateien mit der entsprechenden Endung aus dem Ordner eingelesen. Damit lassen sich auch umfangreiche Datenmengen in EINE Multipunkt FeatureClass importieren und können dann in einem Arbeitsschritt weiter verarbeitet werden.

Damit hab ich mal 6,8GB Daten aus 50 Blattschnitten Laserscan-Punktwolken (jeweils mehrere XYZ-Dateien pro Blattschnitt), 253 Millionen Punkte in 20 Minuten eingelesen. Das Umwandeln in ein Terrain und daraus Umwandlung in ein Raster dauerte jeweils weitere 20 Minuten. Da lohnte sich keine Model Builder- oder Python-Bastelei.
Wie umfangreich sind deine Daten? Ich denke es würde auch mit größeren Datenmengen noch funktionieren, eine FileGDB soll doch bis 1TB Daten speichern können, oder?
Grüße
Rena
Hallo Rena,

danke für den Tipp. Ich hab das ausprobiert aber GIS streikt. Es wird zwar 1 Ausgabe-Shapefile erstellt, aber das ist dann nur 1 kb groß (also jedes Teilshape wie *.prj, *.sbn usw). Es enthält beim Import dann eine leere Attributtabelle. Laut Task-Manager tut sich da auch nichts im Hintergrund...
Wenn ich jedoch GIS beenden möchte, kommt eine "Pending Geoprocessing Process Warning", braucht der Prozess evtl. eine Art Warmlaufzeit?

Beste Grüße
Jochen
Hallo Jochen,
Laserscan- oder DGM-Punkte NIE in ein Shape schreiben, immer in eine File-Geodatabase (bei einem Testlauf hier hats aber auch mit shape funktioniert, allerdings nur bei kleinen Datenmengen). Und wenn du danach ein Terrain erstellen willst (sehr empfehlenswert, bei großen Datenmengen ist TIN zu schwerfällig) am besten gleich in eine FeatureClass in einem Dataset in einer FileGeodatabase.

Wenn kein Ausgabefile erstellt wird, prüfe nochmal genau deine Daten (eine txt in WordPad öffnen) und die Einstellungen: ist das Dezimaltrennzeichen korrekt angegeben (Voreinstellung ist der Punkt, also amerikanische Schreibweise)? Ist das Dateisuffix richtig?
X , Y, Z –Werte müssen in den ersten drei Spalten stehen, Spaltenüberschriften sind nicht notwendig aber auch nicht störend. Als Eingabeformat "xyz" auswählen. Auch die durchschnittliche Punktdichte wird benötigt, die kannst du mit dem 3DAnalyst Tool "Punktdateiinformationen" ermitteln. Das ist auch gleich ein Test dafür, ob deine Einstellungen stimmen. Wenn das Tool Punktdateiinformationen die Dateien korrekt einlesen kann sollte auch das Einlesen der Punkte funktionieren.
Grüße
Rena
PS: mir ist auch nicht ganz klar, warum du in deinem ModelBuilder das Tool "add XY-Coordinates" benötigst. Die txt-Files enthalten doch bereits die Koordinaten, die musst du doch nicht neu berechnen??
Hallo Rena,

stimmt, das muss ich mal ausprobieren, ob ich das brauch. In der Attributtabelle der shapesfiles (wenn ich meinen Modelbuilder-Ansatz verwende) tauchen bislang keine x-y-z Spalten auf. Daher schreibe ich in einem Schritt x-y Koordinaten rein und im nächsten dann auch noch die z-Koordinaten. Die z-Koordinaten brauche ich für das Aufrastern (conversion tools feature to raster), um daraus ein DGM zu erstellen.
und bei der Funktion "ASCII 3D to Feature Class" habe ich nur die Ausgabefunktion *shp, ich kann da gar nicht auf Geodatabase umstellen, das bessert GIS automatisch aus und schreibt als Endung wieder *shp hin
Hallo Jochen,

1. Welche ArcGIS-Version hast du? Ich hab hier ArcGIS 10, die FileGeodatabase gibts aber schon mindestens seit 9.3?

2. Wenn du die Daten mit dem 3D-Analysttool einliest sollte ein 3D-Shape rauskommen (in der Attributspalte Shape* steht dann PunktZM bzw. MultipunktZM). Der Rechtswert ist dann im Shape intern gespeichert und muss nicht unbedingt nochmal in der Attributtabelle auftauchen. Beim Generieren des Rasters kann der im 3D-Shape enthaltene Z-Wert genommen werden (im DropDownMenü Feld bei 3D-Daten "Shape.Z" auswählen).

3. Wennich hier mit ArcGIS 10.0 und dem Tool "3DASCII zu FeatureClass" die Punkte einlese wird eine MULTIPOINT-Datei erstellt (Voreinstellung bei "Typ der Ausgabe"). d.h. es hat nicht mehr jeder Einzelpunkt eine Attributzeile. Das ist der eigentliche Grund, warum das Einlesen vieler Millionen Punkte recht schnell geht: es werden immer Multipunktgruppen von bis zu 5000 Punkten gebildet. Für die kann man dann nicht mehr individuelle Z-Werte berechnen. Muss man ja auch nicht, Siehe 2.

4. In ArcGIS 10.0 wird als Ausgabedatensatz immer zunächst die "Default.gdb" vorgeschlagen, also eine FileGDB. Das kann natürlich in älteren Programmversionen anders sein, da gabs noch keine Default.gdb. Wenn man eine andere FileGDB nimmt, muss diese bereits vorhanden sein. Die wird nicht beim Import erzeugt!
Hallo Rena,

1. ich benutze hier Version 10.1 SP 1

2. also doch 1 Shapefile, oder? das mit der Geodatabase verwirrt mich etwas...

3. ok, das klingt gut, ich probiere das Multipoint file

4. das wurde dann wohl bei Version 10.1 wieder umgestellt, hier ist es ja wie bei 3.) das Multipoint Shapefile

Beste Grüße
Jochen
zu 2: Das war von mir missverständlich ausgedrückt.
Auch in der Attributtabelle einer FeatureClass einer FileGeodatabase gibts eine Spalte namens Shape*, in der der eingetragen ist, um welche Art von features es sich handelt. Also beispielsweise Polylinie, Polygon, Punkt oder eben MultipointZM. Dabei steht Z für Höhenwerte, also 3D-Daten.

Wenn du aus dem „3D-ASCII zu FeatureClass" Tool als Ablageort einfach einen Dateiordner wählst und einen Namen eintippst, wird halt eine Shape-Datei gespeichert. Eine FeatureClass braucht als Speicherort eine VORHANDENE FileGDB.
Und nochwas: mit Laserscandaten besser nicht auf verknüpften Netzlaufwerken arbeiten, sondern in einem Dateiordner direkt auf C:/.


Nochmal zur besten Vorgehensweise:
a) ZUERST eine FileGeodatabase erzeugen. In dieser ein Dataset anlegen.
Die FileGDB aus ArcMap /ArcCAtalog als StandardGeodatabase wählen (rechtsklick, Kontextmenü "Als Standard Geodatabase definieren") und das ArcMap-Dokument abspeichern.

b) Dann im Tool "3D ASCII zu FeatureClass" als Ausgabe-Feature-Class das Dataset in der FileGDB wählen. Wenn du wie oben angegeben diese FileGDB als Standard definiert hat, wird die auch gleich als Ablage vorgeschlagen. Das erspart nochmal ein paar Mausklicks.

c) Aus dieser Multipoint-FeatureClass kannst du jetzt entweder ein Terrain erzeugen (wenn es sich um unregelmäßige Laserscanpunktwolken handelt)
oder
direkt über das Tool Feature zu Raster ein Rastergrid erzeugen (wenn es sich um Punkte in festgelegten, regelmäßigen Abständen handelt)

d) Falls du aus unregelmäßige Punktwolken ein Terrain erzeugt hast, musst du das noch in ein Raster umwandeln.

Ich mach Einlesen und Umwandelt von Höhendaten nie im ModelBuilder, weil der mir hier keine Zeit spart, eher im Gegenteil.
Grüße
Rena
Hallo Rena,

vielen Dank für deine ausführliche Beschreibung.
Ich habe den Prozess nach deiner Anleitung gestartet und werde berichten, wie es gegangen ist.

Das arbeiten mit einer Geodatabase ist mir noch recht neu, aber es sieht sehr komfortabel aus.

Beste Grüße
Jochen
Hallo Rena,

das hat ja super funktioniert und wesentlich schneller als mein Ansatz.
Danke nochmal für die ausführliche Beschreibung.

Beste Grüße
Jochen
Hallo Jochen, hallo Rena,

mit Interesse habe ich gelesen, wie schnell und einfach das mit dem Tool "3D-ASCII to FeatureClass" aus den 3D-Analyst-Tools geht. Ich könnte richtig neidisch werden, denn leider habe ich kein 3D-Analyst, sondern nur Spatial Analyst. Mir wird daher wohl nichts anderes übrig bleiben, als auch weiterhin für jede xyz-Datei die Vorgehensweise

- Make XY Event Layer (Toolbox Data Management) und
- Feature To Raster Conversion (Toolbox Spatial Analyst)

anzuwenden, um zu Grids zu kommen. Und das braucht seine Zeit.

Viele Grüße
Gigi