Example on geographic plotting and constraint variationΒΆ

In this example we show how to plot wind fields on a map and change the default constraint coefficients using PyDDA.

This shows how important it is to have the proper intitial state and constraints when you derive your wind fields. In the first figure, the sounding was used as the initial state, but for the latter two examples we use a zero initial state which provides for more questionable winds at the edges of the Dual Doppler Lobes.

  • ../../_images/sphx_glr_plot_fun_with_constraints_001.png
  • ../../_images/sphx_glr_plot_fun_with_constraints_002.png
  • ../../_images/sphx_glr_plot_fun_with_constraints_003.png

Out:

/home/travis/miniconda3/envs/testenv/lib/python3.6/site-packages/pyart/io/cfradial.py:384: UserWarning: WARNING: valid_min not used since it
cannot be safely cast to variable data type
  data = self.ncvar[:]
/home/travis/miniconda3/envs/testenv/lib/python3.6/site-packages/pyart/io/cfradial.py:384: UserWarning: WARNING: valid_max not used since it
cannot be safely cast to variable data type
  data = self.ncvar[:]
/home/travis/build/openradar/PyDDA/pydda/vis/streamline_plot.py:396: RuntimeWarning: invalid value encountered in less
  w_filled = np.ma.masked_where(w[level, :, :] < np.min(w_vel_contours),
/home/travis/miniconda3/envs/testenv/lib/python3.6/site-packages/cartopy/mpl/geoaxes.py:1400: UserWarning: linewidths is ignored by contourf
  result = matplotlib.axes.Axes.contourf(self, *args, **kwargs)
/home/travis/build/openradar/PyDDA/pydda/retrieval/wind_retrieve.py:544: RuntimeWarning: invalid value encountered in arccos
  theta_2 = np.arccos((x-rad2[1])/b)
/home/travis/miniconda3/envs/testenv/lib/python3.6/site-packages/cartopy/mpl/geoaxes.py:1366: UserWarning: The following kwargs were not used by contour: 'color'
  result = matplotlib.axes.Axes.contour(self, *args, **kwargs)
/home/travis/build/openradar/PyDDA/pydda/retrieval/wind_retrieve.py:544: RuntimeWarning: invalid value encountered in arccos
  theta_2 = np.arccos((x-rad2[1])/b)
/home/travis/miniconda3/envs/testenv/lib/python3.6/site-packages/cartopy/mpl/geoaxes.py:1366: UserWarning: The following kwargs were not used by contour: 'color'
  result = matplotlib.axes.Axes.contour(self, *args, **kwargs)
<string>:6: DeprecationWarning: object of type <class 'numpy.float64'> cannot be safely interpreted as an integer.
<string>:6: DeprecationWarning: object of type <class 'numpy.float64'> cannot be safely interpreted as an integer.
/home/travis/build/openradar/PyDDA/pydda/retrieval/angles.py:24: RuntimeWarning: invalid value encountered in arccos
  elev = np.arccos((Re**2 + slantrsq - rh**2)/(2 * Re * slantr))
Calculating weights for radars 0 and 1
/home/travis/build/openradar/PyDDA/pydda/retrieval/wind_retrieve.py:544: RuntimeWarning: invalid value encountered in arccos
  theta_2 = np.arccos((x-rad2[1])/b)
Calculating weights for models...
Starting solver
rmsVR = 6.776753068494561
Total points:92092.0
| Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Max w
|  17.4621| 117.0106|   0.0000|   0.0000|   0.0000|   0.0000|  11.8748
Norm of gradient: 0.03349348061301913
Iterations before filter: 10
| Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Max w
|   3.3575|  44.2463|   0.0000|   0.0000|   0.0000|   0.0000|  18.9282
Norm of gradient: 0.03588228276176865
Iterations before filter: 20
| Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Max w
|   1.6812|  26.5519|   0.0000|   0.0000|   0.0000|   0.0000|  20.4316
Norm of gradient: 0.01949543388468692
Iterations before filter: 30
| Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Max w
|   0.9545|  18.2288|   0.0000|   0.0000|   0.0000|   0.0000|  23.1132
Norm of gradient: 0.011626804623717037
Iterations before filter: 40
| Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Max w
|   0.5441|  14.1966|   0.0000|   0.0000|   0.0000|   0.0000|  25.9112
Norm of gradient: 0.004176457573828959
Iterations before filter: 50
| Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Max w
|   0.3880|  11.8996|   0.0000|   0.0000|   0.0000|   0.0000|  27.1279
Norm of gradient: 0.0009014047516688181
Iterations before filter: 60
| Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Max w
|   0.3880|  11.8996|   0.0000|   0.0000|   0.0000|   0.0000|  27.1279
Norm of gradient: 0.0009014047516688181
Iterations before filter: 70
Applying low pass filter to wind field...
Iterations after filter: 1
Iterations after filter: 2
Done! Time = 130.1
/home/travis/build/openradar/PyDDA/pydda/vis/streamline_plot.py:396: RuntimeWarning: invalid value encountered in less
  w_filled = np.ma.masked_where(w[level, :, :] < np.min(w_vel_contours),
/home/travis/miniconda3/envs/testenv/lib/python3.6/site-packages/cartopy/mpl/geoaxes.py:1400: UserWarning: linewidths is ignored by contourf
  result = matplotlib.axes.Axes.contourf(self, *args, **kwargs)
/home/travis/build/openradar/PyDDA/pydda/retrieval/wind_retrieve.py:544: RuntimeWarning: invalid value encountered in arccos
  theta_2 = np.arccos((x-rad2[1])/b)
/home/travis/miniconda3/envs/testenv/lib/python3.6/site-packages/cartopy/mpl/geoaxes.py:1366: UserWarning: The following kwargs were not used by contour: 'color'
  result = matplotlib.axes.Axes.contour(self, *args, **kwargs)
/home/travis/build/openradar/PyDDA/pydda/retrieval/wind_retrieve.py:544: RuntimeWarning: invalid value encountered in arccos
  theta_2 = np.arccos((x-rad2[1])/b)
/home/travis/miniconda3/envs/testenv/lib/python3.6/site-packages/cartopy/mpl/geoaxes.py:1366: UserWarning: The following kwargs were not used by contour: 'color'
  result = matplotlib.axes.Axes.contour(self, *args, **kwargs)
<string>:6: DeprecationWarning: object of type <class 'numpy.float64'> cannot be safely interpreted as an integer.
<string>:6: DeprecationWarning: object of type <class 'numpy.float64'> cannot be safely interpreted as an integer.
/home/travis/build/openradar/PyDDA/pydda/retrieval/angles.py:24: RuntimeWarning: invalid value encountered in arccos
  elev = np.arccos((Re**2 + slantrsq - rh**2)/(2 * Re * slantr))
Calculating weights for radars 0 and 1
/home/travis/build/openradar/PyDDA/pydda/retrieval/wind_retrieve.py:544: RuntimeWarning: invalid value encountered in arccos
  theta_2 = np.arccos((x-rad2[1])/b)
Calculating weights for models...
Starting solver
rmsVR = 6.776753068494561
Total points:92092.0
| Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Max w
|  17.4621| 117.0106|   0.0000|   0.0000|   0.0000|   0.0000|  11.8748
Norm of gradient: 0.03349348061301913
Iterations before filter: 10
| Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Max w
|   3.3575|  44.2463|   0.0000|   0.0000|   0.0000|   0.0000|  18.9282
Norm of gradient: 0.03588228276176865
Iterations before filter: 20
| Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Max w
|   1.6812|  26.5519|   0.0000|   0.0000|   0.0000|   0.0000|  20.4316
Norm of gradient: 0.01949543388468692
Iterations before filter: 30
| Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Max w
|   0.9545|  18.2288|   0.0000|   0.0000|   0.0000|   0.0000|  23.1132
Norm of gradient: 0.011626804623717037
Iterations before filter: 40
| Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Max w
|   0.5441|  14.1966|   0.0000|   0.0000|   0.0000|   0.0000|  25.9112
Norm of gradient: 0.004176457573828959
Iterations before filter: 50
| Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Max w
|   0.3880|  11.8996|   0.0000|   0.0000|   0.0000|   0.0000|  27.1279
Norm of gradient: 0.0009014047516688181
Iterations before filter: 60
| Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Max w
|   0.3880|  11.8996|   0.0000|   0.0000|   0.0000|   0.0000|  27.1279
Norm of gradient: 0.0009014047516688181
Iterations before filter: 70
Applying low pass filter to wind field...
Iterations after filter: 1
Iterations after filter: 2
Done! Time = 130.9
/home/travis/build/openradar/PyDDA/pydda/vis/streamline_plot.py:396: RuntimeWarning: invalid value encountered in less
  w_filled = np.ma.masked_where(w[level, :, :] < np.min(w_vel_contours),
/home/travis/miniconda3/envs/testenv/lib/python3.6/site-packages/cartopy/mpl/geoaxes.py:1400: UserWarning: linewidths is ignored by contourf
  result = matplotlib.axes.Axes.contourf(self, *args, **kwargs)
/home/travis/build/openradar/PyDDA/pydda/retrieval/wind_retrieve.py:544: RuntimeWarning: invalid value encountered in arccos
  theta_2 = np.arccos((x-rad2[1])/b)
/home/travis/miniconda3/envs/testenv/lib/python3.6/site-packages/cartopy/mpl/geoaxes.py:1366: UserWarning: The following kwargs were not used by contour: 'color'
  result = matplotlib.axes.Axes.contour(self, *args, **kwargs)
/home/travis/build/openradar/PyDDA/pydda/retrieval/wind_retrieve.py:544: RuntimeWarning: invalid value encountered in arccos
  theta_2 = np.arccos((x-rad2[1])/b)
/home/travis/miniconda3/envs/testenv/lib/python3.6/site-packages/cartopy/mpl/geoaxes.py:1366: UserWarning: The following kwargs were not used by contour: 'color'
  result = matplotlib.axes.Axes.contour(self, *args, **kwargs)
<string>:6: DeprecationWarning: object of type <class 'numpy.float64'> cannot be safely interpreted as an integer.
<string>:6: DeprecationWarning: object of type <class 'numpy.float64'> cannot be safely interpreted as an integer.

import pydda
import pyart
import cartopy.crs as ccrs
import matplotlib.pyplot as plt


berr_grid = pyart.io.read_grid(pydda.tests.EXAMPLE_RADAR0)
cpol_grid = pyart.io.read_grid(pydda.tests.EXAMPLE_RADAR1)

# Load our radar data
sounding = pyart.io.read_arm_sonde(
    pydda.tests.SOUNDING_PATH)
u_init, v_init, w_init = pydda.initialization.make_constant_wind_field(
    berr_grid, (0.0, 0.0, 0.0))

# Let's make a plot on a map
fig = plt.figure(figsize=(7, 7))
ax = plt.axes(projection=ccrs.PlateCarree())

pydda.vis.plot_horiz_xsection_streamlines_map(
    [cpol_grid, berr_grid], ax=ax, bg_grid_no=-1, level=7, w_vel_contours=[3, 5, 8])
plt.show()

# Let's see what happens when we use a zero initialization
new_grids = pydda.retrieval.get_dd_wind_field([cpol_grid, berr_grid],
                                    u_init, v_init, w_init,
                                    Co=1.0, Cm=1500.0, frz=5000.0,
                                    mask_outside_opt=False)

fig = plt.figure(figsize=(7, 7))
ax = plt.axes(projection=ccrs.PlateCarree())

pydda.vis.plot_horiz_xsection_streamlines_map(
    new_grids, ax=ax, bg_grid_no=-1, level=7, w_vel_contours=[3, 5, 8])
plt.show()

# Or, let's make the radar data more important!
new_grids = pydda.retrieval.get_dd_wind_field([cpol_grid, berr_grid],
                                    u_init, v_init, w_init,
                                    Co=1.0, Cm=1500.0, frz=5000.0,
                                    mask_outside_opt=False)
fig = plt.figure(figsize=(7, 7))
ax = plt.axes(projection=ccrs.PlateCarree())

pydda.vis.plot_horiz_xsection_streamlines_map(
    new_grids, ax=ax, bg_grid_no=-1, level=7, w_vel_contours=[3, 5, 8])
plt.show()

Total running time of the script: ( 4 minutes 30.481 seconds)

Gallery generated by Sphinx-Gallery