[GIS] Creating nanoscale DEM with GDAL

coordinate systemdempython

A bit of a strange question perhaps, but let me give you a breif explanation of the background before my actual questions:

Atomic force microscopy (AFM) is a method that, in short (and to my limited knowledge) lets researchers scan areas at the micro- and nanoscale. It works by "scanning" an area using a probe of sorts. More that is difficult for me to explain, as I have no real understanding of it. What I do know, and what triggered my curiosity was that the result is in fact a "grid" of "height" values (a matrix of say 512×512 values describing the height of the probe at that point).

I then thought: Well, apart from the scale, this is in fact a digital elevation model! And, this means that if I can manage to create a DEM file as understood by GIS tools I could apply GIS analysis to it!

As it is, my significat other works at a lab that has a AFM-machine, and is using it in one of her projects. I got some scan files from her, and have managed, using Python (struct and numpy), to parse these binary files and what I have now is a numpy array of size 512×512 filled with int16-values.

What I'm planning on next, and what I need some help with is the "mapping to a proper DEM"-part. I have some knowledge about DEMS, but when it comes to actual generation of them I am quite new.

What I'm thinking is that I have to georeference my data in some way, and for this I need a custom (planar) coordinate system. I envision my coordinate system would use micro- or nano-meters as units. Then it's just a matter of finding the size of the area scanned with the AFM (this I believe is somewhere in the binary file, assume this known).

update: I also have several scans at different resolutions, but of the same area. For example i have this information about two scans:

larger image:

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

smaller (detail) image:

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

What I', thinking is that my custom coordinate system should have an origin in 0,0 and for the larger image I cam assign pixel 0,0 the coordinate value of (0,0) and pixel 512,512 the coordinate value (51443.5, 51443.5) (Guess you get the picture for the other points needed).

Then, the larger image would map pixel (0,0) to (8776.47, 1486.78) and (512,512) to (8776.47 + 5907.44, 1486.78 + 5907.44)

1st Question is then: How do I create a proj4 def for such a coordinate system? I.e: how do I assign these "real world coordinates" to my custom coordinate system (or, if I follow whubers suggestion and using a local coordinate system and lying about units (i.e. treating my nanometers as kilometers)

Then I have to transfer my numpy 2-dimensional Array to a Georeferenced DEM file format. I was thinking of using GDAL (or, rather the Python bindings).

2nd Question is then: How do I create a georeferenced DEM from "arbitrary" data such as mine? Preferably in Python and using open source libraries.

The rest should then be fairly easy going, just a matter of using the right analysis tools. The problem is, this task is driven by my own curiosity, so I'm not quite sure what I actually should do with a nanoscale DEM. This begs the

3rd question: What to do with a nanoscale DEM? What kind of analysis can be done, what are appropriate tools for DEM analysis and finally: is it feasible to make a map with hillshade and contour lines from this data?:)

I welcome all suggestions and pointers, but keep in mind that I am looking for free alternativies, as this is a strictly hobby-based project with no budget or funding (and I do not have access to any lisenced GIS-applications). In addition, I do know that Bruker, the company that sells these AFM machines, ships some software, but using that would not be any fun.

Best Answer

Well, it seems like I solved problems 1 & 2 at least. Full source code at github, but some explaination here:

For a custom CRS I decided (per Whubers suggestion) to "cheat" and use meters as unit. I found a "local crs" at 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]]

Using Python and GDAL this is fairly easy to read: 

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

Also, greating a DEM with GDAL was actually rather straightforward (I ended up with a single-band geo-tiff). The line parser.read_layer(0) returns my previously described 512x512 matrix.

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

The triockiest part was figuring out how to properly "georeference" my file, i ended up using SetGeoTransform, getting the parameters as so:

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]

This last part is probably the one I am the most unsure about, what I really was looking for was something line *gdal_transform -ullr*, but I couldn't find way to do this programatically.

I am able to open my GeoTIFF in Qgis and view it (and visually comparing it with the result from the Bruker program it looks right), but I haven't really answered my question 3; what to do with this data. So, here I'm open for suggestions!

Related Question