Erstellen von nanoskaligem DEM mit GDAL

14

Eine etwas seltsame Frage, aber lassen Sie mich vor meinen eigentlichen Fragen eine kurze Erklärung des Hintergrunds geben:

Die Rasterkraftmikroskopie (AFM) ist eine Methode, mit der Forscher (und meines Wissens nach) Bereiche im Mikro- und Nanobereich scannen können. Es funktioniert, indem ein Bereich mit einer Art Sonde "gescannt" wird. Mehr ist für mich schwer zu erklären, da ich kein wirkliches Verständnis dafür habe. Was ich weiß und was meine Neugierde auslöste, war, dass das Ergebnis tatsächlich ein "Gitter" von "Höhen" -Werten ist (eine Matrix von beispielsweise 512 x 512 Werten, die die Höhe der Sonde an diesem Punkt beschreiben).

Ich dachte dann: Naja, abgesehen von der Skala ist dies tatsächlich ein digitales Höhenmodell! Das bedeutet, dass ich eine GIS-Analyse darauf anwenden kann, wenn ich eine DEM-Datei erstellen kann, wie sie von GIS-Tools verstanden wird!

So wie es ist, arbeitet mein Lebensgefährte in einem Labor, das eine AFM-Maschine hat, und verwendet sie in einem ihrer Projekte. Ich habe einige Scandateien von ihr erhalten und habe es mit Python (struct und numpy) geschafft, diese Binärdateien zu analysieren. Was ich jetzt habe, ist ein numpy-Array der Größe 512x512, das mit int16-Werten gefüllt ist.

Was ich als nächstes vorhabe und bei dem ich Hilfe brauche, ist die "Zuordnung zu einem richtigen DEM" -Teil. Ich habe einige Kenntnisse über DEMS, aber wenn es um die tatsächliche Erzeugung von DEMS geht, bin ich ziemlich neu.

Was ich denke ist, dass ich meine Daten in irgendeiner Weise georeferenzieren muss und dafür ein benutzerdefiniertes (planares) Koordinatensystem benötige. Ich stelle mir vor, mein Koordinatensystem würde Mikro- oder Nanometer als Einheiten verwenden. Dann geht es nur noch darum, die Größe des mit dem AFM gescannten Bereichs zu ermitteln (dies ist meines Erachtens irgendwo in der Binärdatei enthalten, vorausgesetzt, dies ist bekannt).

Update : Ich habe auch mehrere Scans mit unterschiedlichen Auflösungen, aber aus dem gleichen Bereich. Zum Beispiel habe ich diese Informationen über zwei Scans:

größeres Bild:

Scan Size: 51443.5 nm
X Offset: 0 nm
Y Offset: 0 nm

kleineres (detail) bild:

Scan Size: 5907.44 nm
X Offset: 8776.47 nm
Y Offset: 1486.78 nm

Ich denke, mein benutzerdefiniertes Koordinatensystem sollte einen Ursprung in 0,0 haben, und für das größere Bild kann ich Pixel 0,0 den Koordinatenwert (0,0) und Pixel 512.512 den Koordinatenwert (51443,5, 51443,5) zuweisen ) (Vermutlich bekommst du das Bild für die anderen benötigten Punkte).

Dann würde das größere Bild Pixel (0,0) auf (8776,47, 1486,78) und (512,512) auf (8776,47 + 5907,44, 1486,78 + 5907,44) abbilden.

