This notebook shows some (very) brief examples of current PyDSD functionality and how to interact with it. PyDSD is used to read in disdrometer and particle probe data, plot the distributions, calculate parameterizations, radar equivalent parameters, Convective/Stratiform partitioning, and plot the results.
Let's start by reading in parsivel data which has been stored in a netCDF format by the Department of Energy ARM Program for the CACTI field campaign in Argentina. This data can be retrieved from the ARM archive.
In [1]:
%matplotlib inline
import pydsd
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import numpy as np
import pydsd.plot
import pydsd.partition
In [2]:
filename = './corldM1.b1.20181112.000000.cdf.v0' #Parsivel 05, March 2nd
Reading in files is quite easy. There are several readers in pydsd in either pydsd.io or pydsd.io.aux_readers depending on their level of support. In this case we will use the read_parsivel_arm_netcdf function. Generally we design readers for each new format we encounter so if you find a format that we cannot read, feel free to forward it to us, or even better submit a reader yourself. They are very easy to implement.
In [3]:
dsd = pydsd.read_parsivel_arm_netcdf(filename)
So at this point we have the drop size distribution read in. Let's start by looking at the format of the dsd object. A full listing of the features can be found at http://josephhardinee.github.io/PyDSD/pydsd.html#module-pydsd.DropSizeDistribution .
Generally data is stored in the fields dictionary. This is a dictionary of dictionaries where each key corresponds to a variable. Inside each of these dictionaries the data is stored in data while other metadata is attached at the top level.
In [4]:
dsd.fields.keys()
Out[4]:
In [5]:
dsd.fields['Nd']
Out[5]:
We can now start to plot and visualize some of this data. PyDSD has several built in plotting functions, and you can always pull the data out yourself for more custom plotting routines. Many built-in plotting routines are available in the pydsd.plot module.
In [6]:
pydsd.plot.plot_dsd(dsd)
plt.title('Drop Size Distribution - November 12, 2018')
Out[6]:
Plotting routines usually take in extra arguments to help customize the plot. Please check the documentation for the full list of arguments. We'll revisit plotting in a bit once we've calculated a few more interesting parameters.
Depending on the type of Disdrometer, the drop spectra may be stored. PyDSD has some capabilities for working with spectra data including filtering, and reducing this back down to a DSD measurement.
TODO: Examples of spectra processing
In [7]:
dsd.calculate_dsd_parameterization()
Now PyDSD will calculate various parameters and store these on the dsd object in the fields dictionary.
In [8]:
dsd.fields.keys()
Out[8]:
And we can dig down further and see what is in one of these objects. Similar to the defaults read in, this will store the values in the data member, and metadata attached to the object.
In [9]:
dsd.fields['D0']
Out[9]:
Now we can use some built in plotting functions to examine these. The first plot type is a time series plot of arbitrary 1D parameters. For example we can plot the D0 variable corresponding to median drop diameter. The plots know how to handle time and various pieces of metadata. Note that you can pass through arguments accepted by plt.plot and we will pass these on so you can customize the plot.
In [10]:
pydsd.plot.plot_ts(dsd, 'D0', marker='.')
Out[10]:
These plots accept an axis argument as well incase you want to use this to make side by side comparison plots.
In [11]:
plt.figure(figsize=(12,6))
ax = plt.subplot(1,2,1)
pydsd.plot.plot_ts(dsd, 'D0', marker='.', ax =ax)
ax = plt.subplot(1,2,2)
pydsd.plot.plot_ts(dsd, 'Nw', marker='.', ax =ax)
plt.tight_layout()
We have other standard types of plots built in to make life easier. For instance a normal thing to do is compare the relationship of the median drop diameter, with the normalized intercept parameter (Nw).
In [12]:
pydsd.plot.plot_NwD0(dsd)
Out[12]:
As before we can compose these by passing in an ax argument.
In [13]:
fig = plt.figure(figsize=(12,6))
ax = plt.subplot(1,2,1)
pydsd.plot.plot_dsd(dsd, ax=ax)
ax = plt.subplot(1,2,2)
pydsd.plot.plot_NwD0(dsd, ax=ax)
plt.tight_layout()
Finally let's visualize a few more of the calculated fields. We can also look at what all new fields have appeared.
In [14]:
dsd.fields.keys()
Out[14]:
In [15]:
plt.figure(figsize=(12,12))
plt.subplot(2,2,1)
pydsd.plot.plot_ts(dsd, 'D0', x_min_tick_format='hour')
plt.xlabel('Time(hrs)')
plt.ylabel('$D_0$')
# plt.xlim(5,24)
plt.subplot(2,2,2)
pydsd.plot.plot_ts(dsd, 'Nw', x_min_tick_format='hour')
plt.xlabel('Time(hrs)')
plt.ylabel('$log_{10}(N_w)$')
# plt.xlim(5,24)
plt.subplot(2,2,3)
pydsd.plot.plot_ts(dsd, 'Dmax', x_min_tick_format='hour')
plt.xlabel('Time(hrs)')
plt.ylabel('Maximum Drop Size')
# plt.xlim(5,24)
plt.subplot(2,2,4)
pydsd.plot.plot_ts(dsd, 'mu', x_min_tick_format='hour')
plt.xlabel('Time(hrs)')
plt.ylabel('$\mu$')
plt.tight_layout()
plt.show()
Note the fit submodule has alternative algorithms for calculating various DSD parameter fits.
We can calculate radar equivalent parameters as well. We use the PyTMatrix library under the hood for this. Let's look at what these measurements would look like if we did T-Matrix scatttering at X-band, which is the default.
In [16]:
dsd.calculate_radar_parameters()
This assumes BC shape relationship, X band, 10C. All of the scattering options are fully configurable.
In [17]:
dsd.set_scattering_temperature_and_frequency(scattering_temp=10, scattering_freq=9700000000.0)
dsd.set_canting_angle(7)
Note this updates the parameters, but does not re-scatter the fields until you ask it to for computational reasons. Let's do that now while also changing the DSR we are using, and the maximum diameter we will scatter for.
In [18]:
dsd.calculate_radar_parameters(dsr_func = pydsd.DSR.bc, max_diameter=7)
Similar to before these new fields will be added to the DropSizeDistribution object in the fields dictionary.
In [19]:
dsd.fields.keys()
Out[19]:
Now we can plot these variables up using our pydsd.plot.plot_ts function as before.
In [20]:
plt.figure(figsize=(12,12))
plt.subplot(2,2,1)
pydsd.plot.plot_ts(dsd, 'Zh', x_min_tick_format='hour')
plt.xlabel('Time(hrs)')
plt.ylabel('Reflectivity(dBZ)')
# plt.xlim(5,24)
plt.subplot(2,2,2)
pydsd.plot.plot_ts(dsd, 'Zdr', x_min_tick_format='hour')
plt.xlabel('Time(minutes)')
plt.ylabel('Differential Reflectivity(dB)')
# plt.xlim(5,24)
plt.subplot(2,2,3)
pydsd.plot.plot_ts(dsd, 'Kdp', x_min_tick_format='hour')
plt.xlabel('Time(hrs)')
plt.ylabel('Specific Differential Phase(deg/km)')
# plt.xlim(5,24)
plt.subplot(2,2,4)
pydsd.plot.plot_ts(dsd, 'Ai', x_min_tick_format='hour')
plt.xlabel('Time(hrs)')
plt.ylabel('Specific Attenuation')
# plt.xlim(5,24)
plt.tight_layout()
plt.show()
In [21]:
(r_z_a, r_z_b), opt = dsd.calculate_R_Zh_relationship()
print(f'RR(Zh) = {r_z_a} Zh **{r_z_b}')
(r_kdp_a, r_kdp_b), opt = dsd.calculate_R_Kdp_relationship()
print(f'RR(KDP) = {r_kdp_a} KDP **{r_kdp_b}')
(r_zk_a, r_zk_b1, r_zk_b2), opt = dsd.calculate_R_Zh_Kdp_relationship()
print(f'RR(Zh, KDP) = {r_zk_a} Zh **{r_zk_b1} * KDP ** {r_zk_b2}')
Let's visualize how good of a fit each of these estimators is. We have the measured rain rate from the disdrometer in the fields variable.
In [22]:
rr_z = r_z_a * np.power(dsd._idb(dsd.fields['Zh']['data']), r_z_b)
rr_kdp = r_kdp_a * np.power(dsd.fields['Kdp']['data'], r_kdp_b)
rr_zk = r_zk_a * np.power(dsd._idb(dsd.fields['Zh']['data']), r_zk_b1)* np.power(dsd.fields['Kdp']['data'], r_zk_b2)
plt.figure(figsize=(12,4))
plt.subplot(1,3,1)
plt.scatter(rr_z, dsd.fields['rain_rate']['data'])
plt.plot(np.arange(0,80))
plt.xlabel('Predicted rain rate')
plt.ylabel('Measured Rain Rate')
plt.title('R(Z)')
plt.subplot(1,3,2)
plt.scatter(rr_kdp, dsd.fields['rain_rate']['data'])
plt.plot(np.arange(0,80))
plt.xlabel('Predicted rain rate')
plt.ylabel('Measured Rain Rate')
plt.title('R(KDP)')
plt.subplot(1,3,3)
plt.scatter(rr_zk, dsd.fields['rain_rate']['data'])
plt.plot(np.arange(0,80))
plt.xlabel('Predicted rain rate')
plt.ylabel('Measured Rain Rate')
plt.title('R(Zh,KDP)')
Out[22]:
As expected, estimators that use polarimetry tend to do much better.
In [23]:
cs = pydsd.partition.cs_partition.cs_partition_bringi_2010(dsd.fields['Nw']['data'], dsd.fields['D0']['data'])
We have a few ways we can choose to visualize this. One is to just look at the ouput where 0-unclassified, 1-Stratiform, 2-convective, 3-transition
In [24]:
plt.plot(cs)
Out[24]:
We can also color code the D0,Nw points to get a better visual understanding of this algorithm.
In [25]:
plt.scatter(dsd.fields['D0']['data'], np.log10(dsd.fields['Nw']['data']), c=cs)
Out[25]:
In [ ]: