Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

use s2geometry for fundamental operations #284

Draft
wants to merge 2 commits into
base: master
Choose a base branch
from
Draft
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
1 change: 1 addition & 0 deletions pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@ keywords = [
dependencies = [
"numpy>=1.23",
"astropy>=5.0.4",
"s2geometry @ git+https://github.com/zacharyburnett/s2geometry@ci/build",
]

[project.optional-dependencies]
Expand Down
31 changes: 19 additions & 12 deletions spherical_geometry/polygon.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@

# THIRD-PARTY
import numpy as np
import s2geometry

# LOCAL
from spherical_geometry import great_circle_arc, vector
Expand Down Expand Up @@ -64,6 +65,11 @@ def __init__(self, points, inside=None):
raise ValueError("Polygon made of too few points")

self._points = points = np.asanyarray(points)

s2loop = s2geometry.S2Loop()
s2loop.Init([s2geometry.S2Point(*point) for point in points])
self._s2polygon = s2geometry.S2Polygon(s2loop)

new_inside = self._find_new_inside()

if inside is None:
Expand Down Expand Up @@ -381,7 +387,7 @@ def convex_hull(cls, points):

return cls(hull, inside)

def invert_polygon(self):
def invert_polygon(self) -> "SingleSphericalPolygon":
"""
Compute the inverse (complement) of a single polygon
"""
Expand All @@ -390,15 +396,6 @@ def invert_polygon(self):
poly._inside = np.asanyarray(self._find_new_outside())
return poly

def _contains_point(self, point, P, r):
point = np.asanyarray(point)
if np.array_equal(r, point):
return True

intersects = great_circle_arc.intersects(P[:-1], P[1:], r, point)
crossings = np.sum(intersects)
return (crossings % 2) == 0

def contains_point(self, point):
r"""
Determines if this `SingleSphericalPolygon` contains a given point.
Expand All @@ -413,7 +410,7 @@ def contains_point(self, point):
contains : bool
Returns `True` if the polygon contains the given *point*.
"""
return self._contains_point(point, self._points, self._inside)
return polygon_contains_point(point, self)

def contains_lonlat(self, lon, lat, degrees=True):
r"""
Expand All @@ -435,7 +432,7 @@ def contains_lonlat(self, lon, lat, degrees=True):
Returns `True` if the polygon contains the given *point*.
"""
point = vector.lonlat_to_vector(lon, lat, degrees=degrees)
return self._contains_point(point, self._points, self._inside)
return polygon_contains_point(point, self)

# Alias for contains_lonlat
contains_radec = contains_lonlat
Expand Down Expand Up @@ -1286,3 +1283,13 @@ def draw(self, m, **plot_args):
"""
for polygon in self._polygons:
polygon.draw(m, **plot_args)

def polygon_contains_point(point: tuple[float, float, float], polygon: SingleSphericalPolygon) -> bool:
point = np.asanyarray(point)
if np.array_equal(polygon._inside, point):
return True

intersects = great_circle_arc.intersects(polygon._points[:-1], polygon._points[1:], polygon._inside, point)
crossings = np.sum(intersects)
return (crossings % 2) == 0

Loading