Skip to content

Commit

Permalink
Merge branch 'main' into ahijevyc/neighborhood_filter
Browse files Browse the repository at this point in the history
  • Loading branch information
ahijevyc committed Sep 17, 2024
2 parents d6d8a33 + 92f38db commit 47b9cda
Show file tree
Hide file tree
Showing 32 changed files with 1,364 additions and 730 deletions.
2 changes: 1 addition & 1 deletion benchmarks/asv.conf.json
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@
"setuptools_scm": [""],
"xarray": [""],
"netcdf4": [""],
"pip+pyfma": [""]
"pip+git+https://github.com/philipc2/pyfma.git": [""]
},


Expand Down
10 changes: 10 additions & 0 deletions benchmarks/face_bounds.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,3 +28,13 @@ def time_face_bounds(self, grid_path):
def peakmem_face_bounds(self, grid_path):
"""Peak memory usage obtain ``Grid.face_bounds."""
face_bounds = self.uxgrid.bounds

class Bounds:
def setup(self):
self.uxgrid = ux.open_grid(r"C:\Users\chmie\PycharmProjects\ncar-uxarray\uxarray-hongyu\benchmarks\oQU480.grid.nc")

def teardown(self):
del self.uxgrid

def time_bounds(self):
self.uxgrid.bounds
96 changes: 38 additions & 58 deletions benchmarks/mpas_ocean.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,106 +27,81 @@
"120km": [current_path / grid_filename_120, current_path / data_filename_120]}


class Gradient:

param_names = ['resolution']
params = ['480km', '120km']
class DatasetBenchmark:
"""Class used as a template for benchmarks requiring a ``UxDataset`` in
this module across both resolutions."""
param_names = ['resolution',]
params = [['480km', '120km'],]


def setup(self, resolution, *args, **kwargs):


def setup(self, resolution):
self.uxds = ux.open_dataset(file_path_dict[resolution][0], file_path_dict[resolution][1])

def teardown(self, resolution):
def teardown(self, resolution, *args, **kwargs):
del self.uxds

def time_gradient(self, resolution):
self.uxds[data_var].gradient()
class GridBenchmark:
"""Class used as a template for benchmarks requiring a ``Grid`` in this
module across both resolutions."""
param_names = ['resolution', ]
params = [['480km', '120km'], ]

def peakmem_gradient(self, resolution):
grad = self.uxds[data_var].gradient()
def setup(self, resolution, *args, **kwargs):

class Integrate:
self.uxgrid = ux.open_grid(file_path_dict[resolution][0])

param_names = ['resolution']
params = ['480km', '120km']
def teardown(self, resolution, *args, **kwargs):
del self.uxgrid


def setup(self, resolution):
self.uxds = ux.open_dataset(file_path_dict[resolution][0], file_path_dict[resolution][1])

def teardown(self, resolution):
del self.uxds
class Gradient(DatasetBenchmark):
def time_gradient(self, resolution):
self.uxds[data_var].gradient()

def peakmem_gradient(self, resolution):
grad = self.uxds[data_var].gradient()

class Integrate(DatasetBenchmark):
def time_integrate(self, resolution):
self.uxds[data_var].integrate()

def peakmem_integrate(self, resolution):
integral = self.uxds[data_var].integrate()

class GeoDataFrame:

param_names = ['resolution', 'exclude_antimeridian']
params = [['480km', '120km'],
[True, False]]
class GeoDataFrame(DatasetBenchmark):
param_names = DatasetBenchmark.param_names + ['exclude_antimeridian']
params = DatasetBenchmark.params + [[True, False]]


def setup(self, resolution, exclude_antimeridian):
self.uxds = ux.open_dataset(file_path_dict[resolution][0], file_path_dict[resolution][1])

def teardown(self, resolution, exclude_antimeridian):
del self.uxds

def time_to_geodataframe(self, resolution, exclude_antimeridian):
self.uxds[data_var].to_geodataframe(exclude_antimeridian=exclude_antimeridian)

def peakmem_to_geodataframe(self, resolution, exclude_antimeridian):
gdf = self.uxds[data_var].to_geodataframe(exclude_antimeridian=exclude_antimeridian)


class ConnectivityConstruction:

param_names = ['resolution']
params = ['480km', '120km']


def setup(self, resolution):
self.uxds = ux.open_dataset(file_path_dict[resolution][0], file_path_dict[resolution][1])


def teardown(self, resolution):
del self.uxds

class ConnectivityConstruction(DatasetBenchmark):
def time_n_nodes_per_face(self, resolution):
self.uxds.uxgrid.n_nodes_per_face

def time_face_face_connectivity(self, resolution):
ux.grid.connectivity._populate_face_face_connectivity(self.uxds.uxgrid)


class MatplotlibConversion:
param_names = ['resolution', 'periodic_elements']
params = (['480km', '120km'], ['include', 'exclude', 'split'])

def setup(self, resolution, periodic_elements):
self.uxds = ux.open_dataset(file_path_dict[resolution][0], file_path_dict[resolution][1])

