Creando DEM a nanoescala con GDAL

14

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

Microscopía de fuerza atómica (AFM) es un método que, en resumen (y para mi limitado conocimiento) permite a los investigadores escanear Zonas a micro y nanoescala. Funciona al "escanear" un área usando una especie de sonda. Más me resulta difícil explicarlo, 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 valores de 512x512, que describen la altura de la sonda en ese punto).

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

Tal como está, mis otros trabajos significativos en un laboratorio que tiene una máquina AFM, y la está utilizando en uno de sus proyectos. Conseguí algunos archivos escaneados de ella, y he manejado, usando Python (struct y numpy), para analizar estos archivos binarios y lo que tengo ahora es una gran variedad de tamaños 512x512 llenos de valores int16.

Lo que planeo a continuación, y con lo que necesito ayuda es con la parte "mapeando a un DEM adecuado". 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 (planar). Imagino que mi sistema de coordenadas usaría micrómetros o nanómetros como unidades. Entonces solo es cuestión de encontrar el tamaño del área escaneada con el AFM (creo que esto está en algún lugar del archivo binario, supongamos que se sabe).

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

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 yo ', pensando es que mi sistema de coordenadas personalizado debe tener un origen en 0,0 y para la imagen más grande asigno al píxel 0,0 el valor de la coordenada de (0,0) y al píxel 512,512 el valor de la coordenada ( 51443.5, 51443.5) (Supongo que obtendrá la imagen para los otros puntos necesarios).

Luego, la imagen más grande asignaría el píxel (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 una definición de proj4 para este 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 sobre las unidades (es decir, tratar mis nanómetros como kilómetros)

Luego, tengo que transferir mi matriz numpy bidimensional 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 utilizando bibliotecas de código abierto.

El resto debería ser bastante sencillo, solo es cuestión de utilizar 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 realmente debería hacer con un DEM a nanoescala. Esto ruega a la

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

Doy la bienvenida a 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 fondos (y no tengo acceso a ninguna aplicación GIS autorizada) . Además, sí sé que Bruker, la compañía que vende estas máquinas AFM, entrega algún software, pero usarlo no sería divertido.

    
pregunta atlefren 01.06.2013 - 12:00

1 respuesta

4

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

Para un CRS personalizado, decidí (por sugerencia de Whubers) "hacer trampa" y usar los 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]]

El uso de Python y GDAL 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, crear 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 interesante fue descubrir cómo "georreferenciar" correctamente mi archivo. Terminé usando SetGeoTransform , obteniendo los parámetros como tal:

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 la que más me preocupa, lo que realmente estaba buscando era algo en la línea * gdal_transform -ullr *, pero no pude encontrar la manera de hacerlo programáticamente.

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

    
respondido por el atlefren 04.06.2013 - 10:30

Lea otras preguntas en las etiquetas