[GIS] Using GDAL in python to read elevation at a point from dted files

dtedelevationgdalpython

How do I get the elevation at a certain point using gdal to read in a dted file? I don't see much useful information out on the web.

I'm sorry I should have been more specific. How would I accomplish this in python?

Best Answer

You can use the gdal.Band ReadAsArray method. See also the GDAL API tutorial and the example below:

from osgeo import gdal,ogr
import numpy

def map2pixel(mx,my,gt):
    """
    Convert from map to pixel coordinates.
    Only works for geotransforms with no rotation.
    """

    px = int((mx - gt[0]) / gt[1]) #x pixel
    py = int((my - gt[3]) / gt[5]) #y pixel

    return px,py

src_filename = '/path/to/dted/file'
src_ds=gdal.Open(src_filename) 
gt=src_ds.GetGeoTransform()
rb=src_ds.GetRasterBand(1)

#For a single XY coordinate (must be same projection as the raster)
x=123.456
y=12.345
px,py=map2pixel(x,y,gt)
val = rb.ReadAsArray(px,py,1,1) #Read a 1x1 array from the raster at px py

#If you want to do this for every feature in a point dataset:
shp_filename = '/tmp/test.shp'

ds=ogr.Open(shp_filename)
lyr=ds.GetLayer()
for feat in lyr:
    geom = feat.GetGeometryRef()
    mx,my=geom.GetX(), geom.GetY()  #coord in map units, must be same projection as the raster
    px,py=map2pixel(mx,my,gt)
    val=rb.ReadAsArray(px,py,1,1) #Read a 1x1 array from the raster at px py
Related Question