WS30: Improved road generation

ws30/README.TXT

@@ -1,4 +1,4 @@
-A set of trivial python scripts to generate World Scener 3.0 (ws30) scenery
+A set of trivial python scripts to generate World Scenery 3.0 (ws30) scenery
 features from openstreetmap data, using the overpass API to retrieve data at
 runtime.
 

ws30/calc_tile.py

@@ -0,0 +1,206 @@
+# -*- coding: utf-8 -*-
+"""
+shamelessly translated from calc-tile.pl
+"""
+import logging
+from math import floor
+import os
+from typing import List, Tuple
+import unittest
+
+import numpy as np
+
+
+def bucket_span(lat: float) -> float:
+    """Latitude Range -> Tile Width (deg)"""
+    abs_lat = abs(lat)
+    if abs_lat >= 89:
+        return 360
+    elif abs_lat >= 88:
+        return 8
+    elif abs_lat >= 86:
+        return 4
+    elif abs_lat >= 83:
+        return 2
+    elif abs_lat >= 76:
+        return 1
+    elif abs_lat >= 62:
+        return .5
+    elif abs_lat >= 22:
+        return .25
+    elif abs_lat >= 0:
+        return .125
+    return 360
+
+
+def format_lon(lon):
+    """Format longitude as e/w."""
+    if lon < 0.:
+        return "w%03d" % int(0. - lon)
+    else:
+        return "e%03d" % int(lon)
+
+
+def format_lat(lat):
+    """Format latitude as n/s."""
+    if lat < 0.:
+        return "s%02d" % int(0. - lat)
+    else:
+        return "n%02d" % int(lat)
+
+
+def directory_name(lon_lat: Tuple[float, float], separator: str = None) -> str:
+    """Generate the directory name for a location."""
+    (lon, lat) = lon_lat
+    lon_floor = floor(lon)
+    lat_floor = floor(lat)
+    lon_chunk = floor(lon/10.0) * 10
+    lat_chunk = floor(lat/10.0) * 10
+    if separator:
+        return '{}{}{}{}{}'.format(format_lon(lon_chunk), format_lat(lat_chunk), separator,
+                                   format_lon(lon_floor), format_lat(lat_floor))
+    return os.path.join(format_lon(lon_chunk) + format_lat(lat_chunk), format_lon(lon_floor) + format_lat(lat_floor))
+
+
+def calc_tile_index(lon_lat: Tuple[float, float], x: int = 0, y: int = 0) -> int:
+    """See  http://wiki.flightgear.org/Tile_Index_Scheme"""
+    (lon, lat) = lon_lat
+    if x == 0 and y == 0:
+        y = calc_y(lat)
+        x = calc_x(lon, lat)
+
+    index = (int(floor(lon)) + 180) << 14
+    index += (int(floor(lat)) + 90) << 6
+    index += y << 3
+    index += x
+    return index
+
+
+def calc_tile_location(tile_index: int) -> Tuple[Tuple[float, float], Tuple[float, float]]:
+    """The lon/lat as well as x/y tuples for the location of a tile by its index"""
+    lon = tile_index >> 14
+    index = tile_index - (lon << 14)
+    lon = lon - 180
+
+    lat = index >> 6
+    index = index - (lat << 6)
+    lat = lat - 90
+
+    y = index >> 3
+    index = index - (y << 3)
+    x = index
+
+    return (lon, lat), (x, y)
+
+
+def log_tile_info(tile_index: int) -> None:
+    """logs information about the tile to logging.debug."""
+    (lon, lat), (x, y) = calc_tile_location(tile_index)
+    center_lat = lat + y / 8.0 + 0.0625
+    center_lon = bucket_span(center_lat)
+    if center_lon >= 1.0:
+        center_lon = lon + center_lon / 2.0
+    else:
+        center_lon = lon + x * center_lon + center_lon / 2.0
+
+    tile_info = 'Tile location: {}; lon: {}, lat: {}; x: {}, y: {}; tile span degrees: {}; center lon/lat: {}/{}.'
+    logging.debug(tile_info.format(tile_index, lon, lat, x, y, bucket_span(lat), center_lon, center_lat))
+
+
+def construct_path_to_files(base_directory: str, scenery_type: str, center_global: Tuple[float, float]) -> str:
+    """Returns the path to the stg-files in a FG scenery directory hierarchy at a given global lat/lon location.
+    The scenery type is e.g. 'Terrain', 'Object', 'Buildings'."""
+    return os.path.join(base_directory, scenery_type, directory_name(center_global))
+
+
+def construct_stg_file_name(center_global: Tuple[float, float]) -> str:
+    """Returns the file name of the stg-file at a given global lat/lon location"""
+    return construct_stg_file_name_from_tile_index(calc_tile_index(center_global))
+
+
+def construct_stg_file_name_from_tile_index(tile_idx: int) -> str:
+    return str(tile_idx) + '.stg'
+
+
+def construct_btg_file_name_from_tile_index(tile_idx: int) -> str:
+    return str(tile_idx) + '.btg.gz'
+
+
+def construct_btg_file_name_from_airport_code(airport_code: str) -> str:
+    return airport_code + '.btg.gz'
+
+
+def get_north_lat(lat, y):
+    return float(floor(lat)) + y / 8.0 + .125
+
+
+def get_south_lat(lat, y):
+    return float(floor(lat)) + y / 8.0
+
+
+def get_west_lon(lon, lat, x):
+    if x == 0:
+        return float(floor(lon))
+    else: 
+        return float(floor(lon)) + x * (bucket_span(lat))
+
+
+def get_east_lon(lon, lat, x):
+    if x == 0:
+        return float(floor(lon)) + (bucket_span(lat))
+    else: 
+        return float(floor(lon)) + x * (bucket_span(lat)) + (bucket_span(lat))
+
+
+def calc_x(lon: float, lat: float) -> int:
+    """
+    FIXME: is this correct? Also: some returns do not take calculations into account.
+    """
+    epsilon = 0.0000001
+    span = bucket_span(lat)
+    if span < epsilon:
+        lon = 0
+        return 0
+    elif span <= 1.0:
+        return int((lon - floor(lon)) / span)
+    else:
+        if lon >= 0:
+            lon = int(int(lon/span) * span)
+        else:
+            lon = int(int((lon+1)/span) * span - span)
+            if lon < -180:
+                lon = -180
+        return 0
+
+
+def calc_y(lat: float) -> int:
+    return int((lat - floor(lat)) * 8)
+    
+
+def get_stg_files_in_boundary(boundary_west: float, boundary_south: float, boundary_east: float, boundary_north: float,
+                              path_to_scenery: str, scenery_type: str) -> List[str]:
+    """Based on boundary rectangle returns a list of stg-files (incl. full path) to be found within the boundary of
+    the scenery"""
+    stg_files = []
+    for my_lat in np.arange(boundary_south, boundary_north, 0.125):  # latitude; 0.125 as there are always 8 tiles
+        for my_lon in np.arange(boundary_west, boundary_east, bucket_span(my_lat)):  # longitude
+            coords = (my_lon, my_lat)
+            stg_files.append(os.path.join(construct_path_to_files(path_to_scenery, scenery_type, coords),
+                             construct_stg_file_name(coords)))
+    return stg_files
+
+
+# ================ UNITTESTS =======================
+
+
+class TestCalcTiles(unittest.TestCase):
+    def test_calc_tiles(self):
+        self.assertEqual(5760, calc_tile_index((-179.9, 0.1)))
+        self.assertEqual(5752, calc_tile_index((-179.9, -0.1)))
+        self.assertEqual(5887623, calc_tile_index((179.9, 0.1)))
+        self.assertEqual(5887615, calc_tile_index((179.9, -0.1)))
+        self.assertEqual(2954880, calc_tile_index((0.0, 0.0)))
+        self.assertEqual(2938495, calc_tile_index((-0.1, -0.1)))
+
+    def test_file_name(self):
+        self.assertEqual("3088961.stg", construct_stg_file_name((8.29, 47.08)))

