{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Computing indicators on DAP subsets\n", "\n", "\n", "In a typical programming environment, the standard way to compute an indicator for a given location would be to extract the time series at the given location, then run the computation on this subset. When interacting with a remote server, things are a bit more complicated. One option would be to first call a subsetting process to extract the data at the desired location, then run the climate indicator process on that subsetted file. The other option showcased here is to pass a DAP url that encodes the [subsetting operation](https://www.unidata.ucar.edu/software/tds/current/tutorial/Subset.html).\n", "\n", "This tutorial shows how to get the index for the desired location and pass them as a DAP link to a Finch indicator process." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import os\n", "\n", "import xarray as xr\n", "from birdy import WPSClient\n", "\n", "# Link to file storing precipitation\n", "pr = \"https://pavics.ouranos.ca/twitcher/ows/proxy/thredds/dodsC/birdhouse/testdata/flyingpigeon/cmip3/pr.sresa2.miub_echo_g.run1.atm.da.nc\"\n", "\n", "# Open connection to Finch WPS server\n", "pavics_url = \"https://pavics.ouranos.ca/twitcher/ows/proxy/finch/wps\"\n", "url = os.environ.get(\"WPS_URL\", pavics_url)\n", "wps = WPSClient(url)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
              "Dimensions:    (lat: 6, bnds: 2, lon: 7, time: 7200)\n",
              "Coordinates:\n",
              "  * lat        (lat) float64 42.68 46.39 50.1 53.81 57.52 61.23\n",
              "  * lon        (lon) float64 281.2 285.0 288.8 292.5 296.2 300.0 303.8\n",
              "  * time       (time) object 2046-01-01 12:00:00 ... 2065-12-30 12:00:00\n",
              "Dimensions without coordinates: bnds\n",
              "Data variables:\n",
              "    lat_bnds   (lat, bnds) float64 ...\n",
              "    lon_bnds   (lon, bnds) float64 ...\n",
              "    time_bnds  (time, bnds) object ...\n",
              "    pr         (time, lat, lon) float32 ...\n",
              "Attributes: (12/17)\n",
              "    comment:        Spinup: restart files from end of experiment 20C3M (corre...\n",
              "    title:          MIUB  model output prepared for IPCC Fourth Assessment SR...\n",
              "    cmor_version:   0.96\n",
              "    institution:    MIUB (University of Bonn, Bonn, Germany)\n",
              "    source:         ECHO-G(1999): atmosphere: ECHAM4 (T30L19) with partial se...\n",
              "    contact:        Stephanie Legutke (legutke@dkrz.de), Seung-Ki Min(skmin@u...\n",
              "    ...             ...\n",
              "    calendar:       360_day\n",
              "    project_id:     IPCC Fourth Assessment\n",
              "    Conventions:    CF-1.0\n",
              "    id:             pcmdi.ipcc4.miub_echo_g.sresa2.run1.atm.da\n",
              "    history:        Mon Aug  1 11:42:37 2011: ncks -4 -L 7 -d lat,42.0,64.0 -...\n",
              "    NCO:            4.0.9
" ], "text/plain": [ "\n", "Dimensions: (lat: 6, bnds: 2, lon: 7, time: 7200)\n", "Coordinates:\n", " * lat (lat) float64 42.68 46.39 50.1 53.81 57.52 61.23\n", " * lon (lon) float64 281.2 285.0 288.8 292.5 296.2 300.0 303.8\n", " * time (time) object 2046-01-01 12:00:00 ... 2065-12-30 12:00:00\n", "Dimensions without coordinates: bnds\n", "Data variables:\n", " lat_bnds (lat, bnds) float64 ...\n", " lon_bnds (lon, bnds) float64 ...\n", " time_bnds (time, bnds) object ...\n", " pr (time, lat, lon) float32 ...\n", "Attributes: (12/17)\n", " comment: Spinup: restart files from end of experiment 20C3M (corre...\n", " title: MIUB model output prepared for IPCC Fourth Assessment SR...\n", " cmor_version: 0.96\n", " institution: MIUB (University of Bonn, Bonn, Germany)\n", " source: ECHO-G(1999): atmosphere: ECHAM4 (T30L19) with partial se...\n", " contact: Stephanie Legutke (legutke@dkrz.de), Seung-Ki Min(skmin@u...\n", " ... ...\n", " calendar: 360_day\n", " project_id: IPCC Fourth Assessment\n", " Conventions: CF-1.0\n", " id: pcmdi.ipcc4.miub_echo_g.sresa2.run1.atm.da\n", " history: Mon Aug 1 11:42:37 2011: ncks -4 -L 7 -d lat,42.0,64.0 -...\n", " NCO: 4.0.9" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "# Open remote dataset and extract location indices\n", "ds = xr.open_dataset(pr)\n", "ds" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [ "nbval-ignore-output" ] }, "outputs": [ { "data": { "text/plain": [ "{'lat': array(1), 'lon': array(2), 'time': slice(5040, 6840, None)}" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Use the `remap_label_indexers` function to convert coordinates to *positional* indexes.\n", "import datetime as dt\n", "\n", "coords = dict(lat=45, lon=290)\n", "index = xr.core.indexing.map_index_queries(ds, coords, method=\"nearest\").dim_indexers\n", "\n", "# The `nearest` method cannot be used with slices, so we do another call for the time period.\n", "ti = xr.core.indexing.map_index_queries(\n", " ds, dict(time=slice(\"2060-01-01\", \"2064-12-30\"))\n", ").dim_indexers\n", "\n", "# Merge the spatial and temporal indices\n", "index.update(ti)\n", "index" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Subsetting URLs\n", "\n", "The subset syntax consists in a `?` followed by comma separated list of variable names, each followed by a slice `[start, step, stop]` for each dimension. So for example, to get the very first time step of the precipitation time series over the entire grid, we'd write\n", "\n", "`?pr[0:1:0][0:1:5][0:1:6]`\n", "\n", "Note that this uses a 0-based indexing system, so `[0:1:1]` is a slice including the first and second elements.\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
              "Dimensions:  (time: 1, lat: 6, lon: 7)\n",
              "Dimensions without coordinates: time, lat, lon\n",
              "Data variables:\n",
              "    pr       (time, lat, lon) float32 ...\n",
              "Attributes: (12/17)\n",
              "    comment:        Spinup: restart files from end of experiment 20C3M (corre...\n",
              "    title:          MIUB  model output prepared for IPCC Fourth Assessment SR...\n",
              "    cmor_version:   0.96\n",
              "    institution:    MIUB (University of Bonn, Bonn, Germany)\n",
              "    source:         ECHO-G(1999): atmosphere: ECHAM4 (T30L19) with partial se...\n",
              "    contact:        Stephanie Legutke (legutke@dkrz.de), Seung-Ki Min(skmin@u...\n",
              "    ...             ...\n",
              "    calendar:       360_day\n",
              "    project_id:     IPCC Fourth Assessment\n",
              "    Conventions:    CF-1.0\n",
              "    id:             pcmdi.ipcc4.miub_echo_g.sresa2.run1.atm.da\n",
              "    history:        Mon Aug  1 11:42:37 2011: ncks -4 -L 7 -d lat,42.0,64.0 -...\n",
              "    NCO:            4.0.9
" ], "text/plain": [ "\n", "Dimensions: (time: 1, lat: 6, lon: 7)\n", "Dimensions without coordinates: time, lat, lon\n", "Data variables:\n", " pr (time, lat, lon) float32 ...\n", "Attributes: (12/17)\n", " comment: Spinup: restart files from end of experiment 20C3M (corre...\n", " title: MIUB model output prepared for IPCC Fourth Assessment SR...\n", " cmor_version: 0.96\n", " institution: MIUB (University of Bonn, Bonn, Germany)\n", " source: ECHO-G(1999): atmosphere: ECHAM4 (T30L19) with partial se...\n", " contact: Stephanie Legutke (legutke@dkrz.de), Seung-Ki Min(skmin@u...\n", " ... ...\n", " calendar: 360_day\n", " project_id: IPCC Fourth Assessment\n", " Conventions: CF-1.0\n", " id: pcmdi.ipcc4.miub_echo_g.sresa2.run1.atm.da\n", " history: Mon Aug 1 11:42:37 2011: ncks -4 -L 7 -d lat,42.0,64.0 -...\n", " NCO: 4.0.9" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "xr.open_dataset(pr + \"?pr[0:1:0][0:1:5][0:1:6]\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the returned array has no `time`, `lat` or `lon` variables. We only requested the `pr` variable, not these other coordinate variables. To remedy the situation, we add these coordinate variables to the request." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
              "Dimensions:  (lat: 6, lon: 7, time: 1)\n",
              "Coordinates:\n",
              "  * lat      (lat) float64 42.68 46.39 50.1 53.81 57.52 61.23\n",
              "  * lon      (lon) float64 281.2 285.0 288.8 292.5 296.2 300.0 303.8\n",
              "  * time     (time) object 2046-01-01 12:00:00\n",
              "Data variables:\n",
              "    pr       (time, lat, lon) float32 ...\n",
              "Attributes: (12/17)\n",
              "    comment:        Spinup: restart files from end of experiment 20C3M (corre...\n",
              "    title:          MIUB  model output prepared for IPCC Fourth Assessment SR...\n",
              "    cmor_version:   0.96\n",
              "    institution:    MIUB (University of Bonn, Bonn, Germany)\n",
              "    source:         ECHO-G(1999): atmosphere: ECHAM4 (T30L19) with partial se...\n",
              "    contact:        Stephanie Legutke (legutke@dkrz.de), Seung-Ki Min(skmin@u...\n",
              "    ...             ...\n",
              "    calendar:       360_day\n",
              "    project_id:     IPCC Fourth Assessment\n",
              "    Conventions:    CF-1.0\n",
              "    id:             pcmdi.ipcc4.miub_echo_g.sresa2.run1.atm.da\n",
              "    history:        Mon Aug  1 11:42:37 2011: ncks -4 -L 7 -d lat,42.0,64.0 -...\n",
              "    NCO:            4.0.9
" ], "text/plain": [ "\n", "Dimensions: (lat: 6, lon: 7, time: 1)\n", "Coordinates:\n", " * lat (lat) float64 42.68 46.39 50.1 53.81 57.52 61.23\n", " * lon (lon) float64 281.2 285.0 288.8 292.5 296.2 300.0 303.8\n", " * time (time) object 2046-01-01 12:00:00\n", "Data variables:\n", " pr (time, lat, lon) float32 ...\n", "Attributes: (12/17)\n", " comment: Spinup: restart files from end of experiment 20C3M (corre...\n", " title: MIUB model output prepared for IPCC Fourth Assessment SR...\n", " cmor_version: 0.96\n", " institution: MIUB (University of Bonn, Bonn, Germany)\n", " source: ECHO-G(1999): atmosphere: ECHAM4 (T30L19) with partial se...\n", " contact: Stephanie Legutke (legutke@dkrz.de), Seung-Ki Min(skmin@u...\n", " ... ...\n", " calendar: 360_day\n", " project_id: IPCC Fourth Assessment\n", " Conventions: CF-1.0\n", " id: pcmdi.ipcc4.miub_echo_g.sresa2.run1.atm.da\n", " history: Mon Aug 1 11:42:37 2011: ncks -4 -L 7 -d lat,42.0,64.0 -...\n", " NCO: 4.0.9" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "xr.open_dataset(pr + \"?pr[0:1:0][0:1:5][0:1:6],time[0:1:0],lat,lon\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's go back to our original index and convert it into a DAP subset URL." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
              "Dimensions:  (lat: 1, lon: 1, time: 1800, time_1: 1801)\n",
              "Coordinates:\n",
              "  * lat      (lat) float64 46.39\n",
              "  * lon      (lon) float64 288.8\n",
              "  * time     (time) object 2060-01-01 12:00:00 ... 2064-12-30 12:00:00\n",
              "Dimensions without coordinates: time_1\n",
              "Data variables:\n",
              "    pr       (time_1, lat, lon) float32 ...\n",
              "Attributes: (12/17)\n",
              "    comment:        Spinup: restart files from end of experiment 20C3M (corre...\n",
              "    title:          MIUB  model output prepared for IPCC Fourth Assessment SR...\n",
              "    cmor_version:   0.96\n",
              "    institution:    MIUB (University of Bonn, Bonn, Germany)\n",
              "    source:         ECHO-G(1999): atmosphere: ECHAM4 (T30L19) with partial se...\n",
              "    contact:        Stephanie Legutke (legutke@dkrz.de), Seung-Ki Min(skmin@u...\n",
              "    ...             ...\n",
              "    calendar:       360_day\n",
              "    project_id:     IPCC Fourth Assessment\n",
              "    Conventions:    CF-1.0\n",
              "    id:             pcmdi.ipcc4.miub_echo_g.sresa2.run1.atm.da\n",
              "    history:        Mon Aug  1 11:42:37 2011: ncks -4 -L 7 -d lat,42.0,64.0 -...\n",
              "    NCO:            4.0.9
" ], "text/plain": [ "\n", "Dimensions: (lat: 1, lon: 1, time: 1800, time_1: 1801)\n", "Coordinates:\n", " * lat (lat) float64 46.39\n", " * lon (lon) float64 288.8\n", " * time (time) object 2060-01-01 12:00:00 ... 2064-12-30 12:00:00\n", "Dimensions without coordinates: time_1\n", "Data variables:\n", " pr (time_1, lat, lon) float32 ...\n", "Attributes: (12/17)\n", " comment: Spinup: restart files from end of experiment 20C3M (corre...\n", " title: MIUB model output prepared for IPCC Fourth Assessment SR...\n", " cmor_version: 0.96\n", " institution: MIUB (University of Bonn, Bonn, Germany)\n", " source: ECHO-G(1999): atmosphere: ECHAM4 (T30L19) with partial se...\n", " contact: Stephanie Legutke (legutke@dkrz.de), Seung-Ki Min(skmin@u...\n", " ... ...\n", " calendar: 360_day\n", " project_id: IPCC Fourth Assessment\n", " Conventions: CF-1.0\n", " id: pcmdi.ipcc4.miub_echo_g.sresa2.run1.atm.da\n", " history: Mon Aug 1 11:42:37 2011: ncks -4 -L 7 -d lat,42.0,64.0 -...\n", " NCO: 4.0.9" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "xr.open_dataset(\n", " pr + \"?pr[5040:1:6840][1:1:1][2:1:2],lat[1:1:1],lon[2:1:2],time[5040:1:6839]\"\n", ")" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "tags": [ "nbval-ignore-output" ] }, "outputs": [ { "data": { "text/plain": [ "{'lat': array(1), 'lon': array(2), 'time': slice(5040, 6840, None)}" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "index" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "?pr[5040:1:6839][1:1:1][2:1:2],lat[1:1:1],lon[2:1:2],time[5040:1:6839]\n" ] } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "import numpy as np\n", "\n", "\n", "def dap_slice(index):\n", " \"\"\"Convert python index dictionary into DAP subset index dictionary.\"\"\"\n", " dap = {}\n", " for key, val in index.items():\n", " if isinstance(val, slice):\n", " dap[key] = f\"[{val.start}:{val.step or 1}:{val.stop - 1}]\"\n", " elif isinstance(val, int) or (isinstance(val, np.ndarray) and val.size == 1):\n", " dap[key] = f\"[{val}:1:{val}]\"\n", " else:\n", " raise ValueError(\n", " f\"Index {val} can't be converted to a DAP subset index (for {key}).\"\n", " )\n", " return dap\n", "\n", "\n", "def dap_subset(da, index):\n", " \"\"\"Return DAP subset URL.\"\"\"\n", " s = dap_slice(index)\n", " vs = [\n", " da,\n", " ] + list(da.coords.values())\n", " url = \"?\" + \",\".join([x.name + \"\".join([s[dim] for dim in x.dims]) for x in vs])\n", " return url\n", "\n", "\n", "sub = dap_subset(ds.pr, index)\n", "print(sub)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
              "Dimensions:  (lat: 1, lon: 1, time: 1800)\n",
              "Coordinates:\n",
              "  * lat      (lat) float64 46.39\n",
              "  * lon      (lon) float64 288.8\n",
              "  * time     (time) object 2060-01-01 12:00:00 ... 2064-12-30 12:00:00\n",
              "Data variables:\n",
              "    pr       (time, lat, lon) float32 ...\n",
              "Attributes: (12/17)\n",
              "    comment:        Spinup: restart files from end of experiment 20C3M (corre...\n",
              "    title:          MIUB  model output prepared for IPCC Fourth Assessment SR...\n",
              "    cmor_version:   0.96\n",
              "    institution:    MIUB (University of Bonn, Bonn, Germany)\n",
              "    source:         ECHO-G(1999): atmosphere: ECHAM4 (T30L19) with partial se...\n",
              "    contact:        Stephanie Legutke (legutke@dkrz.de), Seung-Ki Min(skmin@u...\n",
              "    ...             ...\n",
              "    calendar:       360_day\n",
              "    project_id:     IPCC Fourth Assessment\n",
              "    Conventions:    CF-1.0\n",
              "    id:             pcmdi.ipcc4.miub_echo_g.sresa2.run1.atm.da\n",
              "    history:        Mon Aug  1 11:42:37 2011: ncks -4 -L 7 -d lat,42.0,64.0 -...\n",
              "    NCO:            4.0.9
" ], "text/plain": [ "\n", "Dimensions: (lat: 1, lon: 1, time: 1800)\n", "Coordinates:\n", " * lat (lat) float64 46.39\n", " * lon (lon) float64 288.8\n", " * time (time) object 2060-01-01 12:00:00 ... 2064-12-30 12:00:00\n", "Data variables:\n", " pr (time, lat, lon) float32 ...\n", "Attributes: (12/17)\n", " comment: Spinup: restart files from end of experiment 20C3M (corre...\n", " title: MIUB model output prepared for IPCC Fourth Assessment SR...\n", " cmor_version: 0.96\n", " institution: MIUB (University of Bonn, Bonn, Germany)\n", " source: ECHO-G(1999): atmosphere: ECHAM4 (T30L19) with partial se...\n", " contact: Stephanie Legutke (legutke@dkrz.de), Seung-Ki Min(skmin@u...\n", " ... ...\n", " calendar: 360_day\n", " project_id: IPCC Fourth Assessment\n", " Conventions: CF-1.0\n", " id: pcmdi.ipcc4.miub_echo_g.sresa2.run1.atm.da\n", " history: Mon Aug 1 11:42:37 2011: ncks -4 -L 7 -d lat,42.0,64.0 -...\n", " NCO: 4.0.9" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "xr.open_dataset(pr + sub)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using the subset URL in a WPS process\n", "\n", "Now this subset url can be used as a normal netCDF link in WPS processes. Here, let's compute the average precipitation during wet days (sdii) over our subset. As expected, the output is only computed for the five years we requested on a single grid point closest to the coordinates we chose." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'sdii' (time: 5, lat: 1, lon: 1)>\n",
              "[5 values with dtype=float64]\n",
              "Coordinates:\n",
              "  * lat      (lat) float64 46.39\n",
              "  * lon      (lon) float64 288.8\n",
              "  * time     (time) object 2060-01-01 00:00:00 ... 2064-01-01 00:00:00\n",
              "Attributes:\n",
              "    units:          mm d-1\n",
              "    cell_methods:   time: mean (interval: 30 minutes)\n",
              "    history:        pr=max(0,pr) applied to raw data;\\n[2024-03-19 13:28:16] ...\n",
              "    standard_name:  lwe_thickness_of_precipitation_amount\n",
              "    long_name:      Average precipitation during days with daily precipitatio...\n",
              "    description:    Annual simple daily intensity index (sdii) or annual aver...
" ], "text/plain": [ "\n", "[5 values with dtype=float64]\n", "Coordinates:\n", " * lat (lat) float64 46.39\n", " * lon (lon) float64 288.8\n", " * time (time) object 2060-01-01 00:00:00 ... 2064-01-01 00:00:00\n", "Attributes:\n", " units: mm d-1\n", " cell_methods: time: mean (interval: 30 minutes)\n", " history: pr=max(0,pr) applied to raw data;\\n[2024-03-19 13:28:16] ...\n", " standard_name: lwe_thickness_of_precipitation_amount\n", " long_name: Average precipitation during days with daily precipitatio...\n", " description: Annual simple daily intensity index (sdii) or annual aver..." ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "resp = wps.sdii(pr + sub)\n", "out = resp.get(asobj=True)\n", "out.output.sdii" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.4" } }, "nbformat": 4, "nbformat_minor": 4 }