Revamp example notebook (josephhardinee#90)

A revamp of the examples notebook to make it more inline with the modern version of the library.  It could still use some improvements, but for now it's vastly better than the old (broken) one.

Changelog.txt

@@ -1,3 +1,6 @@
+v1.0.5:
+	DOCS: Changed docs to generated on 3.7 branch. 
+	RELEASE: Dropped testing support for Python 2.7
 v0.15.1:
 	BUGFIX: Fixed a bug in time handling for ParsivelReader
 v0.15.0:

pydsd/DropSizeDistribution.py

@@ -91,10 +91,6 @@ def __init__(self, reader, time_start=None, location=None):
             self.rain_rate = reader.fields["rain_rate"]
         except:
             self.rain_rate = None
-        try:
-            self.velocity = reader.fields["terminal_velocity"]
-        except:
-            self.velocity = None
         try:
             self.Z = reader.fields["reflectivity"]
         except:
@@ -131,6 +127,15 @@ def __init__(self, reader, time_start=None, location=None):
         self.set_scattering_temperature_and_frequency()
         self.set_canting_angle()
 
+        try:  # We need to make sure this is a dictionary
+            self.velocity = reader.fields["terminal_velocity"]
+        except:
+            self.velocity = None
+        if self.velocity is None:
+            self.calculate_fall_speed(self.diameter["data"], inplace=True)
+        if type(self.velocity) is np.ndarray:
+            self.velocity = {"data": self.velocity}
+
     def set_scattering_temperature_and_frequency(
         self, scattering_temp=10, scattering_freq=9.7e9
     ):
@@ -266,7 +271,8 @@ def _setup_scattering(self, wavelength, dsr_func, max_diameter):
         self.dsr_func = dsr_func
         self.scatterer.psd_integrator.D_max = max_diameter
         self.scatterer.psd_integrator.geometries = (
-            tmatrix_aux.geom_horiz_back, tmatrix_aux.geom_horiz_forw
+            tmatrix_aux.geom_horiz_back,
+            tmatrix_aux.geom_horiz_forw,
         )
         self.scatterer.or_pdf = orientation.gaussian_pdf(
             self.scattering_params["canting_angle"]
@@ -453,6 +459,43 @@ def _calculate_exponential_params(self, moment_1=2, moment_2=4):
 
         return Lambda, N0
 
+    def set_air_density(self, air_density):
+        """ Set air density at ground level
+        """
+        self.air_density = 1000.0
+
+    def calculate_fall_speed(self, diameter, density=1000, inplace=False):
+        """ Calculate terminal fall velocity for drops.
+
+        Parameters
+        ----------
+        diameter: array_like[float]
+            Array of diameters to calculate fall speed for. 
+        density: float, optional
+            Density of the air in millibars. Defaults to 1000mb. 
+        
+        Returns
+        -------
+        terminal_fall_speed: array_like[float]
+            Array of fall speeds matching size of diameter, adjusted for air density.
+        """
+        self.set_air_density(density)
+        velocity = 9.65 - 10.3 * np.exp(-0.6 * diameter)
+        velocity[0] = 0.0
+
+        speed_adjustment = (density / 1000.0) ** 0.4  # Based on Yu 2016
+        terminal_fall_speed = velocity * speed_adjustment
+        if self.velocity is None:
+            self.velocity = {}
+            self.velocity["data"] = terminal_fall_speed
+            self.velocity["air_density"] = self.air_density
+
+        if inplace is True:
+            self.velocity["data"] = terminal_fall_speed
+            self.velocity["air_density"] = self.air_density
+
+        return terminal_fall_speed
+
     def calculate_RR(self):
         """Calculate instantaneous rain rate.
 
@@ -462,15 +505,13 @@ def calculate_RR(self):
         for t in range(0, self.numt):
             # self.rain_rate['data'][t] = 0.6*3.1415 * 10**(-3) * np.dot(np.multiply(self.rain_rate['data'],np.multiply(self.Nd['data'][t],self.spread['data'] )),
             #    np.array(self.diameter['data'])**3)
-            velocity = 9.65 - 10.3 * np.exp(-0.6 * self.diameter["data"])
-            velocity[0] = 0.5
             self.fields["rain_rate"]["data"][t] = (
                 0.6
                 * np.pi
                 * 1e-03
                 * np.sum(
                     self._mmultiply(
-                        velocity,
+                        self.velocity["data"],
                         self.Nd["data"][t],
                         self.spread["data"],
                         np.array(self.diameter["data"]) ** 3,

pydsd/io/ParsivelReader.py

@@ -141,7 +141,9 @@ def _prep_data(self):
         )
         self.fields["terminal_velocity"] = common.var_to_dict(
             "Terminal Fall Velocity",
-            self.vd,  # np.ndarray(self.vd),
+            np.array(
+                self.vd[0]
+            ),  # Should we do something different here? Don't think we want the time series.
             "m/s",
             "Terminal fall velocity for each bin",
         )
@@ -150,7 +152,7 @@ def _prep_data(self):
             self.time = self._get_epoch_time()
         except:
             self.time = {
-                "data": np.array(self.time),
+                "data": np.array(self.time, dtype=float),
                 "units": None,
                 "title": "Time",
                 "full_name": "Native file time",
@@ -306,31 +308,31 @@ def _get_epoch_time(self):
     )
 
     v_spread = [
-        .1,
-        .1,
-        .1,
-        .1,
-        .1,
-        .1,
-        .1,
-        .1,
-        .1,
-        .1,
-        .2,
-        .2,
-        .2,
-        .2,
-        .2,
-        .4,
-        .4,
-        .4,
-        .4,
-        .4,
-        .8,
-        .8,
-        .8,
-        .8,
-        .8,
+        0.1,
+        0.1,
+        0.1,
+        0.1,
+        0.1,
+        0.1,
+        0.1,
+        0.1,
+        0.1,
+        0.1,
+        0.2,
+        0.2,
+        0.2,
+        0.2,
+        0.2,
+        0.4,
+        0.4,
+        0.4,
+        0.4,
+        0.4,
+        0.8,
+        0.8,
+        0.8,
+        0.8,
+        0.8,
         1.6,
         1.6,
         1.6,

pydsd/plot/plot.py

@@ -9,7 +9,8 @@
 from pylab import cm
 from matplotlib.dates import DateFormatter
 from matplotlib.dates import SecondLocator, MinuteLocator, HourLocator, DayLocator
-
+from datetime import timezone
+import datetime
 
 def plot_dsd(
     dsd,
@@ -69,7 +70,7 @@ def plot_dsd(
         pdb.set_trace()
     plt.pcolormesh(
         dsd.time["data"].filled(),
-        dsd.diameter["data"].filled(),
+        dsd.diameter["data"],
         dsd.fields["Nd"]["data"].T,
         vmin=vmin,
         vmax=vmax,
@@ -380,8 +381,8 @@ def plot_ts(
     dsd,
     varname,
     date_format="%H:%M",
-    tz=None,
-    x_min_tick_format="minute",
+    tz=timezone.utc,
+    x_min_tick_format="hour",
     title=None,
     ax=None,
     fig=None,
@@ -411,8 +412,8 @@ def plot_ts(
     fig = parse_ax(fig)
 
     x_fmt = DateFormatter(date_format, tz=tz)
-
-    ax.plot_date(dsd.time["data"], dsd.fields[varname]["data"], **kwargs)
+    sample_time = [datetime.datetime.fromtimestamp(i) for i in dsd.time['data'].filled()]
+    ax.plot_date(sample_time, dsd.fields[varname]["data"], **kwargs)
 
     ax.xaxis.set_major_formatter(x_fmt)
     if x_min_tick_format == "second":
@@ -426,6 +427,11 @@ def plot_ts(
 
     if title is not None:
         ax.set_title(title)
+    else:
+        plt.title(dsd.fields[varname]['long_name'])
+
+    plt.ylabel(dsd.fields[varname]['units'])
+    plt.xlabel('Time')
     return fig, ax