Tworzenie nanoskali DEM za pomocą GDAL

Być może to trochę dziwne pytanie, ale pozwólcie, że wyjaśnię to krótko, zanim zacznę moje rzeczywiste pytania:

Mikroskopia sił atomowych (AFM) to metoda, która w skrócie (i o mojej ograniczonej wiedzy) pozwala badaczom skanować obszary w mikro- i nanoskali. Działa poprzez „skanowanie” obszaru za pomocą rodzaju sondy. Trudno mi to wyjaśnić, ponieważ nie rozumiem tego. To, co wiem, i to, co wzbudziło moją ciekawość, to fakt, że rezultatem jest w rzeczywistości „siatka” wartości „wysokości” (macierz powiedzmy wartości 512x512 opisująca wysokość sondy w tym punkcie).

Pomyślałem wtedy: cóż, oprócz skali, jest to w rzeczywistości cyfrowy model wysokości! Oznacza to, że jeśli uda mi się utworzyć plik DEM w rozumieniu narzędzi GIS, mogę zastosować do niego analizę GIS!

W tej chwili moje znaczenie ma inne, które pracują w laboratorium, które ma maszynę AFM i wykorzystuje ją w jednym ze swoich projektów. Dostałem od niej kilka plików skanowania i za pomocą Pythona (struct i numpy) udało mi się parsować te pliki binarne, a teraz mam tablicę numpy o wielkości 512x512 wypełnioną wartościami int16.

To, co planuję dalej i potrzebuję pomocy, to „mapowanie do właściwego DEM”. Mam trochę wiedzy na temat DEMS, ale jeśli chodzi o ich generowanie, jestem całkiem nowy.

Myślę, że muszę w jakiś sposób georeferenować swoje dane, a do tego potrzebuję niestandardowego (płaskiego) układu współrzędnych. Przewiduję, że mój układ współrzędnych użyłby mikro- lub nanomierzy jako jednostek. Potem jest tylko kwestia znalezienia rozmiaru obszaru skanowanego za pomocą AFM (myślę, że jest to gdzieś w pliku binarnym, załóżmy, że to znane).

aktualizacja : Mam również kilka skanów w różnych rozdzielczościach, ale tego samego obszaru. Na przykład mam te informacje o dwóch skanach:

większy obraz:

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

mniejszy (szczegółowy) obraz:

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

Myślę, że mój niestandardowy układ współrzędnych powinien mieć początek w 0,0, a dla większego obrazu przypisuję pikselowi 0,0 wartość współrzędnej (0,0), a piksel 512,512 wartość współrzędnej (51443.5, 51443.5 ) (Chyba masz zdjęcie dla innych potrzebnych punktów).

Następnie większy obraz zamapowałby piksele (0,0) na (8776.47, 1486.78) i (512.512) na (8776.47 + 5907.44, 1486,78 + 5907.44)

Pierwsze pytanie brzmi : jak utworzyć definicję proj4 dla takiego układu współrzędnych? Tj .: w jaki sposób przypisać te „rzeczywiste współrzędne świata” do mojego niestandardowego układu współrzędnych (lub, jeśli zastosuję się do sugestii Whubera i użyję lokalnego układu współrzędnych i będę kłamał o jednostkach (tj. Traktując moje nanometry jako kilometry)

Następnie muszę przenieść moją numpy 2-wymiarową macierz do formatu pliku DEM z georeferencją. Myślałem o użyciu GDAL (a raczej powiązań Python).

Drugie pytanie brzmi zatem : jak utworzyć georeferencyjny DEM z „dowolnych” danych, takich jak moje? Najlepiej w Pythonie i przy użyciu bibliotek open source.

Reszta powinna być wtedy dość łatwa, wystarczy kwestia korzystania z odpowiednich narzędzi analitycznych. Problem polega na tym, że to zadanie wynika z mojej własnej ciekawości, więc nie jestem do końca pewien, co właściwie powinienem zrobić z DEM w nanoskali. To błaga

Trzecie pytanie : co zrobić z nanoskalową DEM? Jakiego rodzaju analizę można wykonać, jakie są odpowiednie narzędzia do analizy DEM i na koniec: czy możliwe jest stworzenie mapy z cieniem wzgórza i liniami konturowymi na podstawie tych danych? :)

Z zadowoleniem przyjmuję wszystkie sugestie i wskazówki, ale pamiętam, że szukam darmowych alternatyw, ponieważ jest to projekt ściśle hobbystyczny bez budżetu i finansowania (i nie mam dostępu do żadnych aplikacji GIS pozbawionych wolności). Ponadto wiem, że Bruker, firma sprzedająca te maszyny AFM, dostarcza oprogramowanie, ale korzystanie z niego nie byłoby przyjemnością.

Wygląda na to, że przynajmniej rozwiązałem problemy 1 i 2. Pełny kod źródłowy na github , ale kilka wyjaśnień tutaj:

W przypadku niestandardowego CRS postanowiłem (zgodnie z sugestią Whubersa) „oszukiwać” i używać mierników jako jednostki. Znalazłem „lokalnego crs” na apatialreference.org ( SR-ORG: 6707 ):

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

Korzystanie z Pythona i GDAL jest dość łatwe do odczytania:

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

Ponadto powitanie DEM za pomocą GDAL było właściwie dość proste (skończyłem na jednopasmowym geo-tiffie). Linia parser.read_layer (0) zwraca moją wcześniej opisaną macierz 512x512.

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

Najtrudniejszą częścią było ustalenie, jak poprawnie „georeferencyjnie” mój plik, ostatecznie wykorzystałem SetGeoTransform , otrzymując parametry w następujący sposób:

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]

Ta ostatnia część jest prawdopodobnie tą, której nie jestem najbardziej pewna. Naprawdę szukałem czegoś w stylu * gdal_transform -ullr *, ale nie mogłem znaleźć sposobu na zrobienie tego programowo.

Jestem w stanie otworzyć mój GeoTIFF w Qgis i wyświetlić go (i wizualnie porównując go z wynikiem z programu Bruker, wygląda dobrze), ale tak naprawdę nie odpowiedziałem na moje pytanie 3; co zrobić z tymi danymi. Więc jestem otwarty na sugestie!