apps/Radx/src/RadxVolTimingStats - updating Python plotting scripts

codebase/apps/Radx/src/RadxVolTimingStats/PlotRangeHeight.py

@@ -11,7 +11,7 @@
 import sys
 
 import numpy as np
-import scipy.stats as stats
+#import scipy.stats as stats
 import matplotlib.mlab as mlab
 import matplotlib.pyplot as plt
 import matplotlib.mlab as mlab
@@ -47,7 +47,7 @@ def main():
                       help='File path for timing data')
     parser.add_option('--title',
                       dest='title',
-                      default='Range-Height Plot',
+                      default='Range-Height-Plot',
                       help='Title for plot')
     parser.add_option('--width',
                       dest='figWidthMm',
@@ -61,34 +61,68 @@ def main():
                       dest='maxHtKm',
                       default=20,
                       help='Max height to be plotted in km')
+    parser.add_option('--compute',
+                      dest='compute', default=False,
+                      action="store_true",
+                      help='Compute heights from range and elevations')
+    parser.add_option('--elevations',
+                      dest='elevations',
+                      default='0.0,0.5,1.0,2.0,3.0,4.0,5.0,6.0,8.0,10.0,15.0,20.0,30.0,45.0',
+                      help='Enter elevation list, comma-delimited, if compute is true')
+    parser.add_option('--nGates',
+                      dest='nGates',
+                      default=500,
+                      help='Number of gates, if compute is true')
+    parser.add_option('--gateSpacingM',
+                      dest='gateSpacingM',
+                      default=250,
+                      help='Gate spacing in meters, if compute is true')
     (options, args) = parser.parse_args()
-    
+
     if (options.debug):
         print("Running ", appName, file=sys.stderr)
-        print("  data file: ", options.file, file=sys.stderr)
 
-    # read in headers
+    if (options.compute):
 
-    iret = readColumnHeaders(options.file)
-    if (iret != 0):
-        sys.exit(1)
+        if (options.debug):
+            print("  computing heights from elevation angles", file=sys.stderr)
+            print("  elevations: ", options.elevations, file=sys.stderr)
+            print("  nGates: ", options.nGates, file=sys.stderr)
+            print("  gateSpacingM: ", options.gateSpacingM, file=sys.stderr)
 
-    if (options.debug):
-        print("nGates: ", nGates, file=sys.stderr)
-        print("maxRangeKm: ", maxRangeKm, file=sys.stderr)
-        print("beamWidth: ", beamWidth, file=sys.stderr)
-        print("elevs: ", elevs, file=sys.stderr)
-        print("colHeaders: ", colHeaders, file=sys.stderr)
+        # render the plot, computing on the fly
+
+        doPlotFromCompute()
+        #plot_range_height(options)
 
-    # read in data
+    else:
+
+        if (options.debug):
+            print("  reading heights from file", file=sys.stderr)
+            print("  data file: ", options.file, file=sys.stderr)
 
-    iret = readInputData(options.file)
-    if (iret != 0):
-        sys.exit(1)
+        # read in headers
+
+        iret = readColumnHeaders(options.file)
+        if (iret != 0):
+            sys.exit(1)
+
+        if (options.debug):
+            print("nGates: ", nGates, file=sys.stderr)
+            print("maxRangeKm: ", maxRangeKm, file=sys.stderr)
+            print("beamWidth: ", beamWidth, file=sys.stderr)
+            print("elevs: ", elevs, file=sys.stderr)
+            print("colHeaders: ", colHeaders, file=sys.stderr)
+
+        # read in data
+
+        iret = readInputData(options.file)
+        if (iret != 0):
+            sys.exit(1)
 
-    # render the plot
-
-    doPlot()
+        # render the plot
+            
+        doPlotFromFile()
 
     sys.exit(0)
 
@@ -203,9 +237,9 @@ def readInputData(filePath):
     return 0
 
 ########################################################################
-# Plot
+# Plot data from file
 
-def doPlot():
+def doPlotFromFile():
 
     rangeKm = np.array(colData["rangeKm"]).astype(np.double)
 
@@ -291,6 +325,114 @@ def doPlot():
 
     plt.show()
 
