Skip to content

Commit

Permalink
update bounds subset, add latitude and longitude interval cross section
Browse files Browse the repository at this point in the history
  • Loading branch information
philipc2 committed Dec 19, 2024
1 parent 79275c5 commit 7770e01
Show file tree
Hide file tree
Showing 6 changed files with 173 additions and 97 deletions.
20 changes: 14 additions & 6 deletions uxarray/cross_sections/dataarray_accessor.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
from __future__ import annotations


from typing import TYPE_CHECKING
from typing import TYPE_CHECKING, Tuple

if TYPE_CHECKING:
pass
Expand All @@ -19,6 +19,8 @@ def __repr__(self):

methods_heading += " * constant_latitude(lat)\n"
methods_heading += " * constant_longitude(lon)\n"
methods_heading += " * latitude_interval(lats)\n"
methods_heading += " * longitude_interval(lons)\n"

return prefix + methods_heading

Expand Down Expand Up @@ -104,13 +106,19 @@ def constant_longitude(self, lon: float):

return self.uxda.isel(n_face=faces)

def gca(self, *args, **kwargs):
raise NotImplementedError
def latitude_interval(self, lats: Tuple[float, float], *args, **kwargs):
faces = self.uxda.uxgrid.get_faces_between_latitudes(
lats,
)

def bounded_latitude(self, *args, **kwargs):
raise NotImplementedError
return self.uxda.isel(n_face=faces)

def bounded_longitude(self, *args, **kwargs):
def longitude_interval(self, lons: Tuple[float, float], *args, **kwargs):
faces = self.uxda.uxgrid.get_faces_between_longitudes(lons)

return self.uxda.isel(n_face=faces)

def gca(self, *args, **kwargs):
raise NotImplementedError

def gca_gca(self, *args, **kwargs):
Expand Down
38 changes: 33 additions & 5 deletions uxarray/cross_sections/grid_accessor.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
from __future__ import annotations


from typing import TYPE_CHECKING
from typing import TYPE_CHECKING, Tuple

if TYPE_CHECKING:
from uxarray.grid import Grid
Expand All @@ -19,6 +19,8 @@ def __repr__(self):

methods_heading += " * constant_latitude(lat, return_face_indices)\n"
methods_heading += " * constant_longitude(lon, return_face_indices)\n"
methods_heading += " * latitude_interval(lats, return_face_indices)\n"
methods_heading += " * longitude_interval(lons, return_face_indices)\n"
return prefix + methods_heading

