Skip to content

Commit

Permalink
Revert "Added new selection based on proximity to a boundary"
Browse files Browse the repository at this point in the history
This reverts commit c8c277d.
  • Loading branch information
tomdurrant committed Dec 5, 2023
1 parent c8c277d commit b313f68
Show file tree
Hide file tree
Showing 3 changed files with 36 additions and 115 deletions.
61 changes: 20 additions & 41 deletions tests/core/test_sel.py
Original file line number Diff line number Diff line change
@@ -1,11 +1,10 @@
"""Unit testing for SpecDataset wrapper around DataArray."""
import os

import numpy as np
import pytest
import numpy as np

from wavespectra import read_era5, read_ww3
from wavespectra.core.attributes import attrs
from wavespectra import read_ww3, read_era5
from wavespectra.core.select import Coordinates

FILES_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "../sample_files")
Expand Down Expand Up @@ -49,16 +48,6 @@ def test_dset_sel_bbox(self):
)
dset[attrs.SITENAME].size == self.dset[attrs.SITENAME].size

def test_dset_sel_boundary(self):
"""Assert that sel/bbox method runs."""
dset = self.dset.spec.sel(
lons=self.dset.lon.values,
lats=self.dset.lat.values,
method="boundary",
tolerance=0.0,
)
dset[attrs.SITENAME].size == self.dset[attrs.SITENAME].size

def test_dset_sel_bbox_outside(self):
"""Assert that sel/bbox method raises for bbox outside dataset."""
with pytest.raises(ValueError):
Expand Down Expand Up @@ -134,7 +123,6 @@ def test_weighting(self):
upper = np.array([max(s1, s2) for s1, s2 in zip(site0, site1)])
assert (upper - idw > 0).all() and (idw - lower > 0).all()


class TestSelCoordinatesConventions:
"""Test Sel methods with different coordinates conventions."""

Expand Down Expand Up @@ -164,15 +152,25 @@ def test_nearest_both_180(self):
dset = self.dset.copy(deep=True)
dset["lon"].values = [-10, 10]
dset["lat"].values = [30, 30]
ds = dset.spec.sel(method="nearest", lons=[-9], lats=[31], tolerance=5.0)
ds = dset.spec.sel(
method="nearest",
lons=[-9],
lats=[31],
tolerance=5.0
)
assert ds.lon == -10

def test_nearest_dset_360_slice_180(self):
"""Test nearest with Dataset in [0 <--> 360] and slice in [-180 <--> 180]."""
dset = self.dset.copy(deep=True)
dset["lon"].values = [0, 350]
dset["lat"].values = [30, 30]
ds = dset.spec.sel(method="nearest", lons=[-9], lats=[31], tolerance=5.0)
ds = dset.spec.sel(
method="nearest",
lons=[-9],
lats=[31],
tolerance=5.0
)
assert ds.lon == -10

def test_nearest_dset_180_slice_360(self):
Expand All @@ -187,7 +185,7 @@ def test_nearest_dset_180_slice_360(self):
lats=[31],
tolerance=5.0,
dset_lons=dset.lon.values,
dset_lats=dset.lat.values,
dset_lats=dset.lat.values
)
assert ds.lon == 350

Expand Down Expand Up @@ -249,24 +247,7 @@ def test_bbox_both_180(self):
lats=[-35, 35],
tolerance=0.0,
dset_lons=dset.lon.values,
dset_lats=dset.lat.values,
)
assert np.array_equal(ds.lon.values, dset.lon.values)

def test_bbox_both_180(self):
"""Test bbox with both dataset and slice in [-180 <--> 180]."""
dset = self.dset.copy(deep=True)
dset["lon"].values = [-10, 10]
dset["lat"].values = [30, 30]
lon = -9

ds = dset.spec.sel(
method="boundary",
lons=[-10, 10],
lats=[-35, 35],
tolerance=50,
dset_lons=dset.lon.values,
dset_lats=dset.lat.values,
dset_lats=dset.lat.values
)
assert np.array_equal(ds.lon.values, dset.lon.values)

