From 42e3036f131a1ec6f5b287d7563321bd16bb4bde Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Fri, 4 Oct 2024 00:00:03 -0500 Subject: [PATCH 01/25] initial work on fast cross-sections --- uxarray/grid/grid.py | 45 +++++++++++++++ uxarray/grid/intersections.py | 105 +++++++++++++++++++++++++++++++++- 2 files changed, 149 insertions(+), 1 deletion(-) diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 0341426d0..0258c56a0 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -65,6 +65,11 @@ _populate_edge_node_distances, ) +from uxarray.grid.intersections import ( + constant_lat_intersections, + constant_lon_intersections, +) + from spatialpandas import GeoDataFrame from uxarray.plot.accessor import GridPlotAccessor @@ -1878,3 +1883,43 @@ def isel(self, **dim_kwargs): raise ValueError( "Indexing must be along a grid dimension: ('n_node', 'n_edge', 'n_face')" ) + + def get_edges_at_constant_latitude(self, lat): + """TODO:""" + lat = np.deg2rad(lat) + edge_node_connectivity = self.edge_node_connectivity.values + edge_node_x = self.node_x.values[edge_node_connectivity] + edge_node_y = self.node_y.values[edge_node_connectivity] + edge_node_z = self.node_z.values[edge_node_connectivity] + + edges, _ = constant_lat_intersections( + lat, edge_node_x, edge_node_y, edge_node_z, self.n_edge + ) + + return edges + + def get_edges_at_constant_longitude(self, lon): + """TODO:""" + lon = np.deg2rad(lon) + edge_node_connectivity = self.edge_node_connectivity.values + edge_node_x = self.node_x.values[edge_node_connectivity] + edge_node_y = self.node_y.values[edge_node_connectivity] + edge_node_z = self.node_z.values[edge_node_connectivity] + + edges, _ = constant_lon_intersections( + lon, edge_node_x, edge_node_y, edge_node_z, self.n_edge + ) + + return edges + + def get_faces_at_constant_latitude(self, lat): + """TODO:""" + lat = np.deg2rad(lat) + edges = self.get_edges_at_constant_latitude(lat) + return self.edge_face_connectivity[edges].data.ravel() + + def get_faces_at_constant_longitude(self, lon): + """TODO:""" + lon = np.deg2rad(lon) + edges = self.get_edges_at_constant_longitude(lon) + return self.edge_face_connectivity[edges].data.ravel() diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index 21daab6dd..b37124a0c 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -1,5 +1,5 @@ import numpy as np -from uxarray.constants import MACHINE_EPSILON, ERROR_TOLERANCE +from uxarray.constants import MACHINE_EPSILON, ERROR_TOLERANCE, INT_DTYPE from uxarray.grid.utils import _newton_raphson_solver_for_gca_constLat from uxarray.grid.arcs import point_within_gca, extreme_gca_latitude, in_between import platform @@ -7,6 +7,109 @@ from uxarray.utils.computing import cross_fma, allclose, dot, cross, norm +from numba import njit + + +@njit +def constant_lat_intersections(lat, edge_node_x, edge_node_y, edge_node_z, n_edge): + """Determine which edges intersect a constant line of latitude on a sphere. + + Parameters + ---------- + lat: + Constant latitude value in radians. + edge_node_x: + Array of shape (n_edge, 2) containing x-coordinates of the edge nodes. + edge_node_y: + Array of shape (n_edge, 2) containing y-coordinates of the edge nodes. + edge_node_z: + Array of shape (n_edge, 2) containing z-coordinates of the edge nodes. + n_edge: + Total number of edges to check. + + Returns + ------- + intersecting_edges: + List of indices of edges that intersect the constant latitude. + intersection_points: + Array of shape (n_intersections, 3) containing intersection points in Cartesian coordinates. + """ + intersecting_edges = [] + intersection_points = [] + + # Calculate the constant z-value for the given latitude + z_constant = np.sin(lat) + + # Iterate through each edge and check for intersections + for i in range(n_edge): + # Get the z-coordinates of the edge's nodes + z0 = edge_node_z[i, 0] + z1 = edge_node_z[i, 1] + + # Check if the edge crosses the constant latitude plane + if (z0 - z_constant) * (z1 - z_constant) < 0: + # Calculate the intersection point using linear interpolation + t = (z_constant - z0) / (z1 - z0) + x_intersect = edge_node_x[i, 0] + t * ( + edge_node_x[i, 1] - edge_node_x[i, 0] + ) + y_intersect = edge_node_y[i, 0] + t * ( + edge_node_y[i, 1] - edge_node_y[i, 0] + ) + z_intersect = z_constant + + # Append the intersecting edge index and intersection point + intersecting_edges.append(i) + intersection_points.append((x_intersect, y_intersect, z_intersect)) + + return np.array(intersecting_edges, dtype=INT_DTYPE), np.array(intersection_points) + + +@njit +def constant_lon_intersections(lon, edge_node_x, edge_node_y, edge_node_z, n_edge): + """TODO: Not quite ready, some weird behavior around poles.""" + intersecting_edges = [] + intersection_points = [] + + # Iterate through each edge and check for intersections + for i in range(n_edge): + # Get the Cartesian coordinates of the edge's nodes + x0, x1 = edge_node_x[i, 0], edge_node_x[i, 1] + y0, y1 = edge_node_y[i, 0], edge_node_y[i, 1] + z0, z1 = edge_node_z[i, 0], edge_node_z[i, 1] + + # Convert Cartesian coordinates to spherical coordinates (latitude and longitude) + lat0 = np.arcsin(z0) + lat1 = np.arcsin(z1) + lon0 = np.arctan2(y0, x0) + lon1 = np.arctan2(y1, x1) + + # Check if the edge spans the desired longitude + if (lon0 <= lon and lon1 >= lon) or (lon0 >= lon and lon1 <= lon): + # Calculate the intersection latitude using spherical trigonometry + d_lon = lon1 - lon0 + d_lat = lat1 - lat0 + if d_lon == 0: + continue + + # Calculate the intersection latitude using linear interpolation in the spherical space + t = (lon - lon0) / d_lon + lat_intersect = lat0 + t * d_lat + + # Convert intersection point back to Cartesian coordinates + x_intersect = np.cos(lat_intersect) * np.cos(lon) + y_intersect = np.cos(lat_intersect) * np.sin(lon) + z_intersect = np.sin(lat_intersect) + + # Check if the intersection latitude is within the bounds of the edge's latitudes + if min(lat0, lat1) <= lat_intersect <= max(lat0, lat1): + # Append the intersecting edge index and intersection point + intersecting_edges.append(i) + intersection_points.append((x_intersect, y_intersect, z_intersect)) + + return np.array(intersecting_edges, dtype=INT_DTYPE), np.array(intersection_points) + + def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=True): """Calculate the intersection point(s) of two Great Circle Arcs (GCAs) in a Cartesian coordinate system. From 5c8d9d5e9bd07eb82a6ee6725a24cf067cb45c7f Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Fri, 4 Oct 2024 11:22:05 -0500 Subject: [PATCH 02/25] add API --- uxarray/core/dataarray.py | 5 ++ uxarray/grid/grid.py | 57 ++++++++++++++--------- uxarray/grid/intersections.py | 88 ++++++++++++++++++----------------- 3 files changed, 84 insertions(+), 66 deletions(-) diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py index aa4e9a7b4..36f5ad55e 100644 --- a/uxarray/core/dataarray.py +++ b/uxarray/core/dataarray.py @@ -1121,3 +1121,8 @@ def _slice_from_grid(self, sliced_grid): dims=self.dims, attrs=self.attrs, ) + + def constant_latitude_cross_section(self, lat: float): + faces = self.uxgrid.get_faces_at_constant_latitude(lat) + + return self.isel(n_face=faces) diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 0258c56a0..8bdb1f333 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -67,7 +67,6 @@ from uxarray.grid.intersections import ( constant_lat_intersections, - constant_lon_intersections, ) from spatialpandas import GeoDataFrame @@ -91,6 +90,8 @@ import copy +from uxarray.constants import INT_FILL_VALUE + class Grid: """Represents a two-dimensional unstructured grid encoded following the @@ -1884,9 +1885,19 @@ def isel(self, **dim_kwargs): "Indexing must be along a grid dimension: ('n_node', 'n_edge', 'n_face')" ) + def constant_latitude_cross_section(self, lat: float, return_face_indices=False): + faces = self.get_faces_at_constant_latitude(lat) + + grid_at_constant_lat = self.isel(n_face=faces) + + if return_face_indices: + return grid_at_constant_lat, faces + else: + return grid_at_constant_lat + def get_edges_at_constant_latitude(self, lat): """TODO:""" - lat = np.deg2rad(lat) + edge_node_connectivity = self.edge_node_connectivity.values edge_node_x = self.node_x.values[edge_node_connectivity] edge_node_y = self.node_y.values[edge_node_connectivity] @@ -1895,31 +1906,31 @@ def get_edges_at_constant_latitude(self, lat): edges, _ = constant_lat_intersections( lat, edge_node_x, edge_node_y, edge_node_z, self.n_edge ) - return edges - def get_edges_at_constant_longitude(self, lon): - """TODO:""" - lon = np.deg2rad(lon) - edge_node_connectivity = self.edge_node_connectivity.values - edge_node_x = self.node_x.values[edge_node_connectivity] - edge_node_y = self.node_y.values[edge_node_connectivity] - edge_node_z = self.node_z.values[edge_node_connectivity] - - edges, _ = constant_lon_intersections( - lon, edge_node_x, edge_node_y, edge_node_z, self.n_edge - ) - - return edges + # def get_edges_at_constant_longitude(self, lon): + # """TODO:""" + # lon = np.deg2rad(lon) + # edge_node_connectivity = self.edge_node_connectivity.values + # edge_node_x = self.node_x.values[edge_node_connectivity] + # edge_node_y = self.node_y.values[edge_node_connectivity] + # edge_node_z = self.node_z.values[edge_node_connectivity] + # + # edges, _ = constant_lon_intersections( + # lon, edge_node_x, edge_node_y, edge_node_z, self.n_edge + # ) + # + # return edges def get_faces_at_constant_latitude(self, lat): """TODO:""" - lat = np.deg2rad(lat) edges = self.get_edges_at_constant_latitude(lat) - return self.edge_face_connectivity[edges].data.ravel() + faces = self.edge_face_connectivity[edges].data.ravel() - def get_faces_at_constant_longitude(self, lon): - """TODO:""" - lon = np.deg2rad(lon) - edges = self.get_edges_at_constant_longitude(lon) - return self.edge_face_connectivity[edges].data.ravel() + return faces[faces != INT_FILL_VALUE] + + # def get_faces_at_constant_longitude(self, lon): + # """TODO:""" + # lon = np.deg2rad(lon) + # edges = self.get_edges_at_constant_longitude(lon) + # return self.edge_face_connectivity[edges].data.ravel() diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index b37124a0c..76dfcb704 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -34,6 +34,8 @@ def constant_lat_intersections(lat, edge_node_x, edge_node_y, edge_node_z, n_edg intersection_points: Array of shape (n_intersections, 3) containing intersection points in Cartesian coordinates. """ + lat = np.deg2rad(lat) + intersecting_edges = [] intersection_points = [] @@ -65,49 +67,49 @@ def constant_lat_intersections(lat, edge_node_x, edge_node_y, edge_node_z, n_edg return np.array(intersecting_edges, dtype=INT_DTYPE), np.array(intersection_points) -@njit -def constant_lon_intersections(lon, edge_node_x, edge_node_y, edge_node_z, n_edge): - """TODO: Not quite ready, some weird behavior around poles.""" - intersecting_edges = [] - intersection_points = [] - - # Iterate through each edge and check for intersections - for i in range(n_edge): - # Get the Cartesian coordinates of the edge's nodes - x0, x1 = edge_node_x[i, 0], edge_node_x[i, 1] - y0, y1 = edge_node_y[i, 0], edge_node_y[i, 1] - z0, z1 = edge_node_z[i, 0], edge_node_z[i, 1] - - # Convert Cartesian coordinates to spherical coordinates (latitude and longitude) - lat0 = np.arcsin(z0) - lat1 = np.arcsin(z1) - lon0 = np.arctan2(y0, x0) - lon1 = np.arctan2(y1, x1) - - # Check if the edge spans the desired longitude - if (lon0 <= lon and lon1 >= lon) or (lon0 >= lon and lon1 <= lon): - # Calculate the intersection latitude using spherical trigonometry - d_lon = lon1 - lon0 - d_lat = lat1 - lat0 - if d_lon == 0: - continue - - # Calculate the intersection latitude using linear interpolation in the spherical space - t = (lon - lon0) / d_lon - lat_intersect = lat0 + t * d_lat - - # Convert intersection point back to Cartesian coordinates - x_intersect = np.cos(lat_intersect) * np.cos(lon) - y_intersect = np.cos(lat_intersect) * np.sin(lon) - z_intersect = np.sin(lat_intersect) - - # Check if the intersection latitude is within the bounds of the edge's latitudes - if min(lat0, lat1) <= lat_intersect <= max(lat0, lat1): - # Append the intersecting edge index and intersection point - intersecting_edges.append(i) - intersection_points.append((x_intersect, y_intersect, z_intersect)) - - return np.array(intersecting_edges, dtype=INT_DTYPE), np.array(intersection_points) +# @njit +# def constant_lon_intersections(lon, edge_node_x, edge_node_y, edge_node_z, n_edge): +# """TODO: Not quite ready, some weird behavior around poles.""" +# intersecting_edges = [] +# intersection_points = [] +# +# # Iterate through each edge and check for intersections +# for i in range(n_edge): +# # Get the Cartesian coordinates of the edge's nodes +# x0, x1 = edge_node_x[i, 0], edge_node_x[i, 1] +# y0, y1 = edge_node_y[i, 0], edge_node_y[i, 1] +# z0, z1 = edge_node_z[i, 0], edge_node_z[i, 1] +# +# # Convert Cartesian coordinates to spherical coordinates (latitude and longitude) +# lat0 = np.arcsin(z0) +# lat1 = np.arcsin(z1) +# lon0 = np.arctan2(y0, x0) +# lon1 = np.arctan2(y1, x1) +# +# # Check if the edge spans the desired longitude +# if (lon0 <= lon and lon1 >= lon) or (lon0 >= lon and lon1 <= lon): +# # Calculate the intersection latitude using spherical trigonometry +# d_lon = lon1 - lon0 +# d_lat = lat1 - lat0 +# if d_lon == 0: +# continue +# +# # Calculate the intersection latitude using linear interpolation in the spherical space +# t = (lon - lon0) / d_lon +# lat_intersect = lat0 + t * d_lat +# +# # Convert intersection point back to Cartesian coordinates +# x_intersect = np.cos(lat_intersect) * np.cos(lon) +# y_intersect = np.cos(lat_intersect) * np.sin(lon) +# z_intersect = np.sin(lat_intersect) +# +# # Check if the intersection latitude is within the bounds of the edge's latitudes +# if min(lat0, lat1) <= lat_intersect <= max(lat0, lat1): +# # Append the intersecting edge index and intersection point +# intersecting_edges.append(i) +# intersection_points.append((x_intersect, y_intersect, z_intersect)) +# +# return np.array(intersecting_edges, dtype=INT_DTYPE), np.array(intersection_points) def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=True): From aeb544b23c28e58bed7bb22a7fc977f942f69ab0 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Fri, 4 Oct 2024 11:31:08 -0500 Subject: [PATCH 03/25] add section for cross sections --- docs/user-guide/cross-sections.ipynb | 156 +++++++++++++++++++++++++++ docs/userguide.rst | 4 + 2 files changed, 160 insertions(+) create mode 100644 docs/user-guide/cross-sections.ipynb diff --git a/docs/user-guide/cross-sections.ipynb b/docs/user-guide/cross-sections.ipynb new file mode 100644 index 000000000..20cdcaa4a --- /dev/null +++ b/docs/user-guide/cross-sections.ipynb @@ -0,0 +1,156 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4a432a8bf95d9cdb", + "metadata": {}, + "source": [ + "# Cross-Sections" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b35ba4a2c30750e4", + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-04T14:54:06.168033Z", + "start_time": "2024-10-04T14:54:03.916424Z" + } + }, + "outputs": [], + "source": [ + "import uxarray as ux\n", + "import geoviews.feature as gf\n", + "import cartopy.crs as ccrs\n", + "import geoviews as gv\n", + "\n", + "projection = ccrs.Robinson()" + ] + }, + { + "cell_type": "markdown", + "id": "395a3db7-495c-4cff-b733-06bbe522a604", + "metadata": {}, + "source": [ + "## Data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b4160275c09fe6b0", + "metadata": {}, + "outputs": [], + "source": [ + "base_path = \"../../test/meshfiles/ugrid/outCSne30/\"\n", + "grid_path = base_path + \"outCSne30.ug\"\n", + "data_path = base_path + \"outCSne30_vortex.nc\"\n", + "\n", + "uxds = ux.open_dataset(grid_path, data_path)\n", + "\n", + "uxds" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "350df379-f496-46b7-ab86-5a2e630e3cea", + "metadata": {}, + "outputs": [], + "source": [ + "uxds.uxgrid.normalize_cartesian_coordinates()" + ] + }, + { + "cell_type": "markdown", + "id": "a7a40958-0a4d-47e4-9e38-31925261a892", + "metadata": {}, + "source": [ + "## Constant Latitude" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9facc02e-e4d4-4be3-9473-1a47e512c274", + "metadata": {}, + "outputs": [], + "source": [ + "lat = 1.0\n", + "\n", + "uxda_constant_lat = uxds[\"psi\"].constant_latitude_cross_section(lat)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bb5132fc-0d49-4484-95d4-8197c955a954", + "metadata": {}, + "outputs": [], + "source": [ + "constant_lat_line = gv.Path([(-180, lat), (180, lat)]).opts(\n", + " projection=projection, color=\"grey\", line_width=2, alpha=0.5\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "484b77a6-86da-4395-9e63-f5ac56e37deb", + "metadata": {}, + "outputs": [], + "source": [ + "(\n", + " uxda_constant_lat.plot(\n", + " rasterize=False,\n", + " backend=\"bokeh\",\n", + " cmap=\"inferno\",\n", + " projection=projection,\n", + " global_extent=True,\n", + " coastline=True,\n", + " )\n", + " * gf.grid(projection=projection)\n", + " * constant_lat_line\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "c4a7ee25-0b60-470f-bab7-92ff70563076", + "metadata": {}, + "source": [ + "## Constant Longitude" + ] + }, + { + "cell_type": "markdown", + "id": "54d9eff1-67f1-4691-a3b0-1ee0c874c98f", + "metadata": {}, + "source": [ + "## Arbitrary Great Circle Arc (GCA)" + ] + } + ], + "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.11.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/userguide.rst b/docs/userguide.rst index 98618cfe4..24e2ca61b 100644 --- a/docs/userguide.rst +++ b/docs/userguide.rst @@ -43,6 +43,9 @@ These user guides provide detailed explanations of the core functionality in UXa `Subsetting `_ Select specific regions of a grid +`Cross-Sections `_ + Select cross-sections of a grid + `Remapping `_ Remap (a.k.a Regrid) between unstructured grids @@ -79,6 +82,7 @@ These user guides provide additional detail about specific features in UXarray. user-guide/mpl.ipynb user-guide/advanced-plotting.ipynb user-guide/subset.ipynb + user-guide/cross-sections.ipynb user-guide/topological-aggregations.ipynb user-guide/area_calc.ipynb user-guide/holoviz.ipynb From a2b247f606a979fb0742310aedb4f933bda629d4 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Fri, 4 Oct 2024 17:33:01 -0500 Subject: [PATCH 04/25] parallel numba --- docs/user-guide/cross-sections.ipynb | 1863 +++++++++++++++++++++++++- uxarray/grid/grid.py | 50 +- uxarray/grid/intersections.py | 36 +- 3 files changed, 1888 insertions(+), 61 deletions(-) diff --git a/docs/user-guide/cross-sections.ipynb b/docs/user-guide/cross-sections.ipynb index 20cdcaa4a..e7e6f5ebc 100644 --- a/docs/user-guide/cross-sections.ipynb +++ b/docs/user-guide/cross-sections.ipynb @@ -10,7 +10,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "b35ba4a2c30750e4", "metadata": { "ExecuteTime": { @@ -18,7 +18,606 @@ "start_time": "2024-10-04T14:54:03.916424Z" } }, - "outputs": [], + "outputs": [ + { + "data": { + "application/javascript": [ + "(function(root) {\n", + " function now() {\n", + " return new Date();\n", + " }\n", + "\n", + " var force = true;\n", + " var py_version = '3.4.2'.replace('rc', '-rc.').replace('.dev', '-dev.');\n", + " var reloading = false;\n", + " var Bokeh = root.Bokeh;\n", + "\n", + " if (typeof (root._bokeh_timeout) === \"undefined\" || force) {\n", + " root._bokeh_timeout = Date.now() + 5000;\n", + " root._bokeh_failed_load = false;\n", + " }\n", + "\n", + " function run_callbacks() {\n", + " try {\n", + " root._bokeh_onload_callbacks.forEach(function(callback) {\n", + " if (callback != null)\n", + " callback();\n", + " });\n", + " } finally {\n", + " delete root._bokeh_onload_callbacks;\n", + " }\n", + " console.debug(\"Bokeh: all callbacks have finished\");\n", + " }\n", + "\n", + " function load_libs(css_urls, js_urls, js_modules, js_exports, callback) {\n", + " if (css_urls == null) css_urls = [];\n", + " if (js_urls == null) js_urls = [];\n", + " if (js_modules == null) js_modules = [];\n", + " if (js_exports == null) js_exports = {};\n", + "\n", + " root._bokeh_onload_callbacks.push(callback);\n", + "\n", + " if (root._bokeh_is_loading > 0) {\n", + " console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n", + " return null;\n", + " }\n", + " if (js_urls.length === 0 && js_modules.length === 0 && Object.keys(js_exports).length === 0) {\n", + " run_callbacks();\n", + " return null;\n", + " }\n", + " if (!reloading) {\n", + " console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n", + " }\n", + "\n", + " function on_load() {\n", + " root._bokeh_is_loading--;\n", + " if (root._bokeh_is_loading === 0) {\n", + " console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n", + " run_callbacks()\n", + " }\n", + " }\n", + " window._bokeh_on_load = on_load\n", + "\n", + " function on_error() {\n", + " console.error(\"failed to load \" + url);\n", + " }\n", + "\n", + " var skip = [];\n", + " if (window.requirejs) {\n", + " window.requirejs.config({'packages': {}, 'paths': {}, 'shim': {}});\n", + " root._bokeh_is_loading = css_urls.length + 0;\n", + " } else {\n", + " root._bokeh_is_loading = css_urls.length + js_urls.length + js_modules.length + Object.keys(js_exports).length;\n", + " }\n", + "\n", + " var existing_stylesheets = []\n", + " var links = document.getElementsByTagName('link')\n", + " for (var i = 0; i < links.length; i++) {\n", + " var link = links[i]\n", + " if (link.href != null) {\n", + "\texisting_stylesheets.push(link.href)\n", + " }\n", + " }\n", + " for (var i = 0; i < css_urls.length; i++) {\n", + " var url = css_urls[i];\n", + " if (existing_stylesheets.indexOf(url) !== -1) {\n", + "\ton_load()\n", + "\tcontinue;\n", + " }\n", + " const element = document.createElement(\"link\");\n", + " element.onload = on_load;\n", + " element.onerror = on_error;\n", + " element.rel = \"stylesheet\";\n", + " element.type = \"text/css\";\n", + " element.href = url;\n", + " console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n", + " document.body.appendChild(element);\n", + " } var existing_scripts = []\n", + " var scripts = document.getElementsByTagName('script')\n", + " for (var i = 0; i < scripts.length; i++) {\n", + " var script = scripts[i]\n", + " if (script.src != null) {\n", + "\texisting_scripts.push(script.src)\n", + " }\n", + " }\n", + " for (var i = 0; i < js_urls.length; i++) {\n", + " var url = js_urls[i];\n", + " if (skip.indexOf(url) !== -1 || existing_scripts.indexOf(url) !== -1) {\n", + "\tif (!window.requirejs) {\n", + "\t on_load();\n", + "\t}\n", + "\tcontinue;\n", + " }\n", + " var element = document.createElement('script');\n", + " element.onload = on_load;\n", + " element.onerror = on_error;\n", + " element.async = false;\n", + " element.src = url;\n", + " console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n", + " document.head.appendChild(element);\n", + " }\n", + " for (var i = 0; i < js_modules.length; i++) {\n", + " var url = js_modules[i];\n", + " if (skip.indexOf(url) !== -1 || existing_scripts.indexOf(url) !== -1) {\n", + "\tif (!window.requirejs) {\n", + "\t on_load();\n", + "\t}\n", + "\tcontinue;\n", + " }\n", + " var element = document.createElement('script');\n", + " element.onload = on_load;\n", + " element.onerror = on_error;\n", + " element.async = false;\n", + " element.src = url;\n", + " element.type = \"module\";\n", + " console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n", + " document.head.appendChild(element);\n", + " }\n", + " for (const name in js_exports) {\n", + " var url = js_exports[name];\n", + " if (skip.indexOf(url) >= 0 || root[name] != null) {\n", + "\tif (!window.requirejs) {\n", + "\t on_load();\n", + "\t}\n", + "\tcontinue;\n", + " }\n", + " var element = document.createElement('script');\n", + " element.onerror = on_error;\n", + " element.async = false;\n", + " element.type = \"module\";\n", + " console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n", + " element.textContent = `\n", + " import ${name} from \"${url}\"\n", + " window.${name} = ${name}\n", + " window._bokeh_on_load()\n", + " `\n", + " document.head.appendChild(element);\n", + " }\n", + " if (!js_urls.length && !js_modules.length) {\n", + " on_load()\n", + " }\n", + " };\n", + "\n", + " function inject_raw_css(css) {\n", + " const element = document.createElement(\"style\");\n", + " element.appendChild(document.createTextNode(css));\n", + " document.body.appendChild(element);\n", + " }\n", + "\n", + " var js_urls = [\"https://cdn.bokeh.org/bokeh/release/bokeh-3.4.2.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-gl-3.4.2.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-widgets-3.4.2.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-tables-3.4.2.min.js\", \"https://cdn.holoviz.org/panel/1.4.4/dist/panel.min.js\"];\n", + " var js_modules = [];\n", + " var js_exports = {};\n", + " var css_urls = [];\n", + " var inline_js = [ function(Bokeh) {\n", + " Bokeh.set_log_level(\"info\");\n", + " },\n", + "function(Bokeh) {} // ensure no trailing comma for IE\n", + " ];\n", + "\n", + " function run_inline_js() {\n", + " if ((root.Bokeh !== undefined) || (force === true)) {\n", + " for (var i = 0; i < inline_js.length; i++) {\n", + "\ttry {\n", + " inline_js[i].call(root, root.Bokeh);\n", + "\t} catch(e) {\n", + "\t if (!reloading) {\n", + "\t throw e;\n", + "\t }\n", + "\t}\n", + " }\n", + " // Cache old bokeh versions\n", + " if (Bokeh != undefined && !reloading) {\n", + "\tvar NewBokeh = root.Bokeh;\n", + "\tif (Bokeh.versions === undefined) {\n", + "\t Bokeh.versions = new Map();\n", + "\t}\n", + "\tif (NewBokeh.version !== Bokeh.version) {\n", + "\t Bokeh.versions.set(NewBokeh.version, NewBokeh)\n", + "\t}\n", + "\troot.Bokeh = Bokeh;\n", + " }} else if (Date.now() < root._bokeh_timeout) {\n", + " setTimeout(run_inline_js, 100);\n", + " } else if (!root._bokeh_failed_load) {\n", + " console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n", + " root._bokeh_failed_load = true;\n", + " }\n", + " root._bokeh_is_initializing = false\n", + " }\n", + "\n", + " function load_or_wait() {\n", + " // Implement a backoff loop that tries to ensure we do not load multiple\n", + " // versions of Bokeh and its dependencies at the same time.\n", + " // In recent versions we use the root._bokeh_is_initializing flag\n", + " // to determine whether there is an ongoing attempt to initialize\n", + " // bokeh, however for backward compatibility we also try to ensure\n", + " // that we do not start loading a newer (Panel>=1.0 and Bokeh>3) version\n", + " // before older versions are fully initialized.\n", + " if (root._bokeh_is_initializing && Date.now() > root._bokeh_timeout) {\n", + " root._bokeh_is_initializing = false;\n", + " root._bokeh_onload_callbacks = undefined;\n", + " console.log(\"Bokeh: BokehJS was loaded multiple times but one version failed to initialize.\");\n", + " load_or_wait();\n", + " } else if (root._bokeh_is_initializing || (typeof root._bokeh_is_initializing === \"undefined\" && root._bokeh_onload_callbacks !== undefined)) {\n", + " setTimeout(load_or_wait, 100);\n", + " } else {\n", + " root._bokeh_is_initializing = true\n", + " root._bokeh_onload_callbacks = []\n", + " var bokeh_loaded = Bokeh != null && (Bokeh.version === py_version || (Bokeh.versions !== undefined && Bokeh.versions.has(py_version)));\n", + " if (!reloading && !bokeh_loaded) {\n", + "\troot.Bokeh = undefined;\n", + " }\n", + " load_libs(css_urls, js_urls, js_modules, js_exports, function() {\n", + "\tconsole.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n", + "\trun_inline_js();\n", + " });\n", + " }\n", + " }\n", + " // Give older versions of the autoload script a head-start to ensure\n", + " // they initialize before we start loading newer version.\n", + " setTimeout(load_or_wait, 100)\n", + "}(window));" + ], + "application/vnd.holoviews_load.v0+json": "(function(root) {\n function now() {\n return new Date();\n }\n\n var force = true;\n var py_version = '3.4.2'.replace('rc', '-rc.').replace('.dev', '-dev.');\n var reloading = false;\n var Bokeh = root.Bokeh;\n\n if (typeof (root._bokeh_timeout) === \"undefined\" || force) {\n root._bokeh_timeout = Date.now() + 5000;\n root._bokeh_failed_load = false;\n }\n\n function run_callbacks() {\n try {\n root._bokeh_onload_callbacks.forEach(function(callback) {\n if (callback != null)\n callback();\n });\n } finally {\n delete root._bokeh_onload_callbacks;\n }\n console.debug(\"Bokeh: all callbacks have finished\");\n }\n\n function load_libs(css_urls, js_urls, js_modules, js_exports, callback) {\n if (css_urls == null) css_urls = [];\n if (js_urls == null) js_urls = [];\n if (js_modules == null) js_modules = [];\n if (js_exports == null) js_exports = {};\n\n root._bokeh_onload_callbacks.push(callback);\n\n if (root._bokeh_is_loading > 0) {\n console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n return null;\n }\n if (js_urls.length === 0 && js_modules.length === 0 && Object.keys(js_exports).length === 0) {\n run_callbacks();\n return null;\n }\n if (!reloading) {\n console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n }\n\n function on_load() {\n root._bokeh_is_loading--;\n if (root._bokeh_is_loading === 0) {\n console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n run_callbacks()\n }\n }\n window._bokeh_on_load = on_load\n\n function on_error() {\n console.error(\"failed to load \" + url);\n }\n\n var skip = [];\n if (window.requirejs) {\n window.requirejs.config({'packages': {}, 'paths': {}, 'shim': {}});\n root._bokeh_is_loading = css_urls.length + 0;\n } else {\n root._bokeh_is_loading = css_urls.length + js_urls.length + js_modules.length + Object.keys(js_exports).length;\n }\n\n var existing_stylesheets = []\n var links = document.getElementsByTagName('link')\n for (var i = 0; i < links.length; i++) {\n var link = links[i]\n if (link.href != null) {\n\texisting_stylesheets.push(link.href)\n }\n }\n for (var i = 0; i < css_urls.length; i++) {\n var url = css_urls[i];\n if (existing_stylesheets.indexOf(url) !== -1) {\n\ton_load()\n\tcontinue;\n }\n const element = document.createElement(\"link\");\n element.onload = on_load;\n element.onerror = on_error;\n element.rel = \"stylesheet\";\n element.type = \"text/css\";\n element.href = url;\n console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n document.body.appendChild(element);\n } var existing_scripts = []\n var scripts = document.getElementsByTagName('script')\n for (var i = 0; i < scripts.length; i++) {\n var script = scripts[i]\n if (script.src != null) {\n\texisting_scripts.push(script.src)\n }\n }\n for (var i = 0; i < js_urls.length; i++) {\n var url = js_urls[i];\n if (skip.indexOf(url) !== -1 || existing_scripts.indexOf(url) !== -1) {\n\tif (!window.requirejs) {\n\t on_load();\n\t}\n\tcontinue;\n }\n var element = document.createElement('script');\n element.onload = on_load;\n element.onerror = on_error;\n element.async = false;\n element.src = url;\n console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n document.head.appendChild(element);\n }\n for (var i = 0; i < js_modules.length; i++) {\n var url = js_modules[i];\n if (skip.indexOf(url) !== -1 || existing_scripts.indexOf(url) !== -1) {\n\tif (!window.requirejs) {\n\t on_load();\n\t}\n\tcontinue;\n }\n var element = document.createElement('script');\n element.onload = on_load;\n element.onerror = on_error;\n element.async = false;\n element.src = url;\n element.type = \"module\";\n console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n document.head.appendChild(element);\n }\n for (const name in js_exports) {\n var url = js_exports[name];\n if (skip.indexOf(url) >= 0 || root[name] != null) {\n\tif (!window.requirejs) {\n\t on_load();\n\t}\n\tcontinue;\n }\n var element = document.createElement('script');\n element.onerror = on_error;\n element.async = false;\n element.type = \"module\";\n console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n element.textContent = `\n import ${name} from \"${url}\"\n window.${name} = ${name}\n window._bokeh_on_load()\n `\n document.head.appendChild(element);\n }\n if (!js_urls.length && !js_modules.length) {\n on_load()\n }\n };\n\n function inject_raw_css(css) {\n const element = document.createElement(\"style\");\n element.appendChild(document.createTextNode(css));\n document.body.appendChild(element);\n }\n\n var js_urls = [\"https://cdn.bokeh.org/bokeh/release/bokeh-3.4.2.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-gl-3.4.2.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-widgets-3.4.2.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-tables-3.4.2.min.js\", \"https://cdn.holoviz.org/panel/1.4.4/dist/panel.min.js\"];\n var js_modules = [];\n var js_exports = {};\n var css_urls = [];\n var inline_js = [ function(Bokeh) {\n Bokeh.set_log_level(\"info\");\n },\nfunction(Bokeh) {} // ensure no trailing comma for IE\n ];\n\n function run_inline_js() {\n if ((root.Bokeh !== undefined) || (force === true)) {\n for (var i = 0; i < inline_js.length; i++) {\n\ttry {\n inline_js[i].call(root, root.Bokeh);\n\t} catch(e) {\n\t if (!reloading) {\n\t throw e;\n\t }\n\t}\n }\n // Cache old bokeh versions\n if (Bokeh != undefined && !reloading) {\n\tvar NewBokeh = root.Bokeh;\n\tif (Bokeh.versions === undefined) {\n\t Bokeh.versions = new Map();\n\t}\n\tif (NewBokeh.version !== Bokeh.version) {\n\t Bokeh.versions.set(NewBokeh.version, NewBokeh)\n\t}\n\troot.Bokeh = Bokeh;\n }} else if (Date.now() < root._bokeh_timeout) {\n setTimeout(run_inline_js, 100);\n } else if (!root._bokeh_failed_load) {\n console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n root._bokeh_failed_load = true;\n }\n root._bokeh_is_initializing = false\n }\n\n function load_or_wait() {\n // Implement a backoff loop that tries to ensure we do not load multiple\n // versions of Bokeh and its dependencies at the same time.\n // In recent versions we use the root._bokeh_is_initializing flag\n // to determine whether there is an ongoing attempt to initialize\n // bokeh, however for backward compatibility we also try to ensure\n // that we do not start loading a newer (Panel>=1.0 and Bokeh>3) version\n // before older versions are fully initialized.\n if (root._bokeh_is_initializing && Date.now() > root._bokeh_timeout) {\n root._bokeh_is_initializing = false;\n root._bokeh_onload_callbacks = undefined;\n console.log(\"Bokeh: BokehJS was loaded multiple times but one version failed to initialize.\");\n load_or_wait();\n } else if (root._bokeh_is_initializing || (typeof root._bokeh_is_initializing === \"undefined\" && root._bokeh_onload_callbacks !== undefined)) {\n setTimeout(load_or_wait, 100);\n } else {\n root._bokeh_is_initializing = true\n root._bokeh_onload_callbacks = []\n var bokeh_loaded = Bokeh != null && (Bokeh.version === py_version || (Bokeh.versions !== undefined && Bokeh.versions.has(py_version)));\n if (!reloading && !bokeh_loaded) {\n\troot.Bokeh = undefined;\n }\n load_libs(css_urls, js_urls, js_modules, js_exports, function() {\n\tconsole.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n\trun_inline_js();\n });\n }\n }\n // Give older versions of the autoload script a head-start to ensure\n // they initialize before we start loading newer version.\n setTimeout(load_or_wait, 100)\n}(window));" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/javascript": [ + "\n", + "if ((window.PyViz === undefined) || (window.PyViz instanceof HTMLElement)) {\n", + " window.PyViz = {comms: {}, comm_status:{}, kernels:{}, receivers: {}, plot_index: []}\n", + "}\n", + "\n", + "\n", + " function JupyterCommManager() {\n", + " }\n", + "\n", + " JupyterCommManager.prototype.register_target = function(plot_id, comm_id, msg_handler) {\n", + " if (window.comm_manager || ((window.Jupyter !== undefined) && (Jupyter.notebook.kernel != null))) {\n", + " var comm_manager = window.comm_manager || Jupyter.notebook.kernel.comm_manager;\n", + " comm_manager.register_target(comm_id, function(comm) {\n", + " comm.on_msg(msg_handler);\n", + " });\n", + " } else if ((plot_id in window.PyViz.kernels) && (window.PyViz.kernels[plot_id])) {\n", + " window.PyViz.kernels[plot_id].registerCommTarget(comm_id, function(comm) {\n", + " comm.onMsg = msg_handler;\n", + " });\n", + " } else if (typeof google != 'undefined' && google.colab.kernel != null) {\n", + " google.colab.kernel.comms.registerTarget(comm_id, (comm) => {\n", + " var messages = comm.messages[Symbol.asyncIterator]();\n", + " function processIteratorResult(result) {\n", + " var message = result.value;\n", + " console.log(message)\n", + " var content = {data: message.data, comm_id};\n", + " var buffers = []\n", + " for (var buffer of message.buffers || []) {\n", + " buffers.push(new DataView(buffer))\n", + " }\n", + " var metadata = message.metadata || {};\n", + " var msg = {content, buffers, metadata}\n", + " msg_handler(msg);\n", + " return messages.next().then(processIteratorResult);\n", + " }\n", + " return messages.next().then(processIteratorResult);\n", + " })\n", + " }\n", + " }\n", + "\n", + " JupyterCommManager.prototype.get_client_comm = function(plot_id, comm_id, msg_handler) {\n", + " if (comm_id in window.PyViz.comms) {\n", + " return window.PyViz.comms[comm_id];\n", + " } else if (window.comm_manager || ((window.Jupyter !== undefined) && (Jupyter.notebook.kernel != null))) {\n", + " var comm_manager = window.comm_manager || Jupyter.notebook.kernel.comm_manager;\n", + " var comm = comm_manager.new_comm(comm_id, {}, {}, {}, comm_id);\n", + " if (msg_handler) {\n", + " comm.on_msg(msg_handler);\n", + " }\n", + " } else if ((plot_id in window.PyViz.kernels) && (window.PyViz.kernels[plot_id])) {\n", + " var comm = window.PyViz.kernels[plot_id].connectToComm(comm_id);\n", + " comm.open();\n", + " if (msg_handler) {\n", + " comm.onMsg = msg_handler;\n", + " }\n", + " } else if (typeof google != 'undefined' && google.colab.kernel != null) {\n", + " var comm_promise = google.colab.kernel.comms.open(comm_id)\n", + " comm_promise.then((comm) => {\n", + " window.PyViz.comms[comm_id] = comm;\n", + " if (msg_handler) {\n", + " var messages = comm.messages[Symbol.asyncIterator]();\n", + " function processIteratorResult(result) {\n", + " var message = result.value;\n", + " var content = {data: message.data};\n", + " var metadata = message.metadata || {comm_id};\n", + " var msg = {content, metadata}\n", + " msg_handler(msg);\n", + " return messages.next().then(processIteratorResult);\n", + " }\n", + " return messages.next().then(processIteratorResult);\n", + " }\n", + " }) \n", + " var sendClosure = (data, metadata, buffers, disposeOnDone) => {\n", + " return comm_promise.then((comm) => {\n", + " comm.send(data, metadata, buffers, disposeOnDone);\n", + " });\n", + " };\n", + " var comm = {\n", + " send: sendClosure\n", + " };\n", + " }\n", + " window.PyViz.comms[comm_id] = comm;\n", + " return comm;\n", + " }\n", + " window.PyViz.comm_manager = new JupyterCommManager();\n", + " \n", + "\n", + "\n", + "var JS_MIME_TYPE = 'application/javascript';\n", + "var HTML_MIME_TYPE = 'text/html';\n", + "var EXEC_MIME_TYPE = 'application/vnd.holoviews_exec.v0+json';\n", + "var CLASS_NAME = 'output';\n", + "\n", + "/**\n", + " * Render data to the DOM node\n", + " */\n", + "function render(props, node) {\n", + " var div = document.createElement(\"div\");\n", + " var script = document.createElement(\"script\");\n", + " node.appendChild(div);\n", + " node.appendChild(script);\n", + "}\n", + "\n", + "/**\n", + " * Handle when a new output is added\n", + " */\n", + "function handle_add_output(event, handle) {\n", + " var output_area = handle.output_area;\n", + " var output = handle.output;\n", + " if ((output.data == undefined) || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n", + " return\n", + " }\n", + " var id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n", + " var toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n", + " if (id !== undefined) {\n", + " var nchildren = toinsert.length;\n", + " var html_node = toinsert[nchildren-1].children[0];\n", + " html_node.innerHTML = output.data[HTML_MIME_TYPE];\n", + " var scripts = [];\n", + " var nodelist = html_node.querySelectorAll(\"script\");\n", + " for (var i in nodelist) {\n", + " if (nodelist.hasOwnProperty(i)) {\n", + " scripts.push(nodelist[i])\n", + " }\n", + " }\n", + "\n", + " scripts.forEach( function (oldScript) {\n", + " var newScript = document.createElement(\"script\");\n", + " var attrs = [];\n", + " var nodemap = oldScript.attributes;\n", + " for (var j in nodemap) {\n", + " if (nodemap.hasOwnProperty(j)) {\n", + " attrs.push(nodemap[j])\n", + " }\n", + " }\n", + " attrs.forEach(function(attr) { newScript.setAttribute(attr.name, attr.value) });\n", + " newScript.appendChild(document.createTextNode(oldScript.innerHTML));\n", + " oldScript.parentNode.replaceChild(newScript, oldScript);\n", + " });\n", + " if (JS_MIME_TYPE in output.data) {\n", + " toinsert[nchildren-1].children[1].textContent = output.data[JS_MIME_TYPE];\n", + " }\n", + " output_area._hv_plot_id = id;\n", + " if ((window.Bokeh !== undefined) && (id in Bokeh.index)) {\n", + " window.PyViz.plot_index[id] = Bokeh.index[id];\n", + " } else {\n", + " window.PyViz.plot_index[id] = null;\n", + " }\n", + " } else if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n", + " var bk_div = document.createElement(\"div\");\n", + " bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n", + " var script_attrs = bk_div.children[0].attributes;\n", + " for (var i = 0; i < script_attrs.length; i++) {\n", + " toinsert[toinsert.length - 1].childNodes[1].setAttribute(script_attrs[i].name, script_attrs[i].value);\n", + " }\n", + " // store reference to server id on output_area\n", + " output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n", + " }\n", + "}\n", + "\n", + "/**\n", + " * Handle when an output is cleared or removed\n", + " */\n", + "function handle_clear_output(event, handle) {\n", + " var id = handle.cell.output_area._hv_plot_id;\n", + " var server_id = handle.cell.output_area._bokeh_server_id;\n", + " if (((id === undefined) || !(id in PyViz.plot_index)) && (server_id !== undefined)) { return; }\n", + " var comm = window.PyViz.comm_manager.get_client_comm(\"hv-extension-comm\", \"hv-extension-comm\", function () {});\n", + " if (server_id !== null) {\n", + " comm.send({event_type: 'server_delete', 'id': server_id});\n", + " return;\n", + " } else if (comm !== null) {\n", + " comm.send({event_type: 'delete', 'id': id});\n", + " }\n", + " delete PyViz.plot_index[id];\n", + " if ((window.Bokeh !== undefined) & (id in window.Bokeh.index)) {\n", + " var doc = window.Bokeh.index[id].model.document\n", + " doc.clear();\n", + " const i = window.Bokeh.documents.indexOf(doc);\n", + " if (i > -1) {\n", + " window.Bokeh.documents.splice(i, 1);\n", + " }\n", + " }\n", + "}\n", + "\n", + "/**\n", + " * Handle kernel restart event\n", + " */\n", + "function handle_kernel_cleanup(event, handle) {\n", + " delete PyViz.comms[\"hv-extension-comm\"];\n", + " window.PyViz.plot_index = {}\n", + "}\n", + "\n", + "/**\n", + " * Handle update_display_data messages\n", + " */\n", + "function handle_update_output(event, handle) {\n", + " handle_clear_output(event, {cell: {output_area: handle.output_area}})\n", + " handle_add_output(event, handle)\n", + "}\n", + "\n", + "function register_renderer(events, OutputArea) {\n", + " function append_mime(data, metadata, element) {\n", + " // create a DOM node to render to\n", + " var toinsert = this.create_output_subarea(\n", + " metadata,\n", + " CLASS_NAME,\n", + " EXEC_MIME_TYPE\n", + " );\n", + " this.keyboard_manager.register_events(toinsert);\n", + " // Render to node\n", + " var props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n", + " render(props, toinsert[0]);\n", + " element.append(toinsert);\n", + " return toinsert\n", + " }\n", + "\n", + " events.on('output_added.OutputArea', handle_add_output);\n", + " events.on('output_updated.OutputArea', handle_update_output);\n", + " events.on('clear_output.CodeCell', handle_clear_output);\n", + " events.on('delete.Cell', handle_clear_output);\n", + " events.on('kernel_ready.Kernel', handle_kernel_cleanup);\n", + "\n", + " OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n", + " safe: true,\n", + " index: 0\n", + " });\n", + "}\n", + "\n", + "if (window.Jupyter !== undefined) {\n", + " try {\n", + " var events = require('base/js/events');\n", + " var OutputArea = require('notebook/js/outputarea').OutputArea;\n", + " if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n", + " register_renderer(events, OutputArea);\n", + " }\n", + " } catch(err) {\n", + " }\n", + "}\n" + ], + "application/vnd.holoviews_load.v0+json": "\nif ((window.PyViz === undefined) || (window.PyViz instanceof HTMLElement)) {\n window.PyViz = {comms: {}, comm_status:{}, kernels:{}, receivers: {}, plot_index: []}\n}\n\n\n function JupyterCommManager() {\n }\n\n JupyterCommManager.prototype.register_target = function(plot_id, comm_id, msg_handler) {\n if (window.comm_manager || ((window.Jupyter !== undefined) && (Jupyter.notebook.kernel != null))) {\n var comm_manager = window.comm_manager || Jupyter.notebook.kernel.comm_manager;\n comm_manager.register_target(comm_id, function(comm) {\n comm.on_msg(msg_handler);\n });\n } else if ((plot_id in window.PyViz.kernels) && (window.PyViz.kernels[plot_id])) {\n window.PyViz.kernels[plot_id].registerCommTarget(comm_id, function(comm) {\n comm.onMsg = msg_handler;\n });\n } else if (typeof google != 'undefined' && google.colab.kernel != null) {\n google.colab.kernel.comms.registerTarget(comm_id, (comm) => {\n var messages = comm.messages[Symbol.asyncIterator]();\n function processIteratorResult(result) {\n var message = result.value;\n console.log(message)\n var content = {data: message.data, comm_id};\n var buffers = []\n for (var buffer of message.buffers || []) {\n buffers.push(new DataView(buffer))\n }\n var metadata = message.metadata || {};\n var msg = {content, buffers, metadata}\n msg_handler(msg);\n return messages.next().then(processIteratorResult);\n }\n return messages.next().then(processIteratorResult);\n })\n }\n }\n\n JupyterCommManager.prototype.get_client_comm = function(plot_id, comm_id, msg_handler) {\n if (comm_id in window.PyViz.comms) {\n return window.PyViz.comms[comm_id];\n } else if (window.comm_manager || ((window.Jupyter !== undefined) && (Jupyter.notebook.kernel != null))) {\n var comm_manager = window.comm_manager || Jupyter.notebook.kernel.comm_manager;\n var comm = comm_manager.new_comm(comm_id, {}, {}, {}, comm_id);\n if (msg_handler) {\n comm.on_msg(msg_handler);\n }\n } else if ((plot_id in window.PyViz.kernels) && (window.PyViz.kernels[plot_id])) {\n var comm = window.PyViz.kernels[plot_id].connectToComm(comm_id);\n comm.open();\n if (msg_handler) {\n comm.onMsg = msg_handler;\n }\n } else if (typeof google != 'undefined' && google.colab.kernel != null) {\n var comm_promise = google.colab.kernel.comms.open(comm_id)\n comm_promise.then((comm) => {\n window.PyViz.comms[comm_id] = comm;\n if (msg_handler) {\n var messages = comm.messages[Symbol.asyncIterator]();\n function processIteratorResult(result) {\n var message = result.value;\n var content = {data: message.data};\n var metadata = message.metadata || {comm_id};\n var msg = {content, metadata}\n msg_handler(msg);\n return messages.next().then(processIteratorResult);\n }\n return messages.next().then(processIteratorResult);\n }\n }) \n var sendClosure = (data, metadata, buffers, disposeOnDone) => {\n return comm_promise.then((comm) => {\n comm.send(data, metadata, buffers, disposeOnDone);\n });\n };\n var comm = {\n send: sendClosure\n };\n }\n window.PyViz.comms[comm_id] = comm;\n return comm;\n }\n window.PyViz.comm_manager = new JupyterCommManager();\n \n\n\nvar JS_MIME_TYPE = 'application/javascript';\nvar HTML_MIME_TYPE = 'text/html';\nvar EXEC_MIME_TYPE = 'application/vnd.holoviews_exec.v0+json';\nvar CLASS_NAME = 'output';\n\n/**\n * Render data to the DOM node\n */\nfunction render(props, node) {\n var div = document.createElement(\"div\");\n var script = document.createElement(\"script\");\n node.appendChild(div);\n node.appendChild(script);\n}\n\n/**\n * Handle when a new output is added\n */\nfunction handle_add_output(event, handle) {\n var output_area = handle.output_area;\n var output = handle.output;\n if ((output.data == undefined) || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n return\n }\n var id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n var toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n if (id !== undefined) {\n var nchildren = toinsert.length;\n var html_node = toinsert[nchildren-1].children[0];\n html_node.innerHTML = output.data[HTML_MIME_TYPE];\n var scripts = [];\n var nodelist = html_node.querySelectorAll(\"script\");\n for (var i in nodelist) {\n if (nodelist.hasOwnProperty(i)) {\n scripts.push(nodelist[i])\n }\n }\n\n scripts.forEach( function (oldScript) {\n var newScript = document.createElement(\"script\");\n var attrs = [];\n var nodemap = oldScript.attributes;\n for (var j in nodemap) {\n if (nodemap.hasOwnProperty(j)) {\n attrs.push(nodemap[j])\n }\n }\n attrs.forEach(function(attr) { newScript.setAttribute(attr.name, attr.value) });\n newScript.appendChild(document.createTextNode(oldScript.innerHTML));\n oldScript.parentNode.replaceChild(newScript, oldScript);\n });\n if (JS_MIME_TYPE in output.data) {\n toinsert[nchildren-1].children[1].textContent = output.data[JS_MIME_TYPE];\n }\n output_area._hv_plot_id = id;\n if ((window.Bokeh !== undefined) && (id in Bokeh.index)) {\n window.PyViz.plot_index[id] = Bokeh.index[id];\n } else {\n window.PyViz.plot_index[id] = null;\n }\n } else if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n var bk_div = document.createElement(\"div\");\n bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n var script_attrs = bk_div.children[0].attributes;\n for (var i = 0; i < script_attrs.length; i++) {\n toinsert[toinsert.length - 1].childNodes[1].setAttribute(script_attrs[i].name, script_attrs[i].value);\n }\n // store reference to server id on output_area\n output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n }\n}\n\n/**\n * Handle when an output is cleared or removed\n */\nfunction handle_clear_output(event, handle) {\n var id = handle.cell.output_area._hv_plot_id;\n var server_id = handle.cell.output_area._bokeh_server_id;\n if (((id === undefined) || !(id in PyViz.plot_index)) && (server_id !== undefined)) { return; }\n var comm = window.PyViz.comm_manager.get_client_comm(\"hv-extension-comm\", \"hv-extension-comm\", function () {});\n if (server_id !== null) {\n comm.send({event_type: 'server_delete', 'id': server_id});\n return;\n } else if (comm !== null) {\n comm.send({event_type: 'delete', 'id': id});\n }\n delete PyViz.plot_index[id];\n if ((window.Bokeh !== undefined) & (id in window.Bokeh.index)) {\n var doc = window.Bokeh.index[id].model.document\n doc.clear();\n const i = window.Bokeh.documents.indexOf(doc);\n if (i > -1) {\n window.Bokeh.documents.splice(i, 1);\n }\n }\n}\n\n/**\n * Handle kernel restart event\n */\nfunction handle_kernel_cleanup(event, handle) {\n delete PyViz.comms[\"hv-extension-comm\"];\n window.PyViz.plot_index = {}\n}\n\n/**\n * Handle update_display_data messages\n */\nfunction handle_update_output(event, handle) {\n handle_clear_output(event, {cell: {output_area: handle.output_area}})\n handle_add_output(event, handle)\n}\n\nfunction register_renderer(events, OutputArea) {\n function append_mime(data, metadata, element) {\n // create a DOM node to render to\n var toinsert = this.create_output_subarea(\n metadata,\n CLASS_NAME,\n EXEC_MIME_TYPE\n );\n this.keyboard_manager.register_events(toinsert);\n // Render to node\n var props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n render(props, toinsert[0]);\n element.append(toinsert);\n return toinsert\n }\n\n events.on('output_added.OutputArea', handle_add_output);\n events.on('output_updated.OutputArea', handle_update_output);\n events.on('clear_output.CodeCell', handle_clear_output);\n events.on('delete.Cell', handle_clear_output);\n events.on('kernel_ready.Kernel', handle_kernel_cleanup);\n\n OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n safe: true,\n index: 0\n });\n}\n\nif (window.Jupyter !== undefined) {\n try {\n var events = require('base/js/events');\n var OutputArea = require('notebook/js/outputarea').OutputArea;\n if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n register_renderer(events, OutputArea);\n }\n } catch(err) {\n }\n}\n" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.holoviews_exec.v0+json": "", + "text/html": [ + "
\n", + "
\n", + "
\n", + "" + ] + }, + "metadata": { + "application/vnd.holoviews_exec.v0+json": { + "id": "p1002" + } + }, + "output_type": "display_data" + } + ], "source": [ "import uxarray as ux\n", "import geoviews.feature as gf\n", @@ -38,10 +637,1147 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "id": "b4160275c09fe6b0", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.UxDataset> Size: 43kB\n",
+       "Dimensions:  (n_face: 5400)\n",
+       "Dimensions without coordinates: n_face\n",
+       "Data variables:\n",
+       "    psi      (n_face) float64 43kB ...
" + ], + "text/plain": [ + " Size: 43kB\n", + "Dimensions: (n_face: 5400)\n", + "Dimensions without coordinates: n_face\n", + "Data variables:\n", + " psi (n_face) float64 43kB ..." + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "base_path = \"../../test/meshfiles/ugrid/outCSne30/\"\n", "grid_path = base_path + \"outCSne30.ug\"\n", @@ -54,7 +1790,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "id": "350df379-f496-46b7-ab86-5a2e630e3cea", "metadata": {}, "outputs": [], @@ -72,11 +1808,21 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "id": "9facc02e-e4d4-4be3-9473-1a47e512c274", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 5.89 ms, sys: 2.29 ms, total: 8.18 ms\n", + "Wall time: 6.02 ms\n" + ] + } + ], "source": [ + "%%time\n", "lat = 1.0\n", "\n", "uxda_constant_lat = uxds[\"psi\"].constant_latitude_cross_section(lat)" @@ -84,7 +1830,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 8, "id": "bb5132fc-0d49-4484-95d4-8197c955a954", "metadata": {}, "outputs": [], @@ -96,10 +1842,107 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "id": "484b77a6-86da-4395-9e63-f5ac56e37deb", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": {}, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.holoviews_exec.v0+json": "", + "text/html": [ + "
\n", + "
\n", + "
\n", + "" + ], + "text/plain": [ + ":Overlay\n", + " .Polygons.I :Polygons [Longitude,Latitude] (psi)\n", + " .Coastline.I :Feature [Longitude,Latitude]\n", + " .Grid.I :Feature [Longitude,Latitude]\n", + " .Path.I :Path [Longitude,Latitude]" + ] + }, + "execution_count": 9, + "metadata": { + "application/vnd.holoviews_exec.v0+json": { + "id": "p1971" + } + }, + "output_type": "execute_result" + } + ], "source": [ "(\n", " uxda_constant_lat.plot(\n", diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 8bdb1f333..c3ade1f63 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -950,7 +950,7 @@ def face_node_connectivity(self, value): def edge_node_connectivity(self) -> xr.DataArray: """Indices of the two nodes that make up each edge. - Dimensions: ``(n_edge, n_max_edge_nodes)`` + Dimensions: ``(n_edge, two)`` Nodes are in arbitrary order. """ @@ -965,6 +965,23 @@ def edge_node_connectivity(self, value): assert isinstance(value, xr.DataArray) self._ds["edge_node_connectivity"] = value + @property + def edge_node_z(self) -> xr.DataArray: + """Cartesian z location for the two nodes that make up every edge. + + Dimensions: ``(n_edge, two)`` + """ + + if "edge_node_z" not in self._ds: + _edge_node_z = self.node_z.values[self.edge_node_connectivity.values] + + self._ds["edge_node_z"] = xr.DataArray( + data=_edge_node_z, + dims=["n_edge", "two"], + ) + + return self._ds["edge_node_z"] + @property def node_node_connectivity(self) -> xr.DataArray: """Indices of the nodes that surround each node.""" @@ -1898,29 +1915,8 @@ def constant_latitude_cross_section(self, lat: float, return_face_indices=False) def get_edges_at_constant_latitude(self, lat): """TODO:""" - edge_node_connectivity = self.edge_node_connectivity.values - edge_node_x = self.node_x.values[edge_node_connectivity] - edge_node_y = self.node_y.values[edge_node_connectivity] - edge_node_z = self.node_z.values[edge_node_connectivity] - - edges, _ = constant_lat_intersections( - lat, edge_node_x, edge_node_y, edge_node_z, self.n_edge - ) - return edges - - # def get_edges_at_constant_longitude(self, lon): - # """TODO:""" - # lon = np.deg2rad(lon) - # edge_node_connectivity = self.edge_node_connectivity.values - # edge_node_x = self.node_x.values[edge_node_connectivity] - # edge_node_y = self.node_y.values[edge_node_connectivity] - # edge_node_z = self.node_z.values[edge_node_connectivity] - # - # edges, _ = constant_lon_intersections( - # lon, edge_node_x, edge_node_y, edge_node_z, self.n_edge - # ) - # - # return edges + edges = constant_lat_intersections(lat, self.edge_node_z.values, self.n_edge) + return edges.squeeze() def get_faces_at_constant_latitude(self, lat): """TODO:""" @@ -1928,9 +1924,3 @@ def get_faces_at_constant_latitude(self, lat): faces = self.edge_face_connectivity[edges].data.ravel() return faces[faces != INT_FILL_VALUE] - - # def get_faces_at_constant_longitude(self, lon): - # """TODO:""" - # lon = np.deg2rad(lon) - # edges = self.get_edges_at_constant_longitude(lon) - # return self.edge_face_connectivity[edges].data.ravel() diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index 76dfcb704..7ca8a07d9 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -1,5 +1,5 @@ import numpy as np -from uxarray.constants import MACHINE_EPSILON, ERROR_TOLERANCE, INT_DTYPE +from uxarray.constants import MACHINE_EPSILON, ERROR_TOLERANCE from uxarray.grid.utils import _newton_raphson_solver_for_gca_constLat from uxarray.grid.arcs import point_within_gca, extreme_gca_latitude, in_between import platform @@ -7,11 +7,11 @@ from uxarray.utils.computing import cross_fma, allclose, dot, cross, norm -from numba import njit +from numba import njit, prange -@njit -def constant_lat_intersections(lat, edge_node_x, edge_node_y, edge_node_z, n_edge): +@njit(parallel=True, nogil=True) +def constant_lat_intersections(lat, edge_node_z, n_edge): """Determine which edges intersect a constant line of latitude on a sphere. Parameters @@ -36,35 +36,29 @@ def constant_lat_intersections(lat, edge_node_x, edge_node_y, edge_node_z, n_edg """ lat = np.deg2rad(lat) - intersecting_edges = [] - intersection_points = [] + # intersecting_edges = [] + intersecting_edges_mask = np.zeros(n_edge, dtype=np.int32) # Calculate the constant z-value for the given latitude z_constant = np.sin(lat) # Iterate through each edge and check for intersections - for i in range(n_edge): + for i in prange(n_edge): # Get the z-coordinates of the edge's nodes z0 = edge_node_z[i, 0] z1 = edge_node_z[i, 1] # Check if the edge crosses the constant latitude plane - if (z0 - z_constant) * (z1 - z_constant) < 0: - # Calculate the intersection point using linear interpolation - t = (z_constant - z0) / (z1 - z0) - x_intersect = edge_node_x[i, 0] + t * ( - edge_node_x[i, 1] - edge_node_x[i, 0] - ) - y_intersect = edge_node_y[i, 0] + t * ( - edge_node_y[i, 1] - edge_node_y[i, 0] - ) - z_intersect = z_constant + # if (z0 - z_constant) * (z1 - z_constant) < 0: + # intersecting_edges.append(i) + if (z0 - z_constant) * (z1 - z_constant) < 0.0: + intersecting_edges_mask[i] = 1 + + intersecting_edges = np.argwhere(intersecting_edges_mask) - # Append the intersecting edge index and intersection point - intersecting_edges.append(i) - intersection_points.append((x_intersect, y_intersect, z_intersect)) + return intersecting_edges - return np.array(intersecting_edges, dtype=INT_DTYPE), np.array(intersection_points) + # return np.array(intersecting_edges, dtype=INT_DTYPE) # @njit From 82e49eacbd97c2f185ac6a8a6d2134dc2009be51 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Mon, 7 Oct 2024 11:07:00 -0500 Subject: [PATCH 05/25] update notebook --- docs/user-guide/cross-sections.ipynb | 1310 +++----------------------- uxarray/grid/grid.py | 15 +- uxarray/grid/intersections.py | 22 +- 3 files changed, 153 insertions(+), 1194 deletions(-) diff --git a/docs/user-guide/cross-sections.ipynb b/docs/user-guide/cross-sections.ipynb index e7e6f5ebc..ea258f5f7 100644 --- a/docs/user-guide/cross-sections.ipynb +++ b/docs/user-guide/cross-sections.ipynb @@ -14,8 +14,8 @@ "id": "b35ba4a2c30750e4", "metadata": { "ExecuteTime": { - "end_time": "2024-10-04T14:54:06.168033Z", - "start_time": "2024-10-04T14:54:03.916424Z" + "end_time": "2024-10-05T15:19:46.812562Z", + "start_time": "2024-10-05T15:19:41.546324Z" } }, "outputs": [ @@ -539,11 +539,11 @@ "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", - "
\n", + "
\n", "
\n", "" ], "text/plain": [ - " Size: 43kB\n", - "Dimensions: (n_face: 5400)\n", - "Dimensions without coordinates: n_face\n", - "Data variables:\n", - " psi (n_face) float64 43kB ..." + ":DynamicMap []\n", + " :Image [Longitude,Latitude] (Longitude_Latitude psi)" ] }, "execution_count": 2, - "metadata": {}, + "metadata": { + "application/vnd.holoviews_exec.v0+json": { + "id": "p1011" + } + }, "output_type": "execute_result" } ], @@ -1784,18 +748,7 @@ "data_path = base_path + \"outCSne30_vortex.nc\"\n", "\n", "uxds = ux.open_dataset(grid_path, data_path)\n", - "\n", - "uxds" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "350df379-f496-46b7-ab86-5a2e630e3cea", - "metadata": {}, - "outputs": [], - "source": [ - "uxds.uxgrid.normalize_cartesian_coordinates()" + "uxds[\"psi\"].plot(cmap=\"inferno\", periodic_elements=\"split\", projection=projection)" ] }, { @@ -1808,41 +761,19 @@ }, { "cell_type": "code", - "execution_count": 7, - "id": "9facc02e-e4d4-4be3-9473-1a47e512c274", + "execution_count": 14, + "id": "3775daa1-2f1d-4738-bab5-2b69ebd689d9", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "CPU times: user 5.89 ms, sys: 2.29 ms, total: 8.18 ms\n", - "Wall time: 6.02 ms\n" - ] - } - ], + "outputs": [], "source": [ - "%%time\n", - "lat = 1.0\n", + "lat = 30\n", "\n", "uxda_constant_lat = uxds[\"psi\"].constant_latitude_cross_section(lat)" ] }, { "cell_type": "code", - "execution_count": 8, - "id": "bb5132fc-0d49-4484-95d4-8197c955a954", - "metadata": {}, - "outputs": [], - "source": [ - "constant_lat_line = gv.Path([(-180, lat), (180, lat)]).opts(\n", - " projection=projection, color=\"grey\", line_width=2, alpha=0.5\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": 9, + "execution_count": 17, "id": "484b77a6-86da-4395-9e63-f5ac56e37deb", "metadata": {}, "outputs": [ @@ -1855,12 +786,12 @@ "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ - "
\n", - "
\n", + "
\n", + "
\n", "
\n", "" - ] - }, - "metadata": { - "application/vnd.holoviews_exec.v0+json": { - "id": "p1002" - } - }, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "import uxarray as ux\n", "import geoviews.feature as gf\n", @@ -625,7 +28,6 @@ "import cartopy.crs as ccrs\n", "import geoviews as gv\n", "\n", - "\n", "projection = ccrs.Robinson()" ] }, @@ -639,7 +41,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "id": "b4160275c09fe6b0", "metadata": { "ExecuteTime": { @@ -647,108 +49,19 @@ "start_time": "2024-10-05T15:21:17.200324Z" } }, - "outputs": [ - { - "data": {}, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "application/vnd.holoviews_exec.v0+json": "", - "text/html": [ - "
\n", - "
\n", - "
\n", - "" - ], - "text/plain": [ - ":DynamicMap []\n", - " :Image [Longitude,Latitude] (Longitude_Latitude psi)" - ] - }, - "execution_count": 2, - "metadata": { - "application/vnd.holoviews_exec.v0+json": { - "id": "p1011" - } - }, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "base_path = \"../../test/meshfiles/ugrid/outCSne30/\"\n", "grid_path = base_path + \"outCSne30.ug\"\n", "data_path = base_path + \"outCSne30_vortex.nc\"\n", "\n", "uxds = ux.open_dataset(grid_path, data_path)\n", - "uxds[\"psi\"].plot(cmap=\"inferno\", periodic_elements=\"split\", projection=projection)" + "uxds[\"psi\"].plot(\n", + " cmap=\"inferno\",\n", + " periodic_elements=\"split\",\n", + " projection=projection,\n", + " title=\"Global Plot\",\n", + ")" ] }, { @@ -756,12 +69,15 @@ "id": "a7a40958-0a4d-47e4-9e38-31925261a892", "metadata": {}, "source": [ - "## Constant Latitude" + "## Constant Latitude\n", + "\n", + "Cross sections at constant latitudes can be performed using the ``constant_latitude_cross_section`` methods on either a ``ux.Grid`` or ``ux.DataArray``\n", + "\n" ] }, { "cell_type": "code", - "execution_count": 14, + "execution_count": null, "id": "3775daa1-2f1d-4738-bab5-2b69ebd689d9", "metadata": {}, "outputs": [], @@ -773,106 +89,10 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": null, "id": "484b77a6-86da-4395-9e63-f5ac56e37deb", "metadata": {}, - "outputs": [ - { - "data": {}, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "application/vnd.holoviews_exec.v0+json": "", - "text/html": [ - "
\n", - "
\n", - "
\n", - "" - ], - "text/plain": [ - ":Overlay\n", - " .Polygons.I :Polygons [Longitude,Latitude] (psi)\n", - " .Coastline.I :Feature [Longitude,Latitude]\n", - " .Grid.I :Feature [Longitude,Latitude]" - ] - }, - "execution_count": 17, - "metadata": { - "application/vnd.holoviews_exec.v0+json": { - "id": "p4894" - } - }, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "(\n", " uxda_constant_lat.plot(\n", @@ -890,19 +110,10 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": null, "id": "1cbee722-34a4-4e67-8e22-f393d7d36c99", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Global Mean: 1.0000000015332655\n", - "Mean at 30 degrees lat: 0.9797216779775056\n" - ] - } - ], + "outputs": [], "source": [ "print(f\"Global Mean: {uxds['psi'].data.mean()}\")\n", "print(f\"Mean at {lat} degrees lat: {uxda_constant_lat.data.mean()}\")" @@ -916,6 +127,16 @@ "## Constant Longitude" ] }, + { + "cell_type": "markdown", + "id": "9fcc8ec5-c6a8-4bde-a33d-7f37f9116ee2", + "metadata": {}, + "source": [ + "```{warning}\n", + "Constant longitude cross sections are not yet supported.\n", + "```" + ] + }, { "cell_type": "markdown", "id": "54d9eff1-67f1-4691-a3b0-1ee0c874c98f", @@ -923,6 +144,16 @@ "source": [ "## Arbitrary Great Circle Arc (GCA)" ] + }, + { + "cell_type": "markdown", + "id": "ea94ff9f-fe86-470d-813b-45f32a633ffc", + "metadata": {}, + "source": [ + "```{warning}\n", + "Arbitrary great circle arc cross sections are not yet supported.\n", + "```" + ] } ], "metadata": { From b974436aed251d92b01c513e3091920b8be911d2 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Mon, 7 Oct 2024 12:40:22 -0500 Subject: [PATCH 16/25] add docstring for uxda --- uxarray/core/dataarray.py | 32 ++++++++++++++++++++++++++++++-- 1 file changed, 30 insertions(+), 2 deletions(-) diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py index 36f5ad55e..e1d075b01 100644 --- a/uxarray/core/dataarray.py +++ b/uxarray/core/dataarray.py @@ -1122,7 +1122,35 @@ def _slice_from_grid(self, sliced_grid): attrs=self.attrs, ) - def constant_latitude_cross_section(self, lat: float): - faces = self.uxgrid.get_faces_at_constant_latitude(lat) + def constant_latitude_cross_section(self, lat: float, method="fast"): + """Extracts a cross-section of the data array at a specified constant + latitude. + + Parameters + ---------- + lat : float + The latitude at which to extract the cross-section, in degrees. + method : str, optional + The internal method to use when identifying faces at the constant latitude. + Options are: + - 'fast': Uses a faster but potentially less accurate method for face identification. + - 'accurate': Uses a slower but more accurate method. + Default is 'fast'. + + Raises + ------ + ValueError + If no intersections are found at the specified latitude, a ValueError is raised. + + Examples + -------- + >>> uxda.constant_latitude_cross_section(lat=-15.5) + + Notes + ----- + The accuracy and performance of the function can be controlled using the `method` parameter. + For higher precision requreiments, consider using method='acurate'. + """ + faces = self.uxgrid.get_faces_at_constant_latitude(lat, method) return self.isel(n_face=faces) From c1d05382b90a1927ea962380ca4d873b1598b7ae Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Mon, 7 Oct 2024 12:41:27 -0500 Subject: [PATCH 17/25] update docstring of intersections --- uxarray/grid/intersections.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index ddb3d3919..a846c04c0 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -18,10 +18,6 @@ def fast_constant_lat_intersections(lat, edge_node_z, n_edge): ---------- lat: Constant latitude value in radians. - edge_node_x: - Array of shape (n_edge, 2) containing x-coordinates of the edge nodes. - edge_node_y: - Array of shape (n_edge, 2) containing y-coordinates of the edge nodes. edge_node_z: Array of shape (n_edge, 2) containing z-coordinates of the edge nodes. n_edge: @@ -30,7 +26,7 @@ def fast_constant_lat_intersections(lat, edge_node_z, n_edge): Returns ------- intersecting_edges: - List of indices of edges that intersect the constant latitude. + array of indices of edges that intersect the constant latitude. """ lat = np.deg2rad(lat) From d2ab10f51a1046c808433b675c0f15f9b4cf5e04 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Mon, 7 Oct 2024 12:43:28 -0500 Subject: [PATCH 18/25] fix benchmark --- benchmarks/mpas_ocean.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py index d78a5a815..f70adb11a 100644 --- a/benchmarks/mpas_ocean.py +++ b/benchmarks/mpas_ocean.py @@ -159,4 +159,4 @@ class CrossSections(DatasetBenchmark): 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_lazitude_cross_section(lat, method='fast') + self.uxds.uxgrid.constant_latitude_cross_section(lat, method='fast') From 32dde5fb6a78723a45bf280a6bb3b654c0a81b50 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Mon, 7 Oct 2024 15:00:04 -0500 Subject: [PATCH 19/25] add tests for north and south pole --- test/test_cross_sections.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/test/test_cross_sections.py b/test/test_cross_sections.py index 83761e031..d6625b37d 100644 --- a/test/test_cross_sections.py +++ b/test/test_cross_sections.py @@ -10,6 +10,8 @@ quad_hex_grid_path = current_path / 'meshfiles' / "ugrid" / "quad-hexagon" / 'grid.nc' quad_hex_data_path = current_path / 'meshfiles' / "ugrid" / "quad-hexagon" / 'data.nc' +cube_sphere_grid = current_path / "meshfiles" / "geos-cs" / "c12" / "test-c12.native.nc4" + class TestQuadHex: @@ -68,3 +70,25 @@ def test_constant_lat_cross_section_uxds(self): with pytest.raises(ValueError): # no intersections found at this line uxds['t2m'].constant_latitude_cross_section(lat=10.0) + + +class TestGeosCubeSphere: + def test_north_pole(self): + uxgrid = ux.open_grid(cube_sphere_grid) + + lats = [89.85, 89.9, 89.95, 89.99] + + for lat in lats: + cross_grid = uxgrid.constant_latitude_cross_section(lat=lat) + # Cube sphere grid should have 4 faces centered around the pole + assert cross_grid.n_face == 4 + + def test_south_pole(self): + uxgrid = ux.open_grid(cube_sphere_grid) + + lats = [-89.85, -89.9, -89.95, -89.99] + + for lat in lats: + cross_grid = uxgrid.constant_latitude_cross_section(lat=lat) + # Cube sphere grid should have 4 faces centered around the pole + assert cross_grid.n_face == 4 From d7c9714cf816615d9af76b30931397f30342a791 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Mon, 7 Oct 2024 15:49:22 -0500 Subject: [PATCH 20/25] Update intersections.py --- uxarray/grid/intersections.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index a846c04c0..5159b6167 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -17,7 +17,7 @@ def fast_constant_lat_intersections(lat, edge_node_z, n_edge): Parameters ---------- lat: - Constant latitude value in radians. + Constant latitude value in degrees. edge_node_z: Array of shape (n_edge, 2) containing z-coordinates of the edge nodes. n_edge: From 07db497ace0f3e76901390c3d82ccd40cf5c5226 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Mon, 7 Oct 2024 20:53:58 -0500 Subject: [PATCH 21/25] update user guide --- docs/user-guide/cross-sections.ipynb | 33 +++++++++++++++++++++------- 1 file changed, 25 insertions(+), 8 deletions(-) diff --git a/docs/user-guide/cross-sections.ipynb b/docs/user-guide/cross-sections.ipynb index b3886efe9..7d3468608 100644 --- a/docs/user-guide/cross-sections.ipynb +++ b/docs/user-guide/cross-sections.ipynb @@ -7,20 +7,18 @@ "source": [ "# Cross-Sections\n", "\n", - "This section showcases how to obtain cross-sections of an unstructured grid using UXarray.\n" + "This section demonstrates how to extract cross-sections from an unstructured grid using UXarray, which allows the analysis and visualization across slices of grids.\n" ] }, { "cell_type": "code", - "execution_count": null, "id": "b35ba4a2c30750e4", "metadata": { "ExecuteTime": { - "end_time": "2024-10-05T15:19:46.812562Z", - "start_time": "2024-10-05T15:19:41.546324Z" + "end_time": "2024-10-08T01:49:10.528008Z", + "start_time": "2024-10-08T01:49:10.522526Z" } }, - "outputs": [], "source": [ "import uxarray as ux\n", "import geoviews.feature as gf\n", @@ -29,7 +27,9 @@ "import geoviews as gv\n", "\n", "projection = ccrs.Robinson()" - ] + ], + "outputs": [], + "execution_count": 9 }, { "cell_type": "markdown", @@ -71,10 +71,15 @@ "source": [ "## Constant Latitude\n", "\n", - "Cross sections at constant latitudes can be performed using the ``constant_latitude_cross_section`` methods on either a ``ux.Grid`` or ``ux.DataArray``\n", - "\n" + "Cross-sections along constant latitude lines can be obtained using the ``constant_latitude_cross_section`` method, available for both ``ux.Grid`` and ``ux.DataArray`` objects. This functionality allows users to extract and analyze slices of data at specified latitudes, providing insights into variations along horizontal sections of the grid.\n" ] }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "For example, we can obtain a cross-section at 30 degrees latitude by doing the following:", + "id": "2fbe9f6e5bb59a17" + }, { "cell_type": "code", "execution_count": null, @@ -87,6 +92,12 @@ "uxda_constant_lat = uxds[\"psi\"].constant_latitude_cross_section(lat)" ] }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "Since the result is a new ``UxDataArray``, we can directly plot the result to see the cross-section.", + "id": "dcec0b96b92e7f4" + }, { "cell_type": "code", "execution_count": null, @@ -108,6 +119,12 @@ ")" ] }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "You can also perform operations on the cross-section, such as taking the mean.", + "id": "c7cca7de4722c121" + }, { "cell_type": "code", "execution_count": null, From 61f6383df316c1c0fe29efab9e71b4289f8c3b77 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Wed, 9 Oct 2024 11:34:13 -0500 Subject: [PATCH 22/25] Delete docs/user-guide/grid.nc --- docs/user-guide/grid.nc | Bin 26574 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 docs/user-guide/grid.nc diff --git a/docs/user-guide/grid.nc b/docs/user-guide/grid.nc deleted file mode 100644 index 9e208fcdbc2a0276789da865d8a060a12e619977..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 26574 zcmeG^31C#!wP!L3>##&6$f7=$T9#z8Kq7=>2xMUk1PMqh_@6H`FEb;PdE?AW2(~O% ztcZ2_Q2RgqXt7_uOKbgV7cC|qk+$x2DT-LBbwjlJ+q!(U6?)FS@661bBqUK3>7NUk z^WJ^u-h1x3=bn4+x#ztr=T_AW8ZdG|M!$Z*>8O>S@XMS6>%lFD%C4)MTeT!(&=Zb~ zu_Dts#Co=@MqV7U;GjA3en18|=xW$#$VA{g2QWvZ?H1f)L1+<(@?po{J|h77cFg9S z?jGlKUAoV>uAg!B0`?iz4dQ1081GtrSh4Z%x ztp0>GNMa4-IvA_5WK(Sbpq0JQP-1l@7!F?Ui?A!yR}5lO+=m?b)9V2MA*}z;Fap z54DYg$NJ#aBmfjz82f`^W1xiCmj`eLg7zEUp9mkh@TwZ11VL!Zi+NB?tg8c{314^S zvMSg}@_>~hDRAJt4GSQLGG+pdL9k--s#+LE-5(1;RULA@u@JUUD@fQ`2wL91@LX6< zbj%09`3w7EtCm9!aUD5IlJm~o!BM!0FtEr-1eGuEh{Arnfbjs-rP>vr#Nck?SQUU9 zLGDv~zX>~tLC7&;?IX9#M}xo+dH={l%COts<4uFA_=A*+_y-l=Br-fCQ5V&(>M^ zJ3X!hAdc9i02-_w(*W4I4+rVVAO{eS;4NYRjym<;HaIpxwe=ltC8CgYq05s;2kMAt z$oFijZ7ba@05}D~+uv(^g!29$r`b#AEyDSqLz0gyLP3J}o_O*4QZd3f?Mg{DF^w8~ z{SUt`FVnyOL>#d6Kk00q|Dj8MvJY%4%9g2VomkHQgBE=W4+DR^o4Cr?8$=@ir^sGMP3|FeDv zS5OAt;W7d<_I{7evwQyx50k>f zIaq|?+TZVc7nTtQ5`k+I{AP8{$KWN7VaAyV>hHVn6Pw!W`#&p=YE$o|^`)mt%>Vwj zWu=rBEKH}C4w4t+UeC!<*~~PXRcABLY%{=SAK6SHn+0c68a7kT zCfGLf%cdl3rkibc*sMI;%(5v2+eEg_NSn&DnY%XiVACcxBiN=^Y{s5Vo!IPCn}KJW z{L@V09}HiQ z7+GRh>>Mx@QJwI;*`gTyP=fegISY9Al?q*?kg;SRG(gZQtE`<>UoQ+>YB~e|s^%@M zUQ|D?ZjoAFIjgoBT<4Qf>O_lL6j(HCVYR@ou<(}#f-wdEYf4x*E30Z0U$99J$AU&! zX=pZ;uwf=1wP?|x))3O&ij@jWfC$RK`Um4i00~kT>897$IJGL+6!gJRRK=KyNfE?} zOW60GZ(TW=^TK7m0=V9NIgon!1-tmxZet@60|kR}T5(VxBga7P}>XL2$TR?m0{um#OHN?Rm|^Dy%y z;h|F0J08m2&J|}pva2Hx#U>W{z_6!0e4th`@1*5nFZH3+c~?&GmPj6c*!`NKbh4%l zpw-4XN!GS|eZNzPi%!Q@G?OV+h~I=K_87Bz7NI)PgXJPG7Ax~g!-)EVVa?QIc9o7x zpHxTeHj^#6D&(3o5G$c-YJ#EAd0MDhhwM(srvTQJ)&t znslMzXlp>THDrYSL9^MXyU?HI#D4?iT1WW_Y;up;uoa_GY0_f>T;qIt%!Si86U$l) zYjZeg#zZT$Ribh}-5=HUn5t_r(}i-Hg~&s2#m5dN*4)~8^-G1nira^mV-`=4HCW_uMutx`qfC(@TxxDYcxgF#$ZUt481YfD%>vA9ddN|qgo{34u%^I zcSMVtZYS0t>JSC+jDegAY~cUQ3G2sgXlFQBF9=3SPhmJSb}#F`Al45+c@ zXrtz}CK>4>(ga_xu6JT&rtm*ot@$o#!iIY4f&TusZ5}5E)1}bcEuj=-;=~A-WI|1^ zij&?9dbN=1k6Jp*Mrn=>UvD*0>#15WYU)i9p*-j$h58i=hIK9KgkiP<=o`n{)nFLe zDe(!H{ub3a(I1^m47%}yExPJo?NenjQGYW2={B|zqE`=v)EJT<4EuY`L4XZ2Qwi;c z4X4;BZPLx;F4-8CZNm(PBh987i?qdUo!SD&JH`%)|lLyA>NL--?aehLGW1)z$hq(i-5_LU&1V@syGxcgL0TbX$+G zHw}7JH5ye>aZkH6+D4P!8VLoxK{I*nPGLF{2KM<+2N19R`+*Of7*3ijiZRoOs3yi6 zh#F{(Brz;1C@jeKDh zMs*6FC&?;$suxhIQ7>fM$RRliXd$C}mQw$m*xTgJglSTQ^cFqTefb5Qm5&*CzSgty zG!h-CZ}j)D^x2~HFgx{){$96Fu_9>1T6($32%CZKyH2M-G=|9n>c}?wrT%KsLfwk!hI1)O=QJBpAw2PHj_IO@qH|2vvido=yxAuc z&Lty($!DW;7A~!>!WmVNS~t}ji|F1!bCa~8sk`(?D+UcuVt_8es(Y^WD@M;ro&Q)X zR2+KKr1exVy$Q|eP~vEoyVVrYV)4^kuMyP~D0K8m)CKWeR@KZVBNhqhQQe7e28kMD zMsw6s>Rj8x#j}t}F8RFDz-bmXF|0`%cYZ3hIqd0Wv^;x5SZ;cC1lUY`+G zHRK=BAW&HfgxxV%vQIu%qna<+9P2RyON~AA6!oGX!$tR^aMb zFoA2`2}DsBm(s$h_o#r9SVKv?++nZMXc!0vG)02Z=CG0%Rtl#n?)((tQ-v7E1s1ob zC}FB)-qH4$MyntYXu2p?RALqkh+|UJftY+Y(J64_B(*OrqCqci>AO`Qv@cIJw$NVF zDP%%5*F@|Q6Bx_6|y>exH;F{w0feTJ= zd%XD?&l`8n@jUqVD?JAf6nob8ef$2>?7a5D&Y7NyKMXJbeN+FP4@eVosMB3v(*Hb;{X)Kk2}Y#0({S5XI!4|xjfg; z^m4yb`Ns67w$Jr({ahd8^F2Rjx|nY6C%4c1=5b;^r}D{O|A*_Dj;|`<{I{E z$0K1hBDEaXO=5~JKgEmtD?Lhgy53$*3+4?&S%Op)eM&^ZTlvI1S3V(+JtcRmTw!!~?03U{jG3d!H;>#H$U{q%IZIC5u=u`cj+u==t``-E^x5I^anFIP|BYy!tM3@NP{Rdxw z8RBw&#qi(5PX%sZ^D|UyIvr0Ai1Fv`8|r97PR-0i0$jU;z$RnnSkt-kVYsS{`_lEnWmaP!#COhuxb(&ry8 zHFLRa3Du-X+iw2kQR`BG8)mo%i&IekJ{sk2e`TASf~REa6IO|t{92bIlk36s5|SZeeIxZPHprU81EJzd7Z|WIfT3C{WSM&Ze&; zlXZsmNHf_Ac4V^t=ttyHscDt6D~?RoCDtR&oL9|{Sg-UW+NqScjvuj(u^wsWHP&M3 z$Ygz!BT5M}@xU;wdog?=Xs(Het7Wn#$`M3$K#UJmChMp{6oaACd0EKtubF`_$WwA; zveH7fQ7z)=MKEuzE>U{>+i|O4y#Xiu{#0AjQTg0gG?EHyS&O+&;EF~Xw`~36j$?m_@n5awVdw>jPKnh#Lzyc-0B2nv3t+K$ zNm(_0h;gmPWuCwPivQ$mg3evnx}=!aW;n)EP6oc&%woz*RO-FcMfNNOAIlY@0ZGNgsnRYqI3c+g+cYwOA6clI_nF$2>E z*A%JUtmAQSqZ&0r^5~UaZ=J7ZE?oAyRB}aFHPO0Ii^LR+^NFregyfor@QUElFJEED zJIz#rredv=F}DkkH0h;wmkhW;?;)P<2E#tR)ixr>=n}YA_V!ixWVRAlXejwrcN+8L zd(!r+Pg!(tXSN)#aIWt#1;;D&od%l{L$*g`E&v-Le9yK9-`_OJ`QLITc#>uUPbRjX zQtx{=6aH2zTib0l+Jklev(E$?8Y24(nF*bnYdu=-lV$?dqTo!x&o!XjN5kkJ6vo?^ zPxi^6O_qHZ-rqWA_s=?$ec^?pf%x2IpDlQ7Oe8kj`JT;ZzQ1_IEq|4=zcxc|WoW*$ z9Xm6w<);pWsYh{wb=3AB#y?Q93%>K5kSoAIsrUcga9h(P)-olAJHJ)=nsDy*gK84lN%AG7p+{gS{DmsT3=;bnhJWQW`6b~hyu9Kn?I7%4k!69vZZROvB@Q{^3Yg*sQrx;K|tsn;N zJ@)PRYb%I-E;dW~H@^7t256$I+46fY!u@-vl|kElasp>Y_CJr8>GbJ1AjaSNPP8?; zwn4BBvuyxv!)qH!+n_=6fZ4>#CNMS$%+71r|LMAuF4ap8h}}N#)*Pr=44TF8wfp8@ zE?8^+ z_!fkW_du@zIZBStm-k37j@b6wnaA9MZKoMZEgmoT0}nC7-fqECk0b8B^|k5$0HVuw AbpQYW From 5dbb557593adcab2b2a0490f3e1c5e87736d7adb Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Wed, 9 Oct 2024 12:49:58 -0500 Subject: [PATCH 23/25] use accessor --- test/test_cross_sections.py | 20 ++--- uxarray/core/dataarray.py | 35 +------- uxarray/cross_sections/__init__.py | 7 ++ uxarray/cross_sections/dataarray_accessor.py | 70 +++++++++++++++ uxarray/cross_sections/grid_accessor.py | 95 ++++++++++++++++++++ uxarray/grid/grid.py | 66 ++------------ 6 files changed, 189 insertions(+), 104 deletions(-) create mode 100644 uxarray/cross_sections/__init__.py create mode 100644 uxarray/cross_sections/dataarray_accessor.py create mode 100644 uxarray/cross_sections/grid_accessor.py diff --git a/test/test_cross_sections.py b/test/test_cross_sections.py index d6625b37d..26bf8777c 100644 --- a/test/test_cross_sections.py +++ b/test/test_cross_sections.py @@ -35,41 +35,41 @@ class TestQuadHex: def test_constant_lat_cross_section_grid(self): uxgrid = ux.open_grid(quad_hex_grid_path) - grid_top_two = uxgrid.constant_latitude_cross_section(lat=0.1) + grid_top_two = uxgrid.cross_section.constant_latitude(lat=0.1) assert grid_top_two.n_face == 2 - grid_bottom_two = uxgrid.constant_latitude_cross_section(lat=-0.1) + grid_bottom_two = uxgrid.cross_section.constant_latitude(lat=-0.1) assert grid_bottom_two.n_face == 2 - grid_all_four = uxgrid.constant_latitude_cross_section(lat=0.0) + grid_all_four = uxgrid.cross_section.constant_latitude(lat=0.0) assert grid_all_four.n_face == 4 with pytest.raises(ValueError): # no intersections found at this line - uxgrid.constant_latitude_cross_section(lat=10.0) + uxgrid.cross_section.constant_latitude(lat=10.0) def test_constant_lat_cross_section_uxds(self): uxds = ux.open_dataset(quad_hex_grid_path, quad_hex_data_path) - da_top_two = uxds['t2m'].constant_latitude_cross_section(lat=0.1) + da_top_two = uxds['t2m'].cross_section.constant_latitude(lat=0.1) nt.assert_array_equal(da_top_two.data, uxds['t2m'].isel(n_face=[1, 2]).data) - da_bottom_two = uxds['t2m'].constant_latitude_cross_section(lat=-0.1) + da_bottom_two = uxds['t2m'].cross_section.constant_latitude(lat=-0.1) nt.assert_array_equal(da_bottom_two.data, uxds['t2m'].isel(n_face=[0, 3]).data) - da_all_four = uxds['t2m'].constant_latitude_cross_section(lat=0.0) + da_all_four = uxds['t2m'].cross_section.constant_latitude(lat=0.0) nt.assert_array_equal(da_all_four.data , uxds['t2m'].data) with pytest.raises(ValueError): # no intersections found at this line - uxds['t2m'].constant_latitude_cross_section(lat=10.0) + uxds['t2m'].cross_section.constant_latitude(lat=10.0) class TestGeosCubeSphere: @@ -79,7 +79,7 @@ def test_north_pole(self): lats = [89.85, 89.9, 89.95, 89.99] for lat in lats: - cross_grid = uxgrid.constant_latitude_cross_section(lat=lat) + cross_grid = uxgrid.cross_section.constant_latitude(lat=lat) # Cube sphere grid should have 4 faces centered around the pole assert cross_grid.n_face == 4 @@ -89,6 +89,6 @@ def test_south_pole(self): lats = [-89.85, -89.9, -89.95, -89.99] for lat in lats: - cross_grid = uxgrid.constant_latitude_cross_section(lat=lat) + cross_grid = uxgrid.cross_section.constant_latitude(lat=lat) # Cube sphere grid should have 4 faces centered around the pole assert cross_grid.n_face == 4 diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py index e1d075b01..beb839856 100644 --- a/uxarray/core/dataarray.py +++ b/uxarray/core/dataarray.py @@ -30,6 +30,7 @@ from uxarray.plot.accessor import UxDataArrayPlotAccessor from uxarray.subset import DataArraySubsetAccessor from uxarray.remap import UxDataArrayRemapAccessor +from uxarray.cross_sections import UxDataArrayCrossSectionAccessor from uxarray.core.aggregation import _uxda_grid_aggregate import warnings @@ -82,6 +83,7 @@ def __init__(self, *args, uxgrid: Grid = None, **kwargs): plot = UncachedAccessor(UxDataArrayPlotAccessor) subset = UncachedAccessor(DataArraySubsetAccessor) remap = UncachedAccessor(UxDataArrayRemapAccessor) + cross_section = UncachedAccessor(UxDataArrayCrossSectionAccessor) def _repr_html_(self) -> str: if OPTIONS["display_style"] == "text": @@ -1121,36 +1123,3 @@ def _slice_from_grid(self, sliced_grid): dims=self.dims, attrs=self.attrs, ) - - def constant_latitude_cross_section(self, lat: float, method="fast"): - """Extracts a cross-section of the data array at a specified constant - latitude. - - Parameters - ---------- - lat : float - The latitude at which to extract the cross-section, in degrees. - method : str, optional - The internal method to use when identifying faces at the constant latitude. - Options are: - - 'fast': Uses a faster but potentially less accurate method for face identification. - - 'accurate': Uses a slower but more accurate method. - Default is 'fast'. - - Raises - ------ - ValueError - If no intersections are found at the specified latitude, a ValueError is raised. - - Examples - -------- - >>> uxda.constant_latitude_cross_section(lat=-15.5) - - Notes - ----- - The accuracy and performance of the function can be controlled using the `method` parameter. - For higher precision requreiments, consider using method='acurate'. - """ - faces = self.uxgrid.get_faces_at_constant_latitude(lat, method) - - return self.isel(n_face=faces) diff --git a/uxarray/cross_sections/__init__.py b/uxarray/cross_sections/__init__.py new file mode 100644 index 000000000..8694cabca --- /dev/null +++ b/uxarray/cross_sections/__init__.py @@ -0,0 +1,7 @@ +from .dataarray_accessor import UxDataArrayCrossSectionAccessor +from .grid_accessor import GridCrossSectionAccessor + +__all__ = ( + "GridCrossSectionAccessor", + "UxDataArrayCrossSectionAccessor", +) diff --git a/uxarray/cross_sections/dataarray_accessor.py b/uxarray/cross_sections/dataarray_accessor.py new file mode 100644 index 000000000..d840892b6 --- /dev/null +++ b/uxarray/cross_sections/dataarray_accessor.py @@ -0,0 +1,70 @@ +from __future__ import annotations + + +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + pass + + +class UxDataArrayCrossSectionAccessor: + """Accessor for cross-section operations on a ``UxDataArray``""" + + def __init__(self, uxda) -> None: + self.uxda = uxda + + def __repr__(self): + prefix = "\n" + methods_heading = "Supported Methods:\n" + + methods_heading += " * constant_latitude(center_coord, k, element, **kwargs)\n" + + return prefix + methods_heading + + def constant_latitude(self, lat: float, method="fast"): + """Extracts a cross-section of the data array at a specified constant + latitude. + + Parameters + ---------- + lat : float + The latitude at which to extract the cross-section, in degrees. + method : str, optional + The internal method to use when identifying faces at the constant latitude. + Options are: + - 'fast': Uses a faster but potentially less accurate method for face identification. + - 'accurate': Uses a slower but more accurate method. + Default is 'fast'. + + Raises + ------ + ValueError + If no intersections are found at the specified latitude, a ValueError is raised. + + Examples + -------- + >>> uxda.constant_latitude_cross_section(lat=-15.5) + + Notes + ----- + The accuracy and performance of the function can be controlled using the `method` parameter. + For higher precision requreiments, consider using method='acurate'. + """ + faces = self.uxda.uxgrid.get_faces_at_constant_latitude(lat, method) + + return self.uxda.isel(n_face=faces) + + def constant_longitude(self, *args, **kwargs): + raise NotImplementedError + + def gca(self, *args, **kwargs): + raise NotImplementedError + + def bounded_latitude(self, *args, **kwargs): + raise NotImplementedError + + def bounded_longitude(self, *args, **kwargs): + raise NotImplementedError + + def gca_gca(self, *args, **kwargs): + raise NotImplementedError diff --git a/uxarray/cross_sections/grid_accessor.py b/uxarray/cross_sections/grid_accessor.py new file mode 100644 index 000000000..067e8f5fb --- /dev/null +++ b/uxarray/cross_sections/grid_accessor.py @@ -0,0 +1,95 @@ +from __future__ import annotations + + +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + from uxarray.grid import Grid + + +class GridCrossSectionAccessor: + """Accessor for cross-section operations on a ``Grid``""" + + def __init__(self, uxgrid: Grid) -> None: + self.uxgrid = uxgrid + + def __repr__(self): + prefix = "\n" + methods_heading = "Supported Methods:\n" + + methods_heading += " * constant_latitude(lat, )\n" + return prefix + methods_heading + + def constant_latitude(self, lat: float, return_face_indices=False, method="fast"): + """Extracts a cross-section of the grid at a specified constant + latitude. + + This method identifies and returns all faces (or grid elements) that intersect + with a given latitude. The returned cross-section can include either just the grid + or both the grid elements and the corresponding face indices, depending + on the `return_face_indices` parameter. + + Parameters + ---------- + lat : float + The latitude at which to extract the cross-section, in degrees. + return_face_indices : bool, optional + If True, returns both the grid at the specified latitude and the indices + of the intersecting faces. If False, only the grid is returned. + Default is False. + method : str, optional + The internal method to use when identifying faces at the constant latitude. + Options are: + - 'fast': Uses a faster but potentially less accurate method for face identification. + - 'accurate': Uses a slower but more accurate method. + Default is 'fast'. + + Returns + ------- + grid_at_constant_lat : Grid + The grid with faces that interesected at a given lattitude + faces : array, optional + The indices of the faces that intersect with the specified latitude. This is only + returned if `return_face_indices` is set to True. + + Raises + ------ + ValueError + If no intersections are found at the specified latitude, a ValueError is raised. + + Examples + -------- + >>> grid, indices = grid.cross_section.constant_latitude(lat=30.0, return_face_indices=True) + >>> grid = grid.cross_section.constant_latitude(lat=-15.5) + + Notes + ----- + The accuracy and performance of the function can be controlled using the `method` parameter. + For higher precision requreiments, consider using method='acurate'. + """ + faces = self.uxgrid.get_faces_at_constant_latitude(lat, method) + + if len(faces) == 0: + raise ValueError(f"No intersections found at lat={lat}.") + + grid_at_constant_lat = self.uxgrid.isel(n_face=faces) + + if return_face_indices: + return grid_at_constant_lat, faces + else: + return grid_at_constant_lat + + def constant_longitude(self, *args, **kwargs): + raise NotImplementedError + + def gca(self, *args, **kwargs): + raise NotImplementedError + + def bounded_latitude(self, *args, **kwargs): + raise NotImplementedError + + def bounded_longitude(self, *args, **kwargs): + raise NotImplementedError + + def gca_gca(self, *args, **kwargs): + raise NotImplementedError diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 9770bfd7c..9dfdeb296 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -76,6 +76,8 @@ from uxarray.subset import GridSubsetAccessor +from uxarray.cross_sections import GridCrossSectionAccessor + from uxarray.grid.validation import ( _check_connectivity, _check_duplicate_nodes, @@ -217,6 +219,9 @@ def __init__( # declare subset accessor subset = UncachedAccessor(GridSubsetAccessor) + # declare cross section accessor + cross_section = UncachedAccessor(GridCrossSectionAccessor) + @classmethod def from_dataset( cls, dataset: xr.Dataset, use_dual: Optional[bool] = False, **kwargs @@ -1936,67 +1941,6 @@ def isel(self, **dim_kwargs): "Indexing must be along a grid dimension: ('n_node', 'n_edge', 'n_face')" ) - def constant_latitude_cross_section( - self, lat: float, return_face_indices=False, method="fast" - ): - """Extracts a cross-section of the grid at a specified constant - latitude. - - This method identifies and returns all faces (or grid elements) that intersect - with a given latitude. The returned cross-section can include either just the grid - or both the grid elements and the corresponding face indices, depending - on the `return_face_indices` parameter. - - Parameters - ---------- - lat : float - The latitude at which to extract the cross-section, in degrees. - return_face_indices : bool, optional - If True, returns both the grid at the specified latitude and the indices - of the intersecting faces. If False, only the grid is returned. - Default is False. - method : str, optional - The internal method to use when identifying faces at the constant latitude. - Options are: - - 'fast': Uses a faster but potentially less accurate method for face identification. - - 'accurate': Uses a slower but more accurate method. - Default is 'fast'. - - Returns - ------- - grid_at_constant_lat : Grid - The grid with faces that interesected at a given lattitude - faces : array, optional - The indices of the faces that intersect with the specified latitude. This is only - returned if `return_face_indices` is set to True. - - Raises - ------ - ValueError - If no intersections are found at the specified latitude, a ValueError is raised. - - Examples - -------- - >>> grid, indices = grid.constant_latitude_cross_section(lat=30.0, return_face_indices=True) - >>> grid = grid.constant_latitude_cross_section(lat=-15.5) - - Notes - ----- - The accuracy and performance of the function can be controlled using the `method` parameter. - For higher precision requreiments, consider using method='acurate'. - """ - faces = self.get_faces_at_constant_latitude(lat, method) - - if len(faces) == 0: - raise ValueError(f"No intersections found at lat={lat}.") - - grid_at_constant_lat = self.isel(n_face=faces) - - if return_face_indices: - return grid_at_constant_lat, faces - else: - return grid_at_constant_lat - def get_edges_at_constant_latitude(self, lat, method="fast"): """Identifies the edges of the grid that intersect with a specified constant latitude. From 7b78daa47aaad1022a1fc986006d77dda12107f9 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec Date: Wed, 9 Oct 2024 12:53:11 -0500 Subject: [PATCH 24/25] update notebook --- docs/user-guide/cross-sections.ipynb | 49 +++++++++++++++++----------- 1 file changed, 30 insertions(+), 19 deletions(-) diff --git a/docs/user-guide/cross-sections.ipynb b/docs/user-guide/cross-sections.ipynb index 7d3468608..9c0ba703c 100644 --- a/docs/user-guide/cross-sections.ipynb +++ b/docs/user-guide/cross-sections.ipynb @@ -12,13 +12,15 @@ }, { "cell_type": "code", + "execution_count": null, "id": "b35ba4a2c30750e4", "metadata": { "ExecuteTime": { - "end_time": "2024-10-08T01:49:10.528008Z", - "start_time": "2024-10-08T01:49:10.522526Z" + "end_time": "2024-10-09T17:50:50.244285Z", + "start_time": "2024-10-09T17:50:50.239653Z" } }, + "outputs": [], "source": [ "import uxarray as ux\n", "import geoviews.feature as gf\n", @@ -27,9 +29,7 @@ "import geoviews as gv\n", "\n", "projection = ccrs.Robinson()" - ], - "outputs": [], - "execution_count": 9 + ] }, { "cell_type": "markdown", @@ -45,8 +45,8 @@ "id": "b4160275c09fe6b0", "metadata": { "ExecuteTime": { - "end_time": "2024-10-05T15:21:17.331358Z", - "start_time": "2024-10-05T15:21:17.200324Z" + "end_time": "2024-10-09T17:50:51.217211Z", + "start_time": "2024-10-09T17:50:50.540946Z" } }, "outputs": [], @@ -71,32 +71,41 @@ "source": [ "## Constant Latitude\n", "\n", - "Cross-sections along constant latitude lines can be obtained using the ``constant_latitude_cross_section`` method, available for both ``ux.Grid`` and ``ux.DataArray`` objects. This functionality allows users to extract and analyze slices of data at specified latitudes, providing insights into variations along horizontal sections of the grid.\n" + "Cross-sections along constant latitude lines can be obtained using the ``.cross_section.constant_latitude`` method, available for both ``ux.Grid`` and ``ux.DataArray`` objects. This functionality allows users to extract and analyze slices of data at specified latitudes, providing insights into variations along horizontal sections of the grid.\n" ] }, { - "metadata": {}, "cell_type": "markdown", - "source": "For example, we can obtain a cross-section at 30 degrees latitude by doing the following:", - "id": "2fbe9f6e5bb59a17" + "id": "2fbe9f6e5bb59a17", + "metadata": {}, + "source": [ + "For example, we can obtain a cross-section at 30 degrees latitude by doing the following:" + ] }, { "cell_type": "code", "execution_count": null, "id": "3775daa1-2f1d-4738-bab5-2b69ebd689d9", - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-09T17:50:53.093314Z", + "start_time": "2024-10-09T17:50:53.077719Z" + } + }, "outputs": [], "source": [ "lat = 30\n", "\n", - "uxda_constant_lat = uxds[\"psi\"].constant_latitude_cross_section(lat)" + "uxda_constant_lat = uxds[\"psi\"].cross_section.constant_latitude(lat)" ] }, { - "metadata": {}, "cell_type": "markdown", - "source": "Since the result is a new ``UxDataArray``, we can directly plot the result to see the cross-section.", - "id": "dcec0b96b92e7f4" + "id": "dcec0b96b92e7f4", + "metadata": {}, + "source": [ + "Since the result is a new ``UxDataArray``, we can directly plot the result to see the cross-section." + ] }, { "cell_type": "code", @@ -120,10 +129,12 @@ ] }, { - "metadata": {}, "cell_type": "markdown", - "source": "You can also perform operations on the cross-section, such as taking the mean.", - "id": "c7cca7de4722c121" + "id": "c7cca7de4722c121", + "metadata": {}, + "source": [ + "You can also perform operations on the cross-section, such as taking the mean." + ] }, { "cell_type": "code", From 9f862e50a02b7a14f159cc0ea34b312dd52b9801 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Wed, 9 Oct 2024 19:24:34 -0500 Subject: [PATCH 25/25] fix docs --- docs/api.rst | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/docs/api.rst b/docs/api.rst index 2ff00c6d5..31f162f80 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -294,9 +294,6 @@ UxDataArray Cross Sections -------------- -.. seealso:: - - `Cross Sections User Guide Section `_ Grid ~~~~ @@ -316,7 +313,7 @@ UxDataArray :template: autosummary/accessor_method.rst UxDataArray.cross_section - UxDataArray.coss_section.constant_latitude + UxDataArray.cross_section.constant_latitude Remapping ---------