Merge pull request #2 from cevi/improve-point-selection-algorithm

* Improves the point selection algorithm
* Clarifies code documentation

.gitignore

@@ -1,4 +1,3 @@
-python_program/Marschzeittabelle_Generiert.xlsx
 python_program/.idea
 .idea
 python_program/res/swisstlm3d_2021-04_2056_5728.gdb

Readme.md

@@ -31,15 +31,16 @@ 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
 
 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
 

python_program/.gitignore

@@ -1,5 +1,6 @@
 imgs/
 testWalks/
+output/*
 
 # Byte-compiled / optimized / DLL files
 __pycache__/

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
+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```.

python_program/transformation.py → python_program/coord_transformation.py

@@ -28,18 +28,22 @@
 # Aaron Schmocker [[email protected]]
 # 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):

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]))

python_program/find_name.py → 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
 
 
 ########################################################################################################################