Expand All @@ -281,11 +262,9 @@ def test_bbox_dset_360_slice_180(self):
lats=[-35, 35],
tolerance=0.0,
dset_lons=dset.lon.values,
dset_lats=dset.lat.values,
)
assert np.array_equal(
sorted(ds.lon.values % 360), sorted(dset.lon.values % 360)
dset_lats=dset.lat.values
)
assert np.array_equal(sorted(ds.lon.values % 360), sorted(dset.lon.values % 360))

def test_bbox_dset_180_slice_360(self):
"""Test bbox with Dataset in [-180 <--> 180] and slice in [0 <--> 360]."""
Expand All @@ -298,6 +277,6 @@ def test_bbox_dset_180_slice_360(self):
lats=[-35, 35],
tolerance=0.0,
dset_lons=dset.lon.values,
dset_lats=dset.lat.values,
dset_lats=dset.lat.values
)
assert int(ds.lon) == 350
assert int(ds.lon) == 350
80 changes: 12 additions & 68 deletions wavespectra/core/select.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,11 @@
"""Interpolate stations."""
import logging

import numpy as np
import xarray as xr
from scipy.spatial import cKDTree
import logging

from wavespectra.core.attributes import attrs, set_spec_attributes


logger = logging.getLogger(__name__)


Expand Down Expand Up @@ -92,10 +91,7 @@ def distance(self, lon, lat):
List of distances between each station and site.
"""
dist = np.sqrt(
(self.dset_lons % 360 - np.array(lon) % 360) ** 2
+ (self.dset_lats - np.array(lat)) ** 2
)
dist = np.sqrt((self.dset_lons % 360 - np.array(lon) % 360) ** 2 + (self.dset_lats - np.array(lat)) ** 2)
if isinstance(dist, xr.DataArray):
dist = dist.values
return np.abs(dist)
Expand Down Expand Up @@ -123,12 +119,12 @@ def nearer(self, lon, lat, tolerance=np.inf, max_sites=None):
def nearest(self, lon, lat):
"""Nearest station in (dset_lons, dset_lats) to site (lon, lat).
Args:
lon (float): Longitude to locate from lons.
lat (float): Latitude to locate from lats.
Args:
lon (float): Longitude to locate from lons.
lat (float): Latitude to locate from lats.
Returns:
Index and distance of closest station.
Returns:
Index and distance of closest station.
"""
dist = self.distance(lon, lat)
Expand Down Expand Up @@ -168,9 +164,7 @@ def sel_nearest(
improve precision if projected coordinates are provided at high latitudes.
"""
coords = Coordinates(
dset, lons=lons, lats=lats, dset_lons=dset_lons, dset_lats=dset_lats
)
coords = Coordinates(dset, lons=lons, lats=lats, dset_lons=dset_lons, dset_lats=dset_lats)

station_ids = []
for lon, lat in zip(coords.lons, coords.lats):
Expand Down Expand Up @@ -223,9 +217,7 @@ def sel_idw(
improve precision if projected coordinates are provided at high latitudes.
"""
coords = Coordinates(
dset, lons=lons, lats=lats, dset_lons=dset_lons, dset_lats=dset_lats
)
coords = Coordinates(dset, lons=lons, lats=lats, dset_lons=dset_lons, dset_lats=dset_lats)

mask = dset.isel(site=0, drop=True) * np.nan
dsout = []
Expand Down Expand Up @@ -302,9 +294,7 @@ def sel_bbox(dset, lons, lats, tolerance=0.0, dset_lons=None, dset_lats=None):
improve precision if projected coordinates are provided at high latitudes.
"""
coords = Coordinates(
dset, lons=lons, lats=lats, dset_lons=dset_lons, dset_lats=dset_lats
)
coords = Coordinates(dset, lons=lons, lats=lats, dset_lons=dset_lons, dset_lats=dset_lats)

minlon = min(coords.lons) - tolerance
minlat = min(coords.lats) - tolerance
Expand All @@ -331,7 +321,7 @@ def sel_bbox(dset, lons, lats, tolerance=0.0, dset_lons=None, dset_lats=None):
& (coords.dset_lats >= minlat)
& (coords.dset_lons <= minlon)
& (coords.dset_lats <= maxlat)
)[0],
)[0]
)