def teardown(self, resolution, periodic_elements):
del self.uxds
class MatplotlibConversion(DatasetBenchmark):
param_names = DatasetBenchmark.param_names + ['periodic_elements']
params = DatasetBenchmark.params + [['include', 'exclude', 'split']]

def time_dataarray_to_polycollection(self, resolution, periodic_elements):
self.uxds[data_var].to_polycollection()


class ConstructTreeStructures:
param_names = ['resolution']
params = ['480km', '120km']

def setup(self, resolution):
self.uxds = ux.open_dataset(file_path_dict[resolution][0], file_path_dict[resolution][1])

def teardown(self, resolution):
del self.uxds

class ConstructTreeStructures(DatasetBenchmark):
def time_kd_tree(self, resolution):
self.uxds.uxgrid.get_kd_tree()

Expand Down Expand Up @@ -164,3 +139,8 @@ def time_nearest_neighbor_remapping(self):

def time_inverse_distance_weighted_remapping(self):
self.uxds_480["bottomDepth"].remap.inverse_distance_weighted(self.uxds_120.uxgrid)


class HoleEdgeIndices(DatasetBenchmark):
def time_construct_hole_edge_indices(self, resolution):
ux.grid.geometry._construct_hole_edge_indices(self.uxds.uxgrid.edge_face_connectivity)
1 change: 1 addition & 0 deletions docs/internal_api/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -141,6 +141,7 @@ Geometry
grid.geometry._insert_pt_in_latlonbox
grid.geometry._populate_face_latlon_bound
grid.geometry._populate_bounds
grid.geometry._construct_hole_edge_indices

Coordinates
-----------
Expand Down
24 changes: 19 additions & 5 deletions docs/user_api/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -54,9 +54,16 @@ Remapping
---------
.. autosummary::
:toctree: generated/
:template: autosummary/accessor.rst

UxDataArray.remap

.. autosummary::
:toctree: generated/
:template: autosummary/accessor.rst

UxDataset.nearest_neighbor_remap
UxDataset.inverse_distance_weighted_remap
UxDataset.remap.nearest_neighbor
UxDataset.remap.inverse_distance_weighted

Plotting
--------
Expand Down Expand Up @@ -132,10 +139,16 @@ Remapping
---------
.. autosummary::
:toctree: generated/
:template: autosummary/accessor.rst

UxDataArray.remap

.. autosummary::
:toctree: generated/
:template: autosummary/accessor_method.rst

UxDataArray.nearest_neighbor_remap
UxDataArray.inverse_distance_weighted_remap
UxDataArray.nodal_average
UxDataArray.remap.nearest_neighbor
UxDataArray.remap.inverse_distance_weighted

Plotting
--------
Expand Down Expand Up @@ -210,6 +223,7 @@ IO
Grid.to_polycollection
Grid.to_linecollection
Grid.validate
Grid.hole_edge_indices


Methods
Expand Down
33 changes: 13 additions & 20 deletions test/test_arcs.py
Original file line number Diff line number Diff line change
Expand Up @@ -30,13 +30,6 @@ class TestIntersectionPoint(TestCase):

def test_pt_within_gcr(self):
# The GCR that's eexactly 180 degrees will have Value Error raised
gcr_180degree_cart = [
_lonlat_rad_to_xyz(0.0, 0.0),
_lonlat_rad_to_xyz(np.pi, 0.0)
]
pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0)
with self.assertRaises(ValueError):
point_within_gca(np.asarray(pt_same_lon_in), np.asarray(gcr_180degree_cart))

gcr_180degree_cart = [
_lonlat_rad_to_xyz(0.0, np.pi / 2.0),
Expand All @@ -45,27 +38,27 @@ def test_pt_within_gcr(self):

pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0)
with self.assertRaises(ValueError):
point_within_gca(np.asarray(pt_same_lon_in), np.asarray(gcr_180degree_cart))
point_within_gca(pt_same_lon_in, gcr_180degree_cart)

# Test when the point and the GCA all have the same longitude
gcr_same_lon_cart = [
_lonlat_rad_to_xyz(0.0, 1.5),
_lonlat_rad_to_xyz(0.0, -1.5)
]
pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0)
self.assertTrue(point_within_gca(np.asarray(pt_same_lon_in), np.asarray(gcr_same_lon_cart)))
self.assertTrue(point_within_gca(pt_same_lon_in, gcr_same_lon_cart))

pt_same_lon_out = _lonlat_rad_to_xyz(0.0, 1.500000000000001)
res = point_within_gca(np.asarray(pt_same_lon_out), np.asarray(gcr_same_lon_cart))
res = point_within_gca(pt_same_lon_out, gcr_same_lon_cart)
self.assertFalse(res)

pt_same_lon_out_2 = _lonlat_rad_to_xyz(0.1, 1.0)
res = point_within_gca(np.asarray(pt_same_lon_out_2), np.asarray(gcr_same_lon_cart))
res = point_within_gca(pt_same_lon_out_2, gcr_same_lon_cart)
self.assertFalse(res)