ws30/genroads.py

@@ -17,7 +17,7 @@
 # with this program; if not, write to the Free Software Foundation, Inc.,
 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
-
+# Simple generate of line feature - roads, railways, rivers
 # You will need Python Overpass API - "pip install overpy"
 
 import xml.etree.ElementTree as etree
@@ -32,11 +32,13 @@ import calc_tile
 import overpy
 
 nodes = {}
+curr_road_count = 0
+curr_river_count = 0
 road_count = 0
 river_count = 0
 
 if (len(sys.argv) != 6):
-    print("Simple generation of ROAD_LIST files")
+    print("Simple generation of LINEAR_FEATURE_LIST files")
     print("")
     print("Usage: " + sys.argv[0] + " <scenery_dir> <lon1> <lat1> <lon2> <lat2>")
     print("  <scenery_dir>  \tScenery directory to write to")
@@ -63,14 +65,14 @@ def add_to_stg(lat, lon, type):
     with open(stg, 'a') as f:
         f.write("LINE_FEATURE_LIST " + feature_file(lat, lon, type) + " " + type + "\n")
 
-def write_feature(lon, lat, road, type, width):
+def write_feature(lon, lat, road, type, width, lit):
     index = calc_tile.calc_tile_index((lon,lat))
     dirname = os.path.join(scenery_prefix, calc_tile.directory_name((lon, lat)))
     os.makedirs(dirname, exist_ok=True)
     txt = os.path.join(scenery_prefix, calc_tile.directory_name((lon, lat)), feature_file(lat, lon, type))
     #print("Writing " + txt)
     with open(txt, 'a') as f:
