diff --git a/__init__.py b/__init__.py new file mode 100644 index 0000000..52508f1 --- /dev/null +++ b/__init__.py @@ -0,0 +1,39 @@ +#!/usr/bin/env python +# -*- coding: UTF-8 -*- +""" +geo (Geographic) +================ + +**Author:** + +* Dirk Alders + +**Description:** + + This Module support functionalities around geographic issues. + +**Submodules:** + +* :mod:`geo.gps` +* :mod:`geo.osm` +* :mod:`geo.sun` + +**Unittest:** + + See also the :download:`unittest <../../geo/_testresults_/unittest.pdf>` documentation. +""" +__DEPENDENCIES__ = [] + +from geo import gps +from geo import osm +from geo import sun + +logger_name = 'GEO' + +__DESCRIPTION__ = """The Module {\\tt %s} is designed to \\ldots. +For more Information read the sphinx documentation.""" % __name__.replace('_', '\_') +"""The Module Description""" +__INTERPRETER__ = (2, ) +"""The Tested Interpreter-Versions""" + +__all__ = ['gps', 'osm', 'sun'] diff --git a/gps.py b/gps.py new file mode 100644 index 0000000..2c497b7 --- /dev/null +++ b/gps.py @@ -0,0 +1,598 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +""" +geo.gps (Geo-Positioning) +========================= + +**Author:** + +* Dirk Alders + +**Description:** + + This module is a submodule of :mod:`geo` and includes functions and classes for geographic issues (e.g. coordinate, area, tracks, ...). + +**Contentlist:** + +* :class:`geo.gps.area` +* :class:`geo.gps.coordinate` +* :class:`geo.gps.tracklist` +* :class:`geo.gps.track` + +**Unittest:** + + See also the :download:`unittest <../../geo/_testresults_/unittest.pdf>` documentation. +""" + +import calendar +import math +import time +import xml.parsers.expat + + +class area(object): + """ + :param coord1: Corner Coordinate 1 + :type coord1: coordinate + :param coord2: Corner Coordinate 2 + :type coord2: coordinate + + Class to store a geographic area and support some calculations. + + **Example:** + + .. code-block:: python + + >>> import geo + + >>> ar = geo.gps.area(...) + """ + def __init__(self, coord1, coord2): + min_lon = min(coord1[coordinate.LONGITUDE], coord2[coordinate.LONGITUDE]) + max_lon = max(coord1[coordinate.LONGITUDE], coord2[coordinate.LONGITUDE]) + min_lat = min(coord1[coordinate.LATITUDE], coord2[coordinate.LATITUDE]) + max_lat = max(coord1[coordinate.LATITUDE], coord2[coordinate.LATITUDE]) + self.coord1 = coordinate(lon=min_lon, lat=min_lat) + self.coord2 = coordinate(lon=max_lon, lat=max_lat) + + def _max_lat(self): + return self.coord2[coordinate.LATITUDE] + + def _max_lon(self): + return self.coord2[coordinate.LONGITUDE] + + def _min_lat(self): + return self.coord1[coordinate.LATITUDE] + + def _min_lon(self): + return self.coord1[coordinate.LONGITUDE] + + def center_pos(self): + """ + .. warning:: Needs sphinx documentation! + """ + clon = (self.coord1[coordinate.LONGITUDE] + self.coord2[coordinate.LONGITUDE]) / 2 + clat = (self.coord1[coordinate.LATITUDE] + self.coord2[coordinate.LATITUDE]) / 2 + return coordinate(lon=clon, lat=clat) + + def coordinate_in_area(self, coord): + """ + .. warning:: Needs sphinx documentation! + """ + lon = coord[coordinate.LONGITUDE] + lat = coord[coordinate.LATITUDE] + return lon >= self._min_lon() and lon <= self._max_lon() and lat >= self._min_lat() and lat <= self._max_lat() + + def corner_coordinates(self): + """ + .. warning:: Needs sphinx documentation! + """ + return self.coord1, self.coord2 + + def extend_area(self, coord): + """ + .. warning:: Needs sphinx documentation! + """ + if coord[coordinate.LONGITUDE] < self.coord1[coordinate.LONGITUDE]: + self.coord1[coordinate.LONGITUDE] = coord[coordinate.LONGITUDE] + elif coord[coordinate.LONGITUDE] > self.coord2[coordinate.LONGITUDE]: + self.coord2[coordinate.LONGITUDE] = coord[coordinate.LONGITUDE] + + if coord[coordinate.LATITUDE] < self.coord1[coordinate.LATITUDE]: + self.coord1[coordinate.LATITUDE] = coord[coordinate.LATITUDE] + elif coord[coordinate.LATITUDE] > self.coord2[coordinate.LATITUDE]: + self.coord2[coordinate.LATITUDE] = coord[coordinate.LATITUDE] + + def osm_map(self, map_code): + """ + :param map_code: Map code as defined in :class:`geo.osm` (e.g. :class:`geo.osm.MAP_STANDARD`) + + This Method returns a :class:`geo.osm.map` instance. + + .. warning:: Needs sphinx documentation! + """ + # TODO: needs to be implemented + pass + + def __str__(self): + return "%s / %s" % self.coordinates() + + +class coordinate(dict): + """ + :param lon: Londitude + :type lon: float + :param lat: Latitude + :type lat: float + :param height: Height + :type height: float + :param time: Time (Seconds since 1970) + :type time: int + + Class to store a geographic coodinate and support some calculations. + + **Example:** + + .. code-block:: python + + >>> import geo + + >>> ab = geo.gps.coordinate(lat=49.976596,lon=9.1481443) + >>> gb = geo.gps.coordinate(lat=53.6908298,lon=12.1583252) + >>> ab.dist_to(gb) / 1000 + 462.3182843470017 + >>> ab.angle_to(gb) / math.pi * 180 + 39.02285256685333 + """ + LATITUDE = 'lat' + LONGITUDE = 'lon' + HIGHT = 'hight' + TIME = 'time' + + def __init__(self, **kwargs): + dict.__init__(self, **kwargs) + + def __str__(self): + def to_string(lon_or_lat, plus_minus=('N', 'S')): + degrees = int(lon_or_lat) + lon_or_lat -= degrees + minutes = lon_or_lat * 60 + pm = 0 if degrees >= 0 else 1 + return "%d°%.4f%s" % (abs(degrees), abs(minutes), plus_minus[pm]) + lon = self.get(self.LONGITUDE) + lat = self.get(self.LATITUDE) + if lon is not None and lat is not None: + return to_string(lat) + ' ' + to_string(lon, ['E', 'W']) + else: + return None + + def angle_to(self, coord): + """ + This Method calculates the geographic direction in radiant from this to the given coordinate. + + .. note:: North is 0 (turning right). That means east is :class:`math.pi`/2. + + :param coord: Target coordinate. + :type coord: corrdinate + :returns: The geographic direction in radiant. + :rtype: int or float + """ + lat1 = coord[self.LATITUDE] + lon1 = coord[self.LONGITUDE] + lat2 = self[self.LATITUDE] + lon2 = self[self.LONGITUDE] + if lat1 is not None and lat2 is not None and lon1 is not None and lon2 is not None: + dlon = lon1 - lon2 + dlat = lat1 - lat2 + + if dlat > 0: + # case (half circle north) + angle = math.atan(dlon / dlat) + pass + elif dlat < 0: + # case (half circle south) + angle = math.pi + math.atan(dlon / dlat) + elif dlon > 0: + # case (east) + angle = math.pi / 2 + elif dlon < 0: + # case (west) + angle = math.pi * 3 / 2 + else: + # same point + return None + + if angle < 0: + angle += 2 * math.pi + return angle + else: + return None + + def dist_to(self, coord): + """ + This Method calcultes the distance from this coordinate to a given coordinate. + + :param coord: Target coordinate. + :type coord: coordinate + :return: The distance between two coordinates in meters. + :rtype: int or float + :raises: - + """ + lat1 = coord[self.LATITUDE] + lon1 = coord[self.LONGITUDE] + lat2 = self[self.LATITUDE] + lon2 = self[self.LONGITUDE] + if lat1 is not None and lat2 is not None and lon1 is not None and lon2 is not None: + R = 6378140 + dLat = math.radians(lat2 - lat1) + dLon = math.radians(lon2 - lon1) + a = math.sin(dLat / 2) * math.sin(dLat / 2) + math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) * math.sin(dLon / 2) * math.sin(dLon / 2) + c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a)) + return R * c + else: + return None + + +class tracklist(list): + """ + Class to store a a list of tracks and parse xml files created by a navigation system like Etrax Vista. + + **Example:** + + .. code-block:: python + + >>> import geo + + >>> ... + """ + def __init__(self): + list.__init__(self) + self.__xml_track_on_read = None + self.__xml_data_on_read = '' + + def __xml_start_element(self, name, attrs): + self.__xml_data_on_read = '' + if name == 'trk': + # new track found in file + self.__xml_track_on_read = track() + elif name == 'trkpt': + # new waypoint to append + if 'lon' in attrs and 'lat' in attrs: + self.__xml_track_on_read.append(coordinate(lon=float(attrs['lon']), lat=float(attrs['lat']))) + + def __xml_end_element(self, name): + if name == 'trk': + if self.__xml_track_on_read is not None and len(self.__xml_track_on_read) > 0: + self.append(self.__xml_track_on_read) + self.__xml_track_on_read = None + elif name == 'name': + if self.__xml_data_on_read != '': + self.__xml_track_on_read.set_name(self.__xml_data_on_read) + elif name == 'ele': + c = self.__xml_track_on_read[len(self.__xml_track_on_read) - 1] + c[c.HIGHT] = float(self.__xml_data_on_read) + elif name == 'time': + c = self.__xml_track_on_read[len(self.__xml_track_on_read) - 1] + c[c.TIME] = int(calendar.timegm(time.strptime(self.__xml_data_on_read, '%Y-%m-%dT%H:%M:%SZ'))) + + def __xml_char_data(self, data): + self.__xml_data_on_read += data + + def load_from_file(self, xmlfilehandle): + """ + .. warning:: Needs to be documented + """ + # TODO: implement either usage of a filename or filehandle + # parse xml-handle + p = xml.parsers.expat.ParserCreate() + p.StartElementHandler = self.__xml_start_element + p.EndElementHandler = self.__xml_end_element + p.CharacterDataHandler = self.__xml_char_data + p.ParseFile(xmlfilehandle) + + +class track(list): + """ + Class to store a a tracks and support some calculations. + + **Example:** + + .. code-block:: python + + >>> import geo + + >>> ... + """ + def __init__(self): + self._name = None + list.__init__(self) + self.__init_state_variables() + + def __init_state_variables(self): + self._area = None + self._average_speed = None + self._hightcharacteristic = None + self._end_date = None + self._optimized_track = None + self._passed_hight = None + self._speedcharacteristic = None + self._start_date = None + self._total_distance = None + self._total_time = None + + def append(self, coord): + """ + .. warning:: Needs to be documented + """ + list.append(self, coord) + self.__init_state_variables() + + def extend(self, *args, **kwargs): + """ + .. warning:: Needs to be documented + """ + self.__init_state_variables() + return list.extend(self, *args, **kwargs) + + def insert(self, index, coord): + """ + .. warning:: Needs to be documented + """ + list.insert(self, index, coord) + self.__init_state_variables() + + def set_name(self, name): + """ + .. warning:: Needs to be documented + """ + self._name = name + + def area(self): + """ + :rtype: geo.gps.area or None + + .. warning:: Needs to be documented + """ + if self._area is None: + if len(self) > 1: + self._area = area() + for c in self: + self._area.extend_area(c) + return self._area + + def average_speed(self): + """ + .. warning:: Needs to be documented + """ + if self._average_speed is None: + self._average_speed = self.total_distance() / self.total_time() + return self._average_speed + + def hightcharacteristic(self): + """ + .. warning:: Needs to be documented + """ + if self._hightcharacteristic is None: + pass # TODO: implement functionality + return self._hightcharacteristic + + def end_date(self): + """ + .. warning:: Needs to be documented + """ + if self._end_date is None: + if len(self) > 0: + self._end_date = self[len(self) - 1].get(coordinate.TIME) + return self._end_date + + def name(self): + """ + .. warning:: Needs to be documented + """ + return self._name + + def optimized_track(self): + """ + .. warning:: Needs to be documented + """ + # TODO: REWORK TRACK OPTIMIZATION + DIST_MOVED = 15. # 15m + ACCELERATION_MAX = 0.5 * 9.81 # 0.5g + ANGLE_DIF_FOR_BACKWARDS = 10 # +/- 5deg + MAX_DELETED_POINTS = 30 + + def backwards_direction(l_angle, t_angle): + dang = l_angle - t_angle + if dang < 0: + dang += math.degrees(math.pi * 2) + if dang > math.degrees(math.pi) - ANGLE_DIF_FOR_BACKWARDS / 2 and dang < math.degrees(math.pi) + ANGLE_DIF_FOR_BACKWARDS / 2: + return True + else: + return False + + if self._optimized_track is None: + self._optimized_track = track() + self._optimized_track.set_name(self._name) + del_lst = [] + for coord in self: + if len(self._optimized_track) == 0: + # first item has to be added always + self._optimized_track.append(coord) + else: + last = self._optimized_track[-1:][0] + #try: + acc = coord.dist_to(last) / ((coord[coordinate.TIME] - last[coordinate.TIME]) ** 2) + #except: + # acc = 0.0 + if len(self._optimized_track) > 1: + # calculate last angle + last_angle = self._optimized_track[-2:-1][0].angle_to(last) + if last_angle is not None: + last_angle = math.degrees(last_angle) + # calculate this angle + this_angle = last.angle_to(coord) + if this_angle is not None: + this_angle = math.degrees(this_angle) + else: + last_angle = this_angle = None + if coord.dist_to(last) > DIST_MOVED: + # distance ok. + if acc < ACCELERATION_MAX: + # acceleration ok. + if this_angle is None or last_angle is None or not backwards_direction(last_angle, this_angle): + # direction ok + del_lst = [] + self._optimized_track.append(coord) + else: + del_lst.append(coord) + #print "this one was in backwards direction (%d)!" % (len(del_lst)) + #print " %.2fkm: last = %.1f\xc2\xb0 - this = %.1f\xc2\xb0" % (self.total_distance()/1000., last_angle, this_angle) + else: + del_lst.append(coord) + #print "this one was with to high acceleration (%d)!" % (len(del_lst)) + #print " %.2fkm: %.1fm/s\xc2\xb2" % (self.total_distance()/1000., acc) + if len(del_lst) >= MAX_DELETED_POINTS: + print("Moeglicherweise ist die Optimierung des Tracks noch nicht ausgereift genug.") + print(" Bei %.1f km gab es %d Koordinaten die aussortiert worden waeren. Optimierung ausgelassen." % (self.get_total_distance() / 1000., MAX_DELETED_POINTS)) + self._optimized_track.extend(del_lst) + del_lst = [] + return self._optimized_track + + def osm_map(self, map_code): + return self.area().osm_map(map_code) + + def passed_hight(self): + """ + .. warning:: Needs to be documented + """ + if self._passed_hight is None: + if len(self) > 0: + self._passed_hight = 0.0 + hight = self[0][coordinate.HIGHT] + if hight is not None: + for c in self: + last_hight = hight + # hysteresis of 1 meter + hightlist = [c[coordinate.HIGHT] - 1, hight, c[coordinate.HIGHT] + 1] + hightlist.sort() + hight = hightlist[1] + if hight > last_hight: + self._passed_hight += hight - last_hight + return self._passed_hight + + def speedcharacteristic(self): + """ + .. warning:: Needs to be documented + """ + if self._speedcharacteristic is None: + pass # TODO: implement functionality + return self._speedcharacteristic + + def start_date(self): + """ + .. warning:: Needs to be documented + """ + if self._start_date is None: + if len(self) > 0: + self._start_date = self[0].get(coordinate.TIME) + return self._start_date + + def total_distance(self): + """ + .. warning:: Needs to be documented + """ + if self._total_distance is None: + if len(self) > 0: + self._total_distance = 0. + for i in range(0, len(self) - 1): + self._total_distance += self[i].dist_to(self[i + 1]) + return self._total_distance + + def total_time(self): + """ + .. warning:: Needs to be documented + """ + if self._total_time is None: + try: + self._total_time = self.end_date() - self.start_date() + except TypeError: + pass # doing nothing will return None as needed, if calculation is not possible + return self._total_time + + +''' +class gpxmanipu(): + debug=False + class trackmanipu(): + def __init__(self): + self.lines=[] + def AddLine(self, line): + self.lines.append(line) + def SetName(self, name): + SEARCHFORSTARTTAG=0 + SEARCHFORENDTAG=1 + state=SEARCHFORSTARTTAG + for i in range(0, len(self.lines)): + if state==SEARCHFORSTARTTAG: + if self.lines[i].find('')>0: + newline=self.lines[i][:self.lines[i].find('')+6] + newline+=name + state=SEARCHFORENDTAG + if state==SEARCHFORENDTAG: + if self.lines[i].find('')>0: + self.lines[i]=newline+self.lines[i][self.lines[i].find(''):] + else: + self.lines.remove(self.lines[i]) + def GetTrack(self): + rv='' + for line in self.lines: + rv+=line + return rv + def __init__(self, fh): + HEADERSEARCH=0 + TRACKSEARCH=1 + state=HEADERSEARCH + self.header='' + self.tracks=[] + for line in fh: + if state==HEADERSEARCH: + if line.lstrip().startswith(''): + state=TRACKSEARCH + track=self.trackmanipu() + track.AddLine(line) + else: + self.header+=line + elif state==TRACKSEARCH: + if line.lstrip().startswith(''): + track.AddLine(line) + self.tracks.append(track) + track=self.trackmanipu() + else: + track.AddLine(line) + self.footer=track.GetTrack() # This was no track + del(track) + if self.debug: + print "Header found:" + print self.header[:20]+' ... '+self.header[-20:] + print str(len(self.tracks))+' tracks found:' + for track in self.tracks: + track=track.GetTrack() + print track[:20]+' ... '+track[-20:] + print "Footer found:" + print self.footer + def GetGpx(self): + rv=self.header + for track in self.tracks: + rv+=track.GetTrack() + rv+=self.footer + return rv + def SetTrackname(self, number, name): + if len(self.tracks)>number-1: + self.tracks[number].SetName(name) + def DeleteTrack(self, number): + if len(self.tracks)>number-1: + return self.tracks.pop(number) + else: + return None +''' diff --git a/osm.py b/osm.py new file mode 100644 index 0000000..e820bf6 --- /dev/null +++ b/osm.py @@ -0,0 +1,628 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +""" +geo.osm (Open Streetmap) +======================== + +**Author:** + +* Dirk Alders + +**Description:** + + This module is a submodule of :mod:`geo` and supports functions and classes for Open Streetmap. + +**Contentlist:** + +* :func:`geo.osm.landmark_link` + +**Unittest:** + + See also the :download:`unittest <../../geo/_testresults_/unittest.pdf>` documentation. +""" + + +MAP_STANDARD = 'N' +"""MAP definition for Standard Map""" +MAP_LOCAL_TRAFIC = 'TN' +"""MAP definition for Local Trafic Map""" +MAP_CYCLEMAP = 'CN' +"""MAP definition for Cyclemap""" +MAP_HUMANITARIAN = 'HN' +"""MAP definition for Humanitarian Map""" + + +def landmark_link(coord, zoom_level=13, map_code=MAP_STANDARD): + """ + :param coord: Target coordinate. + :type coord: geo.gps.coordinate + :param zoom_level: The zoom level of the map (see https://wiki.openstreetmap.org/wiki/Zoom_levels for more information) + :type zoom_level: int + :param map_code: One of the map_codes :class:`MAP_STANDARD`, :class:`MAP_LOCAL_TRAFIC`, :class:`MAP_CYCLEMAP`, :class:`MAP_HUMANITARIAN`. + :type map_code: str + :return: An openstreetmap-url for marking a position in a map. + :rtype: str + + This Method generates an openstreetmap-url for marking a position in a map. + + .. code-block:: python + + >>> import geo + + >>> gb = geo.gps.coordinate(lat=53.6908298,lon=12.1583252) + >>> geo.osm.landmark_link(gb) + 'http://www.openstreetmap.org?mlat=53.690830&mlon=12.158325&zoom=13&layers=N' + """ + lon = coord[coord.LONGITUDE] + lat = coord[coord.LATITUDE] + # + if lon is not None and lat is not None: + link = 'http://www.openstreetmap.org?mlat=%(' + coord.LATITUDE + ')f&mlon=%(' + coord.LONGITUDE + ')f&zoom=%(zoom)d&layers=%(map)s' + return link % {coord.LATITUDE: lat, + coord.LONGITUDE: lon, + 'zoom': zoom_level, + 'map': map_code} + else: + return None + + +class map_spec(object): + def __init__(self, **kwargs): + coord1 = kwargs.get('coord1') + coord2 = kwargs.get('coord2') + zoom_level = kwargs.get('zoom_level') + x = kwargs.get('x') + y = kwargs.get('y') + + def get_resolution(self): + pass + + def get_coord_range(self): + pass + + def get_map(self): + pass + + def point_to_coord(self, xy): + pass + + def coord_to_point(self, coord): + pass + +class osm_map(object): + def __init__(self, **kwargs): + self.__map_spec__ = map_spec(kwargs) + map_code = kwargs.get('map_code') + + +''' +class get_from_web(): + """ + Class to download images from web. + """ + def __init__(self, props=None): + """ + Init routine for class get_from_web. + @param props: myapptools.AppProp instance with proxy information. This has to be a dictionary (see self._set_props). + """ + self.props = props + if props != None: + # install a callback if property 'Proxy' had been changed. it will be called with the new proxy information. + props.InstallPostSetCallback('Proxy', self._props_callback) + + def _props_callback(self, proxy): + """ + Routione which is called, if proxy information had been changed. It will set the new proxy information. + @param proxy: dictionary with proxy information + """ + self._set_props(**proxy) + + def _set_props(self, use_proxy, host, port, use_user, user, use_passwd, passwd=None): + """ + Routine to set the proxy information. + @param host: host to connect to + @param port: port which is used to connect the proxy + @param user: username for proxy + @param passwd: password to be used for user + """ + proxy_str = None + if not use_proxy: + proxy_support = urllib2.ProxyHandler({}) + opener = urllib2.build_opener(proxy_support, urllib2.HTTPHandler) + urllib2.install_opener(opener) + elif not use_user: + proxy_str = "http://%s:%d" % (host, port) + elif not use_passwd or passwd == None: + proxy_str = "http://%s@%s:%d" % (user, host, port) + else: + proxy_str = "http://%s:%s@%s:%d" % (user, passwd, host, port) + if proxy_str != None: + proxy_support = urllib2.ProxyHandler({"http": proxy_str}) + opener = urllib2.build_opener(proxy_support, urllib2.HTTPHandler) + urllib2.install_opener(opener) + + def get_image(self, url): + """ + Routine to download an image from web. + @param url: url to the image + @return: image (type: python Image) + """ + #try: + f = urllib2.urlopen(url) + im = Image.open(StringIO.StringIO(f.read())) + f.close() + return im + #except ???: + # print "exception in get_from_web().get_image('%s')" % url + # return None + + +class tile_handler(dict): + """ + Class to handle some tile_source classes and a default tile_source. + """ + TILE_SIZE = 256 + + class tile_source(): + """ + Class to get tile by cache or url. It also stores the downloaded tiles to the cache directory. + """ + TILEPATH = 'osm_tiles' + + def __init__(self, gfw, name, url, min_zoom, max_zoom): + """ + Init routine for class tile_source. + @param gfw: instance of get_from_web + @param name: name of tile_type + @param url: base url without zoom or tile information + @param min_zoom: minimum existing zoom level for this tile type + @param max_zoom: maximum existing zoom level for this tile type + """ + self.gfw = gfw + self._name = name + self._url = url + self._zooms = range(min_zoom, max_zoom + 1) + + def get_path(self, x, y, zoom_lvl): + """ + Routine to get the tile-path information for a specific tile. + @param x: horizontal tile number + @param y: vertical tile number + @param zoom_lvl: zoom level for the tile + @return: path to tile as string + """ + def intstr(num): + return str(int(num)) + return os.path.join(__basepath__, self.TILEPATH, self._name, intstr(zoom_lvl), intstr(x), intstr(y) + '.png') + + def get_url(self, x, y, zoom_lvl): + """ + Routine to get the url information for a specific tile. + @param x: horizontal tile number + @param y: vertical tile number + @param zoom_lvl: zoom level for the tile + @return: url to tile as string + """ + def intstr(num): + return str(int(num)) + return self._url + '/' + intstr(zoom_lvl) + '/' + intstr(x) + '/' + intstr(y) + '.png' + + def get_zooms(self): + """ + Routine to get a list of available zoom levels for this source. + @return: zoom levels as a list (e.g. [0,1,2,3]). + """ + return self._zooms + + def get_tile(self, x, y, zoom_lvl, max_age, from_cache): + """ + Routine to get a tile. + @param x: horizontal tile number + @param y: vertical tile number + @param zoom_lvl: zoom level for the tile + @param max_age: maximum age where no www-refresh is needed + @param from_cache: if True the tile is from cache + @return: tile as Image + """ + filename = self.get_path(x, y, zoom_lvl) + url = self.get_url(x, y, zoom_lvl) + if from_cache: + try: + return Image.open(filename) + except: + return None + else: + local_time = calendar.timegm(time.gmtime()) + try: + tile_tm = os.path.getmtime(filename) + except: + tile_tm = local_time - max_age - 1 + if local_time - tile_tm > max_age: # age depending refresh. + im = self.gfw.get_image(url) + try: + self.save(im, filename) + except: + print "exception in tile_handler().get_tile" + #TODO: exception handling. + pass + return im + else: + return None + + def save(self, im, filename): + """ + Routine to save the image to cache (directory). + @param im: image to save (type: python Image) + @param filename: name of the file, which will be created + """ + dirname = os.path.dirname(filename) + if not os.path.exists(dirname): + os.makedirs(dirname) + im.save(filename) + + def __init__(self, gfw, props=None): + """ + Init routine for class tile_handler + @param gfw: instance of get_from_web + @param props: myapptools.AppProp instance with tilehandler information. This has to be a dictionary (see self._set_proxy). + """ + dict.__init__(self) + self.gfw = gfw + self._active_tile_source = None + self._max_age = 3600 + self._append_tile_source(u'OSM-Mapnik', u'http://tile.openstreetmap.org', 0, 18) + self._append_tile_source(u'OSM-CycleMap', u'http://c.tile.opencyclemap.org/cycle', 0, 18) + try: + # install a callback if property 'Tilehandler' had been changed. it will be called with the new tilehandler information. + props.InstallPostSetCallback('Tilehandler', self._props_callback) + except: + #TODO: exception handling + pass + + def _props_callback(self, tilehandler): + """ + Routione which is called, if tilehandler information had been changed. It will set the new tilehandler information. + @param tilehandler: dictionary with tilehandler information + """ + self.set_props(**tilehandler) + + def set_props(self, source, max_age): + """ + Routine to set the proxy information. + @param source: source for tiles. + @param max_age: maximum age for a tile till it will be refreshed. + """ + self._set_default(source) + self._max_age = max_age + + def _append_tile_source(self, name, url, min_zoom, max_zoom): + """ + Routine to append a tilesource. + @param name: Name for this tilesource + @param url: URL for this tile source (without tile depending information e.g. zoom level, ...) + @param min_zoom: Minimum zoom level for this tilesource + @param max_zoom: Maximum zoom level for this tilesource + """ + self[name] = self.tile_source(self.gfw, name, url, min_zoom, max_zoom) + if self._active_tile_source == None: + self._set_default(name) + + def _set_default(self, name): + """ + Routine to set the default tilesorce (by name). + @param name: Name for the default tilesource. + @return: True if name was available, False if not. + """ + if name in self.keys(): + self._active_tile_source = name + return True + else: + return False + + def get_active_source(self): + """ + Routine to get the Name of the active tile source. + @return: name of the active tile source + """ + return self._active_tile_source + + def get_max_age(self): + return self._max_age + + def get_choices(self): + """ + Routine to get the names of the possible tile sources. + @return: list of possible tile sources + """ + return self.keys() + + def get_zooms(self): + """ + Routine to get a list of available zoom levels for this source. + @return: zoom levels as a list (e.g. [0,1,2,3]). + """ + return self[self._active_tile_source].get_zooms() + + def get_url(self, x, y, zoom_lvl): + """ + Routine to get the url information for a specific tile. + @param x: horizontal tile number + @param y: vertical tile number + @param zoom_lvl: zoom level for the tile + @return: url to tile as string + """ + return self[self._active_tile_source].get_url(x, y, zoom_lvl) + + def get_tile(self, x, y, zoom_lvl, from_cache): + """ + Routine to get a tile. + @param x: horizontal tile number + @param y: vertical tile number + @param zoom_lvl: zoom level for the tile + @param from_cache: if True the tile is from cache + @return: tile as Image + """ + return self[self._active_tile_source].get_tile(x, y, zoom_lvl, self._max_age, from_cache) + + def tile_num(self, coordinate, zoom): + """ + Routine which calculates the needed tile for coordinates. + @param coordinate: geo.coordinate instance with geographic information. + @param zoom: zoom information for the needed tile + @return: return a tuple of two float values (x- and y-tile) + """ + lat_rad = math.radians(coordinate[pylibs.geo.coordinate.LATITUDE]) + n = 2.0 ** zoom + xtile = (coordinate[pylibs.geo.coordinate.LONGITUDE] + 180.0) / 360.0 * n + ytile = (1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n + return (xtile, ytile) + + def coordinate(self, xtile, ytile, zoom): + """ + Routine which calculates geographic information out of tile information. + @param xtile: number of the tile (x) + @param ytile: number of the tile (y) + @param zoom: zoom level + @return: geo.coordinate instance with the geographic information + """ + n = 2.0 ** zoom + lon_deg = xtile / n * 360.0 - 180.0 + lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n))) + lat_deg = math.degrees(lat_rad) + return pylibs.geo.coordinate(lon=lon_deg, lat=lat_deg) + + +class osm_map(): + """ + This is an Image including an osm map + """ + def __init__(self, th, callback_refresh): + """ + Init routine for osm_map + @param th: tile_handler class needed to get tiles + @param callback_refresh: function to call, if osm_map had been changed (for refresh view) + callback gets two arguments + * the image or None if no image is available yet + * a description text (of the finished issue) + @param border: additional border around specified area + @param vers: refresh behaviour as described in VER_* + """ + self.th = th + self.callback_refresh = callback_refresh + self._del_map_values() + + self.running = False + self.stoprequest = False + + def _del_map_values(self): + """ + routine to reset (to None) map definitions (e.g. map_resolution, zoom_lvl, view range + """ + self._image = None + self._map_resolution = None + self._zoom_lvl = None + self._center_coordinate = None + + def _set_map_values(self, map_resolution, zoom_lvl, center_coordinate): + def zoom_limitation(th, zoom): + if zoom < min(th.get_zooms()): + return min(th.get_zooms()) + if zoom > max(th.get_zooms()): + return max(th.get_zooms()) + return zoom + self._image = Image.new('RGB', map_resolution, 'white') + self._map_resolution = map_resolution + self._zoom_lvl = zoom_limitation(self.th, zoom_lvl) + self._center_coordinate = center_coordinate + + def disable(self): + self.stop_now() + self.callback_refresh = None + + def get_image(self): + return self._image + + def get_map_resolution(self): + return self._map_resolution + + def get_zoom_lvl(self): + return self._zoom_lvl + + def get_center_coordinate(self): + return self._center_coordinate + + def _paste_tile(self, tile, xy): + """ + routine to paste a single tile at xy in the map image. + @param tile: tile to paste in the image + @param xy: position to paste the tile. also negative or too large values are + possible to paste just parts of a tile + """ + try: + self._image.paste(tile, xy) + except: + print "exception in osm_map()._paste_tile" + #TODO: exception handling + + def create_map_by_res_n_centercoord_n_zoomlvl(self, max_x_res, max_y_res, center_coordinate, zoom_lvl, cache_only=False): + """ + routine to ... + @param max_x_res: maximum x resolution + @param max_y_res: maximum y resolution + @param center_coordinate: center coordinates (object of geo.coordinates) + @param zoom_lvl: zoom_level to use for d + """ + if center_coordinate is not None and zoom_lvl is not None: + self._del_map_values() + # + # needed values for further calculations and map creation + # + self._set_map_values((max_x_res, max_y_res), zoom_lvl, center_coordinate) + self._create_map(cache_only) + + def create_map_by_coord_n_zoomlvl(self, tl_coord, br_coord, zoom_lvl, cache_only=False): + """ + @param tl_coord: top left coordinates (object of geo.coordinates) + @param br_coord: bottom right coordinates (object of geo.coordinates) + @param zoom_lvl: zoom_level to use for + """ + center_coordinate = pylibs.geo.area(tl_coord, br_coord).center_pos() + tl_tile = self.th.tile_num(tl_coord, zoom_lvl) + br_tile = self.th.tile_num(br_coord, zoom_lvl) + max_x_res = int((br_tile[0] - tl_tile[0]) * self.th.TILE_SIZE) + max_y_res = int((br_tile[1] - tl_tile[1]) * self.th.TILE_SIZE) + self.create_map_by_res_n_centercoord_n_zoomlvl(max_x_res, max_y_res, center_coordinate, zoom_lvl, cache_only) + + def create_map_by_res_n_coord(self, max_x_res, max_y_res, tl_coord, br_coord, cache_only=False): + """ + @param max_x_res: maximum x resolution + @param max_y_res: maximum y resolution + @param tl_coord: top left coordinates (object of geo.coordinates) + @param br_coord: bottom right coordinates (object of geo.coordinates) + + coord are not the used coordinated for the map corners, cause the zoom_lvl is quatisised + """ + def coordinates_in_map(max_x, max_y, p1, p2, center, zoom_lvl): + tl = self.get_coord_by_xy(0, 0, (max_x, max_y), center, zoom_lvl) + br = self.get_coord_by_xy(max_x, max_y, (max_x, max_y), center, zoom_lvl) + area = pylibs.geo.area(tl, br) + return area.coordinate_in_area(p1) and area.coordinate_in_area(p2) + center_coordinate = pylibs.geo.area(tl_coord, br_coord).center_pos() + zoom_lvl = max(self.th.get_zooms()) + while not coordinates_in_map(max_x_res, max_y_res, tl_coord, br_coord, center_coordinate, zoom_lvl): + zoom_lvl -= 1 + self.create_map_by_res_n_centercoord_n_zoomlvl(max_x_res, max_y_res, center_coordinate, zoom_lvl, cache_only) + + def stop_now(self): + self.stoprequest = True + while self.running: + pass + self.stoprequest = False + + def get_coord_by_xy(self, x, y, map_resolution=None, center_coordinate=None, zoom_lvl=None): + zoom_lvl = zoom_lvl or self._zoom_lvl + tl_tile = self._get_tl_tile_num(map_resolution, center_coordinate, zoom_lvl) + xy_tile = (tl_tile[0] + float(x) / self.th.TILE_SIZE, tl_tile[1] + float(y) / self.th.TILE_SIZE) + return self.th.coordinate(xy_tile[0], xy_tile[1], zoom_lvl) + + def get_xy_by_coord(self, coord, map_resolution=None, center_coordinate=None, zoom_lvl=None): + tl_tile = self._get_tl_tile_num(map_resolution, center_coordinate, zoom_lvl) + xy_tile = self.th.tile_num(coord, self._zoom_lvl) + x = int((xy_tile[0] - tl_tile[0]) * self.th.TILE_SIZE) + y = int((xy_tile[1] - tl_tile[1]) * self.th.TILE_SIZE) + return (x, y) + + def _get_map_res_tiles(self, map_resolution=None): + """ + returns the map resolution in number of tiles + """ + map_resolution = map_resolution or self._map_resolution + if map_resolution: + return (map_resolution[0] / float(self.th.TILE_SIZE), map_resolution[1] / float(self.th.TILE_SIZE)) + else: + return None + + def _get_tl_tile_num(self, map_resolution=None, center_coordinate=None, zoom_lvl=None): + map_resolution = map_resolution or self._map_resolution + center_coordinate = center_coordinate or self._center_coordinate + zoom_lvl = zoom_lvl or self._zoom_lvl + # + if (map_resolution and center_coordinate and zoom_lvl): + center_tile_num = self.th.tile_num(center_coordinate, zoom_lvl) + map_resolution_tiles = self._get_map_res_tiles(map_resolution) + topleft_tile_num = (center_tile_num[0] - map_resolution_tiles[0] / 2, center_tile_num[1] - map_resolution_tiles[1] / 2) + return topleft_tile_num + else: + return None + + def _get_br_tile_num(self, map_resolution=None, center_coordinate=None, zoom_lvl=None): + topleft_tile_num = self._get_tl_tile_num(map_resolution, center_coordinate, zoom_lvl) + map_resolution_tiles = self._get_map_res_tiles(map_resolution) + bottomright_tile_num = (topleft_tile_num[0] + map_resolution_tiles[0], topleft_tile_num[1] + map_resolution_tiles[1]) + return bottomright_tile_num + + def _get_xy_offset(self): + tl_tile = self._get_tl_tile_num() + x_offs = -int(tl_tile[0] % 1 * self.th.TILE_SIZE) + y_offs = -int(tl_tile[1] % 1 * self.th.TILE_SIZE) + return (x_offs, y_offs) + + def _get_tile_list(self): + tl_tile = self._get_tl_tile_num() + br_tile = self._get_br_tile_num() + tile_list = [] + for x in range(int(tl_tile[0]), int(br_tile[0]) + 1): + for y in range(int(tl_tile[1]), int(br_tile[1]) + 1): + tile_list.append((x, y, self._zoom_lvl)) + return tile_list + + def _create_map(self, cache_only): + """ + @param resoultion: map target resolution + @param xy_offset: offset for top left tile (normally <= 0) + @param zoom_lvl: tile zoom_lvl + @param tile_list: list of tiles [[x1, x2, x3], [y1, y2]] + @param description: description text for callback function + """ + def create_map_by_(xy_offset, tile_list, by_path): + # + # create map from already stored tiles (...by_path) + # + num_tiles = len(tile_list) + x0, y0 = tile_list[0][:2] + num = 0 + for x, y, z in tile_list: + num += 1 + if self.stoprequest: + break + tile = self.th.get_tile(x, y, z, by_path) + if tile != None: + # paste tile only if tile was available + pos = (xy_offset[0] + (x - x0) * self.th.TILE_SIZE, xy_offset[1] + (y - y0) * self.th.TILE_SIZE) + self._paste_tile(tile, pos) + if not by_path: + desc = "Tile " + self.th.get_url(x, y, z) + " added to map." + prog = float(num) / num_tiles + if self.callback_refresh: + self.callback_refresh(self._image, desc, prog) + self.running = True + #TODO: ggf. Klasse um Uebersetzungen zu ermoeglichen. + + create_map_by_(self._get_xy_offset(), self._get_tile_list(), by_path=True) + desc = 'Map creation from cache completeled.' + if self.callback_refresh: + self.callback_refresh(self._image, desc, 1.0) + if not cache_only: + create_map_by_(self._get_xy_offset(), self._get_tile_list(), by_path=False) + desc = 'Map creation completeled.' + if self.callback_refresh: + self.callback_refresh(self._image, desc, 1.0) + self.running = False + + +def show_map(image, description, progress): + print description, "%5.1f%%" % (progress * 100.) + if image != None: + image.show() +''' diff --git a/sun.py b/sun.py new file mode 100644 index 0000000..6a0ec06 --- /dev/null +++ b/sun.py @@ -0,0 +1,182 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +""" +geo.sun (Sun) +============= + +**Author:** + +* Dirk Alders +* Formula from Dr. Roland Brodbeck, Calsky (http://lexikon.astronomie.info/zeitgleichung/neu.html) +* First implementation by Alexander Klupp 2014-01-14 + +**Description:** + + This module is a submodule of :mod:`geo` and supports functions and classes for sun position. + +**Contentlist:** + +* :func:`geo.sun.sunset` +* :func:`geo.sun.sunrise` + +**Unittest:** + + See also the :download:`unittest <../../geo/_testresults_/unittest.pdf>` documentation. +""" + + +import calendar +import math +import time + +europe_berlin = {'lat': 52.520008, 'lon': 13.404954} + + +def JulianischesDatum(Jahr, Monat, Tag, Stunde, Minuten, Sekunden): + if (Monat <= 2): + Monat = Monat + 12 + Jahr = Jahr - 1 + Gregor = (Jahr / 400) - (Jahr / 100) + (Jahr / 4) # Gregorianischer Kalender + return 2400000.5 + 365 * Jahr - 679004 + Gregor + math.floor(30.6001 * (Monat + 1)) + Tag + Stunde / 24 + Minuten / 1440 + Sekunden / 86400 + + +def InPi(x): + n = int(x / (2 * math.pi)) + x = x - n * 2 * math.pi + if (x < 0): + x += 2 * math.pi + return x + + +def eps(T): + # Neigung der Erdachse + return math.pi / 180 * (23.43929111 + (-46.8150 * T - 0.00059 * T ** 2 + 0.001813 * T ** 3) / 3600) + + +def BerechneZeitgleichung(T): + RA_Mittel = 18.71506921 + 2400.0513369 * T + (2.5862e-5 - 1.72e-9 * T) * T ** 2 + M = InPi(2 * math.pi * (0.993133 + 99.997361 * T)) + L = InPi(2 * math.pi * (0.7859453 + M / (2 * math.pi) + (6893 * math.sin(M) + 72 * math.sin(2 * M) + 6191.2 * T) / 1296e3)) + e = eps(T) + RA = math.atan(math.tan(L) * math.cos(e)) + if (RA < 0): + RA += math.pi + if (L > math.pi): + RA += math.pi + RA = 24 * RA / (2 * math.pi) + DK = math.asin(math.sin(e) * math.sin(L)) + # Damit 0 <= RA_Mittel < 24 + RA_Mittel = 24.0 * InPi(2 * math.pi * RA_Mittel / 24.0) / (2 * math.pi) + dRA = RA_Mittel - RA + if (dRA < -12.0): + dRA += 24.0 + if (dRA > 12.0): + dRA -= 24.0 + dRA = dRA * 1.0027379 + return dRA, DK + + +def Sonnenauf_untergang(JD, Zeitzone, coord): + # Zeitzone = 0 #Weltzeit + # Zeitzone = 1 #Winterzeit + # Zeitzone = 2 #Sommerzeit + # JD = JulianischesDatum + + B = math.radians(coord.get('lat')) # geographische Breite Erkelenz + GeographischeLaenge = coord.get('lon') # geographische Laenge + + JD2000 = 2451545 + h = -50.0 / 60.0 * math.pi / 180 + T = (JD - JD2000) / 36525 + + Zeitgleichung, DK = BerechneZeitgleichung(T) + + Zeitdifferenz = 12 * math.acos((math.sin(h) - math.sin(B) * math.sin(DK)) / (math.cos(B) * math.cos(DK))) / math.pi + + AufgangOrtszeit = 12 - Zeitdifferenz - Zeitgleichung + UntergangOrtszeit = 12 + Zeitdifferenz - Zeitgleichung + AufgangWeltzeit = AufgangOrtszeit - GeographischeLaenge / 15 + UntergangWeltzeit = UntergangOrtszeit - GeographischeLaenge / 15 + + Aufgang = AufgangWeltzeit + Zeitzone + if (Aufgang < 0): + Aufgang += 24 + elif (Aufgang >= 24): + Aufgang -= 24 + + AM = round(Aufgang * 60) / 60 # minutengenau runden + + Untergang = UntergangWeltzeit + Zeitzone + if (Untergang < 0): + Untergang += 24 + elif (Untergang >= 24): + Untergang -= 24 + + UM = round(Untergang * 60) / 60 # minutengenau runden + + return AM, UM + + +def sunrise(coord=europe_berlin, date=None): + """ + :param coord: Target coordinate or None (default is central europe). + :type coord: geo.gps.coordinate + :param date: The day to calculate with or None (only year, month and day are relevant; default ist today) + :type date: time.struct_time + :return: The date and time information for the sunrise + :rtype: time.struct_time + + This Method calculates the time for sunrise for a given date and coordinate. + + .. code-block:: python + + >>> import geo + + >>> ... + """ + date = date or time.localtime() + + year, month, day = date[0:3] + dst = date[8] # Sommerzeit + + AM = Sonnenauf_untergang(JulianischesDatum(year, month, day, 12, 0, 0), dst + 1, coord)[0] + + tup_list = list(date) + tup_list[3] = 0 + tup_list[4] = 0 + tup_list[5] = 0 + tm = time.mktime(time.struct_time(tup_list)) + return time.localtime(tm + AM * 3600) + + +def sunset(coord=europe_berlin, date=None): + """ + :param coord: Target coordinate or None (default is central europe). + :type coord: geo.gps.coordinate + :param date: The day to calculate with or None (only year, month and day are relevant; default ist today) + :type date: time.struct_time + :return: The date and time information for the sunset + :rtype: time.struct_time + + This Method calculates the time for sunrise for a given date and coordinate. + + .. code-block:: python + + >>> import geo + + >>> ... + """ + date = date or time.localtime() + + year, month, day = date[0:3] + dst = date[8] # Sommerzeit + + UM = Sonnenauf_untergang(JulianischesDatum(year, month, day, 12, 0, 0), dst + 1, coord)[1] + + tup_list = list(date) + tup_list[3] = 0 + tup_list[4] = 0 + tup_list[5] = 0 + tm = time.mktime(time.struct_time(tup_list)) + return time.localtime(tm + UM * 3600)