source: subversion/applications/utils/srtm2postgis/trunk/read_data.py @ 34941

Last change on this file since 34941 was 9393, checked in by sjors, 11 years ago

Move database activity into a seperate class.

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