Merge branch 'testing' into experimental

auto_rx/utils/plot_sonde_log.py

@@ -0,0 +1,234 @@
+#!/usr/bin/env python
+#
+#   Radiosonde Log Plotter
+#
+#   Copyright (C) 2019  Mark Jessop <vk5qi@rfhead.net>
+#   Released under GNU GPL v3 or later 
+#
+#   Note: This script is very much a first pass, and doesn't have any error checking of data.
+#   
+#    Usage: plot_sonde_log.py [-h] [--metric] [--alt-limit ALT_LIMIT]
+#                             [--temp-limit TEMP_LIMIT] [--decimation DECIMATION]
+#                             filename
+#
+#    positional arguments:
+#      filename              Log File name.
+#
+#    optional arguments:
+#      -h, --help            show this help message and exit
+#      --metric              Use metric altitudes. (Default is to use Feet)
+#      --alt-limit ALT_LIMIT
+#                            Limit plot to supplied altitude (feet or metres,
+#                            depending on user selection)
+#      --temp-limit TEMP_LIMIT
+#                            Limit plot to a lower temperature in degrees. (Default
+#                            is no limit, plot will autoscale)
+#      --decimation DECIMATION
+#                            Decimate input data by X times. (Default = 10)
+#
+
+import argparse
+import os.path
+import sys
+import numpy as np
+import matplotlib.pyplot as plt
+from dateutil.parser import parse
+from math import radians, degrees, sin, cos, atan2, sqrt, pi
+
+
+# Earthmaths code by Daniel Richman (thanks!)
+# Copyright 2012 (C) Daniel Richman; GNU GPL 3
+def position_info(listener, balloon):
+    """
+    Calculate and return information from 2 (lat, lon, alt) tuples
+
+    Returns a dict with:
+
+     - angle at centre
+     - great circle distance
+     - distance in a straight line
+     - bearing (azimuth or initial course)
+     - elevation (altitude)
+
+    Input and output latitudes, longitudes, angles, bearings and elevations are
+    in degrees, and input altitudes and output distances are in meters.
+    """
+
+    # Earth:
+    radius = 6371000.0
+
+    (lat1, lon1, alt1) = listener
+    (lat2, lon2, alt2) = balloon
+
+    lat1 = radians(lat1)
+    lat2 = radians(lat2)
+    lon1 = radians(lon1)
+    lon2 = radians(lon2)
+
+    # Calculate the bearing, the angle at the centre, and the great circle
+    # distance using Vincenty's_formulae with f = 0 (a sphere). See
+    # http://en.wikipedia.org/wiki/Great_circle_distance#Formulas and
+    # http://en.wikipedia.org/wiki/Great-circle_navigation and
+    # http://en.wikipedia.org/wiki/Vincenty%27s_formulae
+    d_lon = lon2 - lon1
+    sa = cos(lat2) * sin(d_lon)
+    sb = (cos(lat1) * sin(lat2)) - (sin(lat1) * cos(lat2) * cos(d_lon))
+    bearing = atan2(sa, sb)
+    aa = sqrt((sa ** 2) + (sb ** 2))
+    ab = (sin(lat1) * sin(lat2)) + (cos(lat1) * cos(lat2) * cos(d_lon))
+    angle_at_centre = atan2(aa, ab)
+    great_circle_distance = angle_at_centre * radius
+
+    # Armed with the angle at the centre, calculating the remaining items
+    # is a simple 2D triangley circley problem:
+
+    # Use the triangle with sides (r + alt1), (r + alt2), distance in a
+    # straight line. The angle between (r + alt1) and (r + alt2) is the
+    # angle at the centre. The angle between distance in a straight line and
+    # (r + alt1) is the elevation plus pi/2.
+
+    # Use sum of angle in a triangle to express the third angle in terms
+    # of the other two. Use sine rule on sides (r + alt1) and (r + alt2),
+    # expand with compound angle formulae and solve for tan elevation by
+    # dividing both sides by cos elevation
+    ta = radius + alt1
+    tb = radius + alt2
+    ea = (cos(angle_at_centre) * tb) - ta
+    eb = sin(angle_at_centre) * tb
+    elevation = atan2(ea, eb)
+
+    # Use cosine rule to find unknown side.
+    distance = sqrt((ta ** 2) + (tb ** 2) - 2 * tb * ta * cos(angle_at_centre))
+
+    # Give a bearing in range 0 <= b < 2pi
+    if bearing < 0:
+        bearing += 2 * pi
+
+    return {
+        "listener": listener, "balloon": balloon,
+        "listener_radians": (lat1, lon1, alt1),
+        "balloon_radians": (lat2, lon2, alt2),
+        "angle_at_centre": degrees(angle_at_centre),
+        "angle_at_centre_radians": angle_at_centre,
+        "bearing": degrees(bearing),
+        "bearing_radians": bearing,
+        "great_circle_distance": great_circle_distance,
+        "straight_distance": distance,
+        "elevation": degrees(elevation),
+        "elevation_radians": elevation
+    }
+
+
+if __name__ == "__main__":
+    # Data format:
+    # 2019-04-17T00:40:40.000Z,P4740856,7611,-35.38981,139.47062,12908.1,-67.9,25.0,RS41,402.500,SATS 9,BATT -1.0
+
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument("filename", type=str, help="Log File name.")
+    parser.add_argument("--metric", action="store_true", default=False, help="Use metric altitudes. (Default is to use Feet)")
+    parser.add_argument("--alt-limit", default=20000, type=int, help="Limit plot to supplied altitude (feet or metres, depending on user selection)")
+    parser.add_argument("--temp-limit", default=None, type=float, help="Limit plot to a lower temperature in degrees. (Default is no limit, plot will autoscale)")
+    parser.add_argument("--decimation", default=10, type=int, help="Decimate input data by X times. (Default = 10)")
+    args = parser.parse_args()
+
+    # Load in the file.
+    data = np.genfromtxt(args.filename,delimiter=',', dtype=None)
+
+    decimation = args.decimation
+
+    # Extract fields.
+    times = data['f0']
+    latitude = data['f3']
+    longitude = data['f4']
+    altitude = data['f5']
+    temperature = data['f6']
+    humidity = data['f7']
+
+    _output = [] # Altitude, Wind Speed, Wind Direction, Temperature, Dew Point
+    # First entry, We assume all the values are unknown for now.
+    _output.append([altitude[0], np.NaN, np.NaN, np.NaN, np.NaN])
+
+    i = decimation
+    while i < len(times):
+        # Check if we are descending. If so, break.
+        if altitude[i] < _output[-1][0]:
+            break
+
+        # If we have valid PTU data, calculate the dew point.
+        if temperature[i] != -273:
+            T = temperature[i]
+            RH = humidity[i]
+            DP = 243.04*(np.log(RH/100)+((17.625*T)/(243.04+T)))/(17.625-np.log(RH/100)-((17.625*T)/(243.04+T)))
+        else:
+            # Otherwise we insert NaNs, so data isn't plotted.
+            T = np.NaN
+            DP = np.NaN
+
+        # Calculate time delta between telemetry frames.
+        _time = parse(times[i])
+        _time_old = parse(times[i-decimation])
+        _delta_seconds = (_time - _time_old).total_seconds()
+
+        # Calculate the movement direction and distance, and then calculate the movement speed.
+        _movement = position_info((latitude[i], longitude[i], altitude[i]), (latitude[i-decimation], longitude[i-decimation], altitude[i-decimation]))
+        _heading = _movement['bearing']
+        _speed = (_movement['great_circle_distance']/_delta_seconds)*1.94384 # Convert to knots
+
+
+        _output.append([altitude[i], _speed, _heading, T, DP])
+
+        i += decimation
+
+    # Convert our output data into something we can process easier.
+    data_np = np.array(_output)
+    
+    if args.metric:
+        _alt = data_np[:,0]
+    else:
+        _alt = data_np[:,0]*3.28084 # Convert to feet.
+
+    _speed = data_np[:,1]
+    _direction = data_np[:,2]/10.0
+    _temp = data_np[:,3]
+    _dp = data_np[:,4]
+
+    # Produce a boolean array to limit the plotted data.
+    _data_limit = _alt < args.alt_limit
+
+    # Plot the data...
+    plt.figure()
+    plt.plot(_speed[_data_limit], _alt[_data_limit], label='Speed (kt)', color='g')
+    plt.plot(_direction[_data_limit], _alt[_data_limit], label='Direction (deg/10)', color='m')
+    plt.plot(_temp[_data_limit], _alt[_data_limit], label='Temp (deg C)', color='b')
+    plt.plot(_dp[_data_limit], _alt[_data_limit], label='DP (deg C)', color='r')
+
+    if args.metric:
+        plt.ylabel("Altitude (metres)")
+    else:
+        plt.ylabel("Altitude (feet)")
+
+    # Determine and set plot axis limits
+    _axes = plt.gca()
+    # Y limit is either a default value, or a user specified altitude.
+    _axes.set_ylim(top=args.alt_limit, bottom=0)
+
+    # X limits are based on a combination of data.
+    # The upper limit is based on the maximum speed within our altitude window
+    if args.temp_limit == None:
+        _temp_in_range= _temp[_data_limit]
+        _dp_in_range= _dp[_data_limit]
+        _min_temp = np.min(_temp_in_range[~np.isnan(_temp_in_range)])
+        _min_dp = np.min(_dp_in_range[~np.isnan(_dp_in_range)])
+        _axes.set_xlim(left=min(_min_temp, _min_dp))
+    else:
+        _axes.set_xlim(left=args.temp_limit)
+
+    plt.title("Sounding File: %s" % os.path.basename(args.filename))
+    plt.grid(which='both')
+    plt.legend(loc='upper right')
+    plt.show()
+
+
+
+