From d7abbaa10f2c4f37d09cb37b52e73de984728d42 Mon Sep 17 00:00:00 2001 From: Stuart Buchanan Date: Sun, 5 Dec 2021 20:38:19 +0000 Subject: [PATCH] WS30: Improved road generation --- ws30/README.TXT | 2 +- ws30/calc_tile.py | 206 ++++++++++++++++++++++++++++++++++++++++++++++ ws30/genroads.py | 107 ++++++++++++++++-------- ws30/genwater.py | 6 ++ 4 files changed, 287 insertions(+), 34 deletions(-) create mode 100755 ws30/calc_tile.py diff --git a/ws30/README.TXT b/ws30/README.TXT index 2b979ed..ab6bac2 100644 --- a/ws30/README.TXT +++ b/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. diff --git a/ws30/calc_tile.py b/ws30/calc_tile.py new file mode 100755 index 0000000..88ced67 --- /dev/null +++ b/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))) diff --git a/ws30/genroads.py b/ws30/genroads.py index a02dbd7..8038e35 100755 --- a/ws30/genroads.py +++ b/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] + " ") print(" \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" - road_count = road_count + 1 + feature_type = "ws30Freeway" + 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" - road_count = road_count + 1 + feature_type = "ws30Freeway" + 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" - road_count = road_count + 1 + feature_type = "ws30Road" + curr_road_count = curr_road_count + 1 if (waterway=="river") or (waterway=="canal") : width = 10.0 - feature_type = "Watercourse" - river_count = river_count + 1 + feature_type = "ws30River" + 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]) -api = overpy.Overpass(url="https://lz4.overpass-api.de/api/interpreter") + # Get River data + 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;" -result = api.query(river_query) -writeOSM(result) + railway_query = "(" + railway_types = ["rail", "preserved", "disused"] + for r in railway_types : + railway_query = railway_query + "way[\"railway\"=\"" + r + "\"](" + osm_bbox + ");" -road_query = "(" -#road_types = ["unclassified", "tertiary", "service", "secondary", "primary", "motorway_junction", "motorway"] -#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 + ");" + railway_query = railway_query + ");(._;>;);out;" + result = api.query(railway_query) + writeOSM(result) -road_query = road_query + ");(._;>;);out;" -result = api.query(road_query) -writeOSM(result) + road_query = "(" + #road_types = ["unclassified", "tertiary", "service", "secondary", "primary", "motorway_junction", "motorway"] + #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") -print("Wrote total of " + str(river_count) + " rivers") + road_query = road_query + ");(._;>;);out;" + 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") # diff --git a/ws30/genwater.py b/ws30/genwater.py index e1ca100..2451490 100755 --- a/ws30/genwater.py +++ b/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") :