+########################################################################
+
+def computeHeight(range_km, elev_deg, earth_radius_km=6375.64):
+
+    """Compute the beam height based on straight-line and curved-earth models."""
+    pseudo_r_km = earth_radius_km * (4.0 / 3.0)
+
+    # Straight-line height calculation
+    ht_straight = range_km * np.sin(np.radians(elev_deg))
+
+    # Curved-earth height calculation using numpy vectorized operations
+    ht_curved = (np.sqrt(range_km**2 +
+                         pseudo_r_km**2 +
+                         2 * range_km * pseudo_r_km * np.sin(np.radians(elev_deg)))
+                 - pseudo_r_km)
+
+    return ht_straight, ht_curved
+
+########################################################################
+# Plot data computing heights from range and elevations
+
+def doPlotFromCompute():
+
+    global scanName, elevList, elevs, colHeaders, colData
+
+    # get elevations list
+
+    elevList = options.elevations
+    elevs = elevList.split(',')
+
+    if (options.debug):
+        print("elevs: ", elevs, file=sys.stderr)
+
+    # compute the range array
+
+    nGates = int(options.nGates)
+    gateSpacingKm = float(options.gateSpacingM)/ 1000.0
+    startRangeKm = gateSpacingKm / 2.0
+    endRangeKm = startRangeKm + nGates * gateSpacingKm
+    rangeKm = np.arange(startRangeKm, endRangeKm, gateSpacingKm)
+    # print("rangeKm: ", rangeKm, file=sys.stderr)
+
+    widthIn = float(options.figWidthMm) / 25.4
+    htIn = float(options.figHeightMm) / 25.4
+
+    fig1 = plt.figure(1, (widthIn, htIn))
+    title = (options.title)
+    fig1.suptitle(title, fontsize=11)
+    ax1 = fig1.add_subplot(1,1,1,xmargin=0.0)
+
+    ax1.set_xlim([0.0, float(endRangeKm)])
+    ax1.set_ylim([0.0, float(options.maxHtKm)])
+
+    # plot the heights for each elevation angle
+
+    colors = [ 'red', 'blue', 'orange', 'green' ]
+    count = 0
+    elevsUsed = []
+
+    for elev in elevs:
+
+        if elev in elevsUsed:
+            continue
+
+        elevsUsed.append(elev)
+
+        elev_float = float(elev)
+
+        gndRangeKm = rangeKm * math.cos(math.radians(elev_float))
+
+        htLabel = "htKmStraight[" + elev + "]"
+        htKmStraight = rangeKm * math.sin(math.radians(elev_float))
+
+        earthRKm = 6375.64
+        pseudoRKm = earthRKm * (4.0 / 3.0)
+        pseudoRKmSq = pseudoRKm * pseudoRKm
+
+        htList = []
+        htList2 = []
+
+        # Compute beam heights
+        ht_straight, ht_curved = computeHeight(gndRangeKm, elev_float)
+
+        col = colors[count % 4]
+
+        x = []
+        y = []
+
+        ax1.plot(gndRangeKm, ht_straight, linewidth=2, dashes = [4, 4], color = col)
+        ax1.plot(gndRangeKm, ht_curved, linewidth=2, color = col)
+
+        count = count + 1
+
+    ax1.set_ylabel('Height above radar (km)')
+    ax1.set_xlabel('Ground range (km)')
+    ax1.set_title('Elevs: ' + elevList, fontsize=8)
+    ax1.grid(True)
+
+    homeDir = os.environ['HOME']
+    downloadsDir = os.path.join(homeDir, 'Downloads')
+    savePath = os.path.join(downloadsDir, options.title + ".png")
+    fig1.savefig(savePath)
+
+    # show
+
+    plt.show()
+    return
+
 ########################################################################
 # Annotate a value
 

codebase/apps/Radx/src/RadxVolTimingStats/PlotVolTimingStats.py

@@ -11,7 +11,7 @@
 import sys
 
 import numpy as np
-import scipy.stats as stats
+#import scipy.stats as stats
 import matplotlib.mlab as mlab
 import matplotlib.pyplot as plt
 import matplotlib.mlab as mlab