Creación de nanoescala DEM con GDAL

Question

Un poco rara la pregunta tal vez, pero te voy a dar una obtener una breve explicación de los antecedentes antes de mi pregunta:

Microscopía de fuerza atómica (AFM) es un método que, en definitiva (y a mi limitado conocimiento) permite a los investigadores de áreas de digitalización en el micro y nanoescala. Funciona mediante la "exploración", un área con una sonda de tipo. Más que es difícil para mí explicar, como yo no tienen una comprensión real de la misma. Lo que sí sé, y lo que provocó mi curiosidad fue que el resultado es, de hecho, una "red" de "altura" de los valores (matriz de decir 512x512 valores que describen la altura de la sonda en ese punto).

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

Como es, mi significativas otro trabaja en un laboratorio que tiene un AFM-máquina, y la está utilizando en uno de sus proyectos. Tengo algunos archivos de escaneado de ella, y han logrado, usando Python (struct y numpy), para analizar estos archivos binarios y de lo que tengo ahora es una colección de matriz de tamaño 512x512 lleno de int16-valores.

Lo que estoy pensando en el próximo, y lo que necesito un poco de ayuda, es la "asignación de una adecuada DEM". Tengo un poco de conocimiento acerca de los DEMS, pero cuando se trata de la generación real de ellos, yo soy bastante nuevo.

Lo que estoy pensando es que tengo la georreferencia de mis datos en alguna forma, y para ello necesito una costumbre (planar) sistema de coordenadas. Imagino a mi sistema de coordenadas sería la utilización de micro - o nano-metros como unidades. Entonces es simplemente una cuestión de encontrar el tamaño del área de escaneado con el AFM (esto creo que está en algún lugar en el archivo binario, asumir este conocido).

actualización: también tengo varios escaneos en diferentes resoluciones, pero de la misma área. Por ejemplo yo tengo esta información acerca de dos exámenes:

ampliar la imagen:

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

los más pequeños (detalle) de la imagen):

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

Lo que me', el pensamiento es que mi sistema de coordenadas de personalizado debe tener un origen en 0,0 y para la imagen más grande que yo cam asignar pixel 0,0 el valor de coordenadas (0,0) y pixel 512,512 el valor de coordenadas (51443.5, 51443.5) (Supongo que se obtiene la imagen de los puntos en los que sea necesario).

Luego, la más grande imagen de mapa de pixel (0,0) a (8776.47, 1486.78) y (512,512) a (8776.47 + 5907.44, 1486.78 + 5907.44)

1ª Pregunta entonces es: ¿Cómo puedo crear un proj4 def para un sistema de coordenadas? I. e: ¿cómo puedo asignar estos "mundo real las coordenadas" a mi costumbre el sistema de coordenadas (o, si sigo whubers sugerencia y el uso de un sistema de coordenadas local y mentir acerca de las unidades (es decir, el tratamiento de mi nanómetros como kilómetros)

Entonces tengo que transferir a mi numpy 2-dimensiones de la Matriz a un Georreferenciados DEM formato de archivo. Yo estaba pensando en usar GDAL (o, más bien, los enlaces Python).

2ª Pregunta entonces es: ¿Cómo puedo crear un georreferenciados DEM de "arbitraria" de datos tales como el mío? Preferiblemente en Python y el uso de bibliotecas de código abierto.

El resto debe entonces ser bastante fácil, sólo es cuestión de utilizar el derecho de las herramientas de análisis. El problema es que esta tarea es impulsado por mi curiosidad, así que no estoy muy seguro de lo que realmente debe hacer con un nanoescala DEM. Esto nos lleva a la

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

Doy la bienvenida a todas las sugerencias y consejos, pero ten en cuenta que estoy buscando sin alternativies, ya que esto es estrictamente una afición basado en proyecto sin presupuesto o financiamiento (y no tengo acceso a cualquier lisenced SIG-aplicaciones). Además, yo sé que Bruker, la empresa que vende estos AFM máquinas, barcos de algunos programas, pero el uso que no sería divertido.

Answer 1

Bueno, parece que he solucionado los problemas 1 y 2 al menos. El código fuente completo en github, pero algunos explicación aquí:

Para que una costumbre CRS decidí (por Whubers sugerencia) para "hacer trampa" y el uso de metros como unidad. Me encontré con un local "crs" 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

También, greating un DEM con GDAL era en realidad bastante sencillo (terminé con una sola banda de geo-tiff). La línea del analizador.read_layer(0) devuelve mi descrito previamente 512x512 de la matriz.

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

El triockiest parte fue averiguar cómo correctamente "georreferenciar" mi archivo, que terminó con SetGeoTransform, obteniendo los parámetros así:

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]

Esta última parte es, probablemente, el que yo soy el más seguro acerca de, lo que yo realmente estaba buscando era algo de la línea *gdal_transform -ullr*, pero no podía encontrar la manera de hacer esto a través de programación.

Soy capaz de abrir mi GeoTIFF en Qgis y ver (y visualmente comparando con el resultado de la Bruker programa se ve a la derecha), pero realmente no he contestado a mi pregunta 3; qué hacer con estos datos. Así que, aquí estoy abierto a sugerencias!