1. Frage ist dann : Wie erstelle ich ein proj4 def für ein solches Koordinatensystem? Dh, wie ordne ich diese "realen Weltkoordinaten" meinem benutzerdefinierten Koordinatensystem zu (oder, wenn ich dem Vorschlag von Whubers folge und ein lokales Koordinatensystem verwende und über Einheiten lüge (dh meine Nanometer als Kilometer behandle)

Dann muss ich mein numpy 2-dimensionales Array in ein georeferenziertes DEM-Dateiformat übertragen. Ich habe darüber nachgedacht, GDAL (oder besser gesagt die Python-Bindungen) zu verwenden.

2. Frage ist dann : Wie erstelle ich ein georeferenziertes DEM aus "beliebigen" Daten wie meinen? Am besten in Python und unter Verwendung von Open Source-Bibliotheken.

Der Rest sollte dann ziemlich einfach sein, nur eine Frage der Verwendung der richtigen Analysewerkzeuge. Das Problem ist, dass diese Aufgabe von meiner eigenen Neugier angetrieben wird, und ich bin mir nicht ganz sicher, was ich mit einem nanoskaligen DEM eigentlich tun soll. Das bittet den

3. Frage : Was tun mit einem nanoskaligen DEM? Welche Art von Analyse kann durchgeführt werden, welche Tools eignen sich für die DEM-Analyse und schließlich: Ist es möglich, aus diesen Daten eine Karte mit Schatten und Höhenlinien zu erstellen? :)

Ich begrüße alle Vorschläge und Hinweise, denke aber daran, dass ich nach freien Alternativen suche, da es sich um ein rein hobbybasiertes Projekt ohne Budget oder Finanzierung handelt (und ich keinen Zugriff auf lizenzierte GIS-Anwendungen habe). Außerdem weiß ich, dass Bruker, das Unternehmen, das diese AFM-Geräte verkauft, eine Software ausliefert, die aber keinen Spaß machen würde.

atlefren
quelle
Lustig und interessant! Können Sie einige Beispieldaten posten? Muss die Projektion eine Nanometerskala haben? Ich denke, es ist vielleicht einfacher zu skalieren, obwohl es ein bisschen "schummelt". Übrigens, ich denke, Sie können mit GDAL / ogr einen langen Weg zurücklegen, obwohl das Projektionsproblem noch angegangen werden muss. gdal.org/gdal_grid.html
alexanno
Vielen Dank! Schätze, das ist eher ein Kommentar als eine Antwort. Was die Nanometerskala angeht, würde ich sagen, dass alles, was funktioniert, großartig ist, aber eine echte nanoskalige Peojektion wäre am coolsten. Wenn es um Beispieldaten geht, muss ich prüfen, ob es einige Einschränkungen gibt, aber im Grunde handelt es sich um eine zweidimensionale Matrix aus int16-Werten.
atlefren
3
Warum brauchen Sie proj4-Parameter? Sie werden diese Daten nicht in ein anderes Koordinatensystem (insbesondere geografisch) übertragen. Imho sollten Sie in der Lage sein, alle Ihre Analysen ohne Koordinatensystem durchzuführen. Was verstehen Sie unter einem DEM? Es gibt verschiedene Typen wie triangulierte Flächen (zB Delauny Triangulation) oder Rasterkarten (Sie haben dies bereits). Dies hängt natürlich stark von der Analysesoftware ab. Natürlich können Sie auch andere Maps erstellen, wenn diese zum Verständnis der Sonde als Ausgabe benötigt werden. Vielleicht haben Sie einen Blick auf nehmen code.google.com/p/mtex für Getreide anlysis.
Mistapink
3
Der Grund, warum ich ein CRS benötige, ist folgender: Wenn ich nur ein GeoTIFF erstelle, ohne ein CRS zuzuweisen, ist die Maßeinheit Pixel. Was ist, wenn ich Entfernungen messen möchte? Und was ist, wenn ich zwei AFM-Scans habe und weiß, wie sie sich zueinander verhalten (in Bezug auf Skalierung und Versatz von einem bestimmten Punkt aus). Wenn ein CRS zugewiesen wird, können mehrere Scans gleichzeitig
angezeigt werden
1
Normalerweise bewältige ich solche Daten, indem ich ein lokales Koordinatensystem (mit einem Ursprung, der mit dem Bildursprung übereinstimmt) einstelle und über den Einheiten "lüge". Sie könnten zum Beispiel festlegen, dass die Einheiten Kilometer sind, wenn sie wirklich Nanometer sind, was es mental einfach macht, hin und her zu konvertieren. Natürlich werden Sie keine Neuprojektionen durchführen, das ist also kein Problem. Das Einrichten dieses Koordinatensystems entspricht dem Georeferenzieren eines DEM. Es kann so einfach sein, wie eine nicht gedrehte Weltdatei zu erstellen.
Whuber