-        f.write(str(width) + " 0 1 1 1 1") # Width plus currently unused generic attributes.
+        f.write(str(width) + " " + str(lit) + " 1 1 1 1") # Width plus currently unused generic attributes.
         for pt in road :
             f.write(" "  + str(pt.lon) + " " + str(pt.lat))
         f.write("\n")
@@ -90,34 +92,46 @@ def write_feature(lon, lat, road, type, width):
 
 
 def parse_way(way) :
-    global road_count, river_count, lat1, lon1, lat2, lon2
+    global curr_road_count, curr_river_count, lat1, lon1, lat2, lon2
     pts = []    
     width = 6.0
-    road = 0    
+    lit = 0
+    road = 0
     river = 0
+    rail = 0
 
     highway = way.tags.get("highway")
     waterway = way.tags.get("waterway")
+    railway = way.tags.get("railway")
     feature_type = "None"
 
     if (highway=="motorway_junction") or (highway=="motorway") or (highway=="motorway_link"):
         width = 15.0
-        feature_type = "Road"
+        feature_type = "ws30Freeway"
-        road_count = road_count + 1
+        curr_road_count = curr_road_count + 1
 
-    if (highway=="secondary") or (highway=="primary")  or (highway=="trunk") or (highway=="trunk_link") or (highway=="primary_link") or (highway=="secondary_link") :
+    if (highway=="primary")  or (highway=="trunk") or (highway=="trunk_link") or (highway=="primary_link") :
         width = 12.0
-        feature_type = "Road"
+        feature_type = "ws30Freeway"
-        road_count = road_count + 1
+        curr_road_count = curr_road_count + 1
+
+    if (highway=="secondary") or (highway=="secondary_link") :
+        width = 12.0
+        feature_type = "ws30Road"
+        curr_road_count = curr_road_count + 1
 
     if (highway=="unclassified") or (highway=="tertiary") or (highway=="tertiary_link") or (highway=="service") or (highway=="residential"):
         width = 6.0
-        feature_type = "Road"
+        feature_type = "ws30Road"
-        road_count = road_count + 1
+        curr_road_count = curr_road_count + 1
     if (waterway=="river") or (waterway=="canal") :
         width = 10.0
-        feature_type = "Watercourse"
+        feature_type = "ws30River"
-        river_count = river_count + 1
+        curr_river_count = curr_river_count + 1
+
+    if (railway=="rail") or (railway=="preserved") or (railway=="disused") :
+        width = 4.2  # Standard guage ~ 1.4m, with twice the space either side
+        feature_type = "ws30Railway"
 
     # Use the width if defined and parseable
     if (way.tags.get("width") != None) :
@@ -130,6 +144,16 @@ def parse_way(way) :
         except ValueError :
             print("Unable to parse width " + width_str)
 
+    # Use the lit tag if defined and parseable
+    if (way.tags.get("lit") != None) :
+        lit_str = way.tags.get("lit")
+        if ((lit_str == "no") or (lit == "disused")) :
+            # Specific tags for unlit ways
+            lit = 0
+        else :
+            # Everything else indicates some form of lighting
+            lit = 1
+
     if (feature_type != "None") :       
         # It's a road or river.  Add it to appropriate tile entries.
         tileids = set()
