# Creating nanoscale DEM with GDAL

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 512x512 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 512x512 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.

-
Fun and interesting! Can you post some sample data? Is it a requirement to have a nanometer scale on the projection? Thinking maybe it's easier to scale up, though it's "cheating" a bit. Btw I guess you can come a long way with GDAL/ogr, although the projection issue will still need to be addressed. gdal.org/gdal_grid.html – alexanno Jun 2 '13 at 21:30
Thanks! Guess this is more a comment than an answer. As far as nanometer scale goes I'd say that whatever works is great, but a true nanoscale peojection would be the coolest. When it comes to sample data I'll have to check if there are some restrictions, but it's basically a 2 dim matrix of int16 values. – atlefren Jun 2 '13 at 22:36
Why would you need proj4 parameters? You won't transfer this data into another coordinate system (esp. geographical). Imho you should be able to do all your analysis without any coordinate system. What do you understand by a DEM? There are several types like triangulated surfaces(e.g. do a delauny triangulation) or raster maps(you have this already). This of course depends highly on the analysis software. Of course can you create other maps if they are needed as output to understand the probe. You might take a look at code.google.com/p/mtex for grain anlysis. – mistapink Jun 2 '13 at 23:40
The reason I (think) i need a CRS is this: if i just create, say a GeoTIFF without assigning a CRS, the unit of measure will be pixels. What if I want to measure distances? And, what if I have two AFM scans and know how they relate to each other (in terms of scale and offset from some point). Having a CRS assigned would facilitate the viewing several scans at once – atlefren Jun 3 '13 at 8:25
I usually cope with such data by setting a local coordinate system (with an origin coincident with the image origin) and "lying" about the units. For instance, you could stipulate that the units are kilometers when they are really nanometers, which makes it easy mentally to convert back and forth. Of course you won't be doing any reprojecting, so that's not an issue. Setting up this coordinate system is the same as georeferencing any DEM; it can be as simple as creating an unrotated world file. – whuber Jun 3 '13 at 13:11

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

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!

-
One interesting challenge could be to compare distances on the DEM to distances between places on the globe to give viewers an idea of just how small nanoscale is. E.g. htwins.net/scale2 – blah238 Jun 4 '13 at 9:15