Antworten:

4

Nun, es scheint, als hätte ich zumindest die Probleme 1 und 2 gelöst. Vollständiger Quellcode bei Github , aber einige Erläuterungen hier:

Für ein benutzerdefiniertes CRS entschied ich mich (nach Whubers Vorschlag) zu "schummeln" und Meter als Einheit zu verwenden. Ich habe ein "lokales crs" bei apatialreference.org gefunden ( SR-ORG: 6707 ):

LOCAL_CS["Non-Earth (Meter)",
    LOCAL_DATUM["Local Datum",0],
    UNIT["Meter",1.0],
    AXIS["X",NORTH],
    AXIS["Y",EAST]]

Mit Python und GDAL ist dies ziemlich einfach zu lesen:

def get_coordsys():
    #load our custom crs
    prj_text = open("coordsys.wkt", 'r').read()
    srs = osr.SpatialReference()
    if srs.ImportFromWkt(prj_text):
        raise ValueError("Error importing PRJ information" )

    return srs

Außerdem war es ziemlich unkompliziert, ein DEM mit GDAL zu erstellen (am Ende hatte ich einen Single-Band-Geotiff). Die Zeile parser.read_layer (0) gibt meine zuvor beschriebene 512x512-Matrix zurück.

def create_dem(afmfile, outfile):

    #parse the afm data
    parser = AFMParser(afmfile)

    #flip to account for the fact that the matrix is top-left, but gdal is bottom-left
    data = flipud(parser.read_layer(0))

    driver = gdal.GetDriverByName("GTiff")
    dst_ds = driver.Create(
        outfile,
        data.shape[1],
        data.shape[0],
        1 ,
        gdal.GDT_Int16 ,
    )

    dst_ds.SetGeoTransform(get_transform(parser, data))
    dst_ds.SetProjection(get_coordsys().ExportToWkt())
    dst_ds.GetRasterBand(1).WriteArray(data)
    dst_ds = None

Der schwierigste Teil war, herauszufinden, wie man meine Datei richtig "georeferenziert". Am Ende benutzte ich SetGeoTransform und erhielt die folgenden Parameter:

def get_transform(parser, data):
    #calculate dims
    scan_size, x_offset, y_offset = parser.get_size()
    x_size = data.shape[0]
    y_size = data.shape[1]
    x_res = scan_size / x_size
    y_res = scan_size / y_size

    #set the transform (offsets doesn't seem to work the way I think)
    #top left x, w-e pixel resolution, rotation, 0 if image is "north up", top left y, rotation, 0 if image is "north up", n-s pixel resolution
    return [0, x_res, 0, 0, 0, y_res]

Dieser letzte Teil ist wahrscheinlich der, bei dem ich mir am unsichersten bin. Was ich wirklich gesucht habe, war eine Zeile * gdal_transform -ullr *, aber ich konnte keinen Weg finden, dies programmatisch zu tun.

Ich kann mein GeoTIFF in Qgis öffnen und anzeigen (und es visuell mit dem Ergebnis des Bruker-Programms vergleichen, es sieht richtig aus), aber ich habe meine Frage 3 nicht wirklich beantwortet. Was tun mit diesen Daten? Also, hier bin ich offen für Vorschläge!

atlefren
quelle
Eine interessante Herausforderung könnte darin bestehen, Entfernungen auf dem DEM mit Entfernungen zwischen Orten auf dem Globus zu vergleichen, um den Zuschauern eine Vorstellung davon zu geben, wie klein die Nanoskala ist. ZB htwins.net/scale2
blah238