Python: Convert DEM to 3D array organized as [lat, long, elevation]

Question

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])

Answer 1

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.