diff --git a/docs/notebooks/example_subset.ipynb b/docs/notebooks/example_subset.ipynb index 423e58c..a69445e 100644 --- a/docs/notebooks/example_subset.ipynb +++ b/docs/notebooks/example_subset.ipynb @@ -12,69 +12,11 @@ "cell_type": "code", "execution_count": 1, "metadata": {}, - "outputs": [], - "source": [ - "# Import packages\n", - "from pathlib import Path\n", - "import numpy as np\n", - "import xarray as xr\n", - "import suxarray as sx" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "# Where the test data is located.\n", - "# The test data is a subset from a Bay-Delta SCHISM model at Carquinez Strait.\n", - "dir_data = Path('../../tests/testdata')\n", - "\n", - "# Read the test data\n", - "files_out2d = [str(p) for p in [dir_data / 'out2d_1.nc', dir_data / 'out2d_2.nc']]\n", - "files_zcoords = [str(p) for p in [dir_data / 'zCoordinates_1.nc', dir_data / 'zCoordinates_2.nc']]\n", - "grid = sx.open_grid(files_out2d, files_zcoords)\n", - "\n", - "# We are converting the data to float64 for later calculation.\n", - "files_salinity = [dir_data / 'salinity_1.nc', dir_data / 'salinity_2.nc']\n", - "ds_salinity = xr.open_mfdataset(files_salinity, parallel=True).astype(np.float64)\n", - "sxds = sx.read_schism_nc(grid, ds_salinity)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Visualize the original data\n", - "Let's visualize the surface salinity data at the first time step from the dataset." - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, "outputs": [ - { - "data": { - "application/javascript": "(function(root) {\n function now() {\n return new Date();\n }\n\n var force = true;\n var py_version = '3.3.4'.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': {'jspanel': 'https://cdn.jsdelivr.net/npm/jspanel4@4.12.0/dist/jspanel', 'jspanel-modal': 'https://cdn.jsdelivr.net/npm/jspanel4@4.12.0/dist/extensions/modal/jspanel.modal', 'jspanel-tooltip': 'https://cdn.jsdelivr.net/npm/jspanel4@4.12.0/dist/extensions/tooltip/jspanel.tooltip', 'jspanel-hint': 'https://cdn.jsdelivr.net/npm/jspanel4@4.12.0/dist/extensions/hint/jspanel.hint', 'jspanel-layout': 'https://cdn.jsdelivr.net/npm/jspanel4@4.12.0/dist/extensions/layout/jspanel.layout', 'jspanel-contextmenu': 'https://cdn.jsdelivr.net/npm/jspanel4@4.12.0/dist/extensions/contextmenu/jspanel.contextmenu', 'jspanel-dock': 'https://cdn.jsdelivr.net/npm/jspanel4@4.12.0/dist/extensions/dock/jspanel.dock', 'gridstack': 'https://cdn.jsdelivr.net/npm/gridstack@7.2.3/dist/gridstack-all', 'notyf': 'https://cdn.jsdelivr.net/npm/notyf@3/notyf.min'}, 'shim': {'jspanel': {'exports': 'jsPanel'}, 'gridstack': {'exports': 'GridStack'}}});\n require([\"jspanel\"], function(jsPanel) {\n\twindow.jsPanel = jsPanel\n\ton_load()\n })\n require([\"jspanel-modal\"], function() {\n\ton_load()\n })\n require([\"jspanel-tooltip\"], function() {\n\ton_load()\n })\n require([\"jspanel-hint\"], function() {\n\ton_load()\n })\n require([\"jspanel-layout\"], function() {\n\ton_load()\n })\n require([\"jspanel-contextmenu\"], function() {\n\ton_load()\n })\n require([\"jspanel-dock\"], function() {\n\ton_load()\n })\n require([\"gridstack\"], function(GridStack) {\n\twindow.GridStack = GridStack\n\ton_load()\n })\n require([\"notyf\"], function() {\n\ton_load()\n })\n root._bokeh_is_loading = css_urls.length + 9;\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 } if (((window['jsPanel'] !== undefined) && (!(window['jsPanel'] instanceof HTMLElement))) || window.requirejs) {\n var urls = ['https://cdn.holoviz.org/panel/1.3.8/dist/bundled/floatpanel/jspanel4@4.12.0/dist/jspanel.js', 'https://cdn.holoviz.org/panel/1.3.8/dist/bundled/floatpanel/jspanel4@4.12.0/dist/extensions/modal/jspanel.modal.js', 'https://cdn.holoviz.org/panel/1.3.8/dist/bundled/floatpanel/jspanel4@4.12.0/dist/extensions/tooltip/jspanel.tooltip.js', 'https://cdn.holoviz.org/panel/1.3.8/dist/bundled/floatpanel/jspanel4@4.12.0/dist/extensions/hint/jspanel.hint.js', 'https://cdn.holoviz.org/panel/1.3.8/dist/bundled/floatpanel/jspanel4@4.12.0/dist/extensions/layout/jspanel.layout.js', 'https://cdn.holoviz.org/panel/1.3.8/dist/bundled/floatpanel/jspanel4@4.12.0/dist/extensions/contextmenu/jspanel.contextmenu.js', 'https://cdn.holoviz.org/panel/1.3.8/dist/bundled/floatpanel/jspanel4@4.12.0/dist/extensions/dock/jspanel.dock.js'];\n for (var i = 0; i < urls.length; i++) {\n skip.push(urls[i])\n }\n } if (((window['GridStack'] !== undefined) && (!(window['GridStack'] instanceof HTMLElement))) || window.requirejs) {\n var urls = ['https://cdn.holoviz.org/panel/1.3.8/dist/bundled/gridstack/gridstack@7.2.3/dist/gridstack-all.js'];\n for (var i = 0; i < urls.length; i++) {\n skip.push(urls[i])\n }\n } if (((window['Notyf'] !== undefined) && (!(window['Notyf'] instanceof HTMLElement))) || window.requirejs) {\n var urls = ['https://cdn.holoviz.org/panel/1.3.8/dist/bundled/notificationarea/notyf@3/notyf.min.js'];\n for (var i = 0; i < urls.length; i++) {\n skip.push(urls[i])\n }\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.3.4.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-gl-3.3.4.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-widgets-3.3.4.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-tables-3.3.4.min.js\", \"https://cdn.holoviz.org/panel/1.3.8/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));", - "application/vnd.holoviews_load.v0+json": "" - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "application/javascript": "\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", - "application/vnd.holoviews_load.v0+json": "" - }, - "metadata": {}, - "output_type": "display_data" - }, { "data": { "text/html": [ - "" ] }, "metadata": {}, "output_type": "display_data" }, + { + "data": { + "application/javascript": "(function(root) {\n function now() {\n return new Date();\n }\n\n const force = false;\n const version = '3.8.2'.replace('rc', '-rc.').replace('.dev', '-dev.');\n const reloading = true;\n const Bokeh = root.Bokeh;\n const BK_RE = /^https:\\/\\/cdn\\.bokeh\\.org\\/bokeh\\/(release|dev)\\/bokeh-/;\n const PN_RE = /^https:\\/\\/cdn\\.holoviz\\.org\\/panel\\/[^/]+\\/dist\\/panel/i;\n\n // Set a timeout for this load but only if we are not already initializing\n if (typeof (root._bokeh_timeout) === \"undefined\" || (force || !root._bokeh_is_initializing)) {\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, Bokeh, 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 // Don't load bokeh if it is still initializing\n console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n return null;\n } else if (js_urls.length === 0 && js_modules.length === 0 && Object.keys(js_exports).length === 0) {\n // There is nothing to load\n run_callbacks();\n return null;\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(e) {\n const src_el = e.srcElement\n console.error(\"failed to load \" + (src_el.href || src_el.src));\n }\n\n const 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 const existing_stylesheets = []\n const links = document.getElementsByTagName('link')\n for (let i = 0; i < links.length; i++) {\n const link = links[i]\n if (link.href != null) {\n existing_stylesheets.push(link.href)\n }\n }\n for (let i = 0; i < css_urls.length; i++) {\n const url = css_urls[i];\n const escaped = encodeURI(url)\n if (existing_stylesheets.indexOf(escaped) !== -1) {\n on_load()\n continue;\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 const scripts = document.getElementsByTagName('script')\n for (let i = 0; i < scripts.length; i++) {\n var script = scripts[i]\n if (script.src != null) {\n existing_scripts.push(script.src)\n }\n }\n for (let i = 0; i < js_urls.length; i++) {\n const url = js_urls[i];\n const escaped = encodeURI(url)\n const shouldSkip = skip.includes(escaped) || existing_scripts.includes(escaped)\n const isBokehOrPanel = BK_RE.test(escaped) || PN_RE.test(escaped)\n const missingOrBroken = Bokeh == null || Bokeh.Panel == null || (Bokeh.version != version && !Bokeh.versions?.has(version)) || Bokeh.versions?.get(version)?.Panel == null;\n if (shouldSkip && !(isBokehOrPanel && missingOrBroken)) {\n if (!window.requirejs) {\n on_load();\n }\n continue;\n }\n const 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 (let i = 0; i < js_modules.length; i++) {\n const url = js_modules[i];\n const escaped = encodeURI(url)\n if (skip.indexOf(escaped) !== -1 || existing_scripts.indexOf(escaped) !== -1) {\n if (!window.requirejs) {\n on_load();\n }\n continue;\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 const url = js_exports[name];\n const escaped = encodeURI(url)\n if (skip.indexOf(escaped) >= 0 || root[name] != null) {\n if (!window.requirejs) {\n on_load();\n }\n continue;\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 const js_urls = [\"https://cdn.holoviz.org/panel/1.9.3/dist/bundled/reactiveesm/es-module-shims@^1.10.0/dist/es-module-shims.min.js\"];\n const js_modules = [];\n const js_exports = {};\n const css_urls = [];\n const 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 (let i = 0; i < inline_js.length; i++) {\n try {\n inline_js[i].call(root, root.Bokeh);\n } catch(e) {\n if (!reloading) {\n throw e;\n }\n }\n }\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 // If the timeout and bokeh was not successfully loaded we reset\n // everything and try loading again\n root._bokeh_timeout = Date.now() + 5000;\n root._bokeh_is_initializing = false;\n root._bokeh_onload_callbacks = undefined;\n root._bokeh_is_loading = 0;\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 const bokeh_loaded = Bokeh != null && ((Bokeh.version === version && Bokeh.Panel) || (Bokeh.versions?.has(version) && Bokeh.versions.get(version)?.Panel));\n if (!reloading && !bokeh_loaded) {\n if (root.Bokeh) {\n root.Bokeh = undefined;\n }\n console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n }\n load_libs(css_urls, js_urls, js_modules, js_exports, Bokeh, function() {\n console.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n run_inline_js();\n if (Bokeh != undefined && !reloading) {\n const NewBokeh = root.Bokeh;\n if (Bokeh.versions === undefined) {\n Bokeh.versions = new Map();\n }\n if (NewBokeh.version !== Bokeh.version) {\n Bokeh[NewBokeh.version] = NewBokeh;\n Bokeh.versions.set(NewBokeh.version, NewBokeh);\n }\n root.Bokeh = Bokeh;\n }\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": "" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/javascript": "\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 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 let retries = 0;\n const open = () => {\n if (comm.active) {\n comm.open();\n } else if (retries > 3) {\n console.warn('Comm target never activated')\n } else {\n retries += 1\n setTimeout(open, 500)\n }\n }\n if (comm.active) {\n comm.open();\n } else {\n setTimeout(open, 500)\n }\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", + "application/vnd.holoviews_load.v0+json": "" + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Import packages\n", + "from pathlib import Path\n", + "import numpy as np\n", + "import xarray as xr\n", + "import suxarray as sx" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# Where the test data is located.\n", + "# The test data is a subset from a Bay-Delta SCHISM model at Carquinez Strait.\n", + "dir_data = Path('../../tests/testdata')\n", + "\n", + "# Read the test data\n", + "files_out2d = [str(p) for p in [dir_data / 'out2d_1.nc', dir_data / 'out2d_2.nc']]\n", + "files_zcoords = [str(p) for p in [dir_data / 'zCoordinates_1.nc', dir_data / 'zCoordinates_2.nc']]\n", + "grid = sx.open_grid(files_out2d, files_zcoords)\n", + "\n", + "# We are converting the data to float64 for later calculation.\n", + "files_salinity = [dir_data / 'salinity_1.nc', dir_data / 'salinity_2.nc']\n", + "ds_salinity = xr.open_mfdataset(files_salinity, parallel=True).astype(np.float64)\n", + "\n", + "ds_salinity.attrs = {\n", + " \"Conventions\": \"CF-1.8\",\n", + " \"title\": \"Bay-Delta SCHISM model at Carquinez Strait\",\n", + " \"source\": \"SCHISM model\",\n", + "}\n", + "\n", + "sxds = sx.read_schism_nc(grid, ds_salinity)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'Conventions': 'CF-1.8',\n", + " 'title': 'Bay-Delta SCHISM model at Carquinez Strait',\n", + " 'source': 'SCHISM model'}" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ds_salinity.attrs" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'Conventions': 'CF-1.8',\n", + " 'title': 'Bay-Delta SCHISM model at Carquinez Strait',\n", + " 'source': 'SCHISM model'}" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sxds.attrs" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Visualize the original data\n", + "Let's visualize the surface salinity data at the first time step from the dataset." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ { "data": {}, "metadata": {}, @@ -272,26 +279,31 @@ "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ - "
\n", - "
\n", + "
\n", + "
\n", "
\n", "" + ], + "text/plain": [ + ":Image [x,y] (x_y volume)" + ] + }, + "execution_count": 10, + "metadata": { + "application/vnd.holoviews_exec.v0+json": { + "id": "8378d635-88bf-4be2-9869-90a8f61c1e1f" + } + }, + "output_type": "execute_result" + } + ], + "source": [ + "# Visualize total water column volume per face at first timestep\n", + "sxda_volumes.sum(dim=\"n_layer\").isel(time=0).plot()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Compare planar vs uxarray volume\n", + "We compute volumes using both implementations, then compare per-layer values and total water-column volumes." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.SxDataArray 'volume_ux' (time: 96, n_layer: 22, n_face: 66)> Size: 1MB\n",
+       "array([[[   0.        ,    0.        ,    0.        , ...,\n",
+       "            0.        ,    0.        ,    0.        ],\n",
+       "        [   0.        ,    0.        ,    0.        , ...,\n",
+       "            0.        ,    0.        ,    0.        ],\n",
+       "        [   0.        ,    0.        ,    0.        , ...,\n",
+       "            0.        ,    0.        ,    0.        ],\n",
+       "        ...,\n",
+       "        [6716.0632124 , 5864.38091235, 6530.21158725, ...,\n",
+       "          470.28959784,  710.04427853,  575.72713514],\n",
+       "        [6265.7617329 , 5541.57429678, 6092.52200288, ...,\n",
+       "          470.28171698,  710.0320048 ,  575.71737839],\n",
+       "        [5845.67645761, 5240.50953111, 5684.13009277, ...,\n",
+       "          470.28976202,  710.04450169,  575.72713514]],\n",
+       "\n",
+       "       [[   0.        ,    0.        ,    0.        , ...,\n",
+       "            0.        ,    0.        ,    0.        ],\n",
+       "        [   0.        ,    0.        ,    0.        , ...,\n",
+       "            0.        ,    0.        ,    0.        ],\n",
+       "        [   0.        ,    0.        ,    0.        , ...,\n",
+       "            0.        ,    0.        ,    0.        ],\n",
+       "...\n",
+       "          434.42529919,  660.95692554,  532.88465672],\n",
+       "        [6222.50068794, 5499.71410167, 6051.2940302 , ...,\n",
+       "          434.41807507,  660.94532129,  532.87548537],\n",
+       "        [5805.3152958 , 5200.87567394, 5645.66589242, ...,\n",
+       "          434.42562756,  660.95692554,  532.88485186]],\n",
+       "\n",
+       "       [[   0.        ,    0.        ,    0.        , ...,\n",
+       "            0.        ,    0.        ,    0.        ],\n",
+       "        [   0.        ,    0.        ,    0.        , ...,\n",
+       "            0.        ,    0.        ,    0.        ],\n",
+       "        [   0.        ,    0.        ,    0.        , ...,\n",
+       "            0.        ,    0.        ,    0.        ],\n",
+       "        ...,\n",
+       "        [6692.95996173, 5842.25828605, 6508.63028394, ...,\n",
+       "          452.99850691,  686.66995425,  555.22781132],\n",
+       "        [6244.20752314, 5520.64313315, 6072.38743961, ...,\n",
+       "          452.99079023,  686.65812684,  555.2182497 ],\n",
+       "        [5825.56724483, 5220.69009412, 5665.34477109, ...,\n",
+       "          452.99850691,  686.66995425,  555.22781132]]],\n",
+       "      shape=(96, 22, 66))\n",
+       "Coordinates:\n",
+       "  * time     (time) datetime64[ns] 768B 2014-12-17T16:30:00 ... 2014-12-19T16...\n",
+       "Dimensions without coordinates: n_layer, n_face
" + ], + "text/plain": [ + " Size: 1MB\n", + "array([[[ 0. , 0. , 0. , ...,\n", + " 0. , 0. , 0. ],\n", + " [ 0. , 0. , 0. , ...,\n", + " 0. , 0. , 0. ],\n", + " [ 0. , 0. , 0. , ...,\n", + " 0. , 0. , 0. ],\n", + " ...,\n", + " [6716.0632124 , 5864.38091235, 6530.21158725, ...,\n", + " 470.28959784, 710.04427853, 575.72713514],\n", + " [6265.7617329 , 5541.57429678, 6092.52200288, ...,\n", + " 470.28171698, 710.0320048 , 575.71737839],\n", + " [5845.67645761, 5240.50953111, 5684.13009277, ...,\n", + " 470.28976202, 710.04450169, 575.72713514]],\n", + "\n", + " [[ 0. , 0. , 0. , ...,\n", + " 0. , 0. , 0. ],\n", + " [ 0. , 0. , 0. , ...,\n", + " 0. , 0. , 0. ],\n", + " [ 0. , 0. , 0. , ...,\n", + " 0. , 0. , 0. ],\n", + "...\n", + " 434.42529919, 660.95692554, 532.88465672],\n", + " [6222.50068794, 5499.71410167, 6051.2940302 , ...,\n", + " 434.41807507, 660.94532129, 532.87548537],\n", + " [5805.3152958 , 5200.87567394, 5645.66589242, ...,\n", + " 434.42562756, 660.95692554, 532.88485186]],\n", + "\n", + " [[ 0. , 0. , 0. , ...,\n", + " 0. , 0. , 0. ],\n", + " [ 0. , 0. , 0. , ...,\n", + " 0. , 0. , 0. ],\n", + " [ 0. , 0. , 0. , ...,\n", + " 0. , 0. , 0. ],\n", + " ...,\n", + " [6692.95996173, 5842.25828605, 6508.63028394, ...,\n", + " 452.99850691, 686.66995425, 555.22781132],\n", + " [6244.20752314, 5520.64313315, 6072.38743961, ...,\n", + " 452.99079023, 686.65812684, 555.2182497 ],\n", + " [5825.56724483, 5220.69009412, 5665.34477109, ...,\n", + " 452.99850691, 686.66995425, 555.22781132]]],\n", + " shape=(96, 22, 66))\n", + "Coordinates:\n", + " * time (time) datetime64[ns] 768B 2014-12-17T16:30:00 ... 2014-12-19T16...\n", + "Dimensions without coordinates: n_layer, n_face" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Compute uxarray-based per-layer face volumes (lat/lon face areas scaled by R^2)\n", + "volumes_ux = _calculate_volumes_uxarray(sxda_subset)\n", + "\n", + "sxda_volumes_ux = sx.SxDataArray(\n", + " xr.DataArray(\n", + " volumes_ux,\n", + " dims=[\"time\", \"n_layer\", \"n_face\"],\n", + " coords={\"time\": sxda_subset.time},\n", + " name=\"volume_ux\",\n", + " ),\n", + " sxgrid=sxda_subset.sxgrid,\n", + ")\n", + "\n", + "sxda_volumes_ux" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Per-layer abs diff max: 37.250629562274526\n", + "Per-layer abs diff mean: 1.7644504003320332\n", + "Per-layer rel diff max: 0.029262378558841222\n", + "Per-layer rel diff mean: 0.0009774319691265771\n" + ] + } + ], + "source": [ + "# Per-layer, per-face comparison statistics\n", + "layer_abs_diff = np.abs(sxda_volumes_ux.values - sxda_volumes.values)\n", + "layer_rel_diff = layer_abs_diff / np.maximum(np.abs(sxda_volumes.values), 1.0e-12)\n", + "\n", + "print(\"Per-layer abs diff max:\", layer_abs_diff.max())\n", + "print(\"Per-layer abs diff mean:\", layer_abs_diff.mean())\n", + "print(\"Per-layer rel diff max:\", layer_rel_diff.max())\n", + "print(\"Per-layer rel diff mean:\", layer_rel_diff.mean())" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Pointwise rel diff mean: -0.00022845575726045358\n", + "Pointwise |rel diff| mean: 0.0011783982481104157\n", + "Pointwise |rel diff| max: 0.023980858656040164\n", + "Timestep-normalized |rel diff| mean: 0.00020286449568462282\n", + "Timestep-normalized |rel diff| max: 0.0021381875189729455\n", + "Per-time mean |rel diff| (first 10): [0.00019809 0.00019846 0.00019906 0.00019991 0.00020094 0.00020207\n", + " 0.00020312 0.00020424 0.00020518 0.00020605]\n", + "Per-time max |rel diff| (first 10): [0.00199495 0.00199951 0.00200727 0.00201832 0.00203294 0.00204866\n", + " 0.00206348 0.00207965 0.00209335 0.00210563]\n" + ] + } + ], + "source": [ + "# Compare total water-column volume per face (sum over layers)\n", + "sxda_volumes_total = sxda_volumes.sum(dim=\"n_layer\")\n", + "sxda_volumes_ux_total = sxda_volumes_ux.sum(dim=\"n_layer\")\n", + "\n", + "# Signed difference: uxarray - planar, shape (time, n_face)\n", + "total_diff = sxda_volumes_ux_total.values - sxda_volumes_total.values\n", + "total_abs_diff = np.abs(total_diff)\n", + "\n", + "eps = 1.0e-12\n", + "\n", + "# 1) Pointwise relative error (denominator is each face at each timestep)\n", + "denom_pointwise = np.maximum(np.abs(sxda_volumes_total.values), eps)\n", + "total_rel_diff_pointwise = total_diff / denom_pointwise\n", + "\n", + "# 2) Timestep-normalized relative error (denominator is max face volume at each timestep)\n", + "denom_timestep = np.maximum(\n", + " np.abs(sxda_volumes_total.values).max(axis=1, keepdims=True),\n", + " eps,\n", + " )\n", + "total_rel_diff_timestep = total_diff / denom_timestep\n", + "\n", + "# Per-timestep summaries\n", + "mean_abs_rel_per_time = np.mean(np.abs(total_rel_diff_timestep), axis=1)\n", + "max_abs_rel_per_time = np.max(np.abs(total_rel_diff_timestep), axis=1)\n", + "\n", + "print(\"Pointwise rel diff mean:\", np.mean(total_rel_diff_pointwise))\n", + "print(\"Pointwise |rel diff| mean:\", np.mean(np.abs(total_rel_diff_pointwise)))\n", + "print(\"Pointwise |rel diff| max:\", np.max(np.abs(total_rel_diff_pointwise)))\n", + "\n", + "print(\"Timestep-normalized |rel diff| mean:\", np.mean(np.abs(total_rel_diff_timestep)))\n", + "print(\"Timestep-normalized |rel diff| max:\", np.max(np.abs(total_rel_diff_timestep)))\n", + "\n", + "print(\"Per-time mean |rel diff| (first 10):\", mean_abs_rel_per_time[:10])\n", + "print(\"Per-time max |rel diff| (first 10):\", max_abs_rel_per_time[:10])" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": {}, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.holoviews_exec.v0+json": "", + "text/html": [ + "
\n", + "
\n", + "
\n", + "" + ], + "text/plain": [ + ":Image [x,y] (x_y column_volume_rel_abs_diff_pointwise)" + ] + }, + "execution_count": 14, + "metadata": { + "application/vnd.holoviews_exec.v0+json": { + "id": "dbdb739c-cf12-4144-86b0-a581fff0fad1" + } + }, + "output_type": "execute_result" + } + ], + "source": [ + "from bokeh.models import PrintfTickFormatter\n", + "\n", + "# Visualize only #3: magnitude of relative difference (pointwise denominator)\n", + "time_i = 0\n", + "\n", + "# Fixed color limit for relative magnitude map\n", + "rel_lim = 0.005\n", + "\n", + "# Build a DataArray from |(uxarray - planar) / planar|\n", + "sxda_total_rel_abs_diff_pointwise = sx.SxDataArray(\n", + " xr.DataArray(\n", + " np.abs(total_rel_diff_pointwise),\n", + " dims=sxda_volumes_total.dims,\n", + " coords=sxda_volumes_total.coords,\n", + " name=\"column_volume_rel_abs_diff_pointwise\",\n", + " ),\n", + " sxgrid=sxda_subset.sxgrid,\n", + ")\n", + "\n", + "# Plot relative difference magnitude map with fixed decimal formatting for\n", + "# both colorbar ticks and hover values (avoids scientific notation).\n", + "sxda_total_rel_abs_diff_pointwise.isel(time=time_i).plot(cmap=\"Reds\").opts(\n", + " width=600,\n", + " height=350,\n", + " clim=(0.0, rel_lim),\n", + " colorbar_opts={\"formatter\": PrintfTickFormatter(format=\"%.4f\")},\n", + " cformatter=\"%.4f\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Calculate mass" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "sxda_mass = sx.SxDataArray(\n", + " xr.DataArray(\n", + " _calculate_mass(sxda_subset),\n", + " dims=[\"time\", \"n_layer\", \"n_face\"],\n", + " coords={\"time\": sxda_subset.time},\n", + " name=\"mass\",\n", + " ),\n", + " sxgrid=sxda_subset.sxgrid,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": {}, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.holoviews_exec.v0+json": "", + "text/html": [ + "
\n", + "
\n", + "
\n", + "" + ], + "text/plain": [ + ":Image [x,y] (x_y mass)" + ] + }, + "execution_count": 16, + "metadata": { + "application/vnd.holoviews_exec.v0+json": { + "id": "7685dc78-5f7a-43e0-b238-d4f2989a5d94" + } + }, + "output_type": "execute_result" + } + ], + "source": [ + "# Visualize total water column mass per face at first timestep\n", + "sxda_mass.sum(dim=\"n_layer\").isel(time=0).plot()" + ] } ], "metadata": { "kernelspec": { - "display_name": "suxarray_dev", + "display_name": "suxarray", "language": "python", "name": "python3" }, @@ -483,7 +5291,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.8" + "version": "3.12.12" } }, "nbformat": 4, diff --git a/pyproject.toml b/pyproject.toml index 894898f..ca4dde9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,6 +28,7 @@ classifiers = [ dependencies = [ "uxarray==2025.12.0", "shapely>=2.0", + "gsw>=3.6", # required for density -> mass calculations in `integrate.py` "Click>=7.0" ] diff --git a/requirements_dev.txt b/requirements_dev.txt index 0f3644b..de52887 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -14,3 +14,4 @@ uxarray==2025.12.0 netcdf>=1.6 xarray>=2023.4.1 holoviews>=1.16.0 +gsw>=3.6 # required for density -> mass calculations in `integrate.py` diff --git a/suxarray/grid/integrate.py b/suxarray/grid/integrate.py new file mode 100644 index 0000000..8529c61 --- /dev/null +++ b/suxarray/grid/integrate.py @@ -0,0 +1,288 @@ +import numpy as np +from .area import _integrate_nodal + + +def _earth_radius_at_latitude(latitude_deg): + """Compute geocentric Earth radius (meters) at geodetic latitude. + + Uses the WGS84 ellipsoid semi-major/minor axes. + """ + # WGS84 ellipsoid axes (meters) + a = 6378137.0 + b = 6356752.314245 + + lat = np.deg2rad(latitude_deg) + cos_lat = np.cos(lat) + sin_lat = np.sin(lat) + + a2 = a * a + b2 = b * b + + numerator = (a2 * cos_lat) ** 2 + (b2 * sin_lat) ** 2 + denominator = (a * cos_lat) ** 2 + (b * sin_lat) ** 2 + return np.sqrt(numerator / denominator) + + +def _face_mean_from_nodal(values, connectivity): + """Average nodal values onto faces. + + Parameters + ---------- + values : ndarray, shape (n_time, n_node, n_interface) + Nodal values to average. + connectivity : ndarray, shape (n_face, n_max_face_nodes) + Face-node connectivity with negative fill values. + + Returns + ------- + ndarray, shape (n_time, n_interface, n_face) + Face-averaged values. + """ + n_time, _, n_interface = values.shape + n_face = connectivity.shape[0] + face_mean = np.zeros((n_time, n_interface, n_face), dtype=np.float64) + + for face_i in range(n_face): + node_ids = connectivity[face_i] + node_ids = node_ids[node_ids >= 0] + face_mean[:, :, face_i] = values[:, node_ids, :].mean(axis=1) + + return face_mean + + +def _calc_layers_thickness(z, bottom_indices): + """ + Calculate the thickness of each layer under a set of nodes. + + Parameters + z: ndarray, shape (n_time, n_node, n_layer) + z-coordinates of nodes + bottom_indices: ndarray, shape (n_node,) + 0-based bottom layer index for each node + + Returns + ndarray, shape (n_time, n_layer-1, n_node) + """ + + thickness = z[:, :, 1:] - z[:, :, :-1] + + if np.any(bottom_indices > 0): + layer_idx = np.arange(thickness.shape[2]) + below_bottom = layer_idx[np.newaxis, :] < bottom_indices[:, np.newaxis] + thickness = np.where(below_bottom[np.newaxis, :, :], 0.0, thickness) + + return thickness + + +def _calculate_prism_volumes(node_x, node_y, z, connectivity, bottom_indices): + """ + Calculate the volume of prism under a selection of nodes for each layer in + the vertical direction. + + Parameters + ---------- + node_x : ndarray, shape (n_node,) + x coordinates of nodes + node_y : ndarray, shape (n_node,) + y coordinates of nodes + z : ndarray, shape (n_time, n_node, n_layer) + z-coordinates + connectivity : ndarray, shape (n_face, n_node_per_face) + Node indices for each face + bottom_indices : ndarray, shape (n_element,) + 0-based bottom layer index for each node + + Returns + ------- + ndarray, shape (n_time, n_layer-1, n_face) + Volume of each prism layer under the selection of nodes. + """ + # Enforce float64 for numba + node_x = node_x.astype(np.float64) + node_y = node_y.astype(np.float64) + z = z.astype(np.float64) + + # Extract dimensions + n_time, n_node, n_layer = z.shape + n_face = connectivity.shape[0] + n_interfaces = n_layer - 1 + + thickness = _calc_layers_thickness(z, bottom_indices) + + # Integrate thickness over elements per timestep + result = np.zeros((n_time, n_interfaces, n_face)) + for t_i in range(n_time): + h = thickness[t_i].T + result[t_i] = _integrate_nodal(node_x, node_y, h, connectivity) + + return result + + +def calc_volumes_from_projected(sxda): + """ + Calculate the volume of prism under a selection of nodes for each layer in + the vertical direction. + + Parameters + ---------- + sxda : SxDataArray + SxDataArray with 3D z-coordinates and bottom indices + + Returns + ------- + ndarray, shape (n_time, n_layer-1, n_face) + Volume of each prism layer under the selection of nodes. + """ + grid = sxda.sxgrid + + return _calculate_prism_volumes( + node_x=grid.node_x.values, + node_y=grid.node_y.values, + z=grid.sgrid_info.zCoordinates.values, + connectivity=grid.face_node_connectivity.values, + bottom_indices=grid.sgrid_info.bottom_index_node.values - 1, + ) + +# TODO: Remove time indexing loops in function and have them only calculate at +# one timestep + +def _calculate_volumes_uxarray(sxda, radius_m=None): + """Calculate per-face prism volumes using uxarray spherical face areas. + + This path uses uxarray's unit-sphere face areas from lat/lon geometry, + then scales by a physical sphere radius in meters. + + Parameters + ---------- + sxda : SxDataArray + SxDataArray with 3D z-coordinates and bottom indices. + radius_m : float, optional + Physical sphere radius in meters used to scale uxarray face areas. + If not provided, computes a local geocentric Earth radius from the + mean node latitude. + + Returns + ------- + ndarray, shape (n_time, n_layer-1, n_face) + Volume of each prism layer under each face element. + """ + grid = sxda.sxgrid + + z = grid.sgrid_info.zCoordinates.values.astype(np.float64) + connectivity = grid.face_node_connectivity.values + + # Thickness at nodes: shape (n_time, n_node, n_interface) + bottom_indices = grid.sgrid_info.bottom_index_node.values - 1 + thickness = _calc_layers_thickness(z, bottom_indices) + + # Face-average nodal thickness: shape (n_time, n_interface, n_face) + # NOTE: This is an arithmetic node mean. For strict parity with the planar + # path on quads, we may need to mirror its bilinear quadrature weighting. + thickness_face = _face_mean_from_nodal(thickness, connectivity) + + # Compute unit-sphere face areas from lon/lat directly to avoid mutating + # the source grid's planar cartesian coordinates. + from uxarray.grid.area import get_all_face_area_from_coords + + node_lon = np.asarray(grid.node_lon.values, dtype=np.float64) + node_lat = np.asarray(grid.node_lat.values, dtype=np.float64) + lon_rad = np.deg2rad(node_lon) + lat_rad = np.deg2rad(node_lat) + node_x_unit = np.cos(lat_rad) * np.cos(lon_rad) + node_y_unit = np.cos(lat_rad) * np.sin(lon_rad) + node_z_unit = np.sin(lat_rad) + + n_nodes_per_face = np.sum(connectivity >= 0, axis=1).astype(np.int64) + face_areas_unit, _ = get_all_face_area_from_coords( + node_x_unit, + node_y_unit, + node_z_unit, + connectivity, + n_nodes_per_face, + quadrature_rule="triangular", + order=4, + latitude_adjusted_area=False, + ) + + # Areas are on a unit sphere, so scale to m^2. + if radius_m is None: + mean_lat = np.asarray(grid.node_lat.values, dtype=np.float64).mean() + radius_m = _earth_radius_at_latitude(mean_lat) + else: + radius_m = float(radius_m) + face_areas_m2 = np.asarray(face_areas_unit, dtype=np.float64) * (radius_m ** 2) + + return thickness_face * face_areas_m2[np.newaxis, np.newaxis, :] + + +# TODO: implement generic prism integration of variables (such as energy, +# chlorophyll or salt) to complement mass calculation function + +# TODO: Test impact of gsw method on performance. + +def _calculate_mass_in_ocean(sxda_salinity, sxda_temperature=None): + """ + Calculate the total water mass in the domain using salinity and temperature + to compute in-situ density via the TEOS-10 equations (gsw). + + This function is better applied in a ocean domain and not in shallow/estuary + domains. + + Parameters + ---------- + sxda_salinity : SxDataArray + Practical salinity (PSU), shape (n_time, n_node, n_layer) + sxda_temperature : SxDataArray, optional + In-situ temperature (°C), shape (n_time, n_node, n_layer). + Defaults to 15.0 °C if not provided. + + Returns + ------- + ndarray, shape (n_time, n_layer-1, n_face) + Water mass per prism layer per face element (kg). + """ + import gsw + + grid = sxda_salinity.sxgrid + + # Pressure from z-coordinates: shape (n_time, n_node, n_layer) + z = grid.sgrid_info.zCoordinates.values + # Expand to (1, n_node, 1) to broadcast with z shape (n_time,n_node,n_layer) + lat = grid.node_lat.values[np.newaxis, :, np.newaxis] + lon = grid.node_lon.values[np.newaxis, :, np.newaxis] + pressure = gsw.p_from_z(z, lat) # shape (n_time, n_node, n_layer) + + # Convert salinity to Absolute Salinity (g/kg) + SP = sxda_salinity.values + SA = gsw.SA_from_SP(SP, pressure, lon, lat) # shape (n_time,n_node,n_layer) + + # Convert temperature to Conservative Temperature + t = sxda_temperature.values if sxda_temperature is not None else 15.0 + CT = gsw.CT_from_t(SA, t, pressure) # shape (n_time, n_node, n_layer) + + # density: shape (n_time, n_node, n_layer) + rho = gsw.density.rho(SA, CT, pressure) + + # Average density to layer intervals: shape (n_time, n_node, n_layer-1) + rho_layer = 0.5 * (rho[:, :, :-1] + rho[:, :, 1:]) + + # Thickness per layer interval: shape (n_time, n_node, n_layer-1) + bottom_indices = grid.sgrid_info.bottom_index_node.values - 1 + thickness = _calc_layers_thickness(z, bottom_indices) + + # rho * h at each node per layer: shape (n_time, n_node, n_layer-1) + rho_h = rho_layer * thickness + + # Integrate rho*h over faces: shape (n_time, n_layer-1, n_face) + connectivity = grid.face_node_connectivity.values + node_x = grid.node_x.values.astype(np.float64) + node_y = grid.node_y.values.astype(np.float64) + n_time, n_node, n_interfaces = rho_h.shape + n_face = connectivity.shape[0] + mass = np.zeros((n_time, n_interfaces, n_face)) + for t_i in range(n_time): + mass[t_i] = _integrate_nodal( + node_x, node_y, rho_h[t_i].T.astype(np.float64), connectivity + ) + + return mass \ No newline at end of file diff --git a/tests/test_volume.py b/tests/test_volume.py new file mode 100644 index 0000000..1b2c742 --- /dev/null +++ b/tests/test_volume.py @@ -0,0 +1,84 @@ +import numpy as np +import pytest + +from suxarray.grid.integrate import _calculate_prism_volumes + + +def test_calculate_prism_volume_synthetic_shaved_cells(): + """Validate prism volumes on a tiny synthetic triangular mesh. + + The mesh has 3 triangular faces of equal area (2 m^2 each), and node-wise + shaved-cell bottom indices that mask lower interfaces at selected nodes. + """ + node_x = np.array([0.0, 2.0, 0.0, -2.0, 0.0], dtype=np.float64) + node_y = np.array([0.0, 0.0, 2.0, 0.0, -2.0], dtype=np.float64) + + # 3 triangles with -1 fill values in the 4th slot. + connectivity = np.array( + [ + [0, 1, 2, -1], + [0, 2, 3, -1], + [0, 3, 4, -1], + ], + dtype=np.int64, + ) + # 0-based node bottom indices. k > 0 masks interfaces [:k] at that node. + bottom_indices = np.array([0, 1, 2, 0, 1], dtype=np.int64) + + # z shape: (n_time=2, n_node=5, n_layer=4) -> 3 interfaces. + z = np.array( + [ + [ + [-3.0, -2.0, -1.0, 0.0], + [-3.0, -2.0, -1.0, 0.0], + [-3.0, -2.0, -1.0, 0.0], + [-3.0, -2.0, -1.0, 0.0], + [-3.0, -2.0, -1.0, 0.0], + ], + [ + [-6.0, -4.0, -2.0, 0.0], + [-6.0, -4.0, -2.0, 0.0], + [-6.0, -4.0, -2.0, 0.0], + [-6.0, -4.0, -2.0, 0.0], + [-6.0, -4.0, -2.0, 0.0], + ], + ], + dtype=np.float64, + ) + + volumes = _calculate_prism_volumes( + node_x=node_x, + node_y=node_y, + z=z, + connectivity=connectivity, + bottom_indices=bottom_indices, + ) + + # Expected shape: (n_time, n_interface, n_face) + assert volumes.shape == (2, 3, 3) + + # Make the geometric factor explicit for readers: + # expected volume = triangle_area * face-mean(thickness after masking) + triangle_area_m2 = 2.0 + expected_face_mean_thickness = np.array( + [ + [ + [1.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0], + [2.0 / 3.0, 2.0 / 3.0, 1.0], + [1.0, 1.0, 1.0], + ], + [ + [2.0 / 3.0, 4.0 / 3.0, 4.0 / 3.0], + [4.0 / 3.0, 4.0 / 3.0, 2.0], + [2.0, 2.0, 2.0], + ], + ], + dtype=np.float64, + ) + expected = triangle_area_m2 * expected_face_mean_thickness + + np.testing.assert_allclose(volumes, expected, rtol=0.0, atol=1e-12) + + # Spot-check a couple of entries with pytest.approx for readability. + assert volumes[0, 0, 0] == pytest.approx(2.0 / 3.0) + assert volumes[1, 2, 2] == pytest.approx(4.0) \ No newline at end of file