@@ -140,35 +164,52 @@ def parse_way(way) :
             idx = calc_tile.calc_tile_index([lon, lat])
             if ((float(lon1) <= lon <= float(lon2)) and (float(lat1) <= lat <= float(lat2)) and (idx not in tileids)) :
                 # Write the feature to a bucket provided it's within the lat/lon bounds and if we've not already written it there
-                write_feature(lon, lat, way.nodes, feature_type, width)
+                write_feature(lon, lat, way.nodes, feature_type, width, lit)
                 tileids.add(idx)
 
 def writeOSM(result):
         for child in result.ways:
             parse_way(child)
 
-# Get River data
+for lat in range(int(lat1), int(lat2)):
+    for lon in range(int(lon1), int(lon2)):
+        curr_road_count = 0
+        curr_river_count = 0
+        osm_bbox = ",".join([str(lat), str(lon), str(lat+1), str(lon+1)])
+        #api = overpy.Overpass(url="https://lz4.overpass-api.de/api/interpreter")
+        api = overpy.Overpass(url="https://overpass.kumi.systems/api/interpreter")
 
-osm_bbox = ",".join([lat1, lon1, lat2, lon2])
+        # Get River data
-api = overpy.Overpass(url="https://lz4.overpass-api.de/api/interpreter")
+        river_query = "(way[\"waterway\"=\"river\"](" + osm_bbox + "); way[\"waterway\"=\"canal\"](" + osm_bbox + "););(._;>;);out;"
+        result = api.query(river_query)
+        writeOSM(result)
 
-river_query = "(way[\"waterway\"=\"river\"](" + osm_bbox + "); way[\"waterway\"=\"canal\"](" + osm_bbox + "););(._;>;);out;"
+        railway_query = "("
-result = api.query(river_query)
+        railway_types = ["rail", "preserved", "disused"]
-writeOSM(result)
+        for r in railway_types :
+            railway_query = railway_query + "way[\"railway\"=\"" + r + "\"](" + osm_bbox + ");"
 
-road_query = "("
+        railway_query = railway_query + ");(._;>;);out;"
-#road_types = ["unclassified", "tertiary", "service", "secondary", "primary", "motorway_junction", "motorway"]
+        result = api.query(railway_query)
-#road_types = ["tertiary", "secondary", "primary", "motorway_junction", "motorway"]
+        writeOSM(result)
-road_types = ["motorway", "trunk", "primary", "secondary", "tertiary", "unclassified", "residential", "motorway_link", "trunk_link", "primary_link", "secondary_link", "tertiary_link"]
-for r in road_types :
-    road_query = road_query + "way[\"highway\"=\"" + r + "\"](" + osm_bbox + ");"
 
-road_query = road_query + ");(._;>;);out;"
+        road_query = "("
-result = api.query(road_query)
+        #road_types = ["unclassified", "tertiary", "service", "secondary", "primary", "motorway_junction", "motorway"]
-writeOSM(result)
+        #road_types = ["tertiary", "secondary", "primary", "motorway_junction", "motorway"]
+        road_types = ["motorway", "trunk", "primary", "secondary", "tertiary", "unclassified", "residential", "motorway_link", "trunk_link", "primary_link", "secondary_link", "tertiary_link"]
+        for r in road_types :
+            road_query = road_query + "way[\"highway\"=\"" + r + "\"](" + osm_bbox + ");"
 
-print("Wrote total of " + str(road_count) + " roads")
+        road_query = road_query + ");(._;>;);out;"
-print("Wrote total of " + str(river_count) + " rivers")
+        result = api.query(road_query)
+        writeOSM(result)
+
+        print(str(lat) + "," + str(lon) + ": " + str(curr_road_count) + " roads " + str(curr_river_count) + " rivers")
+        road_count += curr_road_count
+        river_count += curr_river_count
+
+print(str(lat1) + "," + str(lon1) + " " + str(lat2) + "," + str(lon2))
+print("Wrote total of " + str(road_count) + " roads " + str(river_count) + " rivers")
 
 
 #

ws30/genwater.py

@@ -106,6 +106,12 @@ def parse_way(way) :
             write_feature(lon, lat, way.nodes, feature_type)
             tileids.add(idx)
 
+def get_first_node(way) :
+    return way.nodes[0]
+
+def get_last_node(way) :
+    return way.nodes[-1]
+
 def parse_multi(relation):
     for way in relation.members:
         if (way.tags.get("role") == "outer") :