Crear DEM a nanoescala con GDAL

14

Tal vez una pregunta un poco extraña, pero déjame darte una breve explicación de los antecedentes antes de mis preguntas reales:

La microscopía de fuerza atómica (AFM, por sus siglas en inglés ) es un método que, en resumen (y para mi conocimiento limitado) permite a los investigadores escanear áreas a micro y nanoescala. Funciona "escaneando" un área utilizando una especie de sonda. Es más difícil de explicar, ya que no tengo una comprensión real de ello. Lo que sí sé, y lo que provocó mi curiosidad fue que el resultado es, de hecho, una "cuadrícula" de valores de "altura" (una matriz de, digamos, valores de 512x512 que describe la altura de la sonda en ese punto).

Entonces pensé: Bueno, aparte de la escala, ¡este es de hecho un modelo de elevación digital! ¡Y esto significa que si puedo crear un archivo DEM tal como lo entienden las herramientas SIG, podría aplicarle el análisis SIG!

Tal como están las cosas, mi otra persona trabaja en un laboratorio que tiene una máquina AFM y la está usando en uno de sus proyectos. Obtuve algunos archivos de escaneo de ella, y he logrado, usando Python (estructura y numpy), analizar estos archivos binarios y lo que tengo ahora es una matriz numpy de tamaño 512x512 llena de valores int16.

Lo que estoy planeando a continuación, y con lo que necesito ayuda es la parte de "mapeo a un DEM apropiado". Tengo algunos conocimientos sobre DEMS, pero cuando se trata de la generación real de ellos, soy bastante nuevo.

Lo que estoy pensando es que tengo que georreferenciar mis datos de alguna manera, y para esto necesito un sistema de coordenadas personalizado (plano). Imagino que mi sistema de coordenadas usaría micro o nano metros como unidades. Entonces es solo una cuestión de encontrar el tamaño del área escaneada con el AFM (esto creo que está en algún lugar del archivo binario, suponga que esto se conoce).

Actualización : también tengo varios escaneos en diferentes resoluciones, pero de la misma área. Por ejemplo, tengo esta información sobre dos escaneos:

imagen más grande:

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

imagen más pequeña (detalle):

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

Lo que pienso es que mi sistema de coordenadas personalizado debe tener un origen en 0,0 y para la imagen más grande puedo asignar el píxel 0,0 al valor de coordenadas de (0,0) y el píxel 512.512 el valor de coordenadas (51443.5, 51443.5 ) (Supongo que obtienes la imagen de los otros puntos necesarios).

Luego, la imagen más grande mapearía los píxeles (0,0) a (8776.47, 1486.78) y (512,512) a (8776.47 + 5907.44, 1486.78 + 5907.44)

La primera pregunta es entonces : ¿Cómo creo un def de proj4 para dicho sistema de coordenadas? Es decir: ¿cómo asigno estas "coordenadas del mundo real" a mi sistema de coordenadas personalizado (o, si sigo la sugerencia de whubers y uso un sistema de coordenadas local y mentir sobre las unidades (es decir, tratar mis nanómetros como kilómetros)

Luego tengo que transferir mi matriz de 2 dimensiones numpy a un formato de archivo DEM georreferenciado. Estaba pensando en usar GDAL (o, más bien, los enlaces de Python).

La segunda pregunta es entonces : ¿Cómo creo un DEM georreferenciado a partir de datos "arbitrarios" como el mío? Preferiblemente en Python y usando bibliotecas de código abierto.

El resto debería ser bastante fácil, solo es cuestión de usar las herramientas de análisis adecuadas. El problema es que esta tarea está impulsada por mi propia curiosidad, por lo que no estoy muy seguro de lo que debería hacer con un DEM a nanoescala. Esto suplica al

3ª pregunta : ¿Qué hacer con un DEM a nanoescala? ¿Qué tipo de análisis se puede hacer, cuáles son las herramientas apropiadas para el análisis DEM y finalmente: es factible hacer un mapa con sombreado y líneas de contorno a partir de estos datos? :)

