[gdal-dev] SRTM maps and altitude...

[John Donovan]

Here's a little Python script I knocked up for just this purpose.

# Quick and dirty utility to get the value of a pixel in a georeferenced
# image file, given the map coordinates. Uses nearest-neighbour and doesn't
# try to do anything clever.
# By John Donovan 2010, released into the public domain. Use at your own risk.
# There is absolutely no warranty expressed or implied, and no guarantee that
# it is fit for purpose.
from osgeo import gdal
from osgeo.gdalconst import *
from numpy import matrix
import sys
import time

def main():
    filename = sys.argv[1]
    world = matrix([[float(sys.argv[2])], [float(sys.argv[3])]])

    ds = gdal.Open(filename, GA_ReadOnly)

    if ds is None:
        raise IOError

    try:
        # Transform from map space to image space.
        xform = ds.GetGeoTransform()
        offset = matrix([[xform[0]], [xform[3]]])
        A = matrix([[xform[1], xform[2]], [xform[4], xform[5]]])
        pixel = A.I * (world - offset)

        # Is the point within the bounding box of the image?
        if pixel[0] < 0.0 or pixel[0] > ds.RasterXSize or \
          pixel[1] < 0.0 or pixel[1] > ds.RasterYSize:
            raise IndexError('The coordinates specified are outside the '\
              'bounds of the image.')

        print int(pixel[0]), int(pixel[1]), \
          ds.ReadAsArray(int(pixel[0]), int(pixel[1]), 1, 1).ravel()

    finally:
        ds = None

if __name__ == '__main__':
    main()

It doesn't do any form of interpolation, nor does it take geoids into account, but running it on the SRTM3 version 2.1 data for the point you specified gives reasonable results:

$ python value_at.py N48E008.hgt 8.39774 48.02664
477 1168 [795]

Where the first two numbers are the integer pixel coordinates, and the numbers in brackets are the values, one for each band.

Hope this helps,
-JD

-----Original Message-----
From: gdal-dev-bounces at lists.osgeo.org [mailto:gdal-dev-bounces at lists.osgeo.org] On Behalf Of Gianni Olivieri
Sent: 11 June 2010 15:34
To: gdal-dev at lists.osgeo.org
Subject: [gdal-dev] SRTM maps and altitude...

Hallo,

    I need to extract the altitude of a GPS coordinate from a serie of HGT (SRTM or DEM) (I've download the entire Eurasia...

So, for example, If I have a directory with all htg file and I'd like to know the altitude of this coords:

48°01'35.89" N
8°23'51,86" E

how can I preceed?

The elevation in that point is: 761m (google earth result) I'd like to have a similar results.

Best regards



--
Gianni Olivieri
gianni.olivieri at gmail.com
_______________________________________________
gdal-dev mailing list
gdal-dev at lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

______________________________________________________________________
This email has been scanned for Virtalis Ltd by the MessageLabs Email Security System.
______________________________________________________________________

______________________________________________________________________
This email has been scanned for Virtalis Ltd by the MessageLabs Email Security System. ______________________________________________________________________