Improved checks and visualization

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,15 +356,22 @@ 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)
 parser.add_argument("-p", "--point", 
                     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)

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)