Merge pull request #116 from rcjackson/pydda20

PyDDA 2.0!

.github/workflows/pypi-release.yml

@@ -0,0 +1,74 @@
+name: Build and Upload PySP2 Release to PyPI
+on:
+  release:
+    types:
+      - published
+
+jobs:
+  build-artifacts:
+    runs-on: ubuntu-latest
+    if: github.repository == 'openradar/PyDDA'
+    steps:
+      - uses: actions/checkout@v3
+        with:
+          fetch-depth: 0
+      - uses: actions/setup-python@v4
+        name: Install Python
+        with:
+          python-version: 3.11
+
+      - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip
+          python -m pip install setuptools setuptools-scm wheel twine check-manifest
+      - name: Build tarball and wheels
+        run: |
+          git clean -xdf
+          git restore -SW .
+          python -m build --sdist --wheel .
+      - name: Check built artifacts
+        run: |
+          python -m twine check dist/*
+          pwd
+          if [ -f dist/pydda-0.0.0.tar.gz ]; then
+            echo "❌ INVALID VERSION NUMBER"
+            exit 1
+          else
+            echo "✅ Looks good"
+          fi
+      - uses: actions/upload-artifact@v3
+        with:
+          name: releases
+          path: dist
+
+  test-built-dist:
+    needs: build-artifacts
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/setup-python@v4
+        name: Install Python
+        with:
+          python-version: "3.x"
+      - uses: actions/download-artifact@v3
+        with:
+          name: releases
+          path: dist
+      - name: List contents of built dist
+        run: |
+          ls -ltrh
+          ls -ltrh dist
+  upload-to-pypi:
+    needs: test-built-dist
+    if: github.event_name == 'release'
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/download-artifact@v3
+        with:
+          name: releases
+          path: dist
+      - name: Publish package to PyPI
+        uses: pypa/[email protected]
+        with:
+          user: __token__
+          password: ${{ secrets.PYPI_TOKEN }}
+          verbose: true

.github/workflows/python-package-conda.yml

@@ -17,7 +17,7 @@ jobs:
      strategy:
         fail-fast: false
         matrix:
-          python-version: ["3.9", "3.10", "3.11"]
+          python-version: ["3.9", "3.10", "3.11", "3.12"]
           os: [macOS, ubuntu]
           inlcude:
             - os: macos-latest

REQUIREMENTS.txt

@@ -8,3 +8,5 @@ six
 pooch
 cmweather
 cdsapi
+xarray
+datatree

doc/source/index.rst

@@ -139,12 +139,19 @@ just type in the following commands assuming you have the above dependencies ins
 
 Finally, PyDDA now supports using `Jax <jax.readthedocs.io>`_ and `TensorFlow <tensorflow.org>`_
 for solving the three dimensional wind field. PyDDA requries TensorFlow 2.6 and the
-tensorflow-probability package for TensorFlow to be enabled. Both the Jax and
-TensorFlow-based engines now use automatic differentiation to solve for the gradients
-of each cost function. This therefore will create gradients that are less susceptible
-to boundary artifacts and rounding errors. In addition, both of these packages can
-utilize CUDA-enabled GPUs for much faster processing. These two
+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.
+The Jax optimizer uses the same optimizer as SciPy's (`L-BFGS-B <https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fmin_l_bfgs_b.html>`_).
+
+
+Known issues
+============
+
+The TensorFlow engine uses the unbounded version of this optimizer which removes the constraint that the
+the wind magnitudes must be less than 100 m/s. The removal of this constraint can sometimes
+result in numerical instability, so it is recommended that the user test out both Jax and TensorFlow
+if they desire GPU-accelerated retrievals.
 
 Contents:
 

doc/source/user_guide/index.rst

@@ -14,3 +14,6 @@ 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

doc/source/user_guide/optimizing_wind_retrieval.rst

@@ -95,7 +95,8 @@ that we did in :ref:`retrieving_winds`.
                                  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')
 
@@ -200,7 +201,8 @@ to the above retrieval.
                                  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')
 
@@ -302,7 +304,8 @@ Let's see what happens when we increase the level of smoothing.
                                  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')
 
@@ -405,7 +408,8 @@ artifact near the edge of the Dual Doppler lobe re-appears.
                                  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')
 

doc/source/user_guide/retrieving_winds.rst

@@ -38,8 +38,21 @@ point, or model data as a weak constraint in order to increase the chance that P
 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.
 
-Therefore, let's first add the model data into our Grid objects so that PyDDA can use the model
-data as a constraint.
+------------------------
+Using PyDDA's data model
+------------------------
+
+As of PyDDA 2.0, PyDDA uses `xarray Datasets <https://docs.xarray.dev/en/stable/generated/xarray.Dataset.html>`_
+to represent the underlying datastructure. This makes it easier to integrate PyDDA into xarray-based workflows
+that are the standard for use in the Geoscientific Python Community. Therefore, anyone using PyART Grids will need
+to convert their grids to PyDDA Grids using the :meth:`pydda.io.read_from_pyart_grid` helper function.
+
+.. code-block:: python
+
+    grid_ktlx = pydda.io.read_from_pyart_grid(grid_ktlx)
+    grid_kict = pydda.io.read_from_pyart_grid(grid_kict)
+
+In addition, :meth:`pydda.io.read_grid` will read a cf-Compliant radar grid into a PyDDA Grid.
 
 ----------------------------
 Using models for constraints
@@ -50,10 +63,11 @@ provide a constraint on the horizontal winds. This helps to constrain the backgr
 area where there is suboptimal coverage from the radar network outside of the dual
 Doppler lobes. To add either a RUC or HRRR model time period as a constraint, we need
 to have the original model data in GRIB format and then use the following line to
-load the model data into a PyART grid for processing.
+load the model data into a PyDDA grid for processing.
 
 .. code-block:: python
 
+    # Add constraints
     grid_kict = pydda.constraints.add_hrrr_constraint_to_grid(grid_kict,
         pydda.tests.get_sample_file('ruc2anl_130_20110520_0800_001.grb2'), method='linear')
 
@@ -75,7 +89,7 @@ Retrieving your first wind field
 --------------------------------
 
 We will then take a first attempt at retrieving a wind field using PyDDA. This is done using the
-:code:`pydda.retrieval.get_dd_wind_field` function. This function takes in a minimum of one input, a list
+:meth:`pydda.retrieval.get_dd_wind_field` function. This function takes in a minimum of one input, a list
 of input PyART Grids. We will also specify the constants for the constraints. In this case, we are using
 the mass continuity, radar observation, smoothness, and model constraints.
 
@@ -176,7 +190,8 @@ key on the bottom right interior of the plot.
                                  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')
 

examples/plot_examples.py

@@ -9,19 +9,15 @@
 
 """
 
-import pyart
 import pydda
-import numpy as np
 from matplotlib import pyplot as plt
 
-
-berr_grid = pyart.io.read_grid(pydda.tests.EXAMPLE_RADAR0)
-cpol_grid = pyart.io.read_grid(pydda.tests.EXAMPLE_RADAR1)
-
+berr_grid = pydda.io.read_grid(pydda.tests.EXAMPLE_RADAR0)
+cpol_grid = pydda.io.read_grid(pydda.tests.EXAMPLE_RADAR1)
 
 # Load sounding data and insert as an intialization
-cpol_grid = pydda.initialization.make_constant_wind_field(
-    cpol_grid, (0.0, 0.0, 0.0), vel_field="corrected_velocity"
+berr_grid = pydda.initialization.make_constant_wind_field(
+    berr_grid, (0.0, 0.0, 0.0), vel_field="corrected_velocity"
 )
 
 # Start the wind retrieval. This example only uses the mass continuity
@@ -41,6 +37,7 @@
     wind_tol=0.5,
     engine="scipy",
 )
+
 # Plot a horizontal cross section
 plt.figure(figsize=(9, 9))
 pydda.vis.plot_horiz_xsection_barbs(
@@ -50,9 +47,11 @@
     w_vel_contours=[5, 10, 15],
     barb_spacing_x_km=5.0,
     barb_spacing_y_km=15.0,
+    vmin=0,
+    vmax=70,
 )
 plt.show()
-plt.savefig("Darwin_horiz.png")
+
 # Plot a vertical X-Z cross section
 plt.figure(figsize=(9, 9))
 pydda.vis.plot_xz_xsection_barbs(
@@ -62,6 +61,8 @@
     w_vel_contours=[5, 10, 15],
     barb_spacing_x_km=10.0,
     barb_spacing_z_km=2.0,
+    vmin=0,
+    vmax=70,
 )
 plt.show()
 
@@ -73,5 +74,7 @@
     level=40,
     barb_spacing_y_km=10.0,
     barb_spacing_z_km=2.0,
+    vmin=0,
+    vmax=70,
 )
-plt.savefig("Darwin.png")
+plt.show()