From 1f0d0c1c8272baf21bd77c1260fc4576615ddf38 Mon Sep 17 00:00:00 2001 From: "Dr. Tilmann Bubeck" Date: Thu, 12 Oct 2023 17:13:04 +0200 Subject: [PATCH] Improved checks and visualization --- bin/check-consistency.py | 104 ++++++++++++++++++++++++++-- bin/common.py | 142 +++++++++++++++++++++++++++++++++++++-- bin/visualize.py | 90 +++---------------------- 3 files changed, 245 insertions(+), 91 deletions(-) diff --git a/bin/check-consistency.py b/bin/check-consistency.py index 9531918..1920c6b 100644 --- a/bin/check-consistency.py +++ b/bin/check-consistency.py @@ -9,6 +9,8 @@ from datetime import datetime import os import common +from shapely.geometry import Polygon +import re def findLatLon(p): global records @@ -61,6 +63,9 @@ def getAirspaceName(record): h1 = "" return (n1, h1) +def getAirspaceName2(record): + (n,h) = getAirspaceName(record) + return f'{n} {h}' def findingForTwoPoints(message, record1, record2, p1, p2): global records @@ -150,6 +155,22 @@ def checkPoints(): findNearPoints(record, element["start"]) findNearPoints(record, element["end"]) +def checkDB(): + global records + + for record in records: + for element in record["elements"]: + if element["type"] == "arc": + if not "radius" in element: + (dist1_cm, bearing) = common.geo_distance(element["center"][0],element["center"][1], element["start"][0], element["start"][1]) + (dist2_cm, bearing) = common.geo_distance(element["center"][0],element["center"][1], element["end"][0], element["end"][1]) + diff_m = abs(dist1_cm-dist2_cm)/100 + if diff_m > 30: + lineno = None + if "lineno" in element: + lineno = element["lineno"] + problem(1, f'DB has a big difference radius between start and end of {diff_m:.0f}m', lineno) + def getFirstPoint(element): if element["type"] == "point": return element["location"] @@ -201,13 +222,20 @@ def checkOpenAirspaces(): prio_name = [ "fatal", "error", "warning" ] def problem(prio, message, lineno = None): - global problem_count, prio_name + global problem_count, prio_name, args if lineno != None: in_line = f', line {lineno}' else: in_line = "" - print(f'{prio_name[prio]}{in_line}: {message}') + + print_out = True + if args.errors_only: + if prio >= 2: + print_out = False + if print_out: + print(f'{prio_name[prio]}{in_line}: {message}') + problem_count[prio] = problem_count[prio] + 1 def printProblemCounts(): @@ -247,7 +275,60 @@ def checkEncoding(fp): message = message + "^" message = message + "\n" problem(2, message, lineno) - + +def checkInvalidPolygons(records): + for record in records: + if not record["polygon"].is_valid: + problem(1, "Invalid Polygon for " + getAirspaceName2(record)) + +def checkHeightFL(height): + h = height[2:].strip() + if h == "": + return 1 + else: + if h.isnumeric(): + fl = int(h) + if int(fl / 5) * 5 == fl: + return 2 + return 0 + +def checkHeightFT(height): + m = re.match("(\d+)([a-z]+)\s*(\w+)", height) + if m: + if not m.group(1).isdecimal(): + return 1 + h = int(m.group(1)) + if int(h/100)*100 != h: + return 2 + if not m.group(2) in ["ft"]: + return 1 + if m.group(3) == "MSL": + return 2 + if not m.group(3) in ["AGL", "AMSL"]: + return 1 + return 0 + return 1 + +def checkHeight(height): + if height == "GND": + return 0 + if height.startswith("FL"): + return checkHeightFL(height) + if height[0].isdecimal(): + return checkHeightFT(height) + + return 1 + +def checkHeights(records): + for record in records: + for h in ["floor", "ceiling"]: + if h in record: + prio = checkHeight(record[h]) + if prio > 0: + problem(prio, f'Incorrect height "{record[h]}" in {getAirspaceName2(record)}') + else: + problem(1, f'Missing height "{h}" in {getAirspaceName2(record)}') + def fixOpenAirspaces(fp): (root, ext) = os.path.splitext(args.filename) @@ -275,6 +356,8 @@ def fixOpenAirspaces(fp): findings = [] parser = argparse.ArgumentParser(description='Check OpenAir airspace file for consistency') +parser.add_argument("-e", "--errors-only", action="store_true", + help="Print only errors and no warnings") parser.add_argument("-d", "--distance", help="Specify the distance in meters to see two points as close", type=int, default=100) @@ -282,8 +365,13 @@ def fixOpenAirspaces(fp): help="Find all other points near this point.") parser.add_argument("-F", "--fix-closing", action="store_true", help="Fix all open airspaces by inserting a closing point") +parser.add_argument("-n", "--no-arc", action="store_true", + help="Resolve arcs as straight line") +parser.add_argument("-f", "--fast-arc", action="store_true", + help="Resolve arcs with less quality (10 degree steps)") parser.add_argument("filename") args = parser.parse_args() +common.setArgs(args) # read file into a StringIO, as we have to parse it multiple times content = io.StringIO() @@ -312,13 +400,19 @@ def fixOpenAirspaces(fp): if args.fix_closing: fixOpenAirspaces(content) sys.exit(0) - + +common.resolveRecordArcs(records) +common.createPolygons(records) + +checkInvalidPolygons(records) +checkHeights(records) +checkDB() checkEncoding(content) checkOpenAirspaces() checkNameEncoding(records) checkCircles() checkPoints() - + ret = printProblemCounts() sys.exit(ret) diff --git a/bin/common.py b/bin/common.py index 706539d..ed725eb 100644 --- a/bin/common.py +++ b/bin/common.py @@ -1,7 +1,17 @@ # -*- mode: python; python-indent-offset: 4 -*- import math +from pprint import pprint +from shapely.geometry import Polygon +import sys +args = None + +def setArgs(globalArgs): + global args + + args = globalArgs + def nautical_miles_to_km(nm): return nm * 1.852 @@ -81,14 +91,14 @@ def get_kx_ky(lat): ky = (111.13209 - 0.56605 * cos2 + 0.0012 * cos4) return (kx, ky) +def decimal_degrees_to_dms(decimal_degrees): + mnt,sec = divmod(decimal_degrees*3600,60) + deg,mnt = divmod(mnt, 60) + return (round(deg),round(mnt),round(sec)) + def strDegree(v, width): - v_a = abs(v) - degree_i = int(v_a) - minute = (v_a - degree_i) * 60 - minute_i = int(minute) - sec = (minute - minute_i) * 60 - sec_i = int(round(sec)) - return f'{degree_i:0{width}d}:{minute_i:02d}:{sec_i:02d}' + (degrees,minutes,seconds) = decimal_degrees_to_dms(abs(v)) + return f'{degrees:0{width}d}:{minutes:02d}:{seconds:02d}' def strLatLon(P): result = strDegree(P[0], 2) + " " @@ -104,3 +114,121 @@ def strLatLon(P): result = result + "W " return result +def resolveRecordArcs(records): + for record in records: + resolveArcs(record) + +def resolveArcs(record): + global args + + elements_resolved = [] + for element in record["elements"]: + if element["type"] == "point": + element["computed"] = False + elements_resolved.append(element) + #print(f'resolve(point, {strLatLon(element["location"])}') + elif element["type"] == "arc": + if not args.no_arc: + if "radius" in element: + elements_resolved.extend(resolve_DA(element["center"], nautical_miles_to_km(element["radius"]), element["start"], element["end"], element["clockwise"], True)) + else: + elements_resolved.extend(resolve_DB(element["center"], element["start"], element["end"], element["clockwise"])) + else: + elements_resolved.extend(createElementPoint(element["start"][0], element["start"][1])) + elements_resolved.extend(createElementPoint(element["end"][0], element["end"][1])) + elif element["type"] == "circle": + elements_resolved.extend(resolve_circle(element)) + else: + print(f'Unknown element type: {element["type"]}') + sys.exit(1) + + record["elements_resolved"] = elements_resolved + return elements_resolved + +def createElementPoint(lat, lon): + element = {} + element["type"] = "point" + element["location"] = [ lat, lon] + return element + +def resolve_DA(center, radius_km, start_angle, end_angle, clockwise, use_edge): + global args + elements = [] + + if args.fast_arc: + dir = 10 + else: + dir = 1 + + if clockwise: + while start_angle > end_angle: + end_angle = end_angle + 360 + else: + dir = -dir + while start_angle < end_angle: + start_angle = start_angle + 360 + + #print("start_angle:", start_angle) + #print("end_angle:", end_angle) + #print("clockwise:", clockwise) + + if not use_edge: + start_angle = start_angle + dir + end_angle = end_angle - dir + + angle = start_angle + + while True: + # print("loop angle:", angle) + if clockwise: + if angle >= end_angle: + angle = end_angle + break + else: + if angle <= end_angle: + angle = end_angle + break + (lat, lon) = geo_destination(center[0], center[1], angle, radius_km) + # print(f'lat={lat} lon={lon}') + element = createElementPoint(lat,lon) + element["computed"] = True + elements.append(element) + angle = angle + dir + + return elements + +def resolve_circle(element): + center = element["center"] + radius_km = nautical_miles_to_km(element["radius"]) + elements = resolve_DA(center, radius_km, 0, 180, True, True) + elements.extend(resolve_DA(center, radius_km, 180, 0, True, True)) + return elements + +def resolve_DB(center, start, end, clockwise): + + (dist_s, bearing_s) = geo_distance(center[0], center[1], start[0], start[1]) + (dist_e, bearing_e) = geo_distance(center[0], center[1], end[0], end[1]) + + dist_s_km = (dist_s / 100) / 1000 + + elements = resolve_DA(center, dist_s_km, bearing_s, bearing_e, clockwise, False) + element = createElementPoint(start[0], start[1]) + element["computed"] = False + elements.insert(0, element) + element = createElementPoint(end[0], end[1]) + element["computed"] = False + elements.append(element) + + return elements + + +def createPolygons(records): + for record in records: + createPolygonOfRecord(record) + +def createPolygonOfRecord(record): + points = [] + + for element in record["elements_resolved"]: + points.append(element["location"]) + record["polygon"] = Polygon(points) diff --git a/bin/visualize.py b/bin/visualize.py index bfb4b7a..b8d8b64 100644 --- a/bin/visualize.py +++ b/bin/visualize.py @@ -32,73 +32,14 @@ def plot_to(plt, pos, color="black", label=True): global last_pos global args - #print("plot_to", pos) + #print(f'plot_to({common.strLatLon(pos)}') if last_pos != None: plot_line(plt, last_pos[1], last_pos[0], pos[1], pos[0], color) - + #print(f'plot_to({common.strLatLon(last_pos)} -> {common.strLatLon(pos)}') last_pos = pos if label and args.show_coords: plt.annotate(common.strLatLon(pos), (pos[1], pos[0]), color=color) -def plot_circle(plt, element, color="black"): - center = element["center"] - radius_km = common.nautical_miles_to_km(element["radius"]) - plot_arc_DA(plt, center, radius_km, 0, 180, True, True, color, False) - plot_arc_DA(plt, center, radius_km, 180, 0, True, True, color, False) - -def plot_arc_DA(plt, center, radius_km, start_angle, end_angle, clockwise, use_edge, color="black", label=True): - global args - - if args.fast: - dir = 10 - else: - dir = 1 - - if clockwise: - while start_angle > end_angle: - end_angle = end_angle + 360 - else: - dir = -dir - while start_angle < end_angle: - start_angle = start_angle + 360 - - #print("start_angle:", start_angle) - #print("end_angle:", end_angle) - #print("clockwise:", clockwise) - - if not use_edge: - start_angle = start_angle + dir - end_angle = end_angle - dir - - angle = start_angle - - while True: - # print("loop angle:", angle) - if clockwise: - if angle >= end_angle: - angle = end_angle - break - else: - if angle <= end_angle: - angle = end_angle - break - (lat, lon) = common.geo_destination(center[0], center[1], angle, radius_km) - # print(f'lat={lat} lon={lon}') - plot_to(plt, [lat, lon], color, False) - angle = angle + dir - - -def plot_arc(plt, center, start, end, clockwise, color="black", label=True): - - (dist_s, bearing_s) = common.geo_distance(center[0], center[1], start[0], start[1]) - (dist_e, bearing_e) = common.geo_distance(center[0], center[1], end[0], end[1]) - - dist_s_km = (dist_s / 100) / 1000 - - plot_to(plt, start, color, label) - plot_arc_DA(plt, center, dist_s_km, bearing_s, bearing_e, clockwise, False, color, label) - plot_to(plt, end, color, label) - ax = None fig = None zoom = 1 @@ -118,28 +59,16 @@ def plot(records): color_pos = 25 plot_reset() - # plot_arc(plt, [0,0], [1,0], [0,1], True) if True: for record in records: plot_reset() color = cmap(color_pos) color_pos = (color_pos + 7) % color_num - for element in record["elements"]: + for element in record["elements_resolved"]: if element["type"] == "point": - plot_to(plt, element["location"], color) - elif element["type"] == "arc": - if not args.no_arc: - if "radius" in element: - #pprint(element) - plot_arc_DA(plt, element["center"], common.nautical_miles_to_km(element["radius"]), element["start"], element["end"], element["clockwise"], True, color) - else: - plot_arc(plt, element["center"], element["start"], element["end"], element["clockwise"], color) - else: - plot_to(plt, element["start"], color) - plot_to(plt, element["end"], color) - elif element["type"] == "circle": - plot_circle(plt, element, color) + label = not ("computed" in element and element["computed"] == True) + plot_to(plt, element["location"], color, label) else: print(f'Unknown element type: {element["type"]}') sys.exit(1) @@ -182,13 +111,14 @@ def on_press(event): parser = argparse.ArgumentParser(description='Plot OpenAir airspace file') parser.add_argument("-n", "--no-arc", action="store_true", - help="Draw arcs as straight line") -parser.add_argument("-f", "--fast", action="store_true", - help="Draw arcs with less quality (10 degree steps)") + help="Resolve arcs as straight line") +parser.add_argument("-f", "--fast-arc", action="store_true", + help="Resolve arcs with less quality (10 degree steps)") parser.add_argument("-c", "--show-coords", action="store_true", help="Show latitude/longitude of points in plot") parser.add_argument("filename") args = parser.parse_args() +common.setArgs(args) # read file into a StringIO, as we have to parse it multiple times content = io.StringIO() @@ -204,5 +134,7 @@ def on_press(event): raise error records.append(record) +common.resolveRecordArcs(records) + plot(records)