Agradezco todas las sugerencias y sugerencias, pero tenga en cuenta que estoy buscando alternativas gratuitas, ya que este es un proyecto estrictamente basado en pasatiempos sin presupuesto ni financiación (y no tengo acceso a ninguna aplicación GIS con licencia). Además, sé que Bruker, la compañía que vende estas máquinas AFM, envía algún software, pero usarlo no sería divertido.

python coordinate-system dem atlefren
fuente
¡Divertido e interesante! ¿Puedes publicar algunos datos de muestra? ¿Es un requisito tener una escala nanométrica en la proyección? Pensar que tal vez es más fácil de escalar, aunque es un poco "hacer trampa". Por cierto, supongo que puedes recorrer un largo camino con GDAL / ogr, aunque el problema de la proyección aún deberá abordarse. gdal.org/gdal_grid.html
alexanno
¡Gracias! Supongo que esto es más un comentario que una respuesta. En cuanto a la escala nanométrica, diría que cualquier cosa que funcione es genial, pero una verdadera peoyección a nanoescala sería la mejor. Cuando se trata de datos de muestra, tendré que verificar si hay algunas restricciones, pero es básicamente una matriz de 2 dim de valores int16.
atlefren
3
¿Por qué necesitarías parámetros proj4? No transferirá estos datos a otro sistema de coordenadas (especialmente geográfico). En mi opinión, deberías poder hacer todos tus análisis sin ningún sistema de coordenadas. ¿Qué entiendes por un DEM? Existen varios tipos, como las superficies trianguladas (p. Ej., Hacer una triangulación retardada) o mapas ráster (ya lo tiene). Por supuesto, esto depende en gran medida del software de análisis. Por supuesto, puede crear otros mapas si son necesarios como salida para comprender la sonda. Puede echar un vistazo a code.google.com/p/mtex para el análisis de grano.
equivoca el
3
La razón por la que (creo) que necesito un CRS es esta: si solo creo, digamos un GeoTIFF sin asignar un CRS, la unidad de medida será píxeles. ¿Qué pasa si quiero medir distancias? Y, ¿qué pasa si tengo dos escaneos AFM y sé cómo se relacionan entre sí (en términos de escala y compensación desde algún punto). Tener un CRS asignado facilitaría la visualización de varios escaneos a la vez
atlefren
1
Por lo general, hago frente a esos datos estableciendo un sistema de coordenadas local (con un origen coincidente con el origen de la imagen) y "mintiendo" sobre las unidades. Por ejemplo, podría estipular que las unidades son kilómetros cuando en realidad son nanómetros, lo que facilita la conversión mental de un lado a otro. Por supuesto, no hará ninguna reproyección, por lo que no es un problema. Configurar este sistema de coordenadas es lo mismo que georreferenciar cualquier DEM; puede ser tan simple como crear un archivo mundial sin rotar.
whuber

Respuestas:

4

Bueno, parece que resolví los problemas 1 y 2 al menos. Código fuente completo en github , pero alguna explicación aquí:

Para un CRS personalizado, decidí (según la sugerencia de Whubers) "engañar" y usar medidores como unidad. Encontré un "crs local" en 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]]

Usando Python y GDAL esto es bastante fácil de leer: 

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

Además, mejorar un DEM con GDAL fue bastante sencillo (terminé con un geo-tiff de banda única). La línea parser.read_layer (0) devuelve mi matriz 512x512 descrita anteriormente.

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 parte más complicada fue descubrir cómo "georreferenciar" correctamente mi archivo, terminé usando SetGeoTransform , obteniendo los parámetros de la siguiente manera:

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]

Probablemente, esta última parte es la que no estoy más segura, lo que realmente estaba buscando era algo línea * gdal_transform -ullr *, pero no pude encontrar la forma de hacerlo programáticamente.

Puedo abrir mi GeoTIFF en Qgis y verlo (y compararlo visualmente con el resultado del programa Bruker parece correcto), pero realmente no he respondido mi pregunta 3; Qué hacer con estos datos. ¡Entonces, aquí estoy abierto para sugerencias!

atlefren
fuente
Un desafío interesante podría ser comparar distancias en el DEM con distancias entre lugares en el mundo para dar a los espectadores una idea de cuán pequeña es la nanoescala. Por ejemplo, htwins.net/scale2
blah238