public
Description: Imports data from the NASA Shuttle Radar Topography Mission into a PostGIS database (and other databases in the futue).
Homepage: http://wiki.openstreetmap.org/index.php/Route_altitude_profiles_SRTM#SRTM_import
Clone URL: git://github.com/Sjors/srtm2postgis.git
srtm2postgis / read_data.py
100644 73 lines (50 sloc) 2.176 kb
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
# Read srtm data files and put them in the database.
from osgeo import gdal, gdal_array
import os
 
import re
import zipfile
 
from data import util
 
# Main functions
 
def loadTile(continent, filename):
  # Unzip it
  zf = zipfile.ZipFile('data/' + continent + '/' + filename + ".hgt.zip")
  for name in zf.namelist():
    outfile = open('data/' + continent + '/' + name, 'wb')
    outfile.write(zf.read(name))
    outfile.flush()
    outfile.close()
  
  # Read it
  srtm = gdal.Open('data/' + continent + '/' + filename + '.hgt')
  
  # Clean up
  os.remove('data/' + continent + '/' + filename + '.hgt')
  
  return gdal_array.DatasetReadAsArray(srtm)
 
def connectToDatabasePsycopg2(database):
    conn = psycopg2.connect("dbname='" + database.db + "' host='localhost' user='" + database.db_user + "' password='" + database.db_pass + "'")
    return conn.cursor()
 
def posFromLatLon(lat,lon):
  return (lat * 360 + lon) * 1200 * 1200
 
def verify(db, number_of_tiles, files_hashes, continent, north, south, west, east):
    # For every tile, verify the bottom left coordinate.
    for file in files_hashes:
      # Strip .hgt.zip extension:
      file = file[1][0:-8]
    
      [lat,lon] = util.getLatLonFromFileName(file)
      
      # Only a smaller part of Australia (see below):
      if util.inBoundingBox(lat, lon, north, south, west, east):
      
        print "Verify " + file + "..."
    
        # Get top left altitude from file:
        coordinate_file = loadTile(continent, file)[1][0]
    
        # Get top left altitude from database:
        coordinate_db = db.fetchTopLeftAltitude(lat,lon)
        
        if coordinate_db != coordinate_file:
          print "Mismatch tile " + file[1]
          exit()
    
    # Check the total number of points in the database:
    
    print "Check the total number of points in the database..."
    
    sql = db.query("SELECT count(pos) FROM altitude")
    total = int(sql.getresult()[0][0])
    if not total == number_of_tiles * 1200 * 1200:
      print "Not all tiles have been (completely) inserted!"
      exit()
        
    print "All tiles seem to have made it into the database! Enjoy."
    
    exit()