From 620bcf0a8e7934c724bfa471c074889eaef5a48f Mon Sep 17 00:00:00 2001 From: ThomasHaine Date: Thu, 6 Jul 2023 13:47:05 -0400 Subject: [PATCH 1/3] Add sciserver_notebooks/*.md and rename --- docs/_toc.yml | 5 +- docs/get_started.md | 4 +- docs/network_of_object.md | 4 +- docs/sciserver_notebooks/Fjord.md | 223 ---------------------------- docs/sciserver_notebooks/IGP.md | 219 --------------------------- docs/sciserver_notebooks/LLC4320.md | 70 ++++++--- 6 files changed, 54 insertions(+), 471 deletions(-) delete mode 100644 docs/sciserver_notebooks/Fjord.md delete mode 100644 docs/sciserver_notebooks/IGP.md diff --git a/docs/_toc.yml b/docs/_toc.yml index 89f2e316..d89c1036 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -34,9 +34,10 @@ parts: - file: notebook/topology_tutorial - file: sciserver_example sections: - - file: sciserver_notebooks/IGP - - file: sciserver_notebooks/Fjord + - file: sciserver_notebooks/IGPwinter + - file: sciserver_notebooks/KangerFjord - file: sciserver_notebooks/LLC4320 + - file: sciserver_notebooks/ECCO_plot_stations - caption: Contributing chapters: - file: guide_for_developer diff --git a/docs/get_started.md b/docs/get_started.md index c4f39942..d050e7ae 100644 --- a/docs/get_started.md +++ b/docs/get_started.md @@ -9,7 +9,7 @@ If you have questions about any function, you can always go to [API references]( :gutter: 3 :::{grid-item-card} -:link: sciserver_notebooks/IGP +:link: sciserver_notebooks/IGPwinter :link-type: doc :class-header: bg-light @@ -34,7 +34,7 @@ ECCO dataset for now, but I want to have a 4320 example notebook. ::: :::{grid-item-card} -:link: sciserver_notebooks/Fjord +:link: sciserver_notebooks/KangerFjord :link-type: doc :class-header: bg-light diff --git a/docs/network_of_object.md b/docs/network_of_object.md index 9fdc8c61..5f540c34 100644 --- a/docs/network_of_object.md +++ b/docs/network_of_object.md @@ -7,6 +7,6 @@ Another good way to learn about the package is to learn how different objects ar | [`seaduck.OceData`](api_reference/apiref_OceData.rst) | Interface to data, translate between lat-lon and indices. | [AVISO](notebook/AVISO.ipynb) | | [`seaduck.Topology`](api_reference/apiref_topology.rst) | Describe the how the grids are connected. Similar to `xgcm`. | [topology tutorial](notebook/topology_tutorial.ipynb) | | [`seaduck.KnW`](api_reference/apiref_kernelNweight.rst) | Define what interpolation/derivative to perform. | [ECCO](notebook/global_ECCO.ipynb) | -| [`seaduck.Position`](api_reference/apiref_eulerian.rst) | Interpolate at Eulerian positions. | [Interpolate in a fjord](sciserver_notebooks/Fjord.md) | -| [`seaduck.Particle`](api_reference/apiref_lagrangian.rst) | Perform Lagrangian particle simulation. Sub class of `seaduck.Position` | [Regional simulation](sciserver_notebooks/IGP.md) | +| [`seaduck.Position`](api_reference/apiref_eulerian.rst) | Interpolate at Eulerian positions. | [Interpolate in a fjord](sciserver_notebooks/KangerFjord.md) | +| [`seaduck.Particle`](api_reference/apiref_lagrangian.rst) | Perform Lagrangian particle simulation. Sub class of `seaduck.Position` | [Regional simulation](sciserver_notebooks/IGPwinter.md) | | [`seaduck.OceInterp`](api_reference/apiref_OceInterp.rst) | Uniform interface to Lagrangian and Eulerian operations. | [one minute guide](one_min_guide.ipynb), [ECCO](notebook/global_ECCO.ipynb) | diff --git a/docs/sciserver_notebooks/Fjord.md b/docs/sciserver_notebooks/Fjord.md deleted file mode 100644 index b5db2750..00000000 --- a/docs/sciserver_notebooks/Fjord.md +++ /dev/null @@ -1,223 +0,0 @@ ---- -jupytext: - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.14.6 -kernelspec: - display_name: Oceanography - language: python - name: oceanography ---- - -# Demonstrate Position object with Fjord - -Author: Wenrui Jiang, 14 June 2023 -> **Warning**⚠️ : the notebook was last ran on **2023-06-15** with **seaduck 0.1.3.dev44+g8d60702**. You can find the executable version at https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/Fjord.ipynb. The `eulerian.Position` object is really what connect a point and the numerical model. Its `interpolate` method really is the core of this package. We are going to use a rather interesting example to demonstrate the functionalities of `eulerian.Position`. - -```{code-cell} ipython3 -import oceanspy as ospy -import numpy as np -import matplotlib.pyplot as plt -import matplotlib as mpl -import seaduck as sd -from functools import partial -import xarray as xr -import cmocean - -mpl.rcParams["figure.dpi"] = 300 -``` - -Here on sciserver, we have an interesting dataset simulating the interaction between background circulation and Kangerdlugssuaq Fjord. More information can be found below: - -```{code-cell} ipython3 ---- -jupyter: - outputs_hidden: true ---- -fjord = ospy.open_oceandataset.from_catalog("KangerFjord") -``` - -``` -Opening KangerFjord. -A realistic numerical model constructed to simulate the oceanic conditions -and circulation in a large southeast Greenland fjord (Kangerdlugssuaq) and -the adjacent shelf sea region during winter 2007–2008. -Citation: - * Fraser et al., 2018 - JGR. -``` - -+++ - -Let's first look at what the dataset looks like. We are going to use ETOPO dataset to give you an idea where this domain is located. - -```{code-cell} ipython3 -etopo = ospy.open_oceandataset.from_catalog("ETOPO") -etopo = etopo._ds.sel(X=slice(-40.6, -13.1), Y=slice(62.4, 70.5)).where( - etopo._ds.Depth > 0 -) -``` - -```{code-cell} ipython3 -:tags: [hide-input] - -vmax = 0.2 -cmap = cmocean.cm.balance -``` - -```{code-cell} ipython3 -:tags: [hide-input] - -c = plt.pcolormesh( - fjord._ds.XC, fjord._ds.YC, fjord._ds.Eta[263], vmax=vmax, vmin=-vmax, cmap=cmap -) -plt.pcolormesh( - etopo.X, - etopo.Y, - etopo.Depth, - cmap="bone", - vmin=-3000, - vmax=4500, - zorder=20, -) -plt.gca().set_facecolor("lightsteelblue") -plt.title("Model domain and the sea surface height") -plt.xlabel("Longitude") -plt.ylabel("Latitude") -plt.colorbar(c, label="m") -plt.show() -``` -![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/Fjord_files/Fjord_8_0.png?raw=true) - -+++ {"tags": ["mdformat-skip"]} - -## Doing Interpolation -We are going to use the sea surface height field $\eta$ as an example. - -First, we are going to convert the `xarray.Dataset` to `seaduck.OceData`. - -```{code-cell} ipython3 -tub = sd.OceData(fjord._ds) -``` - -For now, the interpolation points are simply the gridpoints of the dataset. - -```{code-cell} ipython3 -:tags: [hide-input] - -x = tub.XC.ravel() -y = tub.YC.ravel() -z = np.ones_like(x) * -5.0 -t = np.ones_like(x) * sd.utils.convert_time("2008-01-01") -``` - -Let's create the `eulerian.Position` object - -```{code-cell} ipython3 -p = sd.Position().from_latlon(x, y, z, t, data=tub) -``` - -Two kernels are defined here, both of them are the default kernel used by the package. However, we are going to "hack" one of them to demonstrate the "cascade" capacity of interpolation. - -```{code-cell} ipython3 -kernel = sd.KnW() -kernel_to_be_hacked = sd.KnW() -``` - -+++ {"tags": ["mdformat-skip"]} - -First, we do the normal interpolation on $\eta$ in the normal way, and plot it - -```{code-cell} ipython3 -eta = p.interpolate("Eta", kernel) -``` - -```{code-cell} ipython3 -:tags: [hide-input] - -plt.scatter(x, y, c=eta, s=0.7, vmax=vmax, vmin=-vmax, cmap=cmap) -plt.title("Interpolated sea surface height") -plt.colorbar(label="m") -plt.xlabel("Longitude") -plt.ylabel("Latitude") -plt.show() -``` -![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/Fjord_files/Fjord_19_0.png?raw=true) - -We are going to "hack" the code to make it return the size of the kernels used. - -The details of the hack is not very important. But if you are interested you can read the inline comments below. - -```{code-cell} ipython3 -:tags: [hide-input] - -# The partials of this function is going to replace all the weight functions. -# For legal weight functions, the sum of contribution (weight) should be one, -# but here every neighboring point has contribution of one -def dummy_weight(rx, ry, ans): - n = len(rx) - return np.ones((n, ans)) - - -# This is the step that overwrites the existing weight function. Highly not recommended. -kernel_to_be_hacked.hfuncs = [partial(dummy_weight, ans=len(k)) for k in kernel.kernels] -# Create a dataset of ones with the same shape as Eta. -tub["ones"] = xr.ones_like(tub["Eta"]) -``` - -```{code-cell} ipython3 -how_many = p.interpolate("ones", kernel_to_be_hacked) -``` - -Now, we can look at what the interpolation package does under the hood. - -```{code-cell} ipython3 -:tags: [hide-input] - -plt.scatter(x, y, c=how_many, s=0.7, cmap="Set1", vmax=9.5, vmin=0.5) -plt.title("The number of grid points used for interpolation") -plt.colorbar(label="Point count") -plt.xlabel("Longitude") -plt.ylabel("Latitude") -plt.show() -``` -![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/Fjord_files/Fjord_24_0.png?raw=true) - -## Filling between - -As you have seen earlier, the grid of this dataset is very uneven. - -Well, the strength of the `eulerian.Position` object is filling in information in between, so let's do that. - -```{code-cell} ipython3 -:tags: [hide-input] - -Nlon = 600 -Nlat = 100 -Ndep = 1 -ax, ay, az, at = sd.utils.easy_3d_cube( - (-34.5, -28.5, Nlon), (66.5, 67.0, Nlat), (-5.0, -5.0, Ndep), "2008-01-01" -) -``` - -Create `eulerian.Position` in between, and plot it! - -```{code-cell} ipython3 -addition = sd.Position().from_latlon(ax, ay, az, at, data=tub) -more_eta = addition.interpolate("Eta", kernel) -``` - -```{code-cell} ipython3 -:tags: [hide-input] - -plt.scatter(x, y, c=eta, s=0.7, vmax=vmax, vmin=-vmax, cmap=cmap) -plt.scatter(ax, ay, c=more_eta, s=0.7, vmax=vmax, vmin=-vmax, cmap=cmap) -plt.title("Interpolated sea surface height (Reinforced)") -plt.xlabel("Longitude") -plt.ylabel("Latitude") -plt.show() -``` -![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/Fjord_files/Fjord_29_0.png?raw=true) - -I'd say the filling-in is done pretty well! diff --git a/docs/sciserver_notebooks/IGP.md b/docs/sciserver_notebooks/IGP.md deleted file mode 100644 index 997af031..00000000 --- a/docs/sciserver_notebooks/IGP.md +++ /dev/null @@ -1,219 +0,0 @@ ---- -jupytext: - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.14.6 -kernelspec: - display_name: Oceanography - language: python - name: oceanography ---- - -# `Particle` in a East Greenland regional simulation - -Author: Wenrui Jiang, Tom Haine Feb '23 -> **Warning**⚠️ : the notebook was last ran on **2023-06-15** with **seaduck 0.1.3.dev44+g8d60702**. You can find the executable version at https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/IGP.ipynb. -```{code-cell} ipython3 -:tags: [hide-input] - -import matplotlib.pyplot as plt -import numpy as np -import oceanspy as ospy - -import seaduck as sd -``` - -```````{admonition} Access IGP -The regional MITgcm run is the IGPwinter simulationT and is publicly available on SciServer. The simulation output could opened using the OceanSpy package using the [`from_catalog`](https://oceanspy.readthedocs.io/en/latest/generated/oceanspy.open_oceandataset.from_catalog.html#oceanspy.open_oceandataset.from_catalog) method. - -`ecco = ospy.open_oceandataset.from_catalog("ECCO")._ds` - -or the daily-mean data ('daily_ecco'). - -`ecco = ospy.open_oceandataset.from_catalog('daily_ecco')._ds` - -Click [here](https://dev-poseidon-ocean.pantheonsite.io/products/datasets/) for a full list of the dataset hosted and [here](https://oceanspy.readthedocs.io/en/latest/datasets.html#igpwinterhttps://oceanspy.readthedocs.io/en/latest/datasets.html#igpwinter) to find out more. -``````{admonition} Access full ECCO dataset -The ECCO dataset is publicly available on SciServer. The simulation output could opened using the OceanSpy package using the [`from_catalog`](https://oceanspy.readthedocs.io/en/latest/generated/oceanspy.open_oceandataset.from_catalog.html#oceanspy.open_oceandataset.from_catalog) method. - -Choose between the monthly-mean data ('ECCO') - -`ecco = ospy.open_oceandataset.from_catalog("ECCO")._ds` - -or the daily-mean data ('daily_ecco'). - -`ecco = ospy.open_oceandataset.from_catalog('daily_ecco')._ds` - -Click [here](https://dev-poseidon-ocean.pantheonsite.io/products/datasets/) for a full list of the dataset hosted. -``````` - -```{code-cell} ipython3 -od = ospy.open_oceandataset.from_catalog("IGPwinter") -``` - -``` -Opening IGPwinter. -High-resolution numerical simulation carried out in parallel to the observational -component of the Iceland Greenland Seas Project (IGP). -Citation: - * Renfrew et al., 2019 - BAMS. -``` - -+++ - -We are going to artificially create an open boundary in depth. This is to demonstrate that there is no weird behavior associated with open vertical boundaries. - -```{code-cell} ipython3 -ds = od._ds.isel(Z=slice(0, 50), Zl=slice(0, 50)) -``` - -## Prepare the simulation - -Since this notebook works with an open domain, it is going to be a little different a few other notebooks. - -We are interested in looking at the coastal current. First, we need to prepare a `seaduck.OceData` object... - -```{code-cell} ipython3 -oce = sd.OceData(ds) -``` - -...Initialize the particles: We put the ducks on the East Greenland continental shelf... - -```{code-cell} ipython3 -Nx = 20 -Nz = 10 -x = np.linspace(-22, -20, Nx) -z = np.linspace(0, -200, Nz) -x, z = np.meshgrid(x, z) -x = x.ravel() -z = z.ravel() -y = np.ones_like(x) * 71.0 -``` - -...and at the beginning of the simulation... - -```{code-cell} ipython3 -start_time = "2018-01-02" -t = np.ones_like(x) * sd.utils.convert_time(start_time) -``` - -...and integrate forward in time for one month. - -```{code-cell} ipython3 -end_time = "2018-02-01" -tf = sd.utils.convert_time(end_time) -``` - -### Here is where the particles start on the map: - -```{code-cell} ipython3 -plt.pcolormesh(ds["XC"], ds["YC"], np.log10(ds["Depth"] + 10), cmap="Blues") -plt.plot(x, y, "r") -cb = plt.colorbar(label="Depth(m)") -cbar_depth = np.concatenate([np.arange(0, 500, 100), np.arange(500, 4000, 500)]) -cb.ax.set_yticks(np.log10(cbar_depth + 10), cbar_depth, fontsize=6) -plt.xlabel("Longitude") -plt.ylabel("Latitude") -plt.title("Bathymetry of model domain and particle initial position") -plt.show() -``` -![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/IGP_files/IGP_16_0.png?raw=true) - -**Fig.1** The initial position of the particles released (red line) on a horizontal map of the surrounding regions. - -+++ - -We are going to follow a cold puff of fresh water near the Greenland shelf. Since we have `OceanSpy` installed here at SciServer, let me quickly show you a vertical section of this region. (It really takes no effort at all!) - -```{code-cell} ipython3 -od_surv = od.subsample.survey_stations( - Xsurv=[-22.0, -19.0], Ysurv=[71.0, 71.0], delta=1 -) -od_surv._ds = od_surv._ds.isel(time=0) -od_surv.plot.vertical_section(varName="Temp", contourName="Sigma0") -plt.ylim([-750, 0]) -plt.show() -``` -![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/IGP_files/IGP_19_1.png?raw=true) - -**Fig.2** Vertical section at 71N. This vertical section goes a little further east than the particle position to include the shelf break. Colors denotes the potential temperature, while the contours is the potential density anomaly calculated using OceanSpy (with the equation of state the model used). - -+++ - -Since we are in an open domain (in all three dimensions), given long enough time, some particles will go out of the domain and jeopardise the entire simulation. We can define a callback function to stop that from happening. - -If provided, the function will be called every time particles cross walls. It can be used to manipulate the particles. But here we just use it as a stop/continue criterion. - -```{code-cell} ipython3 -def continue_criterion(pt): - x_ = np.logical_and(pt.lon < -10, pt.lon > -35) - y_ = np.logical_and(pt.lat < 72, pt.lat > 65) - z_ = pt.dep > -750 - return np.logical_and(np.logical_and(x_, y_), z_) -``` - -This time we are going to use volume flux (transport) to advect the particles. This is usually better than using the velocity field itself. - -```{code-cell} ipython3 -oce["utrans"] = oce["U"] * oce["drF"] * oce["dyG"] -oce["vtrans"] = oce["V"] * oce["drF"] * oce["dxG"] -oce["wtrans"] = oce["W"] * oce["rA"] -``` - -Finally, create the particle object. And that's all of the preparation we need. - -```{code-cell} ipython3 -p = sd.Particle( - x=x, - y=y, - z=z, - t=t, - data=oce, - callback=continue_criterion, - uname="utrans", - vname="vtrans", - wname="wtrans", - # save_raw = True, - transport=True, -) -``` - -## Perform the particle trajectory simulation - -```{code-cell} ipython3 -stops, raw = p.to_list_of_time([t[0], tf]) -``` - -Retrieve the particle positions from the `seaduck.eulerian.position` objects. - -```{code-cell} ipython3 -:tags: [hide-input] - -lons = np.array([pt.lon for pt in raw]) -lats = np.array([pt.lat for pt in raw]) -``` - -### Plot results - -```{code-cell} ipython3 -:tags: [hide-input] - -plt.pcolormesh(od._ds["XC"], od._ds["YC"], np.log10(od._ds["Depth"] + 10), cmap="Blues") -plt.plot(x, y, "r") -plt.plot(lons, lats, "gold", lw=0.5) -plt.xlim([-35, -10]) -plt.ylim([65, 72]) -plt.xlabel("Longitude") -plt.ylabel("Latitude") -plt.title("Particle trajectories overlaid on bathymetry map") -plt.show() -``` -![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/IGP_files/IGP_32_0.png?raw=true) - -**Fig.3** The particle trajectories overlaid on **Fig.1**. The coloring scheme of bathymetry is also the same as **Fig.1** - -```{code-cell} ipython3 - -``` diff --git a/docs/sciserver_notebooks/LLC4320.md b/docs/sciserver_notebooks/LLC4320.md index 026df8f3..b0dae343 100644 --- a/docs/sciserver_notebooks/LLC4320.md +++ b/docs/sciserver_notebooks/LLC4320.md @@ -4,18 +4,18 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.14.6 + jupytext_version: 1.14.7 kernelspec: display_name: Oceanography language: python name: oceanography --- -# Plot the streamlines on LLC4320 +# Plot surface streamlines on LLC4320 Author: Wenrui Jiang 14 June, 2023 -> **Warning**⚠️ : the notebook was last ran on **2023-06-15** with **seaduck 0.1.3.dev44+g8d60702**. You can find the executable version at https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/LLC4320.ipynb. -LLC4320 is a kilometer scale global simulation with complex grid topology. This is a good test for the performance of the package. +> **Warning**⚠️ : the notebook was last ran on **2023-07-06** with **seaduck 0.1.dev433+g99d4d3b**. You can find the executable version at https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/LLC4320.ipynb. +The LLC4320 ocean circulation model solution is a kilometer-scale, global simulation with complex grid topology. This is a good dataset to test the performance of the `seaduck` package. ```{code-cell} ipython3 import matplotlib.pyplot as plt @@ -31,9 +31,19 @@ import seaduck as sd mpl.rcParams["figure.dpi"] = 600 ``` +```{admonition} Access LLC4320 +The global MITgcm run is the LLC4320 simulation and is publicly available on [SciServer](sciserver.org) (from the Oceanography container). The simulation output can be opened using the [OceanSpy](https://github.com/hainegroup/oceanspyhttps://github.com/hainegroup/oceanspy) package using the [`from_catalog`](https://oceanspy.readthedocs.io/en/latest/generated/oceanspy.open_oceandataset.from_catalog.html#oceanspy.open_oceandataset.from_catalog) method. + +`od = ospy.open_oceandataset.from_catalog("LLC4320")` + +Click [here](https://dev-poseidon-ocean.pantheonsite.io/products/datasets/) for a full list of the SciServer datasets. +``` + ++++ + ## Calculate the streamlines online -First we get the dataset into `OceData` object with help from `OceanSpy`. This step require the grid to be loaded into memory as numpy arrays and creating the cKD tree on top of that, so it is going to take some time. +First we get the dataset into `OceData` object with help from `OceanSpy`. This step requires the grid to be loaded into memory as `numpy` arrays and then create the cKD tree on top of that, so it's going to take some time. ```{code-cell} ipython3 od = ospy.open_oceandataset.from_catalog("LLC4320") @@ -41,7 +51,7 @@ ds = od._ds oce = sd.OceData(ds) ``` -Initiate the particles randomly so that they are distributed evenly on the globe. +Initiate the particles randomly, so that they're distributed evenly on the globe. Use `N` = $1.5 \\times 10^5$ particles. ```{code-cell} ipython3 :tags: [hide-input] @@ -64,16 +74,16 @@ ax.plot(x, y, "o", markersize=0.1, transform=ccrs.PlateCarree()) ax.coastlines() plt.show() ``` -![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/LLC4320_files/LLC4320_5_0.png?raw=true) +![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/LLC4320_files/LLC4320_29_2.png?raw=true) -We will be using just one vertical layer of a single snapshot in the output to make the streamlines. +We'll use just one vertical layer (the first one, which is at the surface) of a single snapshot (the first one) to make the streamlines. ```{code-cell} ipython3 oce["u"] = oce["U"][0, 0] oce["v"] = oce["V"][0, 0] ``` -By setting `wname = None`, we create a particle moving in the 2D mode. +By setting `wname = None`, we fix particles to move on a horizontal surface. Think of isobaric floats. ```{code-cell} ipython3 %%time @@ -89,9 +99,9 @@ p = sd.Particle( ) ``` -The package knows that LLC4320 is a large dataset, so it is not prefetching data to not overwhelm itself. +The package knows that LLC4320 is a large dataset, so by default it doesn't prefetch data. Otherwise, `seaduck` is overwhelmed. In general, this is a smart thing to do. -In general, this is a smart thing to do. However, since we are using the 2D velocity, which is much smaller than all the component, we actually have enough memory for prefetching. We can manually change that by: +However, since we are using a static 2D velocity field here, which is much smaller than the 3D time-varying field, we actually have enough memory to prefetch. So tell `seaduck` this bit of LLC4320 isn't overwhelmingly big: ```{code-cell} ipython3 p.too_large = False @@ -102,7 +112,7 @@ p.uarray = np.array(oce["u"]) p.varray = np.array(oce["v"]) ``` -Everything went well. Now, we are going to run it with 3 hour timestep for 30 days. +Everything went well. Now, we are going to compute trajectories with a 3 hour timestep for 30 days. This will take a while: time for a nap... ```{code-cell} ipython3 step = 10800 @@ -117,7 +127,7 @@ stops, raw = p.to_list_of_time(dest, update_stops=[]) +++ {"tags": ["mdformat-skip"]} -It took around 3 hours to run a month of simulation for $1.5\times 10^{5}$ particles. We can now extract the things we need to plot, namely longitude, latitude, and horizontal speed. +It took around 1.5 hours to run a month of simulation for $1.5\times 10^{5}$ particles. We can now extract the things we need to plot, namely longitude, latitude, and horizontal speed. ```{code-cell} ipython3 lons = np.array([i.lon for i in raw]) @@ -143,19 +153,19 @@ spds = np.load('LLC4320spds.npy') ## Plotting Preparation -The ETOPO topography data set is also accessible on sciserver using `OceanSpy`. This will give us absolutely beautiful plots. +The [ETOPO](https://www.ncei.noaa.gov/products/etopo-global-relief-model) topography dataset is also accessible on SciServer using `OceanSpy`. This will give us absolutely beautiful plots. ```{code-cell} ipython3 etopo = ospy.open_oceandataset.from_catalog("ETOPO") ``` -For this notebook, we only need the over sea level part, because there is already so much going on in the water. +For this notebook, we only need the land, because there is already so much going on in the water. So clip the `etopo` data. ```{code-cell} ipython3 etopo = etopo._ds.where(etopo._ds.Depth > 0) ``` -Here is the function for plotting. The lines will have use speed as colormaps. This part is going to take some time, but I hope you find it worth the wait. +Here is the function for plotting. The lines will be colored with speed. This part is going to take some time, but I hope you think it's worth the wait. ```{code-cell} ipython3 def pretty_stream_plot( @@ -225,26 +235,40 @@ def pretty_stream_plot( ## Voila! -### Looking from Arctic +### Looking from the Arctic ```{code-cell} ipython3 north_projection = ccrs.NorthPolarStereo(central_longitude=38.0) -pretty_stream_plot(lats, lons, spds, north_projection, south=6.5) +pretty_stream_plot( + lats, + lons, + spds, + north_projection, + south=6.5, + save_as="LLC4320_files/LLC4320_29_2.png", +) plt.show() ``` -![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/LLC4320_files/LLC4320_29_2.png?raw=true) +![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/LLC4320_files/LLC4320_32_2.png?raw=true) -**Fig.1** Streamlines of 150,000 particles released at 00:00 Apr 25 th, 2012 in LLC4320 simulated for 30 days. The color shading shows the current speed. This is looking from the North Pole. +**Fig.1** Streamlines of 150,000 particles released at 00:00 Apr 25, 2012 in LLC4320 simulated for 30 days. The color shading shows the current speed. This is looking from the North Pole. +++ -### Looking from Antarctica +### Looking from the Antarctic ```{code-cell} ipython3 south_projection = ccrs.SouthPolarStereo(central_longitude=38.0) -pretty_stream_plot(lats, lons, spds, south_projection, north=-6.5) +pretty_stream_plot( + lats, + lons, + spds, + south_projection, + north=-6.5, + save_as="LLC4320_files/LLC4320_32_2.png", +) plt.show() ``` -![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/LLC4320_files/LLC4320_32_2.png?raw=true) +![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/LLC4320_files/LLC4320_5_0.png?raw=true) **Fig.2** Similar to **Fig.1**, but looking from the South Pole. From a885719a03449ca546a3d1ae3bef009570196ee7 Mon Sep 17 00:00:00 2001 From: ThomasHaine Date: Thu, 6 Jul 2023 13:47:56 -0400 Subject: [PATCH 2/3] Add ECCO_plot_stations.md and rename --- .../sciserver_notebooks/ECCO_plot_stations.md | 162 +++++++++++++ docs/sciserver_notebooks/IGPwinter.md | 205 ++++++++++++++++ docs/sciserver_notebooks/KangerFjord.md | 221 ++++++++++++++++++ 3 files changed, 588 insertions(+) create mode 100644 docs/sciserver_notebooks/ECCO_plot_stations.md create mode 100644 docs/sciserver_notebooks/IGPwinter.md create mode 100644 docs/sciserver_notebooks/KangerFjord.md diff --git a/docs/sciserver_notebooks/ECCO_plot_stations.md b/docs/sciserver_notebooks/ECCO_plot_stations.md new file mode 100644 index 00000000..fb0015ab --- /dev/null +++ b/docs/sciserver_notebooks/ECCO_plot_stations.md @@ -0,0 +1,162 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.7 +kernelspec: + display_name: Oceanography + language: python + name: oceanography +--- + +# `OceInterp` interpolation in the ECCO global state estimate. + +Extract temperature/salinity profiles at specified longitudes, latitudes, and times. This process mimics hydrographic stations from in-situ measurements (like a ship CTD station, or Argo float profile). + +This notebook uses [Oceanspy](https://oceanspy.readthedocs.io/en/latest/) and demonstrates the interface to the Poseidon-viewer on SciServer. + +Author: Tom Haine & Wenrui Jiang, Jun '23 +> **Warning**⚠️ : the notebook was last ran on **2023-07-06** with **seaduck 0.1.dev433+g99d4d3b**. You can find the executable version at https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/ECCO_plot_stations.ipynb. +```{code-cell} ipython3 +import seaduck as sd +import oceanspy as ospy +from poseidon_viewer import get_shapes + +import numpy as np + +import matplotlib.pyplot as plt +plt.rcParams['figure.figsize'] = 12, 8 +``` + +```{admonition} Access ECCO +The global MITgcm run is the ECCO state estimate and is publicly available on [SciServer](sciserver.org) (from the Oceanography container). The simulation output can be opened using the [OceanSpy](https://github.com/hainegroup/oceanspyhttps://github.com/hainegroup/oceanspy) package using the [`from_catalog`](https://oceanspy.readthedocs.io/en/latest/generated/oceanspy.open_oceandataset.from_catalog.html#oceanspy.open_oceandataset.from_catalog) method. + +`ecco = ospy.open_oceandataset.from_catalog("ECCO")` + +Click [here](https://dev-poseidon-ocean.pantheonsite.io/products/datasets/) for a full list of the dataset hosted and [here](https://oceanspy.readthedocs.io/en/latest/datasets.html#igpwinterhttps://oceanspy.readthedocs.io/en/latest/datasets.html#igpwinter) to find out more. +``` + +```{code-cell} ipython3 +od = ospy.open_oceandataset.from_catalog('ECCO') + +od._ds = od._ds.drop_vars({'k', 'k_u', 'k_p1', 'k_l'}) +co_list = [var for var in od._ds.variables if "time" not in od._ds[var].dims] +od._ds = od._ds.set_coords(co_list) + +od._ds = od._ds.drop_vars('time_bnds') +od._ds['Temp'] = od._ds['THETA'] +od._ds['S'] = od._ds['SALT'] +this_time = '1992-03-16T12:00:00.000000000' # Select time of interest. This is cut and pasted (with editing to convert format) from Poseidon-viewer, but in future should be returned in the JSON object. + +varList = ['Temp', 'S'] # Select variables of interest. +Nvars = len(varList) +od._ds = od._ds.drop_vars([var for var in od._ds.data_vars if var not in varList]) +``` + +Define coordinates for interpolation from the Poseidon-viewer (or enter them manually): + +```{code-cell} ipython3 +ps = [{"type":"Point","coordinates":[-57.710530333876235,33.062070736515295]},{"type":"Point","coordinates":[-89.55599488497559,25.453020491528747]},{"type":"Point","coordinates":[-74.70605786848537,13.952759739624767]},{"type":"Point","coordinates":[-59.06355926811892,21.999755420813273]},{"type":"Point","coordinates":[-41.69588507736159,22.43955659234939]},{"type":"Point","coordinates":[-45.113026771172656,14.563543828761837]},{"type":"Point","coordinates":[-26.081872200732885,6.414099524482438]},{"type":"Point","coordinates":[-17.396102656758963,-4.381322875209875]},{"type":"Point","coordinates":[-26.702603318403906,-7.125636489486197]},{"type":"Point","coordinates":[-32.51011235240231,-22.847802807885373]}] +lons, lats = ospy.utils.viewer_to_range(ps) +Nstations = len(lons) +``` + +### `seaduck.OceInterp` interpolates the temperature and salinity to the specified coordinates. +This process makes synthetic hydrographic profiles from ECCO. Compute the potential density anomaly from the T/S data too. + +```{code-cell} ipython3 +tub = sd.OceData(od._ds) +times = sd.utils.convert_time(this_time) +Ntimes = times.size + +depths = tub.Z.ravel() +Ndepths = len(depths) + +sd_lons,sd_times,sd_depths = np.meshgrid(lons,times,depths,indexing='ij') +sd_lats,_,_ = np.meshgrid(lats,times,depths,indexing='ij') +sd_lons = np.ravel(sd_lons) +sd_lats = np.ravel(sd_lats) +sd_times = np.ravel(sd_times) +sd_depths = np.ravel(sd_depths) + +[sd_Temp,sd_S] = np.array(sd.OceInterp(tub,varList, sd_lons, sd_lats, sd_depths, sd_times)) +sd_depths2 = sd_depths.reshape(Nstations,Ntimes,Ndepths).squeeze() +sd_S = sd_S.reshape(Nstations,Ntimes,Ndepths).squeeze() +sd_Temp = sd_Temp.reshape(Nstations,Ntimes,Ndepths).squeeze() + +sd_sigma0 = ospy.utils.densjmd95(sd_S,sd_Temp,0) - 1000.0 +``` + +### Plot station locations: + +```{code-cell} ipython3 +cut_od = od.subsample.cutout(XRange = [min(lons)-10,max(lons)+10],YRange = [min(lats)-5,max(lats)+5]) +fig = plt.figure(figsize=(12, 5)) +ax = cut_od.plot.horizontal_section(varName="Depth", cmap='Greys_r') +colors=['r', 'b', 'k', 'orange', 'darkgreen', 'm', 'yellow', 'c', 'indigo', 'magenta', 'green'] +legestations_ds = np.arange(Nstations) +for i in range(Nstations): + plt.plot(lons[i], lats[i], 'o', color=colors[i], label='station ' + str(legestations_ds[i])) +plt.legend() +plt.show() +``` +![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/ECCO_plot_stations_files/ECCO_plot_stations_11_0.png?raw=true) + +### Plot vertical T/S station profiles: + +```{code-cell} ipython3 +for i in range(Nstations): + + plt.subplot(1,3,1) + plt.plot(sd_Temp[i,:],depths,'o', color=colors[i], label='station ' + str(legestations_ds[i])) + plt.xlabel('Temperature [oC]') + plt.ylabel('Depth [m]') + plt.grid() + + plt.subplot(1,3,2) + plt.plot(sd_S[i,:],depths,'o', color=colors[i], label='station ' + str(legestations_ds[i])) + plt.xlabel('Salinity [g/kg]') + plt.grid() + + + plt.subplot(1,3,3) + plt.plot(sd_sigma0[i,:],depths,'o', color=colors[i], label='station ' + str(legestations_ds[i])) + plt.xlabel('Density anomaly [kg/m^3]') + plt.grid() + +plt.subplot(1,3,1) +plt.xlabel('Temperature [oC]') +plt.ylabel('Depth [m]') +plt.grid() + +plt.subplot(1,3,2) +plt.xlabel('Salinity [g/kg]') +plt.grid() + +plt.subplot(1,3,3) +plt.xlabel('Density anomaly [kg/m^3]') +plt.grid() +plt.legend() + +plt.show() +``` +![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/ECCO_plot_stations_files/ECCO_plot_stations_13_0.png?raw=true) + +### Plot T/S diagram. + +```{code-cell} ipython3 +for i in range(Nstations): + plt.plot(sd_S[i,:],sd_Temp[i,:],'o', color=colors[i], label='station ' + str(legestations_ds[i])) +plt.xlabel('Salinity [g/kg]') +plt.ylabel('Temperature [oC]') +plt.legend() +plt.grid() +plt.show() +``` +![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/ECCO_plot_stations_files/ECCO_plot_stations_9_1.png?raw=true) + +```{code-cell} ipython3 + +``` diff --git a/docs/sciserver_notebooks/IGPwinter.md b/docs/sciserver_notebooks/IGPwinter.md new file mode 100644 index 00000000..b254de77 --- /dev/null +++ b/docs/sciserver_notebooks/IGPwinter.md @@ -0,0 +1,205 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.7 +kernelspec: + display_name: Oceanography + language: python + name: oceanography +--- + +# `Particle` in an East Greenland regional simulation + +Author: Wenrui Jiang, Tom Haine Feb '23 +> **Warning**⚠️ : the notebook was last ran on **2023-07-06** with **seaduck 0.1.dev433+g99d4d3b**. You can find the executable version at https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/IGPwinter.ipynb. +```{code-cell} ipython3 +:tags: [hide-input] + +import matplotlib.pyplot as plt +import numpy as np +import oceanspy as ospy + +import seaduck as sd +``` + +```{admonition} Access IGPwinter +The regional MITgcm run is the IGPwinter simulation and is publicly available on [SciServer](sciserver.org) (from the Oceanography container). The simulation output can be opened using the [OceanSpy](https://github.com/hainegroup/oceanspyhttps://github.com/hainegroup/oceanspy) package using the [`from_catalog`](https://oceanspy.readthedocs.io/en/latest/generated/oceanspy.open_oceandataset.from_catalog.html#oceanspy.open_oceandataset.from_catalog) method. + +`ecco = ospy.open_oceandataset.from_catalog("IGPwinter")` + +Click [here](https://dev-poseidon-ocean.pantheonsite.io/products/datasets/) for a full list of the dataset hosted and [here](https://oceanspy.readthedocs.io/en/latest/datasets.html#igpwinterhttps://oceanspy.readthedocs.io/en/latest/datasets.html#igpwinter) to find out more. +``` + +```{code-cell} ipython3 +od = ospy.open_oceandataset.from_catalog("IGPwinter") +``` + +``` +Opening IGPwinter. +High-resolution numerical simulation carried out in parallel to the observational +component of the Iceland Greenland Seas Project (IGP). +Citation: + * Renfrew et al., 2019 - BAMS. +``` + ++++ + +We're going to artificially create an open boundary in depth. This is to demonstrate that there's no weird behavior associated with open vertical boundaries. + +```{code-cell} ipython3 +ds = od._ds.isel(Z=slice(0, 50), Zl=slice(0, 50)) +``` + +## Prepare the simulation + +Since this notebook works with an open domain (regional model, not global), it's going to differ a bit from other `seaduck` notebooks. + +We're interested in looking at the coastal current. First, we need to prepare a `seaduck.OceData` object... + +```{code-cell} ipython3 +oce = sd.OceData(ds) +``` + +...Initialize the particles: We put the ducks on the East Greenland continental shelf... + +```{code-cell} ipython3 +Nx = 20 +Nz = 10 +x = np.linspace(-22, -20, Nx) +z = np.linspace(0, -200, Nz) +x, z = np.meshgrid(x, z) +x = x.ravel() +z = z.ravel() +y = np.ones_like(x) * 71.0 +``` + +...and at the beginning of the simulation... + +```{code-cell} ipython3 +start_time = "2018-01-02" +t = np.ones_like(x) * sd.utils.convert_time(start_time) +``` + +...and integrate forward in time for one month (the actual trajectory calculation is made below; this is just setting parameters). + +```{code-cell} ipython3 +end_time = "2018-02-01" +tf = sd.utils.convert_time(end_time) +``` + +### Here is where the particles start on the map: + +```{code-cell} ipython3 +plt.pcolormesh(ds["XC"], ds["YC"], np.log10(ds["Depth"] + 10), cmap="Blues") +plt.plot(x, y, "r") +cb = plt.colorbar(label="Depth(m)") +cbar_depth = np.concatenate([np.arange(0, 500, 100), np.arange(500, 4000, 500)]) +cb.ax.set_yticks(np.log10(cbar_depth + 10), cbar_depth, fontsize=6) +plt.xlabel("Longitude") +plt.ylabel("Latitude") +plt.title("Bathymetry of model domain and particle initial position") +plt.show() +``` +![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/IGPwinter_files/IGPwinter_16_0.png?raw=true) + +**Fig.1** The initial position of the particles released (red line) on a horizontal map of the surrounding regions. + ++++ + +We are going to follow a cold puff of fresh water near the Greenland shelf. Since we have `OceanSpy` installed here at SciServer, let me quickly show you how to make a vertical section of this region. (It really takes no effort at all!) + +```{code-cell} ipython3 +od_surv = od.subsample.survey_stations( + Xsurv=[-22.0, -19.0], Ysurv=[71.0, 71.0], delta=1 +) +od_surv._ds = od_surv._ds.isel(time=0) +od_surv.plot.vertical_section(varName="Temp", contourName="Sigma0") +plt.ylim([-750, 0]) +plt.show() +``` +![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/IGPwinter_files/IGPwinter_19_1.png?raw=true) + +**Fig.2** Vertical section at 71N. This vertical section goes a little further east than the initial particle positions to include the shelf break. Colors denotes the potential temperature, while the contours are the potential density anomaly calculated using OceanSpy (with the equation of state the model used). + ++++ + +Since we are in an open domain (in all three dimensions), given long enough time, some particles will leave the domain. This will jeopardise the entire simulation! We can define a callback function to stop that from happening. + +If provided, the function will be called every time particles cross walls. It can be used to manipulate the particles and catch the out-of-domain issue. But here we just use it as a stop/continue criterion. + +```{code-cell} ipython3 +def continue_criterion(pt): + x_ = np.logical_and(pt.lon < -10, pt.lon > -35) + y_ = np.logical_and(pt.lat < 72, pt.lat > 65) + z_ = pt.dep > -750 + return np.logical_and(np.logical_and(x_, y_), z_) +``` + +This time we are going to use volume flux (transport) to advect the particles. This is usually better than using the velocity field itself. + +```{code-cell} ipython3 +oce["utrans"] = oce["U"] * oce["drF"] * oce["dyG"] +oce["vtrans"] = oce["V"] * oce["drF"] * oce["dxG"] +oce["wtrans"] = oce["W"] * oce["rA"] +``` + +Finally, create the particle object: + +```{code-cell} ipython3 +p = sd.Particle( + x=x, + y=y, + z=z, + t=t, + data=oce, + callback=continue_criterion, + uname="utrans", + vname="vtrans", + wname="wtrans", + # save_raw = True, + transport=True, +) +``` + +## Perform the particle trajectory simulation + +Run the simulation! This can take some time, so grab a break... + +```{code-cell} ipython3 +stops, raw = p.to_list_of_time([t[0], tf]) +``` + +Retrieve the particle positions from the `seaduck.eulerian.position` objects. + +```{code-cell} ipython3 +:tags: [hide-input] + +lons = np.array([pt.lon for pt in raw]) +lats = np.array([pt.lat for pt in raw]) +``` + +### Plot results + +```{code-cell} ipython3 +:tags: [hide-input] + +plt.pcolormesh(od._ds["XC"], od._ds["YC"], np.log10(od._ds["Depth"] + 10), cmap="Blues") +plt.plot(x, y, "r") +plt.plot(lons, lats, "gold", lw=0.5) +plt.xlim([-35, -10]) +plt.ylim([65, 72]) +plt.xlabel("Longitude") +plt.ylabel("Latitude") +plt.title("Particle trajectories overlaid on bathymetry map") +plt.show() +``` +![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/IGPwinter_files/IGPwinter_32_0.png?raw=true) + +**Fig.3** The particle trajectories overlaid on **Fig.1**. The color scheme for bathymetry is the same as **Fig.1**. + +```{code-cell} ipython3 + +``` diff --git a/docs/sciserver_notebooks/KangerFjord.md b/docs/sciserver_notebooks/KangerFjord.md new file mode 100644 index 00000000..d74c7fdf --- /dev/null +++ b/docs/sciserver_notebooks/KangerFjord.md @@ -0,0 +1,221 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.7 +kernelspec: + display_name: Oceanography + language: python + name: oceanography +--- + +# Demonstrate `eulerian.Position` object with Fjord + +Author: Wenrui Jiang, 14 June 2023 +> **Warning**⚠️ : the notebook was last ran on **2023-07-06** with **seaduck 0.1.dev433+g99d4d3b**. You can find the executable version at https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/KangerFjord.ipynb. +The `eulerian.Position` object is really what connects a point and the numerical model. Its `interpolate` method really is the core of this package. We're going to use a rather interesting example to demonstrate the functionalities of `eulerian.Position`. + +```{code-cell} ipython3 +import oceanspy as ospy +import numpy as np +import matplotlib.pyplot as plt +import matplotlib as mpl +import seaduck as sd +from functools import partial +import xarray as xr +import cmocean + +mpl.rcParams["figure.dpi"] = 300 +``` + +Here on SciServer, we have an interesting dataset simulating the interaction between background circulation and Kangerdlugssuaq Fjord. More information can be found below, and see the paper by [Fraser et al., 2018](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018JC014435https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018JC014435): + +```{code-cell} ipython3 +fjord = ospy.open_oceandataset.from_catalog("KangerFjord") +``` + +``` +Opening KangerFjord. +A realistic numerical model constructed to simulate the oceanic conditions +and circulation in a large southeast Greenland fjord (Kangerdlugssuaq) and +the adjacent shelf sea region during winter 2007–2008. +Citation: + * Fraser et al., 2018 - JGR. +``` + ++++ + +Let's first explore the dataset a bit. We are going to use [ETOPO](https://www.ncei.noaa.gov/products/etopo-global-relief-modelhttps://www.ncei.noaa.gov/products/etopo-global-relief-model) dataset to give you an idea where this domain is located. + +```{code-cell} ipython3 +etopo = ospy.open_oceandataset.from_catalog("ETOPO") +etopo = etopo._ds.sel(X=slice(-40.6, -13.1), Y=slice(62.4, 70.5)).where( + etopo._ds.Depth > 0 +) +``` + +```{code-cell} ipython3 +:tags: [hide-input] + +vmax = 0.2 +cmap = cmocean.cm.balance +``` + +```{code-cell} ipython3 +:tags: [hide-input] + +c = plt.pcolormesh( + fjord._ds.XC, fjord._ds.YC, fjord._ds.Eta[263], vmax=vmax, vmin=-vmax, cmap=cmap +) +plt.pcolormesh( + etopo.X, + etopo.Y, + etopo.Depth, + cmap="bone", + vmin=-3000, + vmax=4500, + zorder=20, +) +plt.gca().set_facecolor("lightsteelblue") +plt.title("Model domain and the sea surface height") +plt.xlabel("Longitude") +plt.ylabel("Latitude") +plt.colorbar(c, label="m") +plt.show() +``` +![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/KangerFjord_files/Fjord_19_0.png?raw=true) + ++++ {"tags": ["mdformat-skip"]} + +## Doing Interpolation +We are going to use the sea surface height field $\eta$ as an example. + +First, we are going to convert the `xarray.Dataset` to `seaduck.OceData`. + +```{code-cell} ipython3 +tub = sd.OceData(fjord._ds) +``` + +For now, the interpolation points are simply the gridpoints of the dataset. + +```{code-cell} ipython3 +:tags: [hide-input] + +x = tub.XC.ravel() +y = tub.YC.ravel() +z = np.ones_like(x) * -5.0 +t = np.ones_like(x) * sd.utils.convert_time("2008-01-01") +``` + +Let's create the `eulerian.Position` object + +```{code-cell} ipython3 +p = sd.Position().from_latlon(x, y, z, t, data=tub) +``` + +Two interpolation kernels are defined here, both of them are default kernels used by the package. However, we are going to "hack" one of them to demonstrate the "cascade" capacity of interpolation. + +```{code-cell} ipython3 +kernel = sd.KnW() +kernel_to_be_hacked = sd.KnW() +``` + ++++ {"tags": ["mdformat-skip"]} + +First, we do the normal interpolation on $\eta$ in the normal way, and plot it: + +```{code-cell} ipython3 +eta = p.interpolate("Eta", kernel) +``` + +```{code-cell} ipython3 +:tags: [hide-input] + +plt.scatter(x, y, c=eta, s=0.7, vmax=vmax, vmin=-vmax, cmap=cmap) +plt.title("Interpolated sea surface height") +plt.colorbar(label="m") +plt.xlabel("Longitude") +plt.ylabel("Latitude") +plt.show() +``` +![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/KangerFjord_files/Fjord_24_0.png?raw=true) + +Now let's "hack" the code to make it return the size of the kernels used. + +The details of the hack is not very important. But if you're interested you can read the inline comments below. + +```{code-cell} ipython3 +:tags: [hide-input] + +# The partials of this function is going to replace all the weight functions. +# For legal weight functions, the sum of contribution (weight) should be one, +# but here every neighboring point has contribution of one +def dummy_weight(rx, ry, ans): + n = len(rx) + return np.ones((n, ans)) + + +# This is the step that overwrites the existing weight function. Highly not recommended. +kernel_to_be_hacked.hfuncs = [partial(dummy_weight, ans=len(k)) for k in kernel.kernels] +# Create a dataset of ones with the same shape as Eta. +tub["ones"] = xr.ones_like(tub["Eta"]) +``` + +```{code-cell} ipython3 +how_many = p.interpolate("ones", kernel_to_be_hacked) +``` + +Now, we can look at what the interpolation package does under the hood. + +```{code-cell} ipython3 +:tags: [hide-input] + +plt.scatter(x, y, c=how_many, s=0.7, cmap="Set1", vmax=9.5, vmin=0.5) +plt.title("The number of grid points used for interpolation") +plt.colorbar(label="Point count") +plt.xlabel("Longitude") +plt.ylabel("Latitude") +plt.show() +``` +![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/KangerFjord_files/Fjord_29_0.png?raw=true) + +## Filling between + +As you have seen earlier, the grid of this dataset has very uneven spacing. + +Well, the strength of the `eulerian.Position` object is filling in information (interpolation) between the grid points. +So let's do that. We specify part of the model domain to fill in using the parameters to the `utils.easy_3d_cube` call. + +```{code-cell} ipython3 +:tags: [hide-input] + +Nlon = 600 +Nlat = 100 +Ndep = 1 +ax, ay, az, at = sd.utils.easy_3d_cube( + (-34.5, -28.5, Nlon), (66.5, 67.0, Nlat), (-5.0, -5.0, Ndep), "2008-01-01" +) +``` + +Create `eulerian.Position` in between, and plot it! + +```{code-cell} ipython3 +addition = sd.Position().from_latlon(ax, ay, az, at, data=tub) +more_eta = addition.interpolate("Eta", kernel) +``` + +```{code-cell} ipython3 +:tags: [hide-input] + +plt.scatter(x, y, c=eta, s=0.7, vmax=vmax, vmin=-vmax, cmap=cmap) +plt.scatter(ax, ay, c=more_eta, s=0.7, vmax=vmax, vmin=-vmax, cmap=cmap) +plt.title("Interpolated sea surface height (Reinforced)") +plt.xlabel("Longitude") +plt.ylabel("Latitude") +plt.show() +``` +![png](https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/KangerFjord_files/Fjord_8_0.png?raw=true) + +I'd say the filling-in is done pretty well! From 4f508bca8b81edd5a26e7b8d374784ddbf1d16a8 Mon Sep 17 00:00:00 2001 From: ThomasHaine Date: Thu, 6 Jul 2023 13:59:45 -0400 Subject: [PATCH 3/3] Run make --- docs/sciserver_notebooks/ECCO_plot_stations.md | 12 ++++++------ docs/sciserver_notebooks/IGPwinter.md | 4 ++-- docs/sciserver_notebooks/KangerFjord.md | 4 ++-- docs/sciserver_notebooks/LLC4320.md | 4 ++-- 4 files changed, 12 insertions(+), 12 deletions(-) diff --git a/docs/sciserver_notebooks/ECCO_plot_stations.md b/docs/sciserver_notebooks/ECCO_plot_stations.md index fb0015ab..a7c46a05 100644 --- a/docs/sciserver_notebooks/ECCO_plot_stations.md +++ b/docs/sciserver_notebooks/ECCO_plot_stations.md @@ -18,7 +18,7 @@ Extract temperature/salinity profiles at specified longitudes, latitudes, and ti This notebook uses [Oceanspy](https://oceanspy.readthedocs.io/en/latest/) and demonstrates the interface to the Poseidon-viewer on SciServer. Author: Tom Haine & Wenrui Jiang, Jun '23 -> **Warning**⚠️ : the notebook was last ran on **2023-07-06** with **seaduck 0.1.dev433+g99d4d3b**. You can find the executable version at https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/ECCO_plot_stations.ipynb. +> **Warning**⚠️ : the notebook was last ran on **2023-07-06** with **seaduck 0.1.dev433+g99d4d3b**. You can find the executable version at https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/ECCO_plot_stations.ipynb. ```{code-cell} ipython3 import seaduck as sd import oceanspy as ospy @@ -31,7 +31,7 @@ plt.rcParams['figure.figsize'] = 12, 8 ``` ```{admonition} Access ECCO -The global MITgcm run is the ECCO state estimate and is publicly available on [SciServer](sciserver.org) (from the Oceanography container). The simulation output can be opened using the [OceanSpy](https://github.com/hainegroup/oceanspyhttps://github.com/hainegroup/oceanspy) package using the [`from_catalog`](https://oceanspy.readthedocs.io/en/latest/generated/oceanspy.open_oceandataset.from_catalog.html#oceanspy.open_oceandataset.from_catalog) method. +The global MITgcm run is the ECCO state estimate and is publicly available on [SciServer](sciserver.org) (from the Oceanography container). The simulation output can be opened using the [OceanSpy](https://github.com/hainegroup/oceanspyhttps://github.com/hainegroup/oceanspy) package using the [`from_catalog`](https://oceanspy.readthedocs.io/en/latest/generated/oceanspy.open_oceandataset.from_catalog.html#oceanspy.open_oceandataset.from_catalog) method. `ecco = ospy.open_oceandataset.from_catalog("ECCO")` @@ -49,7 +49,7 @@ od._ds = od._ds.drop_vars('time_bnds') od._ds['Temp'] = od._ds['THETA'] od._ds['S'] = od._ds['SALT'] this_time = '1992-03-16T12:00:00.000000000' # Select time of interest. This is cut and pasted (with editing to convert format) from Poseidon-viewer, but in future should be returned in the JSON object. - + varList = ['Temp', 'S'] # Select variables of interest. Nvars = len(varList) od._ds = od._ds.drop_vars([var for var in od._ds.data_vars if var not in varList]) @@ -63,7 +63,7 @@ lons, lats = ospy.utils.viewer_to_range(ps) Nstations = len(lons) ``` -### `seaduck.OceInterp` interpolates the temperature and salinity to the specified coordinates. +### `seaduck.OceInterp` interpolates the temperature and salinity to the specified coordinates. This process makes synthetic hydrographic profiles from ECCO. Compute the potential density anomaly from the T/S data too. ```{code-cell} ipython3 @@ -122,10 +122,10 @@ for i in range(Nstations): plt.subplot(1,3,3) - plt.plot(sd_sigma0[i,:],depths,'o', color=colors[i], label='station ' + str(legestations_ds[i])) + plt.plot(sd_sigma0[i,:],depths,'o', color=colors[i], label='station ' + str(legestations_ds[i])) plt.xlabel('Density anomaly [kg/m^3]') plt.grid() - + plt.subplot(1,3,1) plt.xlabel('Temperature [oC]') plt.ylabel('Depth [m]') diff --git a/docs/sciserver_notebooks/IGPwinter.md b/docs/sciserver_notebooks/IGPwinter.md index b254de77..ffa611a1 100644 --- a/docs/sciserver_notebooks/IGPwinter.md +++ b/docs/sciserver_notebooks/IGPwinter.md @@ -14,7 +14,7 @@ kernelspec: # `Particle` in an East Greenland regional simulation Author: Wenrui Jiang, Tom Haine Feb '23 -> **Warning**⚠️ : the notebook was last ran on **2023-07-06** with **seaduck 0.1.dev433+g99d4d3b**. You can find the executable version at https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/IGPwinter.ipynb. +> **Warning**⚠️ : the notebook was last ran on **2023-07-06** with **seaduck 0.1.dev433+g99d4d3b**. You can find the executable version at https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/IGPwinter.ipynb. ```{code-cell} ipython3 :tags: [hide-input] @@ -26,7 +26,7 @@ import seaduck as sd ``` ```{admonition} Access IGPwinter -The regional MITgcm run is the IGPwinter simulation and is publicly available on [SciServer](sciserver.org) (from the Oceanography container). The simulation output can be opened using the [OceanSpy](https://github.com/hainegroup/oceanspyhttps://github.com/hainegroup/oceanspy) package using the [`from_catalog`](https://oceanspy.readthedocs.io/en/latest/generated/oceanspy.open_oceandataset.from_catalog.html#oceanspy.open_oceandataset.from_catalog) method. +The regional MITgcm run is the IGPwinter simulation and is publicly available on [SciServer](sciserver.org) (from the Oceanography container). The simulation output can be opened using the [OceanSpy](https://github.com/hainegroup/oceanspyhttps://github.com/hainegroup/oceanspy) package using the [`from_catalog`](https://oceanspy.readthedocs.io/en/latest/generated/oceanspy.open_oceandataset.from_catalog.html#oceanspy.open_oceandataset.from_catalog) method. `ecco = ospy.open_oceandataset.from_catalog("IGPwinter")` diff --git a/docs/sciserver_notebooks/KangerFjord.md b/docs/sciserver_notebooks/KangerFjord.md index d74c7fdf..b8e1d460 100644 --- a/docs/sciserver_notebooks/KangerFjord.md +++ b/docs/sciserver_notebooks/KangerFjord.md @@ -14,7 +14,7 @@ kernelspec: # Demonstrate `eulerian.Position` object with Fjord Author: Wenrui Jiang, 14 June 2023 -> **Warning**⚠️ : the notebook was last ran on **2023-07-06** with **seaduck 0.1.dev433+g99d4d3b**. You can find the executable version at https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/KangerFjord.ipynb. +> **Warning**⚠️ : the notebook was last ran on **2023-07-06** with **seaduck 0.1.dev433+g99d4d3b**. You can find the executable version at https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/KangerFjord.ipynb. The `eulerian.Position` object is really what connects a point and the numerical model. Its `interpolate` method really is the core of this package. We're going to use a rather interesting example to demonstrate the functionalities of `eulerian.Position`. ```{code-cell} ipython3 @@ -90,7 +90,7 @@ plt.show() +++ {"tags": ["mdformat-skip"]} ## Doing Interpolation -We are going to use the sea surface height field $\eta$ as an example. +We are going to use the sea surface height field $\eta$ as an example. First, we are going to convert the `xarray.Dataset` to `seaduck.OceData`. diff --git a/docs/sciserver_notebooks/LLC4320.md b/docs/sciserver_notebooks/LLC4320.md index b0dae343..8bf6d4cf 100644 --- a/docs/sciserver_notebooks/LLC4320.md +++ b/docs/sciserver_notebooks/LLC4320.md @@ -14,7 +14,7 @@ kernelspec: # Plot surface streamlines on LLC4320 Author: Wenrui Jiang 14 June, 2023 -> **Warning**⚠️ : the notebook was last ran on **2023-07-06** with **seaduck 0.1.dev433+g99d4d3b**. You can find the executable version at https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/LLC4320.ipynb. +> **Warning**⚠️ : the notebook was last ran on **2023-07-06** with **seaduck 0.1.dev433+g99d4d3b**. You can find the executable version at https://github.com/MaceKuailv/seaduck_sciserver_notebook/blob/master/LLC4320.ipynb. The LLC4320 ocean circulation model solution is a kilometer-scale, global simulation with complex grid topology. This is a good dataset to test the performance of the `seaduck` package. ```{code-cell} ipython3 @@ -32,7 +32,7 @@ mpl.rcParams["figure.dpi"] = 600 ``` ```{admonition} Access LLC4320 -The global MITgcm run is the LLC4320 simulation and is publicly available on [SciServer](sciserver.org) (from the Oceanography container). The simulation output can be opened using the [OceanSpy](https://github.com/hainegroup/oceanspyhttps://github.com/hainegroup/oceanspy) package using the [`from_catalog`](https://oceanspy.readthedocs.io/en/latest/generated/oceanspy.open_oceandataset.from_catalog.html#oceanspy.open_oceandataset.from_catalog) method. +The global MITgcm run is the LLC4320 simulation and is publicly available on [SciServer](sciserver.org) (from the Oceanography container). The simulation output can be opened using the [OceanSpy](https://github.com/hainegroup/oceanspyhttps://github.com/hainegroup/oceanspy) package using the [`from_catalog`](https://oceanspy.readthedocs.io/en/latest/generated/oceanspy.open_oceandataset.from_catalog.html#oceanspy.open_oceandataset.from_catalog) method. `od = ospy.open_oceandataset.from_catalog("LLC4320")`