File: altitude_downloader.py

package info (click to toggle)
cyclograph 1.7.1-1
  • links: PTS
  • area: main
  • in suites: jessie, jessie-kfreebsd
  • size: 1,136 kB
  • ctags: 705
  • sloc: python: 9,682; makefile: 5
file content (195 lines) | stat: -rw-r--r-- 6,800 bytes parent folder | download | duplicates (2)
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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""altitude_downloader.py"""

# Copyright (C) 2008, 2009, 2010, 2011, 2012, 2013, 2014 Federico Brega, Pierluigi Villani

# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 3
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

import logging
import math
import threading
import urllib2
import xml.dom.minidom
import re

ALTITUDE_CACHE = {}
SERVICES = frozenset(("geonames.org", "earthtools.org", "usgs.net"))
_ACTIVE_SERVICE = "geonames.org"

def geonames_org(lat, lng):
    """ Use geonames.org STRM3 to download data"""
    url = 'http://api.geonames.org/srtm3?lat='+lat+'&lng='+lng+'&username=cyclograph'
    line = ''
    tries = 5
    # Sometimes geonames retuns an xml header, which is wrong, so try again.
    # the first character might be a minus sign, we keep it simple
    while not line.isdigit() and tries > 0:
        han = urllib2.urlopen(url)
        line = han.readline().strip()
        tries -= 1
        han.close()
    alt = float(line)
    return alt

def earthtools_org(lat, lng):
    """ Use earthtools.org to download data"""
    #http://www.earthtools.org/height/<latitude>/<longitude>
    url = 'http://www.earthtools.org/height/'+lat+'/'+lng
    han = urllib2.urlopen(url)
    alt = re.search(r"<meters>(-?\d+)<\/meters>",
                     han.read()).group(1)
    alt = float(alt)
    han.close()
    return alt

def usgs_net(lat, lng):
    """ Use usgs.net to download data"""
    url = 'http://gisdata.usgs.net/XMLWebServices/' + \
          'TNM_Elevation_Service.asmx/getElevation?X_Value=' + \
          lng + '&Y_Value=' + lat + \
          '&Elevation_Units=METERS&Source_Layer=-1&Elevation_Only=TRUE'
    han = urllib2.urlopen(url)
    match = re.search(r"<.+>(-?\d+[\.\d]*)<\/.+>",
                       han.read()).group(1)
    alt = float(match)
    return alt

SERVER_DOWNLOADER = {"geonames.org" : geonames_org,
                     "earthtools.org" : earthtools_org,
                     "usgs.net" : usgs_net
                     }

#WARNING! this module is not thread safe, maybe you can use a lock to use
#more threads.
def choose_service(newserv):
    """ Choose the service to use for retreiving altitude data"""
    global _ACTIVE_SERVICE
    global ALTITUDE_CACHE
    if newserv not in SERVICES:
        return -1

    if newserv != _ACTIVE_SERVICE:
        #Because different services may give different datas
        #we clear the cache if we change service.
        ALTITUDE_CACHE = {}
    _ACTIVE_SERVICE = newserv
    return 0

class ImporterThread(threading.Thread):
    """ Download altitude"""
    def __init__(self, outqueue, slopefile, num, start_dist, end_dist, seltype):
        threading.Thread.__init__(self)
        self.outq = outqueue
        self._want_abort = False
        self.slopefile = slopefile
        self.service = _ACTIVE_SERVICE
        self.status = 'OK'
        self.num_cps = num
        self.seltype = seltype
        self.start_dist = start_dist
        self.end_dist = end_dist

    def run(self):
        """ Task to execute while running"""
        if not SERVER_DOWNLOADER.has_key(self.service):
            self.status = "Error: service unknown"
            return
        getalt = SERVER_DOWNLOADER[self.service]


        self.slopefile.setcpToLoad(self.num_cps,
                                   self.start_dist,
                                   self.end_dist,
                                   self.seltype)
        cps = self.slopefile.cps
        usedcp = self.slopefile.cpToLoad

        for i in range(len(usedcp)):
            logging.debug("cpi:"+str(i)+"\tcp:"+str(cps[usedcp[i]]))
            (dist, alt, name, lat, lng) = cps[usedcp[i]]
            if ALTITUDE_CACHE.has_key((lat, lng)):
                # altitude has already been downloaded
                alt = ALTITUDE_CACHE[(lat, lng)]
            else:
                # try to download the altitude
                try:
                    alt = getalt(lat, lng)
                except urllib2.URLError:
                    self.status = 'Error: No network'
                    return
            ALTITUDE_CACHE[(lat, lng)] = alt
            progr = 1000 * i / (len(usedcp) - 1)
            self.outq.put((progr, (dist, alt, name)))
            if self._want_abort:
                self.status = 'Aborted'
                return

    def abort(self):
        """abort worker thread.
        Method for use by main thread to signal an abort.
        """
        self._want_abort = True

def point_conversion(lng, lat):
    """ return radians from coordinates"""
    return (math.radians(float(lng)), math.radians(float(lat)))

def distance(coords):
    """Calculates the distance from a list of coordinates."""
    (lng, lat) = coords[0]
    (lng_old, lat_old) = point_conversion(lng, lat)

    dist = 0.0
    cps = []
    for pnt in coords:
        if len(pnt) != 2:
            continue
        (lng, lat) = pnt
        (lng_new, lat_new) = point_conversion(lng, lat)
        dist += lambert_formulae(lat_old, lat_new, lng_old, lng_new)

        cps.append((dist, lng.strip('\n\t'), lat.strip('\n\t')))
        (lng_old, lat_old) = (lng_new, lat_new)
    return cps

R_EARTH = 6367.45 # km
def haversine_angle(lat_old, lat_new, lng_old, lng_new):
    """Haversine formula"""
    alpha = (math.sin((lat_old-lat_new)/2))**2 \
        + math.cos(lat_old) * math.cos(lat_new) * (math.sin((lng_old-lng_new)/2)**2)
    return 2 * math.atan2(math.sqrt(alpha), math.sqrt(1-alpha))

def haversine(lat_old, lat_new, lng_old, lng_new):
    return R_EARTH * haversine_angle(lat_old, lat_new, lng_old, lng_new)

R_EARTH_EQUATOR = 6378.137 #km WGS 84
r = 298.257223563 # inverse of eccentricity WGS 84

def lambert_formulae(lat_old, lat_new, lng_old, lng_new):
    red_lat_old = math.atan( (r - 1) / r * math.tan(lat_old) )
    red_lat_new = math.atan( (r - 1) / r * math.tan(lat_new) )
    
    sigma = haversine_angle(red_lat_old, red_lat_new, lng_old, lng_new)
    
    if sigma == 0.0:
        return 0.0
    
    P = (red_lat_old + red_lat_new) / 2
    Q = (red_lat_new - red_lat_old) / 2
    
    X = (sigma - math.sin(sigma)) * math.sin(P)**2 * math.cos(Q)**2 / ( math.cos(sigma/2)**2 )
    Y = (sigma + math.sin(sigma)) * math.cos(P)**2 * math.sin(Q)**2 / ( math.sin(sigma/2)**2 )
    
    distance = R_EARTH_EQUATOR * (sigma - (X + Y) / (2*r))
    return distance

# vim:sw=4:softtabstop=4:expandtab