Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
190 changes: 99 additions & 91 deletions packages/minisky/minisky/tools/areafilter.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,12 +6,19 @@
the BOX, CIRCLE, POLY, POLYALT, LINE, and POLYLINE stack commands, and is
used by plugins and traffic logic that need to know which aircraft are
inside an area. Each `AreaFilter` stores its defined shapes by name.

Boxes and circles are tested directly in geographic coordinates; polygons
use a planar shapely geometry internally, which is only valid for polygons
that do not cross the antimeridian or enclose a pole - those are rejected
at definition time.
"""

from __future__ import annotations

from typing import Protocol

import numpy as np
from matplotlib.path import Path
import shapely

from minisky.command import AltM, Keyword, LatLonDeg, LatLonDegrees, command
from minisky.result import Err, Ok, Result
Expand All @@ -23,7 +30,7 @@ class AreaFilter:

def __init__(self) -> None:
# Dictionary of all basic shapes (The shape classes defined in this file) by name
self.basic_shapes: dict[str, Shape] = {}
self.basic_shapes: dict[str, HasArea | Line] = {}

def has_area(self, areaname: str) -> bool:
"""Check if area with name 'areaname' exists."""
Expand Down Expand Up @@ -59,16 +66,19 @@ def define_area(
else:
return Err(f"Unknown shape: {areaname}")

if areatype == "BOX":
shape = Box(areaname, coordinates, top, bottom)
elif areatype == "CIRCLE":
shape = Circle(areaname, coordinates, top, bottom)
elif areatype[:4] == "POLY":
shape = Poly(areaname, coordinates, top, bottom)
elif areatype == "LINE":
shape = Line(areaname, coordinates)
else:
return Err(f"Unknown shape type: {areatype}")
try:
if areatype == "BOX":
shape = Box(areaname, coordinates, top, bottom)
elif areatype == "CIRCLE":
shape = Circle(areaname, coordinates, top, bottom)
elif areatype[:4] == "POLY":
shape = Poly(areaname, coordinates, top, bottom)
elif areatype == "LINE":
shape = Line(areaname, coordinates)
else:
return Err(f"Unknown shape type: {areatype}")
except ValueError as e:
return Err(str(e))

self.basic_shapes[areaname] = shape
return Ok(f"Created {areatype} {areaname}")
Expand Down Expand Up @@ -134,7 +144,7 @@ def define_polyline_area(
"""Draw a multi-segment line through position vertices."""
return self.define_area(name, "LINE", self._coordinates((first, *additional)))

def checkInside(
def contains(

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I haven't taken a look at what parts of the code uses AreaFilter but it looks like we are still associating a "line" (which has no area) as part of an "area filter". I would get rid of Line from AreaFilter.shapes and let callers handle it. Maybe that is too big of a refactor so I would just merge it

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll give it a check later today but these shapes are mostly helpers for the user. I don't think they are that important for calculation. I can rename this to shapes or something

I recall that in bluesky there is a plugin which uses areafilter to create an experiment area. So basically anything outside tje shape is not logged and quickly deleted.

self, areaname: str, lat: np.ndarray, lon: np.ndarray, alt: np.ndarray
) -> np.ndarray:
"""Check if points with coordinates lat, lon, alt are inside area with name 'areaname'.
Expand All @@ -147,12 +157,12 @@ def checkInside(

Returns:
Array of booleans, True == Inside. All False when no area with
the given name exists.
the given name exists, or when the named shape is a line.
"""
if areaname not in self.basic_shapes:
area = self.basic_shapes.get(areaname)
if area is None or isinstance(area, Line):
return np.zeros(len(lat), dtype=bool)
area = self.basic_shapes[areaname]
return area.checkInside(lat, lon, alt)
return area.contains(lat, lon, alt)

def reset(self) -> None:
"""Clear all data."""
Expand All @@ -169,70 +179,40 @@ def deleteArea(self, name: str) -> Result[str, str]:
return Err(f"No area found with name {name}.")


class Shape:
"""
Base class of BlueSky shapes

Handles the naming and altitude bounds common to all shape types.
Derived classes implement checkInside() for their specific geometry.

Attributes:
name: Area name.
coordinates: Flat list of lat/lon coordinates in deg defining the
shape (plus radius in nm for circles).
top: Upper altitude bound [m].
bottom: Lower altitude bound [m].
raw: Dictionary with the raw shape definition (name, kind,
coordinates).
"""

def __init__(self, name: str, coordinates, top: float = 1e9, bottom: float = -1e9) -> None:
self.raw = {"name": name, "shape": self.kind(), "coordinates": coordinates}
self.name = name
self.coordinates = coordinates
self.top = np.maximum(bottom, top)
self.bottom = np.minimum(bottom, top)

def checkInside(self, lat: np.ndarray, lon: np.ndarray, alt: np.ndarray) -> np.ndarray:
"""Returns True (or boolean array) if coordinate lat, lon, alt lies
within this shape.
def _vrange_str(top: float, bottom: float) -> str:
"""Describe an altitude range [m] for shape __str__ output."""
if top < 9e8:
if bottom > -9e8:
return f" with altitude between {bottom} and {top}"
else:
return f" with altitude below {top}"
if bottom > -9e8:
return f" with altitude above {bottom}"
return ""

Reimplement this function in the derived shape classes for this to
work.
"""
return np.zeros(len(lat), dtype=bool)

def _str_vrange(self) -> str:
if self.top < 9e8:
if self.bottom > -9e8:
return f" with altitude between {self.bottom} and {self.top}"
else:
return f" with altitude below {self.top}"
if self.bottom > -9e8:
return f" with altitude above {self.bottom}"
return ""
class HasArea(Protocol):
"""An area shape that supports point-inside tests.

def __str__(self) -> str:
return (
f"{self.name} is a {self.raw['shape']} with coordinates "
+ ", ".join(str(c) for c in self.coordinates)
+ self._str_vrange()
)
Lines are deliberately not part of this protocol: a line has zero area,
so asking whether an aircraft is inside one is a category error.
"""

@classmethod
def kind(cls) -> str:
"""Return a string describing what kind of shape this is."""
return cls.__name__.upper()
def contains(self, lat: np.ndarray, lon: np.ndarray, alt: np.ndarray) -> np.ndarray:
"""Return whether points (lat [deg], lon [deg], alt [m]) lie inside
this shape's geometry and altitude bounds."""
...


class Line(Shape):
class Line:
"""A line shape between two lat/lon positions [deg].

Purely graphical: the inherited checkInside() always returns False.
Purely graphical: a line has no inside, and no contains().
"""

def __init__(self, name: str, coordinates) -> None:
super().__init__(name, coordinates)
self.name = name
self.coordinates = coordinates

def __str__(self) -> str:
return (
Expand All @@ -242,76 +222,104 @@ def __str__(self) -> str:
)


class Box(Shape):
class Box(HasArea):
"""A lat/lon-aligned box shape.

Defined by two opposite corner points [deg] (sorted at construction)
and optional altitude bounds [m].
"""

def __init__(self, name: str, coordinates, top: float = 1e9, bottom: float = -1e9) -> None:
super().__init__(name, coordinates, top, bottom)
self.name = name
self.coordinates = coordinates
self.top = np.maximum(bottom, top)
self.bottom = np.minimum(bottom, top)
# Sort the order of the corner points
self.lat0 = min(coordinates[0], coordinates[2])
self.lon0 = min(coordinates[1], coordinates[3])
self.lat1 = max(coordinates[0], coordinates[2])
self.lon1 = max(coordinates[1], coordinates[3])

def checkInside(self, lat: np.ndarray, lon: np.ndarray, alt: np.ndarray):
def contains(self, lat: np.ndarray, lon: np.ndarray, alt: np.ndarray) -> np.ndarray:
"""Return whether points (lat [deg], lon [deg], alt [m]) lie inside this box."""
return (
((self.lat0 <= lat) & (lat <= self.lat1))
& ((self.lon0 <= lon) & (lon <= self.lon1))
& ((self.bottom <= alt) & (alt <= self.top))
)

def __str__(self) -> str:
return (
f"{self.name} is a BOX with coordinates "
+ ", ".join(str(c) for c in self.coordinates)
+ _vrange_str(self.top, self.bottom)
)


class Circle(Shape):
class Circle(HasArea):
"""A circle shape.

Defined by a center position [deg], a radius [nm], and optional
altitude bounds [m].
"""

def __init__(self, name: str, coordinates, top: float = 1e9, bottom: float = -1e9) -> None:
super().__init__(name, coordinates, top, bottom)
self.name = name
self.coordinates = coordinates
self.top = np.maximum(bottom, top)
self.bottom = np.minimum(bottom, top)
self.clat = coordinates[0]
self.clon = coordinates[1]
self.r = coordinates[2]

def checkInside(self, lat: np.ndarray, lon: np.ndarray, alt: np.ndarray):
def contains(self, lat: np.ndarray, lon: np.ndarray, alt: np.ndarray) -> np.ndarray:
"""Return whether points (lat [deg], lon [deg], alt [m]) lie within
the circle radius [nm] and altitude bounds."""
distance = kwikdist(self.clat, self.clon, lat, lon) # [NM]
inside = (distance <= self.r) & (self.bottom <= alt) & (alt <= self.top)
return inside
return (distance <= self.r) & (self.bottom <= alt) & (alt <= self.top)

def __str__(self) -> str:
return (
f"{self.name} is a CIRCLE with "
f"center ({self.clat}, {self.clon}) "
f"and radius {self.r}." + self._str_vrange()
f"and radius {self.r}." + _vrange_str(self.top, self.bottom)
)


class Poly(Shape):
class Poly(HasArea):
"""A polygon shape.

Defined by a sequence of lat/lon vertices [deg] and optional altitude
bounds [m]; the border is stored as a matplotlib Path for fast
point-in-polygon tests.
bounds [m]. Containment is tested against a planar shapely polygon in
(lat, lon) space, which cannot represent polygons that cross the
antimeridian or enclose a pole; those are rejected with ValueError.
"""

def __init__(self, name: str, coordinates, top: float = 1e9, bottom: float = -1e9) -> None:
super().__init__(name, coordinates, top, bottom)
self.border = Path(np.reshape(coordinates, (len(coordinates) // 2, 2)))

def checkInside(self, lat: np.ndarray, lon: np.ndarray, alt: np.ndarray):
self.name = name
self.coordinates = coordinates
self.top = np.maximum(bottom, top)
self.bottom = np.minimum(bottom, top)
vertices = np.reshape(coordinates, (-1, 2))
# A ring with an edge spanning more than 180 deg of longitude either
# crosses the antimeridian or winds around a pole; both are invalid
# in the planar (lat, lon) space the containment test runs in.
lons = np.append(vertices[:, 1], vertices[0, 1])
if np.any(np.abs(np.diff(lons)) > 180.0):
raise ValueError(
f"Polygon {name} crosses the antimeridian or encloses a pole; "
"split it into separate polygons."
)
self._geom = shapely.Polygon(vertices)

def contains(self, lat: np.ndarray, lon: np.ndarray, alt: np.ndarray) -> np.ndarray:
"""Return whether points (lat [deg], lon [deg], alt [m]) lie inside
the polygon border and altitude bounds."""
points = np.vstack((lat, lon)).T
inside = np.all(
(self.border.contains_points(points), self.bottom <= alt, alt <= self.top),
axis=0,
return shapely.contains_xy(self._geom, lat, lon) & (self.bottom <= alt) & (alt <= self.top)

def __str__(self) -> str:
return (
f"{self.name} is a POLY with coordinates "
+ ", ".join(str(c) for c in self.coordinates)
+ _vrange_str(self.top, self.bottom)
)
return inside
2 changes: 1 addition & 1 deletion packages/minisky/minisky/traffic/trafficgroups.py
Original file line number Diff line number Diff line change
Expand Up @@ -133,7 +133,7 @@ def add_to_group(
self.traffic.idx(sole_member) >= 0 or sole_member in self
)
if sole_member is not None and not sole_is_selection and self.areas.has_area(sole_member):
inside = self.areas.checkInside(
inside = self.areas.contains(
sole_member, self.traffic.lat, self.traffic.lon, self.traffic.alt
)
indices = np.flatnonzero(inside)
Expand Down
2 changes: 1 addition & 1 deletion packages/minisky/pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,6 @@ dependencies = [
"typer>=0.15.0",
"colorama>=0.4.6",
"fastapi[standard]>=0.115.7",
"matplotlib>=3.10.0",
"numpy>=2.2.2",
"openap>=2.4",
"pandas>=2.2.3",
Expand All @@ -17,6 +16,7 @@ dependencies = [
"pyarrow>=19.0.1",
"requests>=2.32.3",
"scipy>=1.15.1",
"shapely>=2.0",
]

[project.scripts]
Expand Down
30 changes: 27 additions & 3 deletions packages/minisky/tests/unit/test_areafilter.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ def area_filter() -> AreaFilter:
def check_single(
area_filter: AreaFilter, name: str, lat: float, lon: float, alt: float = 0.0
) -> bool:
return bool(area_filter.checkInside(name, np.array([lat]), np.array([lon]), np.array([alt]))[0])
return bool(area_filter.contains(name, np.array([lat]), np.array([lon]), np.array([alt]))[0])


class TestDefineArea:
Expand All @@ -30,7 +30,7 @@ def test_unknown_area_absent(self, area_filter: AreaFilter) -> None:
assert not area_filter.has_area("NOPE")

def test_checkinside_unknown_area_returns_false(self, area_filter: AreaFilter) -> None:
result = area_filter.checkInside("NOPE", np.array([52.0]), np.array([4.0]), np.array([0.0]))
result = area_filter.contains("NOPE", np.array([52.0]), np.array([4.0]), np.array([0.0]))
assert not result.any()

def test_reset_clears_areas(self, area_filter: AreaFilter) -> None:
Expand Down Expand Up @@ -123,7 +123,7 @@ def test_array_input(self, area_filter: AreaFilter) -> None:
lat = np.array([52.5, 51.0, 52.9])
lon = np.array([4.5, 4.5, 4.1])
alt = np.zeros(3)
inside = area_filter.checkInside("B", lat, lon, alt)
inside = area_filter.contains("B", lat, lon, alt)
assert inside.tolist() == [True, False, True]


Expand All @@ -137,10 +137,34 @@ def test_center_inside_far_point_outside(self, area_filter: AreaFilter) -> None:
# 2 deg lat is about 120 NM: outside
assert not check_single(area_filter, "C", 54.0, 4.0)

def test_circle_near_pole(self, area_filter: AreaFilter) -> None:
# 100 NM radius centred close to the north pole must remain valid
result = area_filter.define_area("C", "CIRCLE", [89.9, 0.0, 100.0])
assert result.is_ok()
assert check_single(area_filter, "C", 89.9, 0.0)
# 2 deg of latitude south is about 120 NM: outside
assert not check_single(area_filter, "C", 87.9, 0.0)


class TestPoly:
def test_triangle_centroid_inside(self, area_filter: AreaFilter) -> None:
# Triangle (52,4) (53,4) (52.5,5)
area_filter.define_area("P", "POLY", [52.0, 4.0, 53.0, 4.0, 52.5, 5.0])
assert check_single(area_filter, "P", 52.5, 4.3)
assert not check_single(area_filter, "P", 52.5, 5.5)

def test_antimeridian_polygon_rejected(self, area_filter: AreaFilter) -> None:
# Quad spanning lon 170 to -170 across the antimeridian
result = area_filter.define_area(
"P", "POLY", [10.0, 170.0, 10.0, -170.0, -10.0, -170.0, -10.0, 170.0]
)
assert result.is_err()
assert not area_filter.has_area("P")


class TestLine:
def test_line_has_no_inside(self, area_filter: AreaFilter) -> None:
area_filter.define_area("L", "LINE", [52.0, 4.0, 53.0, 5.0])
assert area_filter.has_area("L")
# A point exactly on the line is still not "inside" it
assert not check_single(area_filter, "L", 52.5, 4.5)
Loading