Tutorial

ci/environment-docs.yml

@@ -10,6 +10,7 @@ dependencies:
   - metpy
   - myst-nb
   - netcdf4
+  - pint-xarray
   - pip
   - pooch
   - pydantic

docs/source/_static/custom.js

@@ -0,0 +1,30 @@
+console.log('RUNNING CUSTOM JAVASCRIPT')
+
+function getBodyThemeName() {
+  var bodyTheme = document.body.dataset.theme
+  const prefersDark = window.matchMedia('(prefers-color-scheme: dark)').matches
+
+  if (bodyTheme === 'auto') {
+    if (prefersDark) {
+      bodyTheme = 'dark'
+    } else {
+      bodyTheme = 'light'
+    }
+  }
+
+  return bodyTheme
+}
+
+document.documentElement.setAttribute('theme', getBodyThemeName())
+
+var observer = new MutationObserver(function (mutations) {
+  mutations.forEach(function (mutation) {
+    if (mutation.type == 'attributes') {
+      var bodyTheme = getBodyThemeName()
+      console.log('document.body.dataset.theme changed to ' + bodyTheme)
+      document.documentElement.setAttribute('theme', bodyTheme)
+    }
+  })
+})
+
+observer.observe(document.body, { attributeFilter: ['data-theme'] })

docs/source/conf.py

@@ -96,6 +96,7 @@
     'doc_path': 'docs',
 }
 html_static_path = ['_static']
