Création de DEM à l'échelle nanométrique avec GDAL

14

Une question étrange peut-être, mais permettez-moi de vous donner une brève explication du contexte avant mes vraies questions:

La microscopie à force atomique (AFM) est une méthode qui, en bref (et à ma connaissance limitée) permet aux chercheurs de scanner des zones à l'échelle micro et nanométrique. Il fonctionne en "scannant" une zone à l'aide d'une sorte de sonde. C'est plus difficile à expliquer pour moi, car je ne le comprends pas vraiment. Ce que je sais, et ce qui a déclenché ma curiosité, c'est que le résultat est en fait une "grille" de valeurs de "hauteur" (une matrice de valeurs de 512x512, par exemple, décrivant la hauteur de la sonde à ce point).

J'ai alors pensé: Eh bien, à part l'échelle, c'est en fait un modèle numérique d'élévation! Et cela signifie que si je parviens à créer un fichier DEM tel que compris par les outils SIG, je pourrais lui appliquer une analyse SIG!

Dans l'état actuel des choses, mes autres travaux importants dans un laboratoire qui a une machine AFM, et l'utilise dans l'un de ses projets. J'ai obtenu des fichiers d'analyse d'elle et j'ai réussi, en utilisant Python (struct et numpy), à analyser ces fichiers binaires et ce que j'ai maintenant est un tableau numpy de taille 512x512 rempli de valeurs int16.

Ce que je prévois ensuite, et ce dont j'ai besoin d'aide, c'est la partie "mappage vers un DEM approprié". J'ai quelques connaissances sur DEMS, mais en ce qui concerne leur génération réelle, je suis assez nouveau.

Ce que je pense, c'est que je dois géoréférencer mes données d'une manière ou d'une autre, et pour cela j'ai besoin d'un système de coordonnées (planaire) personnalisé. J'imagine que mon système de coordonnées utiliserait des micro ou nano-mètres comme unités. Ensuite, il s'agit simplement de trouver la taille de la zone numérisée avec l'AFM (cela, je crois, se trouve quelque part dans le fichier binaire, supposons que cela soit connu).

mise à jour : J'ai également plusieurs scans à différentes résolutions, mais de la même zone. Par exemple, j'ai ces informations sur deux scans:

image plus grande:

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

image plus petite (détail):

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

Ce que je pense, c'est que mon système de coordonnées personnalisé devrait avoir une origine en 0,0 et pour l'image plus grande, je dois attribuer au pixel 0,0 la valeur de coordonnées de (0,0) et au pixel 512,512 la valeur de coordonnées (51443.5, 51443.5 ) (Je suppose que vous obtenez l'image pour les autres points nécessaires).

Ensuite, l'image plus grande mapperait le pixel (0,0) à (8776.47, 1486.78) et (512.512) à (8776.47 + 5907.44, 1486.78 + 5907.44)

La 1ère question est alors : comment créer une définition de projet pour un tel système de coordonnées? C'est-à-dire: comment attribuer ces "coordonnées du monde réel" à mon système de coordonnées personnalisé (ou, si je suis la suggestion des whubers et en utilisant un système de coordonnées local et en mentant sur les unités (c'est-à-dire en traitant mes nanomètres comme des kilomètres)

Ensuite, je dois transférer mon tableau bidimensionnel numpy vers un format de fichier DEM géoréférencé. Je pensais utiliser GDAL (ou plutôt les liaisons Python).

La deuxième question est alors : comment créer un DEM géoréférencé à partir de données "arbitraires" comme la mienne? De préférence en Python et en utilisant des bibliothèques open source.

Le reste devrait alors être assez simple, il suffit d'utiliser les bons outils d'analyse. Le problème est que cette tâche est motivée par ma propre curiosité, donc je ne suis pas sûr de ce que je dois réellement faire avec un DEM à l'échelle nanométrique. Cela supplie le

3ème question : que faire avec un MNT à l'échelle nanométrique? Quel type d'analyse peut être fait, quels sont les outils appropriés pour l'analyse DEM et enfin: est-il possible de faire une carte avec des ombrages et des courbes de niveau à partir de ces données? :)

Je me réjouis de toutes les suggestions et conseils, mais gardez à l'esprit que je recherche des alternatives gratuites, car il s'agit d'un projet strictement basé sur les loisirs, sans budget ni financement (et je n'ai accès à aucune application SIG répertoriée). De plus, je sais que Bruker, la société qui vend ces machines AFM, expédie des logiciels, mais ce ne serait pas amusant.

