diff --git a/.gitignore b/.gitignore index 09f38f1a..1e387bd8 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,3 @@ -python_program/Marschzeittabelle_Generiert.xlsx python_program/.idea .idea python_program/res/swisstlm3d_2021-04_2056_5728.gdb diff --git a/Readme.md b/Readme.md index 1ba66988..9b412b1c 100644 --- a/Readme.md +++ b/Readme.md @@ -31,8 +31,8 @@ Weitere Funktionalitäten sind in Planung und können gerne auch per Enhancement **Wichtig:** Die manuelle Planung bleibt ein grundlegender Bestandteil der Vorbereitung auf eine Wanderung. Dieses Projekt zielt lediglich auf die Beschleunigung mechanische, sich wiederholender Prozesse wie die Erstellung einer -Marschzeittabelle auf der Grundlage einer bestehenden Route. Dabei wird das sorgfältige Planen und -Durchdenken einer Aktivität in keinerlei Hinsicht ersetzt! +Marschzeittabelle auf der Grundlage einer bestehenden Route. Dabei wird das sorgfältige Planen und Durchdenken einer +Aktivität in keinerlei Hinsicht ersetzt! ## Glossar @@ -40,6 +40,7 @@ Begriff | Erklärung -------- | -------- | J+S | Jugend+Sport: Ist das Sportförder-Programm des Bundes für Kinder und Jugendsport in der Schweiz. Es unterstützt Sportkurse und Lager in rund 70 Sportarten und Disziplinen. | Marschzeittabelle | Tabelle zum Abschätzen der Marschzeit einer Wanderung. Enthält oft auch Angaben zu Pausen und ein Höhenprofil. J+S empfehlt für jede Wanderung eine Marschzeittabelle zu erstellen. | +SchweizMobil | Ermöglicht das Planen von Wanderungen und Velo-Touren der digitalen Landeskarte. Die magnetischen Wege beschleunigen den Prozess des Einzeichnens einer Route. # Generating J+S-Walk-Time-Tables diff --git a/python_program/.gitattributes b/python_program/.gitattributes deleted file mode 100644 index dfe07704..00000000 --- a/python_program/.gitattributes +++ /dev/null @@ -1,2 +0,0 @@ -# Auto detect text files and perform LF normalization -* text=auto diff --git a/python_program/.gitignore b/python_program/.gitignore index 8a07824e..382747fd 100644 --- a/python_program/.gitignore +++ b/python_program/.gitignore @@ -1,5 +1,6 @@ imgs/ testWalks/ +output/* # Byte-compiled / optimized / DLL files __pycache__/ diff --git a/python_program/Readme.md b/python_program/Readme.md index dcae842b..7fecd8cb 100644 --- a/python_program/Readme.md +++ b/python_program/Readme.md @@ -1,9 +1,18 @@ # How to run it? -Make sure you have installed python 3 and all requirements listed in the requirements.txt file. -Now you can run ```main.py``` to launch the script. +Make sure you have installed python 3 and all requirements listed in the requirements.txt file. Now you can +run ```main.py``` to launch the script. The produced files get saved in the ./output directory. +In the ```main.py``` you can specify the ```DEPARTURE_TIME```, ```GPX_FILE_PATH```, and ```VELOCITY```. Note: the script +is only tested with GPX-files exported form SchweizMobil, but it should work with arbitrary GPX-files. -## About SwissTopo Services +## About swisstopo Services -Request limits, see https://www.geo.admin.ch/de/geo-dienstleistungen/geodienste/terms-of-use.html \ No newline at end of file +This script highly depends on the services form swisstopo (Federal Office of Topography swisstopo). Those services are +free-of-charge and do not require registration (since 01.03.2021). However, we should follow there "Fair Use" rules. +More details see [API documentation](https://api3.geo.admin.ch/services/sdiservices.html) +and [Terms of use](https://www.geo.admin.ch/de/geo-dienstleistungen/geodienste/terms-of-use.html). + +Thus, in some cases, e.g. finding a swisstopo name for a coordinate, we will not use the official API but rebuild the +query algorithm locally. Therefore, please make sure, you have downloaded the necessary files in the ./res folder. The +files needed include ```swissNAMES3D_LIN.csv```, ```swissNAMES3D_PKT.csv```, and ```swissNAMES3D_PLY.csv```. \ No newline at end of file diff --git a/python_program/calculations.py b/python_program/calculations.py deleted file mode 100644 index 1e21db5a..00000000 --- a/python_program/calculations.py +++ /dev/null @@ -1,15 +0,0 @@ -def calcTime(delta_height, delta_dist, speed): - """ - - Calculates the walking time form one point to another - - for this calculation the basic formula form Jugend+Sport is used for preciser we could use the formula - form SchweizMobil or use more way points. But since we want to create a "normal" walk table as specified by - Jugend+Sport we use there basic formula - - """ - - if delta_height is None or delta_dist is None: - return 0 - - return (delta_dist + (delta_height / 100 if delta_height > 0 else 0)) / speed diff --git a/python_program/transformation.py b/python_program/coord_transformation.py similarity index 96% rename from python_program/transformation.py rename to python_program/coord_transformation.py index afbefd20..15c393df 100644 --- a/python_program/transformation.py +++ b/python_program/coord_transformation.py @@ -28,18 +28,22 @@ # Aaron Schmocker [aaron@duckpond.ch] # vim: tabstop=4 shiftwidth=4 softtabstop=4 expandtab -# Source: http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/projections.html (see PDFs under "Documentation") +# Source: http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/projections.html +# (see PDFs under "Documentation") # Updated 9 dec 2014 -# Please validate your results with NAVREF on-line service: http://www.swisstopo.admin.ch/internet/swisstopo/en/home/apps/calc/navref.html (difference ~ 1-2m) + +# Please validate your results with NAVREF on-line service: +# http://www.swisstopo.admin.ch/internet/swisstopo/en/home/apps/calc/navref.html (difference ~ 1-2m) import math class GPSConverter(object): - ''' + + """ GPS Converter class which is able to perform convertions between the CH1903 and WGS84 system. - ''' + """ # Convert CH y/x/h to WGS height def CHtoWGSheight(self, y, x, h): diff --git a/python_program/create_map.py b/python_program/create_map.py index bb9bd2f9..d7356048 100644 --- a/python_program/create_map.py +++ b/python_program/create_map.py @@ -1,18 +1,43 @@ import math from io import BytesIO +import gpxpy +import grequests import numpy as np -import requests from PIL import Image, ImageDraw +from python_program.coord_transformation import GPSConverter from python_program.gpx_utils import calc_perimeter -from python_program.transformation import GPSConverter TILE_SIZES = [64000, 25600, 12800, 5120, 2560, 1280, 640, 512, 384, 256, 128, 64, 25.6] +""" Tile width m, see https://api3.geo.admin.ch/services/sdiservices.html#wmts """ + ZOOM_LEVELS = [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28] +""" Zoom levels, see https://api3.geo.admin.ch/services/sdiservices.html#wmts """ + +TILE_WIDTH = 256 +""" Width of a tile in pixels """ + +TILE_MATRIX_SET: str = '2056' +""" EPSG code for LV03/CH1903 """ + +def plot_route_on_map(raw_gpx_data: gpxpy.gpx, + way_points: [], + file_name: str, + layer: str = 'ch.swisstopo.pixelkarte-farbe', + tile_format_ext: str = 'jpeg'): + """ + + Creates a map of the route and marking the selected way points on it. + + raw_gpx_data : raw data form imported GPX file + way_points : selected way points of the walk-time table + tile_format_ext : Format of the tile, allowed values jpeg or png, default jpeg + layer : Map layer, see https://wmts.geo.admin.ch/EPSG/2056/1.0.0/WMTSCapabilities.xml for options + + """ -def plot_route_on_map(raw_gpx_data, way_points): lv03_min, lv03_max = calc_perimeter(raw_gpx_data) # zoom level of the map snippets (a value form 0 to 12) @@ -21,47 +46,63 @@ def plot_route_on_map(raw_gpx_data, way_points): x_count = 0 y_count = 0 - while zoom_level < 12 and x_count * y_count < 64: + while zoom_level < 12 and x_count * y_count < 36: zoom_level = zoom_level + 1 # calc the number of the bottom left tile - x_tile = math.floor((lv03_min[0] - 420_000) / TILE_SIZES[zoom_level]) - 1 - y_tile = math.floor((350_000 - lv03_min[1]) / TILE_SIZES[zoom_level]) + 1 + tile_base_row = math.floor((lv03_min[0] - 420_000) / TILE_SIZES[zoom_level]) - 1 + tile_base_col = math.floor((350_000 - lv03_min[1]) / TILE_SIZES[zoom_level]) + 1 # calc number of tiles in each direction x_count = math.ceil((lv03_max[0] - lv03_min[0]) / TILE_SIZES[zoom_level]) + 2 y_count = math.ceil((lv03_max[1] - lv03_min[1]) / TILE_SIZES[zoom_level]) + 2 - lv03_min = (x_tile * TILE_SIZES[zoom_level] + 420_000, 350_000 - y_tile * TILE_SIZES[zoom_level]) - print(zoom_level, ': ', x_count, ' ', y_count) - - # creates the urls of the image tiles as a 3x3 grid around the centered tile - base_url = 'https://wmts.geo.admin.ch/1.0.0/ch.swisstopo.pixelkarte-farbe/default/current/2056/' + str( - ZOOM_LEVELS[zoom_level]) + '/' + lv03_min = (tile_base_row * TILE_SIZES[zoom_level] + 420_000, 350_000 - tile_base_col * TILE_SIZES[zoom_level]) # load tiles and combine map parts - card_snippet_as_image = Image.new("RGB", (254 * x_count, 254 * y_count)) + card_snippet_as_image = Image.new("RGBA", (TILE_WIDTH * x_count, TILE_WIDTH * y_count)) + + base_url = 'http://wmts.geo.admin.ch/1.0.0/' + layer + \ + '/default/current/' + TILE_MATRIX_SET + '/' + str(ZOOM_LEVELS[zoom_level]) + '/' + + # create urls + urls = [] + for i in range(0, x_count): + for j in range(1, y_count + 1): + urls.append(base_url + str(tile_base_row + i) + '/' + str(tile_base_col - j) + '.' + tile_format_ext) + + print('URL of a sample tile: ' + urls[0]) + + # Request images in parallel + rs = (grequests.get(u) for u in urls) + results = grequests.map(rs) + + # Since swisstopo services are free, we must guarantee to use the service fairly + # See https://www.geo.admin.ch/de/geo-dienstleistungen/geodienste/terms-of-use.html + if len(urls) > 300: + raise Exception('Fair use limit exceeded!') + index = 0 for i in range(0, x_count): for j in range(1, y_count + 1): - url = base_url + str(x_tile + i) + '/' + str(y_tile - j) + '.jpeg' - response = requests.get(url) + response = results[index] + index = index + 1 img = Image.open(BytesIO(response.content)) - img.thumbnail((254, 254), Image.ANTIALIAS) + img.thumbnail((TILE_WIDTH, TILE_WIDTH), Image.ANTIALIAS) w, h = img.size card_snippet_as_image.paste(img, (i * w, h * (y_count - j), i * w + w, h * (y_count - j) + h)) # mark point on the map draw = ImageDraw.Draw(card_snippet_as_image) - pixels_per_meter = (254.0 / TILE_SIZES[zoom_level]) + pixels_per_meter = (TILE_WIDTH / TILE_SIZES[zoom_level]) old_coords = None for track in raw_gpx_data.tracks: for segment in track.segments: for point in segment.points: wgs84_point = [point.latitude, point.longitude, point.elevation] - img_x, img_y = calc_img_cood(card_snippet_as_image.size, lv03_min, pixels_per_meter, wgs84_point) + img_x, img_y = calc_img_coord(card_snippet_as_image.size, lv03_min, pixels_per_meter, wgs84_point) if old_coords is not None: draw.line((old_coords, (img_x, img_y)), fill=(255, 165, 0), width=5) @@ -70,18 +111,19 @@ def plot_route_on_map(raw_gpx_data, way_points): for point in way_points: wgs84_point = [point[1].latitude, point[1].longitude, point[1].elevation] - img_x, img_y = calc_img_cood(card_snippet_as_image.size, lv03_min, pixels_per_meter, wgs84_point) + img_x, img_y = calc_img_coord(card_snippet_as_image.size, lv03_min, pixels_per_meter, wgs84_point) - circle_coords = (img_x - 18, img_y - 18, img_x + 18, img_y + 18) + pkt_rad = 10 + circle_coords = (img_x - pkt_rad, img_y - pkt_rad, img_x + pkt_rad, img_y + pkt_rad) draw.ellipse(circle_coords, outline=(255, 0, 0), width=5) # saves the image as '.jpg' - card_snippet_as_image.save('imgs/map.jpg') + card_snippet_as_image.save('output/' + file_name + '_map.png') card_snippet_as_image.show() -def calc_img_cood(image_size, lv03_min, pixels_per_meter, wgs84_point): +def calc_img_coord(image_size, lv03_min, pixels_per_meter, wgs84_point): converter = GPSConverter() lv03_point = np.round(converter.WGS84toLV03(wgs84_point[0], wgs84_point[1], wgs84_point[2])) diff --git a/python_program/find_name.py b/python_program/find_swisstopo_name.py similarity index 75% rename from python_program/find_name.py rename to python_program/find_swisstopo_name.py index bec787df..f2864040 100644 --- a/python_program/find_name.py +++ b/python_program/find_swisstopo_name.py @@ -21,11 +21,22 @@ def add_to_database(file, db, typeIndex, name, x, y): def find_name(coord, dist): + """ + See also https://api3.geo.admin.ch/api/faq/index.html#which-layers-have-a-tooltip + fair use limit 20 Request per minute + (see https://www.geo.admin.ch/de/geo-dienstleistungen/geodienste/terms-of-use.html) + + Therefore this process should be done locally. + + """ + flurname_list = [swiss_name for swiss_name in database if abs(swiss_name.x - coord[0]) < dist * 4 and abs(swiss_name.y - coord[1]) < dist * 4] + # print(list(map(lambda x: x.name, flurname_list))) + # Suche nach der min. Distanz, dabei werden gewisse Objekte bevorzugt: # Turm, Haupthuegel, Huegel, Pass, Strassenpass, Alpiner Gipfel: 2.5 # Kapelle: 2 @@ -36,26 +47,24 @@ def find_name(coord, dist): flurname_list.sort(key=lambda swiss_name: math.sqrt((abs(swiss_name.x - coord[0]) ** 2 + abs(swiss_name.y - coord[1]) ** 2)) / ( - 2 if swiss_name.object_type in ['Haupthuegel', 'Huegel', 'Pass', 'Strassenpass', 'Alpiner Gipfel', 'Gipfel', - 'Gletscher'] else + 2 if swiss_name.object_type in ['Haltestelle Bahn', 'Huegel', 'Pass', 'Strassenpass', 'Alpiner Gipfel', + 'Gipfel', + ] else 1.25 if swiss_name.object_type in ['Kapelle', 'Turm', 'Schwimmbadareal', 'Campingplatzareal', 'Golfplatzareal', 'Zooareal', 'Freizeitanlagenareal', 'Abwasserreinigungsareal', 'Friedhof', - 'Spitalareal', 'Quartierteil', 'Ort', 'See', 'Bach'] else + 'Spitalareal', 'Quartierteil', 'Ort', 'See', 'Bach', + 'Lokalname swisstopo'] else 1.15 if swiss_name.object_type in ['Haltestelle Bus', 'Haltestelle Schiff', 'Uebrige Bahnen', - 'Haltestelle Bahn'] else + 'Haupthuegel'] else 1.05 if swiss_name.object_type in ['Gebaeude', 'Offenes Gebaeude', 'Schul- und Hochschulareal'] else - 1 if swiss_name.object_type in ['Lokalname swisstopo', 'Flurname swisstopo', 'Tal', 'Grat', - 'Graben'] else 0.95)) - - print(list(map(lambda place: place.name + ' ' + str(place.x) + ' ' + str(place.y), flurname_list))) + 1 if swiss_name.object_type in ['Flurname swisstopo', 'Tal', 'Grat', 'Graben', 'Gletscher'] else 0.95)) if len(flurname_list) == 0: return '' return ('bei/m ' if sqrt( abs(flurname_list[0].x - coord[0]) ** 2 + abs(flurname_list[0].y - coord[1]) ** 2) > dist else '') + \ - flurname_list[ - 0].name + flurname_list[0].name ######################################################################################################################## diff --git a/python_program/find_walk_table_points.py b/python_program/find_walk_table_points.py index 7a0b22c8..610537bc 100644 --- a/python_program/find_walk_table_points.py +++ b/python_program/find_walk_table_points.py @@ -1,48 +1,227 @@ +from copy import copy +from typing import List, Tuple + import geopy.distance import gpxpy +from gpxpy.gpx import GPXTrackPoint -def select_waypoints(raw_gpx_data: gpxpy.gpx): +def select_waypoints(raw_gpx_data: gpxpy.gpx, walk_point_limit=21): """ Algorithm that selects suitable points for the Marschzeittabelle. + raw_gpx_data : raw gpx data from imported GPX-File + walk_point_limit : max number of points in the walk-time table, default 21 + + ------------------------------------------------------------------------- + The aim is to choose points that are as evenly distributed as possible - and that cover the topology of the path as well as possible. + and that cover the topology of the path as well as possible. The process + is done in three steps: preselection, remove_unnecessary_points, reduce_number_of_points. """ - total_distance = 0 - way_points_walk_table = [] + total_distance, pts_step_1 = preselection_step(raw_gpx_data) + pts_step_2 = remove_unnecessary_points(pts_step_1) + pts_step_3 = reduce_number_of_points(pts_step_2, walk_point_limit) + + return total_distance, pts_step_2, pts_step_3 + + +def reduce_number_of_points(pts_step_2: List[Tuple[int, GPXTrackPoint]], walk_point_limit: int): + """ + + Final selection: Iteratively reduce the number of points to the maximum specified in walk_point_limit. To + achieve this we iteratively increase a maximum derivation (drv_limit) value until we have dropped enough points. + + For this purpose, three preselected, adjacent points (A, B, and B) are used to construct a secant line + (a straight line between A and C). Point B is now tested against the drv_limit. If C has a derivation + bigger than the limit, C gets dropped or moved to a new location, i.g. if the distance between point C and a + point on secant line at the same location is bigger than drv_limit, C gets dropped or moved. + + Point B gets moved if and only if after the removal of point B a previously dropped point D is now further way + from the secant line than the drv_limit. In this case be gets dropped and D gets added again, thus be gets + replaced with point D in ths case. In all other cases B gets dropped. + + Once we've checked that. Points A, B, and C are shifted. If B gets dropped, A remains the same. + This process is repeated until we have reached the walk_point_limit. If no point gets dropped during a pass + through the selected points, we increase drv_limit by 2 meters and try again. + + """ + + pts_step_3: List[Tuple[int, GPXTrackPoint]] = copy(pts_step_2) + + pt_A: Tuple[int, GPXTrackPoint] = None + pt_B: Tuple[int, GPXTrackPoint] = None + + drv_limit = 0 + + pt_dropped = True + pts_dropped = [] + + while len(pts_step_3) > walk_point_limit: + + # Increase drv_limit if no points as been dropped in the last pass + if not pt_dropped: + drv_limit += 1 + pt_dropped = False + + for pt_C in pts_step_3: + + # after moving point B we skip one iteration + if pt_B == pt_C: + continue + + if pt_A is not None and pt_B is not None: + + m, b = calc_secant_line(pt_A, pt_C) + secant_elev = calc_secant_elevation(m, b, pt_B) + + if abs(secant_elev - pt_B[1].elevation) < drv_limit: + + # Check if B must be replaced by a previously dropped point D + pt_D = None + for pt in list(filter(lambda p: pt_A[0] < p[0] < pt_C[0], pts_dropped)): + secant_elev = calc_secant_elevation(m, b, pt) + if abs(secant_elev - pt[1].elevation) >= drv_limit: + pt_D = pt + break + + if pt_D is not None: # Replace B with point D + pts_step_3.remove(pt_B) + index = next(x for x, val in enumerate(pts_step_3) if val[0] > pt_D[0]) + pts_step_3.insert(index, pt_D) + pts_dropped.append(pt_B) + + else: # remove point B + pt_dropped = True + + pts_step_3.remove(pt_B) + pts_dropped.append(pt_B) + + pt_B = pt_A + + pt_A = pt_B + pt_B = pt_C + + pt_A = None + pt_B = None + return pts_step_3 + + +def remove_unnecessary_points(pts_step_1: []): + """ + + Now we loop through the preselected points and tighten the selection criteria. + Into the list pts_step_2 we include points according to the following rules: + - first or last point of a track segment + - a points (called point C) with a delta elevation to the previous selected point (called point A, i.g. the last + in point in pts_step_2) of at least 20 meters that are also at least 250 meters apart form the point A (i.g. + the distance between A and C is at least 250 meters). + - a point (called point C) that is at least 1'500 meters apart form the the last point (called C, i.g. the + distance between A and C is at least 1'500 meters). + + After we selected a point C, we loop through all points B out of pts_step_1 that lie between point A and C. + Now we construct a secant line (a straight line between A and C). A point B is now tested against the drv_limit. + If B has a derivation bigger than the limit, B gets also included into the list pts_step_2, i.g. if the + distance between point B and a point on secant line at the same location is bigger than drv_limit, C gets + added. + + """ + + drv_limit = 20 + + last_pt: Tuple[int, GPXTrackPoint] = pts_step_1[0] + pts_step_2: List[Tuple[int, GPXTrackPoint]] = [last_pt] + + for pt in pts_step_1[1:]: + + if last_pt is not None: + last_coord = get_coordinates(last_pt[1]) + coord = get_coordinates(pt[1]) + + if (abs(last_pt[1].elevation - pt[1].elevation) > 20 + and geopy.distance.distance(last_coord, coord).m > 250) \ + or geopy.distance.distance(last_coord, coord).m > 1500: + + m, b = calc_secant_line(last_pt, pt) + + # add point with max derivation + for intermediate_point in list(filter(lambda p: last_pt[0] < p[0] < pt[0], pts_step_1)): + + derivation = abs(calc_secant_elevation(m, b, intermediate_point) - intermediate_point[1].elevation) + if drv_limit <= derivation: + pts_step_2.append(intermediate_point) + drv_limit = derivation + m, b = calc_secant_line(intermediate_point, pt) + last_pt = pt + + pts_step_2.append(pt) + last_pt = pt + + last_point = pts_step_1[len(pts_step_1) - 1] + + if last_pt != last_point: + pts_step_2.append(last_point) + + return pts_step_2 + + +def get_coordinates(last_pt: GPXTrackPoint): + return last_pt.latitude, last_pt.longitude + + +def preselection_step(raw_gpx_data: gpxpy.gpx): + """ + + Preselection: Select all points form the tracking file according which satisfies one of the following criteria. + This guarantees that all important points are considered. We call the set of selected points pts_step_1. + + Preselection criteria: + - first or last point of a track segment + - significant local extremum in track elevation + - significant changes in elevation based on the approximated slope + - distance to last selected point bigger than 250 meters + + """ + + cumulated_distance = 0.0 + pts_step_1: List[Tuple[int, GPXTrackPoint]] = [] oldHeight = 0 + slope = None coord = None - points = [] - + # (1) Preselection for track in raw_gpx_data.tracks: for segment in track.segments: # insert first point - way_points_walk_table.append((0, segment.points[0])) + last_pt = segment.points[0] + pts_step_1.append((0, last_pt)) + last_pt = get_coordinates(last_pt) for point in segment.points: - newCoord = (point.latitude, point.longitude) + newCoord = get_coordinates(point) if coord is not None: distDelta = geopy.distance.distance(coord, newCoord).km - total_distance += distDelta - points.append(point) + cumulated_distance += distDelta if distDelta != 0: newSlope = (oldHeight - point.elevation) / distDelta if slope is not None and newSlope != 0 and ( - (slope / newSlope < -0.75 or abs(slope / newSlope) > 2)): - way_points_walk_table.append((total_distance, point)) + (slope / newSlope < -0.75 or abs(slope / newSlope) > 1.5)): + pts_step_1.append((cumulated_distance, point)) + + if geopy.distance.distance(last_pt, newCoord).m > 250: + pts_step_1.append((cumulated_distance, point)) + last_pt = newCoord slope = newSlope @@ -50,111 +229,73 @@ def select_waypoints(raw_gpx_data: gpxpy.gpx): oldHeight = point.elevation # insert last point - if abs(way_points_walk_table[len(way_points_walk_table) - 1][0] - total_distance) < 0.5: - del way_points_walk_table[len(way_points_walk_table) - 1] - way_points_walk_table.append((total_distance, segment.points[len(segment.points) - 1])) - - # remove points - final_way_points_walk_table = [way_points_walk_table[0]] - old_index = 0 - - for i, point in enumerate(way_points_walk_table): - - old_point = final_way_points_walk_table[len(final_way_points_walk_table) - 1] - - if old_point == point: - continue - - x1, y1 = old_point[0], old_point[1].elevation - x2, y2 = point[0], point[1].elevation - - m = (y1 - y2) / (x1 - x2) - b = (x1 * y2 - x2 * y1) / (x1 - x2) - - diff = 0 - max_diff = 0 - for j in range(old_index, i): - diff += abs(m * way_points_walk_table[j][0] + b - way_points_walk_table[j][1].elevation) - - if max_diff < diff: - max_diff = diff - - diff /= float(i - old_index) + if abs(pts_step_1[len(pts_step_1) - 1][0] - cumulated_distance) < 0.5: + del pts_step_1[len(pts_step_1) - 1] + pts_step_1.append((cumulated_distance, segment.points[len(segment.points) - 1])) - if diff > 15 or max_diff > 100 or point[0] - old_point[0] > 1.5: + return cumulated_distance, pts_step_1 - new_index = int(old_index + 2.0 / 3.0 * (i - old_index)) - new_point = way_points_walk_table[new_index] - if new_point[0] - old_point[0] > 0.25: - final_way_points_walk_table.append(new_point) - - else: - new_index = int(old_index + 1.0 / 2.0 * (i - old_index)) - new_point = way_points_walk_table[new_index] - final_way_points_walk_table.remove(old_point) - final_way_points_walk_table.append(new_point) - - old_index = new_index +def calc_secant_elevation(m: float, b: float, pt_B: Tuple[int, GPXTrackPoint]): + """ - final_way_points_walk_table.append(way_points_walk_table[len(way_points_walk_table) - 1]) + Calculates the elevation of a point on the secant line defined by m and b. The point on the + secant line is chosen such that location matches the location of pkr_B. - # check for points on one line - old_point = None - current_point = None - removedAPoint = True + elevation = m * loc_of_pt_B + b - deviation = 0 + """ - while removedAPoint or len(final_way_points_walk_table) > 30: + # pt_B[0] returns the location of point B in km + return m * (pt_B[0] * 1000) + b - if not removedAPoint: - deviation += 5 - removedAPoint = False +def calc_secant_line(pt_A: Tuple[int, GPXTrackPoint], pt_C: Tuple[int, GPXTrackPoint]): + """ - for next_point in final_way_points_walk_table: + Constructs a secant line through points A and C, i.g. a linear function passing through point A and C. + Returns the slope and the intersect of a linear function through A and C. - if old_point is not None: + """ - x1, y1 = old_point[0], old_point[1].elevation - x2, y2 = next_point[0], next_point[1].elevation + # the locations of pt_A and pt_C given in km. + x1, y1 = pt_A[0] * 1000, pt_A[1].elevation + x2, y2 = pt_C[0] * 1000, pt_C[1].elevation - m = (y1 - y2) / (x1 - x2) - b = (x1 * y2 - x2 * y1) / (x1 - x2) + # if the location of A and C is identical, the slope m is defined as 0 + m: float = (y1 - y2) / (x1 - x2) if (x1 - x2) != 0 else 0.0 - if abs(m * current_point[0] + b - current_point[1].elevation) < deviation: - final_way_points_walk_table.remove(current_point) - removedAPoint = True + # if the location of A and C is identical, the intercept is given by y1 (since y1 == y2) + b: float = (x1 * y2 - x2 * y1) / (x1 - x2) if (x1 - x2) != 0 else y1 - old_point = current_point - current_point = next_point + return m, b - print('Anzahl Punkte: ' + str(len(final_way_points_walk_table))) - print(raw_gpx_data.length_3d()) - print(raw_gpx_data.length_2d()) - print(total_distance) +def prepare_for_plot(gpx: gpxpy.gpx): + """ + Prepares a gpx file for plotting. + Returns two list, one with the elevation of all points in the gpx file and one with the associated, + accumulated distance form the starting point. - return total_distance, way_points_walk_table, final_way_points_walk_table + """ + coord: Tuple[float, float] = None + accumulated_distance: float = 0.0 -def prepare_for_plot(gpx): - coord = None - total_distance = 0 - distances = [] - heights = [] + distances: List[float] = [] + heights: List[float] = [] for track in gpx.tracks: for segment in track.segments: for point in segment.points: - newCoord = (point.latitude, point.longitude) + + newCoord = get_coordinates(point) if coord is not None: distDelta = geopy.distance.distance(coord, newCoord).km - total_distance += distDelta + accumulated_distance += distDelta - distances.append(total_distance) + distances.append(accumulated_distance) heights.append(point.elevation) coord = newCoord diff --git a/python_program/gpx_utils.py b/python_program/gpx_utils.py index 2c310682..b3afd1b1 100644 --- a/python_program/gpx_utils.py +++ b/python_program/gpx_utils.py @@ -1,9 +1,10 @@ +import gpxpy import numpy as np -from python_program.transformation import GPSConverter +from python_program.coord_transformation import GPSConverter -def calc_perimeter(raw_gpx_data): +def calc_perimeter(raw_gpx_data: gpxpy.gpx): min_latitude = None max_latitude = None min_longitude = None diff --git a/python_program/main.py b/python_program/main.py index bbeedc11..53cf0cd4 100644 --- a/python_program/main.py +++ b/python_program/main.py @@ -6,22 +6,30 @@ from python_program.create_map import plot_route_on_map from python_program.walk_table import plot_elevation_profile, create_walk_table -# Open GPX-File with the way-points -gpx_file = open('./testWalks/badiglk1.gpx', 'r') -raw_gpx_data = gpxpy.parse(gpx_file) -# define the departure time of the hike -start_time = datetime(year=2021, month=8, day=16, hour=10, minute=00) +def generate_automated_walk_table(departure_date, gpx_file_path, velocity): + # Open GPX-File with the way-points + gpx_file = open(gpx_file_path, 'r') + raw_gpx_data = gpxpy.parse(gpx_file) -# get Meta-Data -name = raw_gpx_data.tracks[0].name + # get Meta-Data + name = raw_gpx_data.tracks[0].name -# calc Points for walk table -total_distance, temp_points, way_points = select_waypoints(raw_gpx_data) + # calc Points for walk table + total_distance, temp_points, way_points = select_waypoints(raw_gpx_data) + plot_elevation_profile(raw_gpx_data, way_points, temp_points, file_name=name) + create_walk_table(departure_date, velocity, way_points, total_distance, file_name=name) + plot_route_on_map(raw_gpx_data, way_points, file_name=name) -# Walk-Speed in km/h -walk_speed = 3.75 -plot_elevation_profile(raw_gpx_data, way_points, file_name=name + '.png') -create_walk_table(start_time, walk_speed, way_points, total_distance) -plot_route_on_map(raw_gpx_data, way_points) +if __name__ == "__main__": + # defines the departure time of the hike + DEPARTURE_TIME = datetime(year=2021, month=8, day=16, hour=10, minute=00) + + # Path to the GPX-File with the pre-planed route + GPX_FILE_PATH = './testWalks/wanderweg_823.gpx' + + # Velocity , walk-speed in km/h + VELOCITY = 3.75 + + generate_automated_walk_table(DEPARTURE_TIME, GPX_FILE_PATH, VELOCITY) diff --git a/python_program/requirements.txt b/python_program/requirements.txt index 0c3b85e4..9438342d 100644 --- a/python_program/requirements.txt +++ b/python_program/requirements.txt @@ -4,3 +4,4 @@ numpy==1.21.1 geopy==2.2.0 openpyxl==3.0.7 matplotlib==3.4.3 +grequests==0.6.0 diff --git a/python_program/walk_table.py b/python_program/walk_table.py index 817109cb..09e3821a 100644 --- a/python_program/walk_table.py +++ b/python_program/walk_table.py @@ -1,17 +1,31 @@ import math from datetime import timedelta +from typing import Tuple, List +import gpxpy import numpy as np import openpyxl +from gpxpy.gpx import GPXTrackPoint from matplotlib import pyplot as plt -from python_program.calculations import calcTime -from python_program.find_name import find_name +from python_program.find_swisstopo_name import find_name from python_program.find_walk_table_points import prepare_for_plot -from python_program.transformation import GPSConverter +from python_program.coord_transformation import GPSConverter -def plot_elevation_profile(raw_data_points, way_points, file_name): +def plot_elevation_profile(raw_data_points: gpxpy.gpx, + way_points: List[Tuple[int, GPXTrackPoint]], + temp_points: List[Tuple[int, GPXTrackPoint]], + file_name: str): + """ + + Plots the elevation profile of the path contained in the GPX-file. In addition the + plot contains the approximated elevation profile by the way_points. + + Saves the plot as an image in the ./output directory as an image called {{file_name}}<.png + + """ + # plot heights of exported data from SchweizMobil distances, heights = prepare_for_plot(raw_data_points) plt.plot(distances, heights, label='Wanderweg') @@ -21,7 +35,7 @@ def plot_elevation_profile(raw_data_points, way_points, file_name): plt.ylim(ymax=max(heights) + additional_space, ymin=min(heights) - additional_space) # add way_points to plot - # plt.scatter([dist[0] for dist in temp_points], [height[1].elevation for height in temp_points], c='gray', ) + plt.scatter([dist[0] for dist in temp_points], [height[1].elevation for height in temp_points], c='lightgray', ) plt.scatter([dist[0] for dist in way_points], [height[1].elevation for height in way_points], c='orange', ) plt.plot([dist[0] for dist in way_points], [height[1].elevation for height in way_points], label='Marschzeittabelle') @@ -36,17 +50,23 @@ def plot_elevation_profile(raw_data_points, way_points, file_name): plt.grid(color='gray', linestyle='dashed', linewidth=0.5) # show the plot and save image - plt.savefig('imgs/' + file_name, dpi=750) + plt.savefig('output/' + file_name + '_elevation_profile.png', dpi=750) plt.show() -def create_walk_table(time_stamp, speed, way_points, total_distance): +def create_walk_table(time_stamp, speed, way_points, total_distance, file_name: str): + """ + + Saves the Excel file as .output/Marschzeittabelle_{{file_name}}.xlsx' + + """ + xfile = openpyxl.load_workbook('res/Marschzeit_Template.xlsx') sheet = xfile.worksheets[0] oldPoint = None time = 0 - print(' Geschwindikeit: ', speed, 'km/h') + print(' Geschwindigkeit: ', speed, 'km/h') print() print('Distanz Höhe Zeit Uhrzeit Ort (Koordinaten und Namen)') @@ -71,14 +91,15 @@ def create_walk_table(time_stamp, speed, way_points, total_distance): time_stamp = time_stamp + timedelta(hours=deltaTime) # print in§fos + name_of_point = find_name((lv03[0] + 2_000_000, lv03[1] + 1_000_000), 50) print( round(abs((oldPoint[0] if oldPoint is not None else 0.0) - point[0]), 1), 'km ', int(lv03[2]), 'm ü. M. ', round(deltaTime, 1), 'h ', time_stamp.strftime('%H:%M'), 'Uhr ', - (int(lv03[0]), int(lv03[1])), find_name((lv03[0] + 2_000_000, lv03[1] + 1_000_000), 75)) + (int(lv03[0]), int(lv03[1])), name_of_point) - sheet['A' + str(8 + i)] = str(find_name((lv03[0] + 2_000_000, lv03[1] + 1_000_000), 75)) + ' (' + str( + sheet['A' + str(8 + i)] = str(name_of_point) + ' (' + str( int(lv03[0])) + ', ' + str(int(lv03[1])) + ')' sheet['C' + str(8 + i)] = int(lv03[2]) if i > 0: @@ -92,4 +113,21 @@ def create_walk_table(time_stamp, speed, way_points, total_distance): print() print() - xfile.save('Marschzeittabelle_Generiert.xlsx') + xfile.save('output/' + file_name + '_Marschzeittabelle.xlsx') + + +def calcTime(delta_height, delta_dist, speed): + """ + + Calculates the walking time form one point to another + + for this calculation the basic formula form Jugend+Sport is used for preciser we could use the formula + form SchweizMobil or use more way points. But since we want to create a "normal" walk table as specified by + Jugend+Sport we use there basic formula + + """ + + if delta_height is None or delta_dist is None: + return 0 + + return (delta_dist + (delta_height / 100 if delta_height > 0 else 0)) / speed