Python – Convert DEM to 3D Array Organized as [Lat, Long, Elevation]

demelevationpythonrasterraster-conversion

I'm trying to convert a .tiff DEM to a numpy array using python. Ideally the array will store latitude, longitude (or some coordinate alternative), elevation).

I'm using Gdal but currently it only saves elevation values and the positional data is lost.

Current code:

import os, sys
import numpy as np
from osgeo import gdal, osr

gdal.AllRegister()
gdal.UseExceptions()

# Open Image
fpath = ('*location*\\DEM.tif')
inRas = gdal.Open(fpath)
if inRas is None:
  print ('Could not open image file')
  sys.exit(1)

# read in the crop data and get info about it
band1 = inRas.GetRasterBand(1)
rows = inRas.RasterYSize
cols = inRas.RasterXSize

cropData = band1.ReadAsArray(0,0,cols,rows)

Sample output from:

print(cropData.shape)
print(cropData.size)
numrows = len(cropData)
numcols = len(cropData[0])
for i in range(numrows):
  for j in range(numcols):
    print(cropData[i,j])

Best Answer

I'm going to assume the DEM you are working with is in a geographic projection. If it's in UTM or something, you can still do what you want, but it'll be a little more complex.

First, you need to get the transformation matrix from the GDAL dataset.

transform = inRas.GetGeoTransform()

xOrigin = transform[0] # x origin
yOrigin = transform[3] # y origin
pixWidth = transform[1] # pixel width. Assuming the source is in degrees, this is how many degrees per pixel
pixHeight = transform[5] # pixel height

So the lat/lon of any point in the array will be:

latitude = yOrigin + (row * pixHeight)
longitude = xOrigin + (col * pixWidth)

If the DEM is not in Geographic coordinates, you will need to read about coordinate transformations in GDAL.

Related Solutions

Python Raster Elevation – How to Get Elevation at Lat/Long Using Python

Here's a more programmatic way of using GDAL than @Aragon's answer. I've not tested it, but it is mostly boiler-plate code that has worked for me in the past. It relies on Numpy and GDAL bindings, but that's about it.

import osgeo.gdal as gdal
import osgeo.osr as osr
import numpy as np
from numpy import ma

def maFromGDAL(filename):
    dataset = gdal.Open(filename, gdal.GA_ReadOnly)

    if dataset is None:
        raise Exception()

    # Get the georeferencing metadata.
    # We don't need to know the CRS unless we want to specify coordinates
    # in a different CRS.
    #projection = dataset.GetProjection()
    geotransform = dataset.GetGeoTransform()

    # We need to know the geographic bounds and resolution of our dataset.
    if geotransform is None:
        dataset = None
        raise Exception()

    # Get the first band.
    band = dataset.GetRasterBand(1)
    # We need to nodata value for our MaskedArray later.
    nodata = band.GetNoDataValue()
    # Load the entire dataset into one numpy array.
    image = band.ReadAsArray(0, 0, band.XSize, band.YSize)
    # Close the dataset.
    dataset = None

    # Create a numpy MaskedArray from our regular numpy array.
    # If we want to be really clever, we could subclass MaskedArray to hold
    # our georeference metadata as well.
    # see here: http://docs.scipy.org/doc/numpy/user/basics.subclassing.html
    # For details.
    masked_image = ma.masked_values(image, nodata, copy=False)
    masked_image.fill_value = nodata

    return masked_image, geotransform

def pixelToMap(gt, pos):
    return (gt[0] + pos[0] * gt[1] + pos[1] * gt[2],
            gt[3] + pos[0] * gt[4] + pos[1] * gt[5])

# Reverses the operation of pixelToMap(), according to:
# https://en.wikipedia.org/wiki/World_file because GDAL's Affine GeoTransform
# uses the same values in the same order as an ESRI world file.
# See: http://www.gdal.org/gdal_datamodel.html
def mapToPixel(gt, pos):
    s = gt[0] * gt[4] - gt[3] * gt[1]
    x = (gt[4] * pos[0] - gt[1] * pos[1] + gt[1] * gt[5] - gt[4] * gt[2]) / s
    y = (-gt[3] * pos[0] + gt[0] * pos[1] + gt[3] * gt[2] - gt[0] * gt[5]) / s
    return (x, y)

def valueAtMapPos(image, gt, pos):
    pp = mapToPixel(gt, pos)
    x = int(pp[0])
    y = int(pp[1])

    if x < 0 or y < 0 or x >= image.shape[1] or y >= image.shape[0]:
        raise Exception()

    # Note how we reference the y column first. This is the way numpy arrays
    # work by default. But GDAL assumes x first.
    return image[y, x]

try:
    image, geotransform = maFromGDAL('myimage.tif')
    val = valueAtMapPos(image, geotransform, (434323.0, 2984745.0))
    print val
except:
    print('Something went wrong.')

[GIS] Creating nanoscale DEM with GDAL

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!

Best Answer

Related Solutions

Python Raster Elevation – How to Get Elevation at Lat/Long Using Python

[GIS] Creating nanoscale DEM with GDAL

Related Question