diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py index ef2dc7692..aadbd45a1 100644 --- a/benchmarks/mpas_ocean.py +++ b/benchmarks/mpas_ocean.py @@ -2,10 +2,10 @@ import urllib.request from pathlib import Path -import uxarray as ux - import numpy as np +import uxarray as ux + current_path = Path(os.path.dirname(os.path.realpath(__file__))) data_var = 'bottomDepth' @@ -33,8 +33,8 @@ class DatasetBenchmark: """Class used as a template for benchmarks requiring a ``UxDataset`` in this module across both resolutions.""" - param_names = ['resolution',] - params = [['480km', '120km'],] + param_names = ['resolution', ] + params = [['480km', '120km'], ] def setup(self, resolution, *args, **kwargs): self.uxds = ux.open_dataset(file_path_dict[resolution][0], file_path_dict[resolution][1]) @@ -42,6 +42,7 @@ def setup(self, resolution, *args, **kwargs): def teardown(self, resolution, *args, **kwargs): del self.uxds + class GridBenchmark: """Class used as a template for benchmarks requiring a ``Grid`` in this module across both resolutions.""" @@ -62,6 +63,7 @@ def time_gradient(self, resolution): def peakmem_gradient(self, resolution): grad = self.uxds[data_var].gradient() + class Integrate(DatasetBenchmark): def time_integrate(self, resolution): self.uxds[data_var].integrate() @@ -74,12 +76,10 @@ class GeoDataFrame(DatasetBenchmark): param_names = DatasetBenchmark.param_names + ['exclude_antimeridian'] params = DatasetBenchmark.params + [[True, False]] - def time_to_geodataframe(self, resolution, exclude_antimeridian): self.uxds[data_var].to_geodataframe(exclude_antimeridian=exclude_antimeridian) - class ConnectivityConstruction(DatasetBenchmark): def time_n_nodes_per_face(self, resolution): self.uxds.uxgrid.n_nodes_per_face @@ -136,6 +136,7 @@ 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_boundary_edge_indices(self.uxds.uxgrid.edge_face_connectivity) @@ -145,6 +146,7 @@ class DualMesh(DatasetBenchmark): def time_dual_mesh_construction(self, resolution): self.uxds.uxgrid.get_dual() + class ConstructFaceLatLon(GridBenchmark): def time_welzl(self, resolution): self.uxgrid.construct_face_centers(method='welzl') @@ -177,9 +179,26 @@ def setup(self, resolution, lat_step): self.lats = np.arange(-45, 45, lat_step) _ = self.uxgrid.bounds +class CrossSections(DatasetBenchmark): + param_names = DatasetBenchmark.param_names + ['n_lat'] + params = DatasetBenchmark.params + [[1, 2, 4, 8]] + + def time_constant_lat_fast(self, resolution, n_lat): + for lat in np.linspace(-89, 89, n_lat): + self.uxds.uxgrid.constant_latitude_cross_section(lat, method='fast') def teardown(self, resolution, lat_step): del self.uxgrid def time_const_lat(self, resolution, lat_step): for lat in self.lats: self.uxgrid.cross_section.constant_latitude(lat) + + +class ZonalAverage(DatasetBenchmark): + def setup(self, resolution, *args, **kwargs): + self.uxds = ux.open_dataset(file_path_dict[resolution][0], file_path_dict[resolution][1]) + bounds = self.uxds.uxgrid.bounds + + def time_zonal_average(self, resolution): + lat_step = 10 + self.uxds['bottomDepth'].zonal_mean(lat=(-45, 45, lat_step)) diff --git a/docs/api.rst b/docs/api.rst index 010f1c96d..a9b20d02e 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -197,6 +197,7 @@ I/O & Conversion UxDataArray.to_dataset UxDataArray.from_xarray + UxDataset ----------- @@ -265,6 +266,8 @@ UxDataArray UxDataArray.plot UxDataArray.plot.polygons UxDataArray.plot.points + UxDataArray.plot.line + UxDataArray.plot.scatter UxDataset ~~~~~~~~~ @@ -422,6 +425,15 @@ on each face. UxDataArray.topological_all UxDataArray.topological_any +Zonal Average +~~~~~~~~~~~~~ +.. autosummary:: + :toctree: generated/ + + UxDataArray.zonal_mean + + + Weighted ~~~~~~~~ .. autosummary:: diff --git a/docs/user-guide/zonal-average.ipynb b/docs/user-guide/zonal-average.ipynb new file mode 100644 index 000000000..d45b90e32 --- /dev/null +++ b/docs/user-guide/zonal-average.ipynb @@ -0,0 +1,180 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "dbade3aaedad9ae7", + "metadata": {}, + "source": [ + "# Zonal Averaging\n", + "\n", + "This section demonstrates how to perform Zonal Averaging using UXarray.\n" + ] + }, + { + "cell_type": "code", + "id": "185e2061bc4c75b9", + "metadata": {}, + "source": [ + "import uxarray as ux\n", + "import numpy as np\n", + "\n", + "uxds = ux.open_dataset(\n", + " \"../../test/meshfiles/ugrid/outCSne30/outCSne30.ug\",\n", + " \"../../test/meshfiles/ugrid/outCSne30/outCSne30_vortex.nc\",\n", + ")\n", + "uxds[\"psi\"].plot(cmap=\"inferno\", periodic_elements=\"split\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "d938d659b89bc9cb", + "metadata": {}, + "source": [ + "## What is a Zonal Average/Mean?\n", + "\n", + "A zonal average (or zonal mean) is a statistical measure that represents the average of a variable along one or more lines of constant latitude. In other words, it's the mean value calculated around the sphere at constant latitudes. \n", + "\n", + "UXarray currently implements a non-conservative Zonal Mean, which weights candidate faces by the length of intersection of a line of constant latitude.\n", + "\n", + "\n", + "```{seealso}\n", + "[NCL Zonal Average](https://www.ncl.ucar.edu/Applications/zonal.shtml)\n", + "```" + ] + }, + { + "cell_type": "code", + "id": "d342f8f449543994", + "metadata": {}, + "source": [ + "zonal_mean_psi = uxds[\"psi\"].zonal_mean()\n", + "zonal_mean_psi" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "65194a38c76e8a62", + "metadata": {}, + "source": [ + "The default latitude range is between -90 and 90 degrees with a step size of 10 degrees. " + ] + }, + { + "cell_type": "code", + "id": "b5933beaf2f598ab", + "metadata": {}, + "source": [ + "(zonal_mean_psi.plot.line() * zonal_mean_psi.plot.scatter(color=\"red\")).opts(\n", + " title=\"Zonal Average Plot (Default)\", xticks=np.arange(-90, 100, 20), xlim=(-95, 95)\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "ba2d641a3076c692", + "metadata": {}, + "source": [ + "The range of latitudes can be modified by using the `lat` parameter. It accepts:\n", + "\n", + "* **Single scalar**: e.g., `lat=45`\n", + "* **List/array**: e.g., `lat=[10, 20]` or `lat=np.array([10, 20])`\n", + "* **Tuple**: e.g., `(min_lat, max_lat, step)`" + ] + }, + { + "cell_type": "code", + "id": "4f665827daac1c46", + "metadata": {}, + "source": [ + "zonal_mean_psi_large = uxds[\"psi\"].zonal_mean(lat=(-90, 90, 1))" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "id": "1998f1a55b67100c", + "metadata": {}, + "source": [ + "(\n", + " zonal_mean_psi_large.plot.line()\n", + " * zonal_mean_psi_large.plot.scatter(color=\"red\", s=1)\n", + ").opts(\n", + " title=\"Zonal Average Plot (Larger Sample)\",\n", + " xticks=np.arange(-90, 100, 20),\n", + " xlim=(-95, 95),\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "id": "4b91ebb8f733a318", + "metadata": {}, + "source": [ + "## Combined Plots\n", + "\n", + "It is often desired to plot the zonal average along side other plots, such as color or contour plots. " + ] + }, + { + "cell_type": "code", + "id": "cb2255761173d53e", + "metadata": {}, + "source": [ + "(\n", + " uxds[\"psi\"].plot(\n", + " cmap=\"inferno\",\n", + " periodic_elements=\"split\",\n", + " height=250,\n", + " width=500,\n", + " colorbar=False,\n", + " ylim=(-90, 90),\n", + " )\n", + " + zonal_mean_psi.plot.line(\n", + " x=\"psi_zonal_mean\",\n", + " y=\"latitudes\",\n", + " height=250,\n", + " width=150,\n", + " ylabel=\"\",\n", + " ylim=(-90, 90),\n", + " xlim=(0.8, 1.2),\n", + " xticks=[0.8, 0.9, 1.0, 1.1, 1.2],\n", + " yticks=[-90, -45, 0, 45, 90],\n", + " grid=True,\n", + " )\n", + ").opts(title=\"Combined Zonal Average & Raster Plot\")" + ], + "outputs": [], + "execution_count": null + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.7" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/userguide.rst b/docs/userguide.rst index 5ebf96de9..708526f60 100644 --- a/docs/userguide.rst +++ b/docs/userguide.rst @@ -52,6 +52,9 @@ These user guides provide detailed explanations of the core functionality in UXa `Cross-Sections `_ Select cross-sections of a grid +`Zonal Means `_ + Compute the zonal averages across lines of constant latitude + `Remapping `_ Remap (a.k.a Regrid) between unstructured grids @@ -104,6 +107,7 @@ These user guides provide additional details about specific features in UXarray. user-guide/advanced-plotting.ipynb user-guide/subset.ipynb user-guide/cross-sections.ipynb + user-guide/zonal-average.ipynb user-guide/remapping.ipynb user-guide/topological-aggregations.ipynb user-guide/weighted_mean.ipynb diff --git a/test/test_cross_sections.py b/test/test_cross_sections.py index 8eaf035da..4942c86fb 100644 --- a/test/test_cross_sections.py +++ b/test/test_cross_sections.py @@ -1,9 +1,10 @@ -from pathlib import Path - -import numpy as np +import uxarray as ux import pytest +import numpy as np +from pathlib import Path +import os -import uxarray as ux +import numpy.testing as nt # Define the current path and file paths for grid and data current_path = Path(__file__).resolve().parent @@ -42,6 +43,7 @@ def test_constant_lon_cross_section_grid(): def test_constant_lat_cross_section_uxds(): uxds = ux.open_dataset(quad_hex_grid_path, quad_hex_data_path) + uxds.uxgrid.normalize_cartesian_coordinates() da_top_two = uxds['t2m'].cross_section.constant_latitude(lat=0.1) np.testing.assert_array_equal(da_top_two.data, uxds['t2m'].isel(n_face=[1, 2]).data) @@ -57,6 +59,7 @@ def test_constant_lat_cross_section_uxds(): def test_constant_lon_cross_section_uxds(): uxds = ux.open_dataset(quad_hex_grid_path, quad_hex_data_path) + uxds.uxgrid.normalize_cartesian_coordinates() da_left_two = uxds['t2m'].cross_section.constant_longitude(lon=-0.1) np.testing.assert_array_equal(da_left_two.data, uxds['t2m'].isel(n_face=[0, 2]).data) @@ -117,54 +120,18 @@ def test_constant_lat_out_of_bounds(): assert len(candidate_faces) == 0 +class TestArcs: + def test_latitude_along_arc(self): + node_lon = np.array([-40, -40, 40, 40]) + node_lat = np.array([-20, 20, 20, -20]) + face_node_connectivity = np.array([[0, 1, 2, 3]], dtype=np.int64) -def test_constant_longitude_interval_grid(): - uxgrid = ux.open_grid(quad_hex_grid_path) - - grid_all_four = uxgrid.cross_section.constant_longitude_interval(lons=(-1, 1)) - - - assert grid_all_four.n_face == 4 - - grid_left_two = uxgrid.cross_section.constant_longitude_interval(lons=(-1, 0.2)) - - assert grid_left_two.n_face == 2 - - grid_right_two = uxgrid.cross_section.constant_longitude_interval(lons=(-0.1, 1.0)) + uxgrid = ux.Grid.from_topology(node_lon, node_lat, face_node_connectivity) - assert grid_right_two.n_face == 2 + # intersection at exactly 20 degrees latitude + out1 = uxgrid.get_faces_at_constant_latitude(lat=20) -def test_constant_latitude_interval_grid(): - uxgrid = ux.open_grid(quad_hex_grid_path) + # intersection at 25.41 degrees latitude (max along the great circle arc) + out2 = uxgrid.get_faces_at_constant_latitude(lat=25.41) - grid_all_four = uxgrid.cross_section.constant_latitude_interval(lats=(-0.3, 0.3)) - - assert grid_all_four.n_face == 4 - - grid_top_two = uxgrid.cross_section.constant_latitude_interval(lats=(-0.1, 0.3)) - assert grid_top_two.n_face == 2 - - grid_bottom_two = uxgrid.cross_section.constant_latitude_interval(lats=(-0.3, 0.1)) - - assert grid_bottom_two.n_face == 2 - -def test_constant_latitude_interval_uxds(): - uxds = ux.open_dataset(quad_hex_grid_path, quad_hex_data_path) - - da_top_two = uxds['t2m'].cross_section.constant_latitude_interval(lats=(-0.1, 0.3)) - np.testing.assert_array_equal(da_top_two.data, uxds['t2m'].isel(n_face=[1, 2]).data) - - da_bottom_two = uxds['t2m'].cross_section.constant_latitude_interval(lats=(-0.3, 0.1)) - np.testing.assert_array_equal(da_bottom_two.data, uxds['t2m'].isel(n_face=[0, 3]).data) - - da_all_four = uxds['t2m'].cross_section.constant_latitude_interval(lats=(-0.3, 0.3)) - np.testing.assert_array_equal(da_all_four.data, uxds['t2m'].data) - -def test_constant_longitude_interval_uxds(): - uxds = ux.open_dataset(quad_hex_grid_path, quad_hex_data_path) - - da_left_two = uxds['t2m'].cross_section.constant_longitude_interval(lons=(-1, 0.2)) - np.testing.assert_array_equal(da_left_two.data, uxds['t2m'].isel(n_face=[0, 2]).data) - - da_right_two = uxds['t2m'].cross_section.constant_longitude_interval(lons=(-0.1, 1.0)) - np.testing.assert_array_equal(da_right_two.data, uxds['t2m'].isel(n_face=[1, 3]).data) + nt.assert_array_equal(out1, out2) diff --git a/test/test_plot.py b/test/test_plot.py index 24fa081a7..51f1ec32a 100644 --- a/test/test_plot.py +++ b/test/test_plot.py @@ -87,3 +87,13 @@ def test_dataarray_methods(): # plot.scatter() is an xarray method assert hasattr(uxds.plot, 'scatter') + +def test_line(): + uxds = ux.open_dataset(gridfile_mpas, gridfile_mpas) + _plot_line = uxds['bottomDepth'].zonal_average().plot.line() + assert isinstance(_plot_line, hv.Curve) + +def test_scatter(): + uxds = ux.open_dataset(gridfile_mpas, gridfile_mpas) + _plot_line = uxds['bottomDepth'].zonal_average().plot.scatter() + assert isinstance(_plot_line, hv.Scatter) diff --git a/test/test_zonal.py b/test/test_zonal.py new file mode 100644 index 000000000..6ee578f72 --- /dev/null +++ b/test/test_zonal.py @@ -0,0 +1,90 @@ +import os +from pathlib import Path + +import dask.array as da +import numpy as np +import pytest + +import numpy.testing as nt + +import uxarray as ux +from uxarray.constants import ERROR_TOLERANCE + +current_path = Path(os.path.dirname(os.path.realpath(__file__))) + + +class TestZonalCSne30: + gridfile_ne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug" + datafile_vortex_ne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30_vortex.nc" + dsfile_var2_ne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30_var2.nc" + test_file_2 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30_test2.nc" + test_file_3 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30_test3.nc" + + def test_non_conservative_zonal_mean_equator(self): + """Tests the zonal mean at the equator. This grid contains points that are exactly """ + grid_path = self.gridfile_ne30 + data_path = self.datafile_vortex_ne30 + uxds = ux.open_dataset(grid_path, data_path) + + res = uxds['psi'].zonal_mean(0) + + assert res.values[0] == pytest.approx(1, abs=ERROR_TOLERANCE) + + def test_non_conservative_zonal_mean(self): + """Tests if the correct number of queries are returned.""" + grid_path = self.gridfile_ne30 + data_path = self.datafile_vortex_ne30 + uxds = ux.open_dataset(grid_path, data_path) + + res = uxds['psi'].zonal_mean((-90.0, 90.0, 1)) + + assert len(res) == 181 + + def test_non_conservative_zonal_mean_at_pole(self): + """Tests the zonal average execution at both poles.""" + grid_path = self.gridfile_ne30 + data_path = self.datafile_vortex_ne30 + uxds = ux.open_dataset(grid_path, data_path) + + # Test at the poles + res_n90 = uxds['psi'].zonal_mean(90) + res_p90 = uxds['psi'].zonal_mean(-90) + + # result should be a scalar + assert len(res_n90.values) == 1 + assert len(res_p90.values) == 1 + + + def test_zonal_mean_dask(self): + """Tests if zonal average returns Dask arrays when appropriate.""" + grid_path = self.gridfile_ne30 + data_path = self.datafile_vortex_ne30 + uxds = ux.open_dataset(grid_path, data_path) + + uxds['psi'] = uxds['psi'].chunk() + + res = uxds['psi'].zonal_mean((-90, 90, 1)) + + assert isinstance(res.data, da.Array) + + res_computed = res.compute() + + assert isinstance(res_computed.data, np.ndarray) + + def test_zonal_weights(self): + grid_path = self.gridfile_ne30 + data_path = self.datafile_vortex_ne30 + uxds = ux.open_dataset(grid_path, data_path) + + za_1 = uxds['psi'].zonal_mean((-90, 90, 1), use_robust_weights=True) + za_2 = uxds['psi'].zonal_mean((-90, 90, 1), use_robust_weights=False) + + nt.assert_almost_equal(za_1.data, za_2.data) + + def test_lat_inputs(self): + grid_path = self.gridfile_ne30 + data_path = self.datafile_vortex_ne30 + uxds = ux.open_dataset(grid_path, data_path) + + assert len(uxds['psi'].zonal_mean(lat=1)) == 1 + assert len(uxds['psi'].zonal_mean(lat=(-90, 90, 1))) == 181 diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py index c2244f934..6d6364c01 100644 --- a/uxarray/core/dataarray.py +++ b/uxarray/core/dataarray.py @@ -1,5 +1,6 @@ from __future__ import annotations + import xarray as xr import numpy as np @@ -36,6 +37,7 @@ from uxarray.cross_sections import UxDataArrayCrossSectionAccessor from uxarray.core.aggregation import _uxda_grid_aggregate from uxarray.core.utils import _map_dims_to_ugrid +from uxarray.core.zonal import _compute_non_conservative_zonal_mean import warnings from warnings import warn @@ -376,6 +378,9 @@ def to_dataset( return uxds + def to_xarray(self): + return xr.DataArray(self) + def integrate( self, quadrature_rule: Optional[str] = "triangular", order: Optional[int] = 4 ) -> UxDataArray: @@ -430,6 +435,80 @@ def integrate( return uxda + def zonal_mean(self, lat=(-90, 90, 10), **kwargs): + """Compute averages along lines of constant latitude. + + Parameters + ---------- + lat : tuple, float, or array-like, default=(-90, 90, 10) + Latitude values in degrees. Can be specified as: + - tuple (start, end, step): Computes means at intervals of `step` in range [start, end] + - float: Computes mean for a single latitude + - array-like: Computes means for each specified latitude + + Returns + ------- + UxDataArray + Contains zonal means with a new 'latitudes' dimension and corresponding coordinates. + Name will be original_name + '_zonal_mean' or 'zonal_mean' if unnamed. + + Examples + -------- + # All latitudes from -90° to 90° at 10° intervals + >>> uxds["var"].zonal_mean() + + # Single latitude at 30° + >>> uxds["var"].zonal_mean(lat=30.0) + + # Range from -60° to 60° at 10° intervals + >>> uxds["var"].zonal_mean(lat=(-60, 60, 10)) + + Notes + ----- + Only supported for face-centered data variables. Candidate faces are determined + using spherical bounding boxes - faces whose bounds contain the target latitude + are included in calculations. + """ + if not self._face_centered(): + raise ValueError( + "Zonal mean computations are currently only supported for face-centered data variables." + ) + + if isinstance(lat, tuple): + # zonal mean over a range of latitudes + latitudes = np.arange(lat[0], lat[1] + lat[2], lat[2]) + latitudes = np.clip(latitudes, -90, 90) + elif isinstance(lat, (float, int)): + # zonal mean over a single latitude + latitudes = [lat] + elif isinstance(lat, (list, np.ndarray)): + # zonal mean over an array of arbitrary latitudes + latitudes = np.asarray(lat) + else: + raise ValueError( + "Invalid value for 'lat' provided. Must either be a single scalar value, tuple (min_lat, max_lat, step), or array-like." + ) + + res = _compute_non_conservative_zonal_mean( + uxda=self, latitudes=latitudes, **kwargs + ) + + dims = list(self.dims[:-1]) + ["latitudes"] + + uxda = UxDataArray( + res, + uxgrid=self.uxgrid, + dims=dims, + coords={"latitudes": latitudes}, + name=self.name + "_zonal_mean" if self.name is not None else "zonal_mean", + attrs={"zonal_mean": True}, + ) + + return uxda + + # Alias for 'zonal_mean', since this name is also commonly used. + zonal_average = zonal_mean + def weighted_mean(self, weights=None): """Computes a weighted mean. diff --git a/uxarray/core/zonal.py b/uxarray/core/zonal.py new file mode 100644 index 000000000..e57975b87 --- /dev/null +++ b/uxarray/core/zonal.py @@ -0,0 +1,60 @@ +import numpy as np +import dask.array as da + + +from uxarray.grid.integrate import _zonal_face_weights, _zonal_face_weights_robust +from uxarray.grid.utils import _get_cartesian_face_edge_nodes + + +def _compute_non_conservative_zonal_mean(uxda, latitudes, use_robust_weights=False): + """Computes the non-conservative zonal mean across one or more latitudes.""" + uxgrid = uxda.uxgrid + n_nodes_per_face = uxgrid.n_nodes_per_face.values + shape = uxda.shape[:-1] + (len(latitudes),) + if isinstance(uxda.data, da.Array): + # Create a Dask array for storing results + result = da.zeros(shape, dtype=uxda.dtype) + else: + # Create a NumPy array for storing results + result = np.zeros(shape, dtype=uxda.dtype) + + faces_edge_nodes_xyz = _get_cartesian_face_edge_nodes( + uxgrid.face_node_connectivity.values, + uxgrid.n_face, + uxgrid.n_max_face_nodes, + uxgrid.node_x.values, + uxgrid.node_y.values, + uxgrid.node_z.values, + ) + + bounds = uxgrid.bounds.values + + for i, lat in enumerate(latitudes): + face_indices = uxda.uxgrid.get_faces_at_constant_latitude(lat) + + z = np.sin(np.deg2rad(lat)) + + faces_edge_nodes_xyz_candidate = faces_edge_nodes_xyz[face_indices, :, :, :] + + n_nodes_per_face_candidate = n_nodes_per_face[face_indices] + + bounds_candidate = bounds[face_indices] + + if use_robust_weights: + weights = _zonal_face_weights_robust( + faces_edge_nodes_xyz_candidate, z, bounds_candidate + )["weight"].to_numpy() + else: + weights = _zonal_face_weights( + faces_edge_nodes_xyz_candidate, + bounds_candidate, + n_nodes_per_face_candidate, + z, + ) + + total_weight = weights.sum() + result[..., i] = ((uxda.data[..., face_indices] * weights) / total_weight).sum( + axis=-1 + ) + + return result diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index aa93192d2..105841339 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -86,6 +86,7 @@ faces_within_lat_bounds, ) + from spatialpandas import GeoDataFrame from uxarray.plot.accessor import GridPlotAccessor @@ -1411,9 +1412,7 @@ def bounds(self): "This initial execution will be significantly longer.", RuntimeWarning, ) - _populate_bounds(self) - return self._ds["bounds"] @bounds.setter @@ -1445,6 +1444,7 @@ def face_bounds_lon(self): @property def face_bounds_lat(self): """Latitude bounds for each face in degrees.""" + if "face_bounds_lat" not in self._ds: bounds = self.bounds.values bounds_lat = np.sort(np.rad2deg(bounds[:, 0, :]), axis=-1) diff --git a/uxarray/plot/accessor.py b/uxarray/plot/accessor.py index 9c0adc9b9..8280d1b38 100644 --- a/uxarray/plot/accessor.py +++ b/uxarray/plot/accessor.py @@ -6,6 +6,7 @@ import cartopy.crs as ccrs import hvplot.pandas +import hvplot.xarray import pandas as pd import uxarray.plot.dataarray_plot as dataarray_plot @@ -550,6 +551,18 @@ def rasterize( **kwargs, ) + def line(self, backend=None, *args, **kwargs): + """Wrapper for ``hvplot.line()``""" + uxarray.plot.utils.backend.assign(backend) + da = self._uxda.to_xarray() + return da.hvplot.line(*args, **kwargs) + + def scatter(self, backend=None, *args, **kwargs): + """Wrapper for ``hvplot.scatter()``""" + uxarray.plot.utils.backend.assign(backend) + da = self._uxda.to_xarray() + return da.hvplot.scatter(*args, **kwargs) + class UxDatasetPlotAccessor: """Plotting accessor for ``UxDataset``."""