diff --git a/docs/conf.py b/docs/conf.py index 88584cc8dc..7fd13ba1f4 100755 --- a/docs/conf.py +++ b/docs/conf.py @@ -521,7 +521,7 @@ def linkcode_resolve(domain, info): # -- Options for MyST parser ---------------------------------------------- myst_heading_anchors = 3 -myst_enable_extensions = ["substitution"] +myst_enable_extensions = ["substitution", "amsmath", "dollarmath"] # -- Options for MyST-nb -------------------------------------------------- nb_execution_mode = "cache" diff --git a/docs/user_guide/examples_v3/tutorial_unitconverters.ipynb b/docs/user_guide/examples/tutorial_unitconverters.ipynb similarity index 52% rename from docs/user_guide/examples_v3/tutorial_unitconverters.ipynb rename to docs/user_guide/examples/tutorial_unitconverters.ipynb index 838c21e5ac..e9847781fa 100644 --- a/docs/user_guide/examples_v3/tutorial_unitconverters.ipynb +++ b/docs/user_guide/examples/tutorial_unitconverters.ipynb @@ -5,7 +5,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Spherical grids and UnitConversion\n" + "# Spherical grids and unit converters\n" ] }, { @@ -15,7 +15,7 @@ "source": [ "In most applications, Parcels works with `spherical` meshes, where longitude and latitude are given in degrees, while depth is given in meters. But it is also possible to use `flat` meshes, where longitude and latitude are given in meters (note that the dimensions are then still called `longitude` and `latitude` for consistency reasons).\n", "\n", - "In all cases, velocities are given in m/s. So Parcels seamlessly converts between meters and degrees, under the hood. For transparency, this tutorial explain how this works.\n" + "In all cases, velocities are given in m/s. Parcels seamlessly converts between meters and degrees, under the hood. For transparency, this guide explains how this works.\n" ] }, { @@ -23,7 +23,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Let's first import the relevant modules, and create dictionaries for the `U`, `V` and `temp` data arrays, with the velocities 1 m/s and the temperature 20C.\n" + "Let's first import the relevant modules, and generate a simple dataset on a 2D spherical mesh, with `U`, `V` and `temperature` data arrays, with the velocities 1 m/s and the temperature 20C.\n" ] }, { @@ -34,6 +34,7 @@ "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", + "import xarray as xr\n", "\n", "import parcels" ] @@ -44,16 +45,15 @@ "metadata": {}, "outputs": [], "source": [ - "xdim, ydim = (10, 20)\n", - "data = {\n", - " \"U\": np.ones((ydim, xdim), dtype=np.float32),\n", - " \"V\": np.ones((ydim, xdim), dtype=np.float32),\n", - " \"temp\": 20 * np.ones((ydim, xdim), dtype=np.float32),\n", - "}\n", - "dims = {\n", - " \"lon\": np.linspace(-15, 5, xdim, dtype=np.float32),\n", - " \"lat\": np.linspace(35, 60, ydim, dtype=np.float32),\n", - "}" + "from parcels._datasets.structured.generated import simple_UV_dataset\n", + "\n", + "nlat = 10\n", + "nlon = 18\n", + "ds = simple_UV_dataset(dims=(1, 1, nlat, nlon), mesh=\"spherical\").isel(time=0, depth=0)\n", + "ds[\"temperature\"] = ds[\"U\"] + 20 # add temperature field of 20 deg\n", + "ds[\"U\"].data[:] = 1.0 # set U to 1 m/s\n", + "ds[\"V\"].data[:] = 1.0 # set V to 1 m/s\n", + "ds" ] }, { @@ -61,9 +61,13 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We can convert these data and dims to a FieldSet object using `FieldSet.from_data`. We add the argument `mesh='spherical'` (this is the default option) to signal that all longitudes and latitudes are in degrees.\n", + "To create a `parcels.FieldSet` object, we define the `parcels.Field`s and the structured grid (`parcels.XGrid`) the fields are defined on. We add the argument `mesh='spherical'` to the `parcels.XGrid` to signal that all longitudes and latitudes are in degrees.\n", + "\n", + "```{note}\n", + "When using a `FieldSet` method for a specific dataset, such as `from_copernicusmarine()`, the grid information is known and parsed by Parcels, so we do not have to add the `mesh` argument.\n", + "```\n", "\n", - "Plotting the `U` field indeed shows a uniform 1m/s eastward flow.\n" + "Plotting the `U` field indeed shows a uniform 1 m/s eastward flow.\n" ] }, { @@ -72,8 +76,25 @@ "metadata": {}, "outputs": [], "source": [ - "fieldset = parcels.FieldSet.from_data(data, dims, mesh=\"spherical\")\n", - "plt.pcolormesh(fieldset.U.lon, fieldset.U.lat, fieldset.U.data[0, :, :], vmin=0, vmax=1)\n", + "grid = parcels.XGrid.from_dataset(ds, mesh=\"spherical\")\n", + "U = parcels.Field(\"U\", ds[\"U\"], grid, interp_method=parcels.interpolators.XLinear)\n", + "V = parcels.Field(\"V\", ds[\"V\"], grid, interp_method=parcels.interpolators.XLinear)\n", + "UV = parcels.VectorField(\"UV\", U, V)\n", + "temperature = parcels.Field(\n", + " \"temperature\", ds[\"temperature\"], grid, interp_method=parcels.interpolators.XLinear\n", + ")\n", + "fieldset = parcels.FieldSet([U, V, UV, temperature])\n", + "\n", + "plt.pcolormesh(\n", + " fieldset.U.grid.lon,\n", + " fieldset.U.grid.lat,\n", + " fieldset.U.data[0, 0, :, :],\n", + " vmin=0,\n", + " vmax=1,\n", + " shading=\"gouraud\",\n", + ")\n", + "plt.ylabel(\"Latitude\")\n", + "plt.xlabel(\"Longitude\")\n", "plt.colorbar()\n", "plt.show()" ] @@ -92,8 +113,12 @@ "metadata": {}, "outputs": [], "source": [ - "print(fieldset.UV[0, 0, 40, -5])\n", - "print(fieldset.temp[0, 0, 40, -5])" + "time = np.array([0])\n", + "z = np.array([0])\n", + "lat = np.array([40])\n", + "lon = np.array([-5])\n", + "print(fieldset.UV[time, z, lat, lon])\n", + "print(fieldset.temperature[time, z, lat, lon])" ] }, { @@ -103,7 +128,7 @@ "source": [ "While the temperature field indeed is 20C, as we defined, these printed velocities are much smaller.\n", "\n", - "This is because Parcels converts under the hood from m/s to degrees/s. This conversion is done with a `UnitConverter` object, which is stored in the `.units` attribute of each Field. Below, we print these\n" + "This is because Parcels converts under the hood from m/s to degrees/s. This conversion is done with a `parcels.converters` object, which is stored in the `.units` attribute of each Field. Below, we print these\n" ] }, { @@ -112,7 +137,7 @@ "metadata": {}, "outputs": [], "source": [ - "for fld in [fieldset.U, fieldset.V, fieldset.temp]:\n", + "for fld in [fieldset.U, fieldset.V, fieldset.temperature]:\n", " print(f\"{fld.name}: {fld.units}\")" ] }, @@ -121,7 +146,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "So the U field has a `GeographicPolar` UnitConverter object, the V field has a `Geographic` Unitconverter and the `temp` field has a `UnitConverter` object.\n", + "So the U field has a `GeographicPolar` UnitConverter object, the V field has a `Geographic` UnitConverter and the `temp` field has a `UnitConverter` object.\n", "\n", "Indeed, if we multiply the value of the V field with 1852 \\* 60 (the number of meters in 1 degree of latitude), we get the expected 1 m/s.\n" ] @@ -132,7 +157,7 @@ "metadata": {}, "outputs": [], "source": [ - "u, v = fieldset.UV[0, 0, 40, -5]\n", + "u, v = fieldset.UV[time, z, lat, lon]\n", "print(v * 1852 * 60)" ] }, @@ -150,7 +175,15 @@ "metadata": {}, "outputs": [], "source": [ - "print(fieldset.UV.eval(0, 0, 40, -5, applyConversion=False))" + "print(\n", + " fieldset.UV.eval(\n", + " time,\n", + " z,\n", + " lat,\n", + " lon,\n", + " applyConversion=False,\n", + " )\n", + ")" ] }, { @@ -166,7 +199,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "If longitudes and latitudes are given in meters, rather than degrees, simply add `mesh='flat'` when creating the FieldSet object.\n" + "If longitudes and latitudes are given in meters, rather than degrees, simply add `mesh='flat'` when creating the XGrid object.\n" ] }, { @@ -175,20 +208,38 @@ "metadata": {}, "outputs": [], "source": [ - "fieldset_flat = parcels.FieldSet.from_data(data, dims, mesh=\"flat\")\n", + "ds_flat = simple_UV_dataset(dims=(1, 1, nlat, nlon), mesh=\"flat\").isel(time=0, depth=0)\n", + "ds_flat[\"temperature\"] = ds_flat[\"U\"] + 20 # add temperature field of 20 deg\n", + "ds_flat[\"U\"].data[:] = 1.0 # set U to 1 m/s\n", + "ds_flat[\"V\"].data[:] = 1.0 # set V to 1 m/s\n", + "grid = parcels.XGrid.from_dataset(ds_flat, mesh=\"flat\")\n", + "U = parcels.Field(\"U\", ds_flat[\"U\"], grid, interp_method=parcels.interpolators.XLinear)\n", + "V = parcels.Field(\"V\", ds_flat[\"V\"], grid, interp_method=parcels.interpolators.XLinear)\n", + "UV = parcels.VectorField(\"UV\", U, V)\n", + "temperature = parcels.Field(\n", + " \"temperature\",\n", + " ds_flat[\"temperature\"],\n", + " grid,\n", + " interp_method=parcels.interpolators.XLinear,\n", + ")\n", + "fieldset_flat = parcels.FieldSet([U, V, UV, temperature])\n", "\n", "plt.pcolormesh(\n", - " fieldset_flat.U.lon,\n", - " fieldset_flat.U.lat,\n", - " fieldset_flat.U.data[0, :, :],\n", + " fieldset_flat.U.grid.lon,\n", + " fieldset_flat.U.grid.lat,\n", + " fieldset_flat.U.data[0, 0, :, :],\n", " vmin=0,\n", " vmax=1,\n", + " shading=\"gouraud\",\n", ")\n", "plt.colorbar()\n", "plt.show()\n", "\n", - "print(\"Velocities:\", fieldset_flat.UV[0, 0, 40, -5])\n", - "for fld in [fieldset_flat.U, fieldset_flat.V, fieldset_flat.temp]:\n", + "print(\n", + " \"Velocities:\",\n", + " fieldset_flat.UV[time, z, lat, lon],\n", + ")\n", + "for fld in [fieldset_flat.U, fieldset_flat.V, fieldset_flat.temperature]:\n", " print(f\"{fld.name}: {fld.units}\")" ] }, @@ -225,23 +276,35 @@ "kh_zonal = 100 # in m^2/s\n", "kh_meridional = 100 # in m^2/s\n", "\n", - "fieldset.add_field(\n", - " parcels.Field(\n", - " \"Kh_zonal\",\n", - " kh_zonal * np.ones((ydim, xdim), dtype=np.float32),\n", - " grid=fieldset.U.grid,\n", - " )\n", + "ds[\"Kh_zonal\"] = xr.DataArray(\n", + " data=kh_zonal * np.ones((nlat, nlon), dtype=np.float32), dims=[\"YG\", \"XG\"]\n", ")\n", - "fieldset.add_field(\n", - " parcels.Field(\n", - " \"Kh_meridional\",\n", - " kh_meridional * np.ones((ydim, xdim), dtype=np.float32),\n", - " grid=fieldset.U.grid,\n", - " )\n", + "\n", + "kh_zonal_field = parcels.Field(\n", + " \"Kh_zonal\",\n", + " ds[\"Kh_zonal\"],\n", + " grid=fieldset.U.grid,\n", + " interp_method=parcels.interpolators.XLinear,\n", + ")\n", + "\n", + "fieldset.add_field(kh_zonal_field)\n", + "\n", + "ds[\"Kh_meridional\"] = xr.DataArray(\n", + " data=kh_meridional * np.ones((nlat, nlon), dtype=np.float32), dims=[\"YG\", \"XG\"]\n", ")\n", "\n", + "kh_meridional_field = parcels.Field(\n", + " \"Kh_meridional\",\n", + " ds[\"Kh_meridional\"],\n", + " grid=grid,\n", + " interp_method=parcels.interpolators.XLinear,\n", + ")\n", + "\n", + "fieldset.add_field(kh_meridional_field)\n", + "\n", "for fld in [fieldset.Kh_zonal, fieldset.Kh_meridional]:\n", - " print(f\"{fld.name}: {fld[0, 0, 40, -5]:e} {fld.units}\")" + " val = fld[time, z, lat, lon]\n", + " print(f\"{fld.name}: {val} {fld.units}\")" ] }, { @@ -261,7 +324,7 @@ "outputs": [], "source": [ "deg_to_m = 1852 * 60\n", - "print(fieldset.Kh_meridional[0, 0, 40, -5] * deg_to_m**2)" + "print(fieldset.Kh_meridional[time, z, lat, lon] * deg_to_m**2)" ] }, { @@ -281,11 +344,11 @@ "\n", "| Field name | Converter object | Conversion for `mesh='spherical'` | Conversion for `mesh='flat'` |\n", "| ---------------- | ----------------------- | --------------------------------------------------------- | ---------------------------- |\n", - "| 'U' | `GeographicPolar` | $1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180})$ | 1 |\n", - "| 'V' | `Geographic` | $1852 \\cdot 60$ | 1 |\n", - "| 'Kh_zonal' | `GeographicPolarSquare` | $(1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180}))^2$ | 1 |\n", - "| 'Kh_meridional' | `GeographicSquare` | $(1852 \\cdot 60)^2$ | 1 |\n", - "| All other fields | `UnitConverter` | 1 | 1 |\n", + "| `\"U\"` | `GeographicPolar` | $1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180})$ | 1 |\n", + "| `\"V\"` | `Geographic` | $1852 \\cdot 60$ | 1 |\n", + "| `\"Kh_zonal\"` | `GeographicPolarSquare` | $(1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180}))^2$ | 1 |\n", + "| `\"Kh_meridional\"` | `GeographicSquare` | $(1852 \\cdot 60)^2$ | 1 |\n", + "| All other fields | `UnitConverter` | 1 | 1 |\n", "\n", "Only four Field names are recognised and assigned an automatic UnitConverter object. This means that things might go very wrong when e.g. a velocity field is not called `U` or `V`.\n", "\n", @@ -298,12 +361,19 @@ "metadata": {}, "outputs": [], "source": [ + "ds[\"Ustokes\"] = xr.DataArray(\n", + " data=np.ones((nlat, nlon), dtype=np.float32), dims=[\"YG\", \"XG\"]\n", + ")\n", + "\n", "fieldset.add_field(\n", " parcels.Field(\n", - " \"Ustokes\", np.ones((ydim, xdim), dtype=np.float32), grid=fieldset.U.grid\n", + " \"Ustokes\",\n", + " ds[\"Ustokes\"],\n", + " grid=fieldset.U.grid,\n", + " interp_method=parcels.interpolators.XLinear,\n", " )\n", ")\n", - "print(fieldset.Ustokes[0, 0, 40, -5])" + "print(fieldset.Ustokes[time, z, lat, lon])" ] }, { @@ -320,17 +390,15 @@ "metadata": {}, "outputs": [], "source": [ - "from parcels.tools.converters import GeographicPolar\n", - "\n", - "fieldset.Ustokes.units = GeographicPolar()\n", - "print(fieldset.Ustokes[0, 0, 40, -5])\n", - "print(fieldset.Ustokes[0, 0, 40, -5] * 1852 * 60 * np.cos(40 * np.pi / 180))" + "fieldset.Ustokes.units = parcels.GeographicPolar()\n", + "print(fieldset.Ustokes[time, z, lat, lon])\n", + "print(fieldset.Ustokes[time, z, lat, lon] * 1852 * 60 * np.cos(40 * np.pi / 180))" ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "test-notebooks", "language": "python", "name": "python3" }, @@ -344,7 +412,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.6" + "version": "3.11.0" } }, "nbformat": 4, diff --git a/docs/user_guide/index.md b/docs/user_guide/index.md index eebe94cc1e..1c456403cb 100644 --- a/docs/user_guide/index.md +++ b/docs/user_guide/index.md @@ -24,6 +24,7 @@ The tutorials written for Parcels v3 are currently being updated for Parcels v4. :titlesonly: examples/tutorial_nemo_curvilinear.ipynb +examples/tutorial_unitconverters.ipynb ``` @@ -32,7 +33,6 @@ examples/tutorial_nemo_curvilinear.ipynb - ```{toctree} :caption: Create ParticleSets