+html_js_files = ['custom.js']
 html_theme_options = dict(
     # analytics_id=''  this is configured in rtfd.io
     # canonical_url="",

docs/source/tutorials/Overview.md

@@ -0,0 +1,108 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.13.8
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
++++ {"tags": []}
+
+# xWRF overview
+
++++
+
+The package xWRF is designed to enable a more pythonic post-processing of [WRF](https://www.mmm.ucar.edu/weather-research-and-forecasting-model) output data. It's aim is to smooth the rough edges around the unique, non CF-compliant output data format and make the data accessible to utilities like [`dask`](https://dask.org/) and the wider [Pangeo](https://pangeo.io/) universe.
+
+It is built as an [Accessor](https://xarray.pydata.org/en/stable/internals/extending-xarray.html) on top of [`xarray`](https://xarray.pydata.org/en/stable/index.html), providing a very simple user interface.
+
++++ {"tags": []}
+
+## Examining the data
+
++++
+
+When opening up a normal [WRF](https://www.mmm.ucar.edu/weather-research-and-forecasting-model) output file with the simple `xarray` netcdf backend, one can see that it does not provide a lot of useful information.
+
+```{code-cell} ipython3
+import xwrf
+
+ds_old = xwrf.tutorial.open_dataset("wrfout")
+ds_old
+```
+
+While all variables are present, e.g. the information about the projection is still in the metadata and also for some fields, there are non-`metpy` compliant units attributes.
+
+So let's try to use the standard `xWRF.postprocess()` function.
+
+```{code-cell} ipython3
+ds_new = xwrf.tutorial.open_dataset("wrfout").xwrf.postprocess()
+ds_new
+```
+
+As you see, xWRF added some coordinate data, reassigned some dimensions and generally increased the amount of information available in the dataset.
+
++++
+
+## Projection treatment
+
++++
+
+As you can see, xWRF has determined the correct projection from the netCDF metadata and has added a `wrf_projection` variable to the dataset storing the [`pyproj`](https://pyproj4.github.io/pyproj/stable/) projection object.
+
+```{code-cell} ipython3
+ds_new['wrf_projection'].item()
+```
+
+Using this projection xWRF has also calculated the regular model grid, which the [WRF](https://www.mmm.ucar.edu/weather-research-and-forecasting-model) simulation is performed on and which can be used for e.g. bilinear interpolation.
+
+```{code-cell} ipython3
+ds_new[['x', 'y']]
+```
+
+## Attribute changes
+
++++
+
+When using xWRF, additional attributes are added to variables in order to make them [CF](https://cfconventions.org/)- and [COMODO](https://web.archive.org/web/20160417032300/http://pycomodo.forge.imag.fr/norm.html)-compliant. It also amends `unit` attributes to work with [`metpy`](https://unidata.github.io/MetPy/latest/index.html) units, enabling a seamless integration with the [Pangeo](https://pangeo.io/) software stack.
+
+Here, for example the x-wind component gets the correct [CF](https://cfconventions.org/) `standard_name` and a [COMODO](https://web.archive.org/web/20160417032300/http://pycomodo.forge.imag.fr/norm.html) `grid_mapping` attribute indicating the respective projection.
+
+```{code-cell} ipython3
+ds_old['U'].attrs, ds_new['U'].attrs
+```
+
+Also, the `units` attribute of the `PSN` variable was cleaned up to conform to [`metpy`](https://unidata.github.io/MetPy/latest/index.html) unit conventions.
+
+:::{note}
+As of now, unit translations are implemented on a manual basis, so please raise an [issue](https://github.com/ncar-xdev/xwrf/issues/new?assignees=&labels=bug%2Ctriage&template=bugreport.yml&title=%5BBug%5D%3A+) with us if you encounter any problems in this regard. In the future, this will be implemented in a more structured manner.
+:::
+
+```{code-cell} ipython3
+ds_old['PSN'].attrs['units'], ds_new['PSN'].attrs['units']
+```
+
+```{code-cell} ipython3
+import metpy
+
+try:
+    ds_old['PSN'].metpy.quantify()
+except metpy.units.UndefinedUnitError as e:
+    print(e)
+ds_new['PSN'].metpy.quantify()
+```
+
+## Diagnostic variables
+
++++
+
+Because some [WRF](https://www.mmm.ucar.edu/weather-research-and-forecasting-model) output fields are quite raw, essential diagnostic variables like the air pressure or air potential temperature are missing. These get added by xWRF by default. Users can choose to keep the fields after the computation of diagnostics is done by using `.xwrf.postprocess(drop_diagnostic_variable_components=False)`.
+
+```{code-cell} ipython3
+ds_new['air_pressure']
+```

docs/source/tutorials/index.md

@@ -3,6 +3,9 @@
 ```{toctree}
 ---
 maxdepth: 2
-caption: Tutorials
+---
+Overview.md
+xgcm.md
+metpy.md
 ---
 ```

docs/source/tutorials/metpy.md

@@ -0,0 +1,113 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.13.8
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
+# Using `metpy` for WRF analysis
+
++++
+
+In this tutorial, we will show you how xWRF enables a seamless integration of [WRF](https://www.mmm.ucar.edu/weather-research-and-forecasting-model) data with [`metpy`](https://unidata.github.io/MetPy/latest/). In the end, we will have a skewT plot at a lat-lon location in the simulation domain. Much of this tutorial was adapted from [`metpy`](https://unidata.github.io/MetPy/latest/tutorials/upperair_soundings.html).
+
++++
+
+## Loading the data
+
++++
+
+First of all, we load the data and use the simple `.xwrf.postprocess()` API.
+
+```{code-cell} ipython3
+import xwrf
+
+ds = xwrf.tutorial.open_dataset("wrfout").xwrf.postprocess()
+ds
+```
+
+## Sampling
+
++++
+
+Then, we sample the dataset at the desired location using `pyproj`. I selected the German Weather Service station at Mannheim, which is close to the middle of the simulation domain.
+
+```{code-cell} ipython3
+import numpy as np
+from pyproj import Transformer, CRS
+
+
+def sample_wrf_ds_at_latlon(ds, lat, long):
+    trf = Transformer.from_crs(CRS.from_epsg(4326), ds.wrf_projection.item(), always_xy=True)
+    x, y = ([var] if np.isscalar(var) else var for var in trf.transform(long, lat))
+    return ds.interp(x=x, y=y, x_stag=x, y_stag=y)
+```
+
+```{code-cell} ipython3
+ds = sample_wrf_ds_at_latlon(ds, 49.5091090731333, 8.553966628049963)
+ds
+```
+
+## Computation of desired quantities
+
++++
+
+Now we have to compute the quantities [`metpy`](https://unidata.github.io/MetPy/latest/) uses for the skewT. For that we have to first quantify the [WRF](https://www.mmm.ucar.edu/weather-research-and-forecasting-model) data and then convert the data to the desired units.
+
+```{code-cell} ipython3
+import metpy
+import metpy.calc as mpcalc
+import pint_xarray
+
+
+ds = ds.metpy.quantify()
+ds['dew_point_temperature'] = mpcalc.dewpoint(
+    mpcalc.vapor_pressure(ds.air_pressure, ds.QVAPOR)
+).pint.to("degC")
+ds['air_temperature'] = mpcalc.temperature_from_potential_temperature(
+    ds.air_pressure, ds.air_potential_temperature
+).pint.to("degC")
+ds['air_pressure'] = ds.air_pressure.pint.to("hPa")
+```
+
+## Plotting
+
++++
+
+Finally, we can create the skew-T plot using the computed quantities.
+
+```{code-cell} ipython3
+import matplotlib.pyplot as plt
+from metpy.plots import SkewT
+
+# Create a new figure. The dimensions here give a good aspect ratio
+fig = plt.figure(figsize=(9, 9))
+skew = SkewT(fig)
+
+# Make the dimensions of the data palatable to metpy
+if len(ds.dims) > 1:
+    ds = ds.isel(x=0, x_stag=0, y=0, y_stag=0)
+
+# Plot the data using normal plotting functions, in this case using
+# log scaling in Y, as dictated by the typical meteorological plot
+skew.plot(ds.air_pressure, ds.air_temperature, 'r', linewidth=2, label='Temperature')
+skew.plot(
+    ds.air_pressure, ds.dew_point_temperature, 'g', linewidth=2, label='Dew Point Temperature'
+)
+skew.plot_barbs(ds.air_pressure, ds.U, ds.V)
+
+# Make the plot a bit prettier
+plt.xlim([-50, 30])
+plt.xlabel('Temperature [C]')
+plt.ylabel('Pressure [hPa]')
+plt.legend().set_zorder(1)
+
+# Show the plot
+plt.show()
+```

docs/source/tutorials/xgcm.md

@@ -0,0 +1,71 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.13.8
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
+# Treating the WRF grid with `xgcm`
+
++++
+
+In this tutorial, we will show you how to leverage the xWRF-provided [COMODO](https://web.archive.org/web/20160417032300/http://pycomodo.forge.imag.fr/norm.html)-compliant attributes with `xgcm` in order to destagger the [WRF](https://www.mmm.ucar.edu/weather-research-and-forecasting-model) output.
+
++++
+
+## Loading the data
+
++++
+
+First of all, we load the data and use the simple `.xwrf.postprocess()` API.
+
+```{code-cell} ipython3
+import xwrf
+
+ds = xwrf.tutorial.open_dataset("wrfout").xwrf.postprocess()
+ds
+```
+
+If we naively try to calculate the wind speed from the `U` and `V` components, we get an error due to them having different shapes.
+
+```{code-cell} ipython3
+:tags: [raises-exception]
+
+from metpy.calc import wind_speed
+
+wind_speed(ds.U, ds.V)
+```
+
+Upon investigating the wind components (here `U`), we can see that they are defined on the [WRF](https://www.mmm.ucar.edu/weather-research-and-forecasting-model)-internal Arakawa-C grid, which causes the shapes to differ.
+
+```{code-cell} ipython3
+ds.U
+```
+
+## Destaggering the WRF grid
+
++++
+
+But, since xWRF prepared the dataset with the appropriate [COMODO](https://web.archive.org/web/20160417032300/http://pycomodo.forge.imag.fr/norm.html) (and `units`) attributes, we can simply use [`xgcm`](https://xgcm.readthedocs.io/en/latest/grids.html) to solve this problem!
+
+```{code-cell} ipython3
+from xgcm import Grid
+
+
+def interp_and_keep_attrs(grid, da, axis):
+    _attrs = da.attrs
+    da = grid.interp(da, axis=axis)
+    da.attrs = _attrs
+    return da
+
+
+grid = Grid(ds)
+ds['U'], ds['V'] = interp_and_keep_attrs(grid, ds.U, 'X'), interp_and_keep_attrs(grid, ds.V, 'Y')
+wind_speed(ds.U, ds.V)
+```

xwrf/tutorial.py

@@ -36,6 +36,7 @@ def _construct_cache_dir(path):
     'mercator': 'data/mercator_sample.nc',
     'tiny': 'data/tiny.nc',
     'met_em_sample': 'data/met_em.d01.2005-08-28_12:00:00.nc',
+    'wrfout': 'data/wrfout_d03_2018-10-16_00:00:00.nc',
 }
 
 
@@ -62,6 +63,7 @@ def open_dataset(
     * ``"lambert_conformal"``
     * ``"mercator"``
     * ``"met_em_sample"``
+    * ``"wrfout"``
 
     Parameters
     ----------