def constant_latitude(
Expand Down Expand Up @@ -127,11 +129,37 @@ def constant_longitude(
def gca(self, *args, **kwargs):
raise NotImplementedError

def bounded_latitude(self, *args, **kwargs):
raise NotImplementedError
def latitude_interval(
self,
lats: Tuple[float, float],
return_face_indices: bool = False,
*args,
**kwargs,
):
faces = self.uxgrid.get_faces_between_latitudes(lats)

def bounded_longitude(self, *args, **kwargs):
raise NotImplementedError
grid_between_lats = self.uxgrid.isel(n_face=faces)

if return_face_indices:
return grid_between_lats, faces
else:
return grid_between_lats

def longitude_interval(
self,
lons: Tuple[float, float],
return_face_indices: bool = False,
*args,
**kwargs,
):
faces = self.uxgrid.get_faces_between_longitudes(lons)

grid_between_lons = self.uxgrid.isel(n_face=faces)

if return_face_indices:
return grid_between_lons, faces
else:
return grid_between_lons

def gca_gca(self, *args, **kwargs):
raise NotImplementedError
9 changes: 9 additions & 0 deletions uxarray/grid/grid.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
from typing import (
Optional,
Union,
Tuple,
)

# reader and writer imports
Expand Down Expand Up @@ -79,6 +80,8 @@
constant_lon_intersections_no_extreme,
constant_lat_intersections_face_bounds,
constant_lon_intersections_face_bounds,
faces_within_lon_bounds,
faces_within_lat_bounds,
)

from spatialpandas import GeoDataFrame
Expand Down Expand Up @@ -2362,3 +2365,9 @@ def get_faces_at_constant_longitude(self, lon: float):

faces = constant_lon_intersections_face_bounds(lon, self.face_bounds_lon.values)
return faces

def get_faces_between_longitudes(self, lons: Tuple[float, float]):
return faces_within_lon_bounds(lons, self.face_bounds_lon.values)

def get_faces_between_latitudes(self, lats: Tuple[float, float]):
return faces_within_lat_bounds(lats, self.face_bounds_lat.values)
135 changes: 110 additions & 25 deletions uxarray/grid/intersections.py
Original file line number Diff line number Diff line change
Expand Up @@ -129,28 +129,26 @@ def constant_lon_intersections_no_extreme(lon, edge_node_x, edge_node_y, n_edge)
return np.unique(intersecting_edges)


# @njit(cache=True)
def constant_lat_intersections_face_bounds(lat, face_bounds_lat):
"""Identifies the candidate faces on a grid that intersect with a given
constant latitude.
This function checks whether the specified latitude, `lat`, in degrees lies within
the latitude bounds of grid faces, defined by `face_min_lat_rad` and `face_max_lat_rad`,
which are given in radians. The function returns the indices of the faces where the
latitude is within these bounds.
@njit(cache=True)
def constant_lat_intersections_face_bounds(lat: float, face_bounds_lat: np.ndarray):
"""
Identify candidate faces that intersect with a given constant latitude line.
Parameters
----------
lat : float
The latitude in degrees for which to find intersecting faces.
TODO:
face_bounds_lat : numpy.ndarray
A 2D array of shape (n_faces, 2), where each row represents the latitude
bounds of a face. The first element of each row is the minimum latitude
and the second element is the maximum latitude of the face.
Returns
-------
candidate_faces : numpy.ndarray
A 1D array containing the indices of the faces that intersect with the given latitude.
A 1D array of integers containing the indices of the faces that intersect
the given latitude.
"""

face_bounds_lat_min = face_bounds_lat[:, 0]
face_bounds_lat_max = face_bounds_lat[:, 1]

Expand All @@ -160,27 +158,25 @@ def constant_lat_intersections_face_bounds(lat, face_bounds_lat):


@njit(cache=True)
def constant_lon_intersections_face_bounds(lon, face_bounds_lon):
"""Identifies the candidate faces on a grid that intersect with a given
constant longitude.
This function checks whether the specified longitude, `lon`, in degrees lies within
the longitude bounds of grid faces, defined by `face_min_lon_rad` and `face_max_lon_rad`,
which are given in radians. The function returns the indices of the faces where the
longitude is within these bounds.
def constant_lon_intersections_face_bounds(lon: float, face_bounds_lon: np.ndarray):
"""
Identify candidate faces that intersect with a given constant longitude line.
Parameters
----------
lon : float
The longitude in degrees for which to find intersecting faces.
TODO:
face_bounds_lon : numpy.ndarray
A 2D array of shape (n_faces, 2), where each row represents the longitude
bounds of a face. The first element of each row is the minimum longitude
and the second element is the maximum longitude of the face.
Returns
-------
candidate_faces : numpy.ndarray
A 1D array containing the indices of the faces that intersect with the given longitude.
A 1D array of integers containing the indices of the faces that intersect
the given longitude.
"""

face_bounds_lon_min = face_bounds_lon[:, 0]
face_bounds_lon_max = face_bounds_lon[:, 1]
n_face = face_bounds_lon.shape[0]
Expand All @@ -191,16 +187,105 @@ def constant_lon_intersections_face_bounds(lon, face_bounds_lon):
cur_face_bounds_lon_max = face_bounds_lon_max[i]

if cur_face_bounds_lon_min < cur_face_bounds_lon_max:
if (lon >= cur_face_bounds_lon_min) & (lon <= cur_face_bounds_lon_max):
if (lon >= cur_face_bounds_lon_min) and (lon <= cur_face_bounds_lon_max):
candidate_faces.append(i)
else:
# antimeridian case
if (lon >= cur_face_bounds_lon_min) | (lon <= cur_face_bounds_lon_max):
if (lon >= cur_face_bounds_lon_min) or (lon <= cur_face_bounds_lon_max):
candidate_faces.append(i)

return np.array(candidate_faces, dtype=INT_DTYPE)


@njit(cache=True)
def faces_within_lon_bounds(lons, face_bounds_lon):
"""
Identify candidate faces that lie within a specified longitudinal interval.
Parameters
----------
lons : tuple or list of length 2
A pair (min_lon, max_lon) specifying the query interval. If `min_lon <= max_lon`,
the interval is [min_lon, max_lon]. If `min_lon > max_lon`, the interval
crosses the antimeridian and should be interpreted as [min_lon, 180] U [-180, max_lon].
face_bounds_lon : numpy.ndarray
A 2D array of shape (n_faces, 2), where each row represents the longitude bounds
of a face. The first element is the minimum longitude and the second is the maximum
longitude for that face. Bounds may cross the antimeridian.
Returns
-------
candidate_faces : numpy.ndarray
A 1D array of integers containing the indices of the faces whose longitude bounds
overlap with the specified interval.
"""
face_bounds_lon_min = face_bounds_lon[:, 0]
face_bounds_lon_max = face_bounds_lon[:, 1]
n_face = face_bounds_lon.shape[0]

min_lon, max_lon = lons

# For example, a query of (160, -160) would cross the antimeridian
antimeridian = min_lon > max_lon

candidate_faces = []
for i in range(n_face):
cur_face_min = face_bounds_lon_min[i]
cur_face_max = face_bounds_lon_max[i]

if not antimeridian:
# Normal case: min_lon <= max_lon
# Check if face overlaps with [min_lon, max_lon]
if (cur_face_max >= min_lon) and (cur_face_min <= max_lon):
candidate_faces.append(i)
else:
# Antimeridian case: interval crosses the -180/180 boundary
# The query interval is effectively [min_lon, 180] U [-180, max_lon]

# Check overlap with [min_lon, 180]
overlap_part1 = (cur_face_max >= min_lon) and (cur_face_min <= 180)

# Check overlap with [-180, max_lon]
overlap_part2 = (cur_face_max >= -180) and (cur_face_min <= max_lon)

if overlap_part1 or overlap_part2:
candidate_faces.append(i)

return np.array(candidate_faces, dtype=INT_DTYPE)


@njit(cache=True)
def faces_within_lat_bounds(lats, face_bounds_lat):
"""
Identify candidate faces that lie within a specified latitudinal interval.
Parameters
----------
lats : tuple or list of length 2
A pair (min_lat, max_lat) specifying the query interval. All returned faces
must be fully contained within this interval.
face_bounds_lat : numpy.ndarray
A 2D array of shape (n_faces, 2), where each row represents the latitude
bounds of a face. The first element is the minimum latitude and the second
is the maximum latitude for that face.
Returns
-------
candidate_faces : numpy.ndarray
A 1D array of integers containing the indices of the faces whose latitude
bounds lie completely within the specified interval.
"""

min_lat, max_lat = lats

face_bounds_lat_min = face_bounds_lat[:, 0]
face_bounds_lat_max = face_bounds_lat[:, 1]

within_bounds = (face_bounds_lat_max <= max_lat) & (face_bounds_lat_min >= min_lat)
candidate_faces = np.where(within_bounds)[0]
return candidate_faces


def _gca_gca_intersection_cartesian(gca_a_xyz, gca_b_xyz):
gca_a_xyz = np.asarray(gca_a_xyz)
gca_b_xyz = np.asarray(gca_b_xyz)
Expand Down
6 changes: 1 addition & 5 deletions uxarray/subset/dataarray_accessor.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,6 @@ def bounding_box(
self,
lon_bounds: Union[Tuple, List, np.ndarray],
lat_bounds: Union[Tuple, List, np.ndarray],
element: Optional[str] = "nodes",
method: Optional[str] = "coords",
**kwargs,
):
"""Subsets an unstructured grid between two latitude and longitude
Expand All @@ -54,9 +52,7 @@ def bounding_box(
element: str
Element for use with `coords` comparison, one of `nodes`, `face centers`, or `edge centers`
"""
grid = self.uxda.uxgrid.subset.bounding_box(
lon_bounds, lat_bounds, element, method
)
grid = self.uxda.uxgrid.subset.bounding_box(lon_bounds, lat_bounds)

return self.uxda._slice_from_grid(grid)

Expand Down
Loading

0 comments on commit 7770e01

Please sign in to comment.