atlefren
la source
Amusant et intéressant! Pouvez-vous publier des exemples de données? Faut-il avoir une échelle nanométrique sur la projection? Penser qu'il est peut-être plus facile de passer à l'échelle, bien que cela "triche" un peu. Btw Je suppose que vous pouvez faire beaucoup de chemin avec GDAL / ogr, bien que le problème de projection devra encore être résolu. gdal.org/gdal_grid.html
alexanno
Merci! Je suppose que c'est plus un commentaire qu'une réponse. En ce qui concerne l'échelle nanométrique, je dirais que tout ce qui fonctionne est génial, mais une véritable peojection à l'échelle nanométrique serait la plus cool. En ce qui concerne les échantillons de données, je devrai vérifier s'il existe des restrictions, mais il s'agit essentiellement d'une matrice à 2 dim de valeurs int16.
atlefren
3
Pourquoi auriez-vous besoin de paramètres proj4? Vous ne transférerez pas ces données dans un autre système de coordonnées (en particulier géographique). À mon humble avis, vous devriez être en mesure de faire toutes vos analyses sans aucun système de coordonnées. Qu'entendez-vous par un DEM? Il existe plusieurs types comme les surfaces triangulées (par exemple faire une triangulation delauny) ou les cartes raster (vous l'avez déjà). Cela dépend bien sûr fortement du logiciel d'analyse. Bien sûr, vous pouvez créer d'autres cartes si elles sont nécessaires en sortie pour comprendre la sonde. Vous pouvez consulter le code.google.com/p/mtex pour l' analyse des céréales.
mistapink
3
La raison pour laquelle j'ai (pense) avoir besoin d'un CRS est la suivante: si je viens de créer, disons un GeoTIFF sans assigner un CRS, l'unité de mesure sera les pixels. Et si je veux mesurer des distances? Et que se passe-t-il si j'ai deux analyses AFM et que je sais comment elles sont liées (en termes d'échelle et de décalage à partir d'un certain point). Un CRS attribué faciliterait la visualisation de plusieurs numérisations à la fois
atlefren
1
Je gère généralement ces données en définissant un système de coordonnées local (avec une origine coïncidant avec l'origine de l'image) et en «mentant» sur les unités. Par exemple, vous pourriez stipuler que les unités sont des kilomètres alors qu'elles sont vraiment des nanomètres, ce qui facilite la conversion mentale d'avant en arrière. Bien sûr, vous ne ferez aucun reprojet, donc ce n'est pas un problème. La configuration de ce système de coordonnées est la même que le géoréférencement de tout DEM; cela peut être aussi simple que de créer un fichier mondial sans rotation.
whuber

Réponses:

4

Eh bien, il semble que j'ai résolu au moins les problèmes 1 et 2. Code source complet sur github , mais quelques explications ici:

Pour un CRS personnalisé, j'ai décidé (par suggestion Whubers) de "tricher" et d'utiliser les compteurs comme unité. J'ai trouvé un "crs local" sur 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]]

En utilisant Python et GDAL, c'est assez facile à lire:

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

De plus, faire un DEM avec GDAL était en fait assez simple (je me suis retrouvé avec un géo-tiff à bande unique). La ligne parser.read_layer (0) renvoie ma matrice 512x512 décrite précédemment.

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

La partie la plus choquante était de savoir comment «géoréférencer» correctement mon fichier, j'ai fini par utiliser SetGeoTransform , obtenant les paramètres comme suit :

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]

Cette dernière partie est probablement celle dont je suis le plus incertain, ce que je cherchais vraiment était quelque chose de ligne * gdal_transform -ullr *, mais je n'ai pas trouvé de moyen de le faire par programme.

Je peux ouvrir mon GeoTIFF dans Qgis et le voir (et le comparer visuellement avec le résultat du programme Bruker, il a l'air correct), mais je n'ai pas vraiment répondu à ma question 3; que faire de ces données. Alors, ici, je suis ouvert aux suggestions!

atlefren
la source
Un défi intéressant pourrait être de comparer les distances sur le MNT aux distances entre des endroits sur le globe pour donner aux téléspectateurs une idée de la petite échelle. Par exemple htwins.net/scale2
blah238