if station_ids.size == 0:
Expand All @@ -350,49 +340,3 @@ def sel_bbox(dset, lons, lats, tolerance=0.0, dset_lons=None, dset_lats=None):
dsout = dsout.assign_coords({attrs.SITENAME: np.arange(len(station_ids))})

return dsout


def sel_boundary(dset, lons, lats, tolerance=0.1, dset_lons=None, dset_lats=None):
"""Select sites from wthin the given dset that are within a given distance of any of the passed points
Primary use case is to extract spectral points from a parent spectral forcing dataset that
are within a certain distance of nest boundary.
Args:
dset (Dataset): Stations SpecDataset to select from.
lons (array): Longitude of sites
lats (array): Latitude of sites
tolerance (float): Minimum distance from any site (degrees)
dset_lons (array): Longitude of stations in dset.
dset_lats (array): Latitude of stations in dset.
Returns:
Selected SpecDataset point within the miniumum_disance any site passed
Note:
Args `dset_lons`, `dset_lats` are not required but can improve performance when
`dset` is chunked with site=1 (expensive to access station coordinates) and
improve precision if projected coordinates are provided at high latitudes.
"""
coords = Coordinates(
dset, lons=lons, lats=lats, dset_lons=dset_lons, dset_lats=dset_lats
)

# Create kdtree of boundary
input_points = list(zip(coords.dset_lons, coords.dset_lats))
boundary = list(zip(coords.lons, coords.lats))
tree = cKDTree(np.array(boundary))

# Find the indices of points in forcing dataset within the distance_threshold from any point in the boundary
within_distance_indices = tree.query_ball_point(np.array(input_points), tolerance)
result_list = [i for i, indices in enumerate(within_distance_indices) if indices]

dsout = dset.isel(**{attrs.SITENAME: result_list})

# Return longitudes in the convention provided
if coords.consistent is False:
dsout.lon.values = coords._swap_longitude_convention(dsout.lon.values)

dsout = dsout.assign_coords({attrs.SITENAME: np.arange(len(result_list))})

return dsout
10 changes: 4 additions & 6 deletions wavespectra/specdataset.py
Original file line number Diff line number Diff line change
@@ -1,13 +1,12 @@
"""Wrapper around the xarray dataset."""
import types
import os
import re
import sys
import types

import xarray as xr

from wavespectra.core.attributes import attrs
from wavespectra.core.select import sel_bbox, sel_boundary, sel_idw, sel_nearest
from wavespectra.core.select import sel_idw, sel_nearest, sel_bbox
from wavespectra.specarray import SpecArray

here = os.path.dirname(os.path.abspath(__file__))
Expand Down Expand Up @@ -112,7 +111,7 @@ def sel(
tolerance=2.0,
dset_lons=None,
dset_lats=None,
**kwargs,
**kwargs
):
"""Select stations near or at locations defined by (lons, lats) vector.
Expand Down Expand Up @@ -149,7 +148,6 @@ def sel(
"idw": sel_idw,
"bbox": sel_bbox,
"nearest": sel_nearest,
"boundary": sel_boundary,
None: sel_nearest,
}
try:
Expand All @@ -172,6 +170,6 @@ def sel(
tolerance=tolerance,
dset_lons=dset_lons,
dset_lats=dset_lats,
**kwargs,
**kwargs
)
return dsout

0 comments on commit b313f68

Please sign in to comment.