diff --git a/README.md b/README.md index 823c6077..ff4d45de 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,9 @@ # PyDDA (Pythonic Direct Data Assimilation) -![alt text](https://github.com/openradar/PyDDA/blob/pydda_devel/pydda%20logo.png "Logo Title Text 1") +![alt text](https://github.com/openradar/PyDDA/blob/master/pydda%20logo.png "Logo Title Text 1") [](https://anaconda.org/conda-forge/pydda) [](https://anaconda.org/conda-forge/pydda/files) -[](https://travis-ci.org/openradar/PyDDA) +[](https://github.com/openradar/PyDDA/actions/workflows/python-package-conda.yml) [](https://openradarscience.org/PyDDA) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3942686.svg)](https://doi.org/10.5281/zenodo.3942686) @@ -19,7 +19,7 @@ This package is a rewrite of the Potvin et al. (2012) and Shapiro et al (2009) w This new package also uses a faster minimization technique, L-BFGS-B, which provides a factor of 2 to 5 speedup versus using the predecessor code, [NASA-Multidop](https://github.com/nasa/MultiDop), as well as a more elegant syntax as well as support for an arbitrary number of radars. The code is also threadsafe and has been tested using - HPC tools such as Dask on large (100+ core) clusers. + HPC tools such as Dask on large (100+ core) clusters. The user has an option to adjust strength of data, mass continuity constraints as well as implement a low pass filter. @@ -28,22 +28,21 @@ field from a grid. Angles.py is from Multidop and was written by Timothy Lang of NASA. -We recommend using Python 3.8+ or better and using anaconda or pip to install +We recommend using Python 3.9+ or better and using anaconda or pip to install the required dependencies of PyDDA. -In addition, in order to use the capability to load HRRR data as a constraint, the [cfgrib](https://github.com/ecmwf/cfgrib) package is needed. Since this does not work on Windows, this is an optional depdenency for those who wish to use HRRR data. To install cfgrib, simply do: +In addition, in order to use the capability to load HRRR data as a constraint, the [cfgrib](https://github.com/ecmwf/cfgrib) package is needed. Since this does not work on Windows, this is an optional dependency for those who wish to use HRRR data. To install cfgrib, simply do: pip install cfgrib Finally, PyDDA now supports Jax and TensorFlow as optional dependencies. PyDDA can be configured to use these two packages to perform the cost function minimization. We highly encourage users to take advantage of these two new engines, as they offer advantages such as better calculation of gradients and faster convergence. In addition, GPUs and TPUs are supported by these packages that can help drastically accelerate the calculations. For Tensorflow type: - conda install -c conda-forge tensorflow=2.6 tensorflow-probability + conda install -c conda-forge tensorflow tensorflow-probability For Jax, type: conda install -c conda-forge jax -======= ## Links to important documentation 1. [Examples](http://openradarscience.org/PyDDA/source/auto_examples/plot_examples.html) @@ -70,23 +69,19 @@ order to install from source, in a bash shell or the Anaconda prompt if you are ``` git clone https://github.com/openradar/PyDDA cd PyDDA -python setup.py install +pip install -e . ``` -======= ## Acknowledgments Core components of the software are adopted from the [Multidop](https://github.com/nasa/MultiDop) package by converting the C code to Python. The development of this software is supported by the Climate Model Development and Validation (CMDV) activity which is funded by the Office of Biological and Environmental Research in the US Department of Energy Office of Science. -======= ## Contributing We have a set of goals that we wish to accomplish using PyDDA, including the assimilation of data from various models in the retrieval, improved visualizations, use of radar data in antenna coordinates, and improved documentation. For more details on what contributions -would be useful to acheiving these goals, see the [PyDDA Roadmap](https://github.com/openradar/PyDDA/blob/master/ROADMAP.md). - -======= +would be useful to achieving these goals, see the [PyDDA Roadmap](https://github.com/openradar/PyDDA/blob/master/ROADMAP.md). ## Further support @@ -95,8 +90,6 @@ relegated to the [openradar Discourse group](https://openradar.discourse.group) enables the entire open radar science community to answer questions related to PyDDA so that both the maintainer and users can answer questions people may have. -======= - ## References You must cite these papers if you use PyDDA: diff --git a/ROADMAP.md b/ROADMAP.md index b1a873b5..e62f7adb 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -1,14 +1,14 @@ # PyDDA Roadmap In this document we show the details of the PyDDA roadmap. The roadmap shows the next steps for PyDDA and what help is needed for -contributions so that we can acheive our goal of providing a quality Direct Data Assmilation framework for the meterorological community +contributions so that we can achieve our goal of providing a quality Direct Data Assimilation framework for the meteorological community to use. Right now, PyDDA currently supports retrievals from an arbitrary number of Doppler radars and can integrate in data from rawinsondes, HRRR and WRF model runs. We would like improve how PyDDA assimilates data into a retrieval. Furthermore, we would like to make PyDDA more accessible to the international meteorology community. Our current goals in this regard are: * Support for a greater number of high resolution (LES) models such as CM1 * Support for integrating in data from the Rapid Refresh - * Coarser resolution reanalyses such as the NCEP reanalysis and ERA Interim would also provide useful information. + * Coarser resolution reanalyses such as the NCEP reanalysis and ERA5 would also provide useful information. Due to PyDDA's capability in customizing the weight each dataset has in the retrieval, using a weak constraint against coarse resolution reanalyses would provide a useful background. * Support for individual point analyses, such as those from wind profilers and METARs @@ -21,5 +21,5 @@ more accessible to the international meteorology community. Our current goals in Contributions can be made from people of all skill levels in Python, even none! If there is anyone who would like to make a contribution to any of these features, we would be more than welcoming of any additions. All we ask is that you follow the [Code of Conduct](https://github.com/openradar/PyDDA/blob/master/CODE_OF_CONDUCT.md) and that your contributions are in accordance with - the [Contibutor's Guide](https://openradarscience.org/PyDDA/contributors_guide/index.html), complete with documentation and unit + the [Contributor's Guide](https://openradarscience.org/PyDDA/contributors_guide/index.html), complete with documentation and unit tests where applicable. diff --git a/doc/source/contributors_guide/index.rst b/doc/source/contributors_guide/index.rst index 4581e833..12c42754 100644 --- a/doc/source/contributors_guide/index.rst +++ b/doc/source/contributors_guide/index.rst @@ -14,8 +14,8 @@ There are just some things that we ask of you. One is that your code be able to be distributed under the BSD 3-clause license, which is available in LICENSE in the main directory. -One, we ask, that when on the GitHub forum or making contributions to PyDDA -that all developers and users follow the PyDDA code of conduct. +We ask that all developers and users follow the PyDDA code of conduct when +participating in the GitHub forum or making contributions. Contributor Covenant Code of Conduct @@ -94,8 +94,8 @@ by other members of the project's leadership. **Attribution** -This Code of Conduct is adapted from the Contributor Covenant, version 1.4, -available at ``_ +This Code of Conduct is adapted from the Contributor Covenant, version 2.1, +available at ``_ Code Style ---------- @@ -121,7 +121,7 @@ Python File Setup ----------------- In a new .py file, the top of the code should have the function, sphinx comments -and the public and private functions within the .py file. Public fuunctions are +and the public and private functions within the .py file. Public functions are listed first and then private functions and classes. Private functions should have an underscore in front of the name. A space is needed between the last function and the closing docstring quotation marks. @@ -134,7 +134,7 @@ standards, modules should be added in the following order: 3. Local application imports Following the main function def line, but before the code within it, a docstring is -needed to explain all arguments, retuns, references, and other information. Please +needed to explain all arguments, returns, references, and other information. Please follow the NumPy documentation style: ``_ @@ -169,7 +169,7 @@ For an example format of the documentation, see this: Returns a 3D float array containing the v component of the wind field. The shape will be the same shape as the fields in Grid. w: 3D float array - Returns a 3D float array containing the u component of the wind field. + Returns a 3D float array containing the w component of the wind field. The shape will be the same shape as the fields in Grid. """ @@ -179,13 +179,13 @@ For an example format of the documentation, see this: Testing ------- -When adding a new function to pyart it is important to add it to the __init__.py +When adding a new function to PyDDA it is important to add it to the __init__.py under the corresponding folder. Create a test function and use assert to test the calculated values against known values. For an example, see: -``_ +``_ Pytest will run this test whenever a pull request is made to the master branch of the openradar/PyDDA repository. This will then allow the maintainers to @@ -235,7 +235,6 @@ forking the repository on GitHub, create your own branch by doing: :: git checkout -b this_branch -git branch this_branch :: diff --git a/doc/source/dev_reference/index.rst b/doc/source/dev_reference/index.rst index f8aedbc9..8384f8e0 100644 --- a/doc/source/dev_reference/index.rst +++ b/doc/source/dev_reference/index.rst @@ -43,9 +43,9 @@ Visualization module for PyDDA. :undoc-members: :show-inheritance: -=========================== -:mod:`initalization` Module -=========================== +============================ +:mod:`initialization` Module +============================ The module for creating custom initial states for the PyDDA retrieval. @@ -64,3 +64,14 @@ The module for creating custom constraints (i.e. from models, satellites) for th :members: :undoc-members: :show-inheritance: + +================ +:mod:`io` Module +================ + +Input/output utilities for reading and writing PyDDA Grids. + +.. automodule:: pydda.io + :members: + :undoc-members: + :show-inheritance: diff --git a/doc/source/index.rst b/doc/source/index.rst index 01901a8e..b1fb0fa8 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -29,7 +29,7 @@ description of each of PyDDA's subroutines. **User Guide** - The cookbook provides in-depth information on how + The user guide provides in-depth information on how to use PyDDA, including how to get started. This is where to look for general conceptual descriptions on how to use parts of PyDDA, including how to make your first wind retrieval and @@ -73,7 +73,7 @@ System Requirements ========================= This works on any modern version of Linux, Mac OS X, and Windows. For Windows, -HRRR data integration is not supported. In addition, since PyDDA takes advtange +HRRR data integration is not supported. In addition, since PyDDA takes advantage of parallelism we recommend: :: An Intel machine with at least 4 cores @@ -134,11 +134,11 @@ just type in the following commands assuming you have the above dependencies ins git clone https://github.com/openradar/PyDDA cd PyDDA - python setup.py install + pip install -e . :: -Finally, PyDDA now supports using `Jax `_ and `TensorFlow `_ -for solving the three dimensional wind field. PyDDA requries TensorFlow 2.6 and the +Finally, PyDDA now supports using `Jax `_ and `TensorFlow `_ +for solving the three dimensional wind field. PyDDA requires TensorFlow 2.x and the tensorflow-probability package for TensorFlow to be enabled. In addition, both of these packages can utilize CUDA-enabled GPUs for much faster processing. These two dependencies are optional as the user can still use PyDDA with the SciPy ecosystem. @@ -148,8 +148,8 @@ The Jax optimizer uses the same optimizer as SciPy's (`L-BFGS-B ` with a 'pydda' tag on it. This +relegated to the `openradar Discourse group `_ with a 'pydda' tag on it. This enables the entire open radar science community to answer questions related to PyDDA so that both the maintainer and users can answer questions people may have. diff --git a/doc/source/user_guide/dealiasing_velocities.rst b/doc/source/user_guide/dealiasing_velocities.rst index 676205e9..c387fc3b 100644 --- a/doc/source/user_guide/dealiasing_velocities.rst +++ b/doc/source/user_guide/dealiasing_velocities.rst @@ -5,14 +5,14 @@ Radar Data Quality Control - Dealiasing In this notebook, we will showcase how to perform quality control of your radar files, specifically dealiasing velocities. By doing this we can provide -PyDDA with quality controlled doppler velocities for dual doppler analysis. +PyDDA with quality controlled Doppler velocities for dual Doppler analysis. ------------- Read the Data ------------- -For this example, we use test data for two NEXRAD radars in -northern Oklahoma. For more information on +For this example, we use test data for two NEXRAD (WSR-88D) radars in +northern Oklahoma: KTLX (Twin Lakes, OK) and KICT (Wichita, KS). For more information on reading the radar data, consult :ref:`reading-radar-data`. Get test data:: @@ -164,7 +164,7 @@ Velocity Texture Displays Let's see what this velocity texture looks like. Additionally, a histogram of velocity texture values will allow for -the determination of a threshold to distingiush the hydrometeor signal +the determination of a threshold to distinguish the hydrometeor signal from artifacts. .. code-block:: python @@ -260,7 +260,7 @@ Filter Doppler Velocity Artifacts +++++++++++++++++++++++++++++++++ Now that we have determined which velocity texture values correspond to -artifacts within the doppler velocity data, we utilize Py-ART's GateFilter +artifacts within the Doppler velocity data, we utilize Py-ART's GateFilter to filter out these artifacts Py-ART's GateFilter function:: @@ -279,7 +279,7 @@ Apply Dealiasing ---------------- Now that we have removed artifacts, we can proceed with dealiasing -the doppler velocity data with Py-ART's Region-Based Dealiasing +the Doppler velocity data with Py-ART's Region-Based Dealiasing Algorithm. The Region-Based Dealiasing finds regions of similar velocities and @@ -449,9 +449,9 @@ Summary ------- Utilizing Py-ART, we read in two radar files within close proximity to each other. -We then corrected the radar doppler velocities to remove artifacts and clutter. -Finally, utilizing Py-ART, we applied a region-based dealiasing alogrithm to -unfold the doppler velocities. +We then corrected the radar Doppler velocities to remove artifacts and clutter. +Finally, utilizing Py-ART, we applied a region-based dealiasing algorithm to +unfold the Doppler velocities. Now that we have corrected velocities, incorporating these radars into PyDDA will be shown in the next notebook. diff --git a/doc/source/user_guide/gridding.rst b/doc/source/user_guide/gridding.rst index 9053949c..464ad051 100644 --- a/doc/source/user_guide/gridding.rst +++ b/doc/source/user_guide/gridding.rst @@ -6,7 +6,7 @@ Converting the radar data to Cartesian coordinates with Py-ART PyDDA expects radar data to be in Cartesian coordinates before it retrieves the wind fields. Most radar data, however, is in the radar's native antenna coordinates. Therefore, the radar data needs to be converted to Cartesian -coordinates. Py-ART's mapping toolbox contains the necessary utilties +coordinates. Py-ART's mapping toolbox contains the necessary utilities We will assume that you have followed the steps outlined in :ref:`reading-radar-data` for reading the radar data in its native coordinates. PyDDA requires dealiased velocities @@ -46,7 +46,7 @@ The :code:`grid_limits` is a 3-tuple of 2-tuples specifying the :math:`z`, :math limits of the grid in meters. The :code:`grid_shape` specifies the shape of the grid in number of points. We then use PyART's `grid_from_radars `_ function to create the grids :code:`grid_ktlx` and :code:`grid_kict`. PyDDA requires that both grids have the same -shape and origin, so we explictly set those in the options for both grids. +shape and origin, so we explicitly set those in the options for both grids. .. code-block:: python diff --git a/doc/source/user_guide/index.rst b/doc/source/user_guide/index.rst index 35a0cd91..b33c2959 100644 --- a/doc/source/user_guide/index.rst +++ b/doc/source/user_guide/index.rst @@ -2,7 +2,9 @@ User Guide ########################## -This is a place to include our user guide. +This user guide covers the full multi-Doppler wind retrieval workflow in PyDDA, +from reading and quality-controlling raw radar data through to retrieving and +visualizing the 3D wind field. .. toctree:: :maxdepth: 3 @@ -14,6 +16,4 @@ This is a place to include our user guide. retrieving_winds.rst optimizing_wind_retrieval.rst visualizing_winds.rst - retrieving_winds.rst - optimizing_wind_retrieval.rst nesting_wind_retrieval.rst diff --git a/doc/source/user_guide/nesting_wind_retrieval.rst b/doc/source/user_guide/nesting_wind_retrieval.rst new file mode 100644 index 00000000..06cc5414 --- /dev/null +++ b/doc/source/user_guide/nesting_wind_retrieval.rst @@ -0,0 +1,60 @@ +.. _nesting-wind-retrieval: + +Nested wind retrievals +====================== + +PyDDA supports nested wind retrievals via :func:`pydda.retrieval.get_dd_wind_field_nested`. +A nested retrieval performs the optimization on a coarse outer grid first, then uses +the result to initialize progressively finer inner grids. This approach can reduce +computational cost while still resolving fine-scale features in a region of interest. + +----------------------------- +Building the nested grid tree +----------------------------- + +Nested grids are organized in an :class:`xarray.DataTree`. Each node of the tree +represents one nesting level and contains one child dataset per radar +(named ``radar_0``, ``radar_1``, …). Per-level retrieval parameters can be stored +as node attributes to override the top-level keyword arguments for that nesting level. + +.. code-block:: python + + import xarray as xr + from xarray import DataTree + import pydda + + # Assume grid_ktlx_coarse, grid_kict_coarse, grid_ktlx_fine, grid_kict_fine + # have already been created with pyart.map.grid_from_radars and converted + # with pydda.io.read_from_pyart_grid. + + tree = DataTree.from_dict( + { + "/nest_0/radar_0": grid_ktlx_coarse, + "/nest_0/radar_1": grid_kict_coarse, + "/nest_1/radar_0": grid_ktlx_fine, + "/nest_1/radar_1": grid_kict_fine, + } + ) + + output_tree = pydda.retrieval.get_dd_wind_field_nested( + tree, + Cm=256.0, + Co=1e-2, + Cx=100, + Cy=100, + Cz=1, + model_fields=["hrrr"], + refl_field="DBZ", + engine="scipy", + ) + +The function returns the same :class:`xarray.DataTree` with an ``output_grids`` +entry added to each node containing the retrieved wind fields for that nesting level. + +.. note:: + + All grids within the same nesting level must share the same grid shape, + origin latitude/longitude, and coordinate arrays. The fine grid does not need + to match the coarse grid. + +For full API details see :func:`pydda.retrieval.get_dd_wind_field_nested`. diff --git a/doc/source/user_guide/optimizing_wind_retrieval.rst b/doc/source/user_guide/optimizing_wind_retrieval.rst index 9b833a49..f73e360a 100644 --- a/doc/source/user_guide/optimizing_wind_retrieval.rst +++ b/doc/source/user_guide/optimizing_wind_retrieval.rst @@ -115,8 +115,8 @@ We can see several potential issues with the wind retrieval. First, there are ar Dual Doppler lobe edges where updrafts are being produced by the optimization code simply because of a discontinuity in the horizontal winds at the edges of the lobes. In addition, there are other discontinuities in the horizontal winds that should be addressed. One thing -we can do to mitigate these discontunities is to increase the weight of the horizontal -smoothnes constraints. Therefore, let's prescribe :code:`Cx = 100.` and :code:`Cy = 100` +we can do to mitigate these discontinuities is to increase the weight of the horizontal +smoothness constraints. Therefore, let's prescribe :code:`Cx = 100.` and :code:`Cy = 100` to the above retrieval. .. code-block:: python diff --git a/doc/source/user_guide/overview.rst b/doc/source/user_guide/overview.rst index 62596305..e0a6f595 100644 --- a/doc/source/user_guide/overview.rst +++ b/doc/source/user_guide/overview.rst @@ -1,13 +1,60 @@ Overview ========================== -Within this userguide, we will detail the entire multi-doppler workflow, including: +PyDDA (Pythonic Direct Data Assimilation) is a Python package for retrieving +three-dimensional wind fields from one or more Doppler weather radars. This +technique is commonly called *multi-Doppler* or *dual-Doppler* wind retrieval +(though PyDDA supports an arbitrary number of radars). -* How to read radar data with Py-ART -* How to clean your radar data -* How to grid your radar data to cartesian coordinates -* How to run PyDDA and tune the parameters -* How to visualize your PyDDA output +Why multi-Doppler retrieval? +---------------------------- -If you are interested in adding a section to the user guide, -please use the template.rst file in the documentation. +A single Doppler radar measures only the component of wind motion *along* the +radar beam (the radial velocity). By combining observations from two or more +radars viewing the same storm volume from different angles, it becomes possible +to decompose the full three-dimensional wind vector (u, v, w) at each grid point. + +The quality of the retrieved vertical velocity **w** depends strongly on the +*beam-crossing angle* (BCA) — the angle between the two radar beams at each +grid point. PyDDA uses a BCA threshold (default: 30°–150°) to mask grid points +where the geometry is too unfavorable for reliable multi-Doppler synthesis. + +How PyDDA works +--------------- + +PyDDA formulates wind retrieval as a variational problem. It minimizes a cost +function **J** that penalizes deviations from: + +* observed radial velocities from each radar (*J*\ :sub:`o`) +* the anelastic mass continuity equation (*J*\ :sub:`m`) +* smoothness of the wind field (*J*\ :sub:`s`) +* optional background/sounding constraints (*J*\ :sub:`b`) +* optional model constraints, e.g., HRRR or WRF (*J*\ :sub:`mod`) + +The minimization uses the L-BFGS-B algorithm (via SciPy, Jax, or TensorFlow) +to find the wind field that best satisfies all active constraints simultaneously. +Each constraint is weighted by a user-supplied coefficient, allowing the +relative importance of each term to be tuned for a given radar network geometry +and scientific application. + +Workflow overview +----------------- + +The full multi-Doppler retrieval workflow consists of the following steps, +each covered in a dedicated section of this user guide: + +1. **Read radar data** — Load raw radar files using Py-ART + (:ref:`reading-radar-data`). +2. **Quality control** — Remove noise, ground clutter, and second-trip echoes, + then dealias the Doppler velocities (:ref:`dealiasing-velocities`). +3. **Grid to Cartesian coordinates** — Project the radar data from native + antenna coordinates onto a common Cartesian grid (:ref:`gridding`). +4. **Retrieve the wind field** — Run PyDDA's variational solver + (:ref:`retrieving_winds`). +5. **Optimize parameters** — Tune the cost function weights to balance + accuracy and smoothness (:ref:`optimizing-wind-retrieval`). +6. **Visualize** — Plot horizontal and vertical cross-sections of the + retrieved wind field (:ref:`visualizing-winds`). + +If you would like to contribute a new section to this user guide, +please use the ``template.rst`` file in the documentation source directory. diff --git a/doc/source/user_guide/read_radar_data.rst b/doc/source/user_guide/read_radar_data.rst index e949e8ce..3ee85e9e 100644 --- a/doc/source/user_guide/read_radar_data.rst +++ b/doc/source/user_guide/read_radar_data.rst @@ -25,9 +25,8 @@ there is convection as well as velocity values for dual doppler analysis. Read in the Data ---------------- -For this example, we use test data found in Py-ART for two X-Band Scanning ARM -Precipitation Radars (X-SAPR) found at the Atmospheric Radiation Measurment -(ARM) Southern Great Plains (SGP) site in Oklahoma. +For this example, we use test data for two NEXRAD (WSR-88D) radars in +northern Oklahoma: KTLX (Twin Lakes, OK) and KICT (Wichita, KS). Get test data:: https://arm-doe.github.io/pyart/API/generated/pyart.testing.get_test_data.html @@ -153,7 +152,7 @@ Plot reflectivity of Both Radars lon_lines=np.arange(-98, -96.75, 0.25), ) -We can see convection on both radar images near eachother with similar timestamps which will be perfect for PyDDA. +We can see convection on both radar images near each other with similar timestamps which will be perfect for PyDDA. ++++++++++++++++++++++++++++ Plot Velocity of Both Radars diff --git a/doc/source/user_guide/retrieving_winds.rst b/doc/source/user_guide/retrieving_winds.rst index e6262fc4..f068d349 100644 --- a/doc/source/user_guide/retrieving_winds.rst +++ b/doc/source/user_guide/retrieving_winds.rst @@ -16,7 +16,7 @@ PyDDA minimizes a cost function :math:`J` that corresponds to various penalties +----------------------------------------------------------------+----------------------------+ | :math:`J_{o} = \sum_{radar} [V_{ar} - \textbf{V}]^2` | Radar winds | +----------------------------------------------------------------+----------------------------+ -| :math:`J_{o} = \sum_{domain} [V_{model} - \textbf{V}]^2` | Model winds | +| :math:`J_{mod} = \sum_{domain} [V_{model} - \textbf{V}]^2` | Model winds | +----------------------------------------------------------------+----------------------------+ | :math:`J_{b} = \sum_{background} [V_{sounding} - \textbf{V}]^2`| Sounding background | +----------------------------------------------------------------+----------------------------+ @@ -36,7 +36,7 @@ retrieval, there are many aspects that must be considered. After the data proces finished, it is now important to constrain the wind field further by adding in either sounding, point, or model data as a weak constraint in order to increase the chance that PyDDA will provide a solution that converges to a physically realistic wind field. For this particular example, -we are lucky enough to have model data from the Rapid Update Cycle that can be used as a constraint. +we are lucky enough to have model data from the Rapid Update Cycle (RUC; the predecessor to the current Rapid Refresh/RAP) that can be used as a constraint. ------------------------ Using PyDDA's data model @@ -210,7 +210,7 @@ We can see in this figure that PyDDA is resolving numerous updrafts in the mid-l is the vertical motion at the edge of the Dual Doppler lobe in the top right corner. This vertical motion is likely caused by the wind source changing from primarily the radar data to the RUC model run outside of the Dual Doppler lobes, causing a slight shift in winds that results in horizontal convergence. This convergence will result in an -updraft in the domain that is an artifiact of this switch in data sources. It is therefore recommended to not +updraft in the domain that is an artifact of this switch in data sources. It is therefore recommended to not use vertical velocity data in updrafts that are touching the Dual Doppler lobe edges to mitigate this issue. In addition, prescribing a stronger background constraint or filtering the data more often may also help mitigate this issue. We will go into this further in :ref:`optimizing-wind-retrieval`. diff --git a/doc/source/user_guide/visualizing_winds.rst b/doc/source/user_guide/visualizing_winds.rst index 4506a720..b208fe41 100644 --- a/doc/source/user_guide/visualizing_winds.rst +++ b/doc/source/user_guide/visualizing_winds.rst @@ -1,5 +1,154 @@ .. _visualizing-winds: - Visualizing the wind retrieval ============================== + +PyDDA's :mod:`pydda.vis` module provides routines for plotting horizontal and +vertical cross-sections of retrieved wind fields overlaid on a gridded radar +background field (e.g., reflectivity). Three wind-vector styles are supported: + +* **Quivers** — arrow length proportional to wind speed. +* **Barbs** — meteorological wind barbs. +* **Streamlines** — continuous flow lines. + +Each style has three geometry variants: + ++---------------------------------------------------+---------------------------------------------+ +| Function | Cross-section | ++===================================================+=============================================+ +| :func:`pydda.vis.plot_horiz_xsection_quiver` | Horizontal (constant altitude) | ++---------------------------------------------------+---------------------------------------------+ +| :func:`pydda.vis.plot_xz_xsection_quiver` | Vertical, east–west (constant y) | ++---------------------------------------------------+---------------------------------------------+ +| :func:`pydda.vis.plot_yz_xsection_quiver` | Vertical, north–south (constant x) | ++---------------------------------------------------+---------------------------------------------+ +| :func:`pydda.vis.plot_horiz_xsection_barbs` | Horizontal (constant altitude) | ++---------------------------------------------------+---------------------------------------------+ +| :func:`pydda.vis.plot_xz_xsection_barbs` | Vertical, east–west (constant y) | ++---------------------------------------------------+---------------------------------------------+ +| :func:`pydda.vis.plot_yz_xsection_barbs` | Vertical, north–south (constant x) | ++---------------------------------------------------+---------------------------------------------+ +| :func:`pydda.vis.plot_horiz_xsection_streamline` | Horizontal (constant altitude) | ++---------------------------------------------------+---------------------------------------------+ + +--------------------------------- +Horizontal cross-section (quiver) +--------------------------------- + +The most common visualization is a horizontal cross-section of horizontal winds +overlaid on reflectivity. The example below plots the wind field at vertical +level 15, with updraft contours at 1, 2, 5, and 10 m/s and quivers spaced +25 km apart. + +.. code-block:: python + + pydda.vis.plot_horiz_xsection_quiver( + grids_out, + level=15, + cmap="ChaseSpectral", + vmin=-10, + vmax=80, + quiverkey_len=10.0, + background_field="DBZ", + bg_grid_no=1, + w_vel_contours=[1, 2, 5, 10], + quiver_spacing_x_km=25.0, + quiver_spacing_y_km=25.0, + quiverkey_loc="bottom_right", + ) + +.. plot:: + + import warnings + + import cartopy.crs as ccrs + import matplotlib.pyplot as plt + import numpy as np + + import pyart + import pydda + + warnings.filterwarnings("ignore") + + ktlx_file = pydda.tests.get_sample_file("cfrad.20110520_081431.542_to_20110520_081813.238_KTLX_SUR.nc") + kict_file = pydda.tests.get_sample_file("cfrad.20110520_081444.871_to_20110520_081914.520_KICT_SUR.nc") + radar_ktlx = pyart.io.read_cfradial(ktlx_file) + radar_kict = pyart.io.read_cfradial(kict_file) + + vel_tex_ktlx = pyart.retrieve.calculate_velocity_texture(radar_ktlx, vel_field='VEL') + vel_tex_kict = pyart.retrieve.calculate_velocity_texture(radar_kict, vel_field='VEL') + radar_ktlx.add_field('velocity_texture', vel_tex_ktlx, replace_existing=True) + radar_kict.add_field('velocity_texture', vel_tex_kict, replace_existing=True) + + gatefilter_ktlx = pyart.filters.GateFilter(radar_ktlx) + gatefilter_ktlx.exclude_above('velocity_texture', 3) + gatefilter_kict = pyart.filters.GateFilter(radar_kict) + gatefilter_kict.exclude_above('velocity_texture', 3) + + vel_dealias_ktlx = pyart.correct.dealias_region_based( + radar_ktlx, vel_field='VEL', centered=True, gatefilter=gatefilter_ktlx) + vel_dealias_kict = pyart.correct.dealias_region_based( + radar_kict, vel_field='VEL', centered=True, gatefilter=gatefilter_kict) + + radar_kict.add_field('corrected_velocity', vel_dealias_kict, replace_existing=True) + radar_ktlx.add_field('corrected_velocity', vel_dealias_ktlx, replace_existing=True) + + grid_limits = ((0., 15000.), (-300000., -100000.), (-250000., 0.)) + grid_shape = (31, 201, 251) + + grid_ktlx = pyart.map.grid_from_radars( + [radar_ktlx], grid_limits=grid_limits, grid_shape=grid_shape, + gatefilter=gatefilter_ktlx, + grid_origin=(radar_kict.latitude['data'].filled(), + radar_kict.longitude['data'].filled())) + grid_kict = pyart.map.grid_from_radars( + [radar_kict], grid_limits=grid_limits, grid_shape=grid_shape, + gatefilter=gatefilter_kict, + grid_origin=(radar_kict.latitude['data'].filled(), + radar_kict.longitude['data'].filled())) + + grid_ktlx = pydda.io.read_from_pyart_grid(grid_ktlx) + grid_kict = pydda.io.read_from_pyart_grid(grid_kict) + grid_kict = pydda.constraints.add_hrrr_constraint_to_grid( + grid_kict, pydda.tests.get_sample_file('ruc2anl_130_20110520_0800_001.grb2'), + method='linear') + grid_kict = pydda.initialization.make_constant_wind_field(grid_kict, (0.0, 0.0, 0.0)) + + grids_out, _ = pydda.retrieval.get_dd_wind_field( + [grid_kict, grid_ktlx], + Cm=256.0, Co=1e-2, Cx=1, Cy=1, Cz=1, Cmod=1e-5, + model_fields=["hrrr"], refl_field='DBZ', wind_tol=0.5, + max_iterations=50, filter_window=15, filter_order=3, engine='scipy') + + pydda.vis.plot_horiz_xsection_quiver( + grids_out, level=15, cmap='ChaseSpectral', vmin=-10, vmax=80, + quiverkey_len=10.0, background_field='DBZ', bg_grid_no=1, + w_vel_contours=[1, 2, 5, 10], quiver_spacing_x_km=25.0, + quiver_spacing_y_km=25.0, quiverkey_loc='bottom_right') + +Key parameters +-------------- + +``level`` + The index of the vertical level (z dimension) to plot for horizontal + cross-sections. + +``background_field`` + The name of the gridded variable to show as the color-filled background + (e.g., ``"DBZ"`` for reflectivity). + +``bg_grid_no`` + Which grid in the input list to use for the background field. Defaults to 0. + +``w_vel_contours`` + List of vertical velocity values (m/s) at which to draw contours. Set to + ``None`` to suppress vertical velocity contours. + +``quiver_spacing_x_km`` / ``quiver_spacing_y_km`` + Horizontal spacing between quiver arrows in kilometers. Larger values + produce a less cluttered plot. + +``quiverkey_len`` + The reference wind speed (m/s) represented by the quiver key arrow. + +For full parameter descriptions see the :ref:`user` reference guide. diff --git a/pydda/__init__.py b/pydda/__init__.py index a7798e83..1b6ca089 100644 --- a/pydda/__init__.py +++ b/pydda/__init__.py @@ -12,7 +12,7 @@ from . import constraints from . import io -__version__ = "2.2" +__version__ = "2.3.0" print("Welcome to PyDDA %s" % __version__) print("If you are using PyDDA in your publications, please cite:") diff --git a/pydda/cost_functions/cost_functions.py b/pydda/cost_functions/cost_functions.py index b537bbdc..1ab64ff1 100644 --- a/pydda/cost_functions/cost_functions.py +++ b/pydda/cost_functions/cost_functions.py @@ -1,6 +1,6 @@ import numpy as np -# Adding jax inport statements +# Adding jax import statements try: import tensorflow as tf diff --git a/pydda/retrieval/wind_retrieve.py b/pydda/retrieval/wind_retrieve.py index c17cb79a..5fb50766 100644 --- a/pydda/retrieval/wind_retrieve.py +++ b/pydda/retrieval/wind_retrieve.py @@ -1,5 +1,5 @@ """ -reated on Mon Aug 7 09:17:40 2017 +Created on Mon Aug 7 09:17:40 2017 @author: rjackson """ @@ -71,7 +71,7 @@ class DDParameters(object): u_back: 1D float array (number of vertical levels) Background u wind v_back: 1D float array (number of vertical levels) - Background u wind + Background v wind u_model: list of 3D float arrays U from each model integrated into the retrieval v_model: list of 3D float arrays @@ -1308,20 +1308,20 @@ def get_dd_wind_field( All grids must have the same shape, x coordinates, y coordinates and z coordinates. u_init: 3D ndarray - The intial guess for the zonal wind field, input as a 3D array + The initial guess for the zonal wind field, input as a 3D array with the same shape as the fields in Grids. If this is None, PyDDA will use the u field in the first Grid as the initalization. v_init: 3D ndarray - The intial guess for the meridional wind field, input as a 3D array + The initial guess for the meridional wind field, input as a 3D array with the same shape as the fields in Grids. If this is None, PyDDA will use the v field in the first Grid as the initalization. w_init: 3D ndarray - The intial guess for the vertical wind field, input as a 3D array + The initial guess for the vertical wind field, input as a 3D array with the same shape as the fields in Grids. If this is None, PyDDA will use the w field in the first Grid as the initalization. engine: str (one of "scipy", "tensorflow", "jax") Setting this flag will use the solver based off of SciPy, TensorFlow, or Jax. - Using Tensorflow or Jax expands PyDDA's capabiability to take advantage of GPU-based systems. + Using TensorFlow or Jax expands PyDDA's capability to take advantage of GPU-based systems. In addition, these two implementations use automatic differentation to calculate the gradient of the cost function in order to optimize the gradient calculation. TensorFlow 2.6 and tensorflow-probability are required for the TensorFlow-based engine. @@ -1402,7 +1402,7 @@ def get_dd_wind_field( Minimum beam crossing angle in degrees between two radars. 30.0 is the typical value used in many publications. max_bca: float - Minimum beam crossing angle in degrees between two radars. 150.0 is the + Maximum beam crossing angle in degrees between two radars. 150.0 is the typical value used in many publications. upper_bc: bool Set this to true to enforce w = 0 at the top of the atmosphere. This is @@ -1428,8 +1428,8 @@ def get_dd_wind_field( Stop iterations after maximum change in winds is less than this value. tolerance: float Tolerance for :math:`L_{2}` norm of gradient before stopping. - max_wind_magnitude: float - Constrain the optimization to have :math:`|u|`, :math:`|w|`, and :math:`|w| < x` m/s. + max_wind_mag: float + Constrain the optimization to have :math:`|u|`, :math:`|v|`, and :math:`|w| < x` m/s. Returns ======= diff --git a/setup.py b/setup.py index 0cc26975..afa6d956 100644 --- a/setup.py +++ b/setup.py @@ -41,7 +41,7 @@ LICENSE = "BSD" PLATFORMS = "Linux, Windows, OSX" MAJOR = 2 -MINOR = 2 +MINOR = 3 MICRO = 0 # SCRIPTS = glob.glob('scripts/*')