# And if we increase the digital place by one, it should be true again
pt_same_lon_out_add_one_place = _lonlat_rad_to_xyz(0.0, 1.5000000000000001)
res = point_within_gca(np.asarray(pt_same_lon_out_add_one_place), np.asarray(gcr_same_lon_cart))
res = point_within_gca(pt_same_lon_out_add_one_place, gcr_same_lon_cart)
self.assertTrue(res)

# Normal case
Expand All @@ -76,7 +69,7 @@ def test_pt_within_gcr(self):
-0.997]])
pt_cart_within = np.array(
[0.25616109352676675, 0.9246590335292105, -0.010021496695000144])
self.assertTrue(point_within_gca(np.asarray(pt_cart_within), np.asarray(gcr_cart_2), True))
self.assertTrue(point_within_gca(pt_cart_within, gcr_cart_2, True))

# Test other more complicate cases : The anti-meridian case

Expand All @@ -87,16 +80,16 @@ def test_pt_within_gcr_antimeridian(self):
gcr_cart = np.array([[0.351, -0.724, 0.593], [0.617, 0.672, 0.410]])
pt_cart = np.array(
[0.9438777657502077, 0.1193199333436068, 0.922714737029319])
self.assertTrue(point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart), is_directed=True))
self.assertTrue(point_within_gca(pt_cart, gcr_cart, is_directed=True))
# If we swap the gcr, it should throw a value error since it's larger than 180 degree
gcr_cart_flip = np.array([[0.617, 0.672, 0.410], [0.351, -0.724,
0.593]])
with self.assertRaises(ValueError):
point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart_flip), is_directed=True)
point_within_gca(pt_cart, gcr_cart_flip, is_directed=True)

# If we flip the gcr in the undirected mode, it should still work
self.assertTrue(
point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart_flip), is_directed=False))
point_within_gca(pt_cart, gcr_cart_flip, is_directed=False))

# 2nd anti-meridian case
# GCR vertex0 in radian : [4.104711496596806, 0.5352983676533828],
Expand All @@ -107,9 +100,9 @@ def test_pt_within_gcr_antimeridian(self):
pt_cart_within = np.array(
[0.6136726305712109, 0.28442243941920053, -0.365605190899831])
self.assertFalse(
point_within_gca(np.asarray(pt_cart_within), np.asarray(gcr_cart_1), is_directed=True))
point_within_gca(pt_cart_within, gcr_cart_1, is_directed=True))
self.assertFalse(
point_within_gca(np.asarray(pt_cart_within), np.asarray(gcr_cart_1), is_directed=False))
point_within_gca(pt_cart_within, gcr_cart_1, is_directed=False))

# The first case should not work and the second should work
v1_rad = [0.1, 0.0]
Expand All @@ -119,10 +112,10 @@ def test_pt_within_gcr_antimeridian(self):
gcr_cart = np.array([v1_cart, v2_cart])
pt_cart = _lonlat_rad_to_xyz(0.01, 0.0)
with self.assertRaises(ValueError):
point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart), is_directed=True)
point_within_gca(pt_cart, gcr_cart, is_directed=True)
gcr_car_flipped = np.array([v2_cart, v1_cart])
self.assertTrue(
point_within_gca(np.asarray(pt_cart), np.asarray(gcr_car_flipped), is_directed=True))
point_within_gca(pt_cart, gcr_car_flipped, is_directed=True))

def test_pt_within_gcr_cross_pole(self):
gcr_cart = np.array([[0.351, 0.0, 0.3], [-0.351, 0.0, 0.3]])
Expand Down
26 changes: 0 additions & 26 deletions test/test_dataarray.py
Original file line number Diff line number Diff line change
Expand Up @@ -113,29 +113,3 @@ def test_geodataframe_caching(self):

# override will recompute the grid
assert gdf_start is not gdf_end

def test_nodal_average(self):

# test on a node-centered dataset
uxds = ux.open_dataset(gridfile_geoflow, dsfile_v1_geoflow)

v1_nodal_average = uxds['v1'].nodal_average()

# final dimension should match number of faces
self.assertEqual(v1_nodal_average.shape[-1], uxds.uxgrid.n_face)

# all other dimensions should remain unchanged
self.assertEqual(uxds['v1'].shape[0:-1], v1_nodal_average.shape[0:-1])

# test on a sample mesh with 4 verts
verts = [[[-170, 40], [180, 30], [165, 25], [-170, 20]]]
data = [1, 2, 3, 4]

uxgrid = ux.open_grid(verts, latlon=True)

uxda = ux.UxDataArray(uxgrid=uxgrid, data=data, dims=('n_node'))

uxda_nodal_average = uxda.nodal_average()

# resulting data should be the mean of the corner nodes of the single face
self.assertEqual(uxda_nodal_average, np.mean(data))
Loading

0 comments on commit 47b9cda

Please sign in to comment.