Script working for plot

api/main.py

@@ -26,7 +26,8 @@ associators = [
 # todo: ellipse conic int (analytic), SX, arc length
 localisations = [
   {"name": "Ellipse Parametric", "id": "ellipse-parametric"},
-  {"name": "Ellipsoid Parametric", "id": "ellipsoid-parametric"}
+  {"name": "Ellipsoid Parametric", "id": "ellipsoid-parametric"},
+  {"name": "Spherical Intersection", "id": "spherical-intersection"}
 ]
 adsbs = [
   {"name": "adsb.30hours.dev", "url": "adsb.30hours.dev"},

event/algorithm/localisation/EllipseParametric.py

@@ -104,7 +104,7 @@ class EllipseParametric:
               samples_intersect[i] = Geometry.ecef2lla(
                 samples_intersect[i][0], 
                 samples_intersect[i][1], 
-                0)
+                samples_intersect[i][2])
               output[target]["points"].append([
                 round(samples_intersect[i][0], 3),
                 round(samples_intersect[i][1], 3),

event/algorithm/localisation/SphericalIntersection.py

@@ -34,4 +34,14 @@ class SphericalIntersection:
 
         output = {}
 
+        # return if no detections
+        if not assoc_detections:
+          return output
+
+        for target in assoc_detections:
+
+            print(assoc_detections)
+            #nDetections = assoc_detections
+            #S = np.zeros()
+
         return output

event/event.py

@@ -12,6 +12,7 @@ import time
 import copy
 import json
 import hashlib
+import os
 
 from algorithm.associator.AdsbAssociator import AdsbAssociator
 from algorithm.localisation.EllipseParametric import EllipseParametric
@@ -19,7 +20,6 @@ from algorithm.localisation.EllipsoidParametric import EllipsoidParametric
 from algorithm.localisation.SphericalIntersection import SphericalIntersection
 from algorithm.truth.AdsbTruth import AdsbTruth
 from common.Message import Message
-
 from data.Ellipsoid import Ellipsoid
 from algorithm.geometry.Geometry import Geometry
 
@@ -195,6 +195,11 @@ async def main():
         await asyncio.sleep(1)
 
 def append_api_to_file(api_object, filename=saveFile):
+
+    if not os.path.exists(filename):
+        with open(filename, 'w') as new_file:
+            pass
+
     with open(filename, 'a') as json_file:
         json.dump(api_object, json_file)
         json_file.write('\n')

script/Dockerfile

@@ -0,0 +1,10 @@
+# use an official Python runtime as a parent image
+FROM python:3.9-slim
+
+# set the working directory to /app
+WORKDIR /app
+
+# install any needed packages specified in requirements.txt
+COPY requirements.txt /app
+RUN pip install --no-cache-dir -r requirements.txt
+

script/README.md

@@ -2,4 +2,14 @@ This folder is for post-processing scripts on API data.
 
 ## Scripts
 
-- TODO
+- **plot_accuracy.py** plots a comparison between ADS-B truth and target localisation data in ENU coordinates.
+
+## Docker Environment
+
+To avoid having to install extra libraries on the host machine (e.g. numpy, matplotlib), a Dockerfile is provided to create this environment.
+
+```bash
+sudo docker build -t 3lips-script .
+sudo docker run -it -v /opt/3lips/save:/app/save -v /opt/3lips/script:/app/script -v /opt/3lips/event/algorithm/geometry:/app/geometry 3lips-script bash
+PYTHONPATH=/app python <script> <args>
+```

script/plot_accuracy.py

@@ -1,6 +1,11 @@
 import argparse
 import json
+import sys
 from datetime import datetime
+import numpy as np
+import matplotlib.pyplot as plt
+
+from geometry.Geometry import Geometry
 
 def parse_posix_time(value):
     try:
@@ -11,23 +16,153 @@ def parse_posix_time(value):
 def parse_command_line_arguments():
     parser = argparse.ArgumentParser(description="Process command line arguments.")
     
-    parser.add_argument("json_string", type=str, help="Input JSON string")
+    parser.add_argument("json_file", type=str, help="Input JSON file path")
     parser.add_argument("target_name", type=str, help="Target name")
     parser.add_argument("--start_time", type=parse_posix_time, help="Optional start time in POSIX seconds")
     parser.add_argument("--stop_time", type=parse_posix_time, help="Optional stop time in POSIX seconds")
 
     return parser.parse_args()
 
+def interpolate_positions(timestamp_vector, truth_timestamp, truth_position):
+    # Convert lists to NumPy arrays for easier manipulation
+    truth_timestamp = np.array(truth_timestamp)
+    truth_position = np.array(truth_position)
+
+    # Interpolate positions for the new timestamp vector
+    interpolated_positions = np.zeros((len(timestamp_vector), truth_position.shape[1]))
+
+    for i in range(truth_position.shape[1]):
+        interpolated_positions[:, i] = np.interp(timestamp_vector, truth_timestamp, truth_position[:, i])
+
+    return interpolated_positions
+
 def main():
   
+    # input handling
     args = parse_command_line_arguments()
-    json_data = json.loads(args.json_string)
+    json_data = []
+    with open(args.json_file, 'r') as json_file:
+        for line in json_file:
+            try:
+                json_object = json.loads(line)
+                json_data.append(json_object)
+            except json.JSONDecodeError:
+                print(f"Error decoding JSON from line: {line}")
+    json_data = [item for item in json_data if item]
     start_time = args.start_time if args.start_time else None
     stop_time = args.stop_time if args.stop_time else None
-    print("JSON String:", json_data)
+    print("JSON String (Last Non-Empty Data):", json_data[-1])
     print("Target Name:", args.target_name)
     print("Start Time:", start_time)
     print("Stop Time:", stop_time)
 
+    # get LLA coords from first radar
+    radar4_lla = [-34.91041, 138.68924, 210]
+
+    # extract data of interest
+    server = json_data[0][0]["server"]
+    timestamp = []
+    position = {}
+    truth_timestamp = []
+    truth_position = []
+    for item in json_data:
+
+        for method in item:
+
+            if method["server"] != server:
+                continue
+
+            if start_time and method["timestamp_event"]/1000 < start_time:
+              continue
+
+            if stop_time and method["timestamp_event"]/1000 > stop_time:
+              continue
+
+            # store target data
+            method_localisation = method["localisation"]
+            if method_localisation not in position:
+              position[method_localisation] = {}
+              position[method_localisation]["timestamp"] = []
+              position[method_localisation]["detections"] = []
+            else:
+              if args.target_name in method["detections_localised"] and \
+              len(method["detections_localised"][args.target_name]["points"]) > 0:
+                position[method_localisation]["timestamp"].append(
+                  method["timestamp_event"]/1000)
+                position[method_localisation]["detections"].append(
+                  method["detections_localised"][args.target_name]["points"][0])
+                # covert to ENU
+                x, y, z = Geometry.lla2ecef(
+                  position[method_localisation]["detections"][-1][0],
+                  position[method_localisation]["detections"][-1][1],
+                  position[method_localisation]["detections"][-1][2])
+                x, y, z = Geometry.ecef2enu(x, y, z, radar4_lla[0], 
+                  radar4_lla[1], radar4_lla[2])
+                if not "detections_enu" in position[method_localisation]:
+                  position[method_localisation]["detections_enu"] = []
+                position[method_localisation]["detections_enu"].append([x, y, z])
+
+            # store truth data
+            if args.target_name in method["truth"]:
+                truth_timestamp.append(
+                  method["truth"][args.target_name]["timestamp"])
+                truth_position.append([
+                  method["truth"][args.target_name]["lat"], 
+                  method["truth"][args.target_name]["lon"], 
+                  method["truth"][args.target_name]["alt"]])
+
+            timestamp.append(method["timestamp_event"])
+
+    # remove duplicates in truth data
+    timestamp = list(dict.fromkeys(timestamp))
+    timestamp = [element/1000 for element in timestamp]
+    truth_timestamp_unique = []
+    truth_position_unique = []
+    for t, p in zip(truth_timestamp, truth_position):
+        if t not in truth_timestamp_unique:
+            truth_timestamp_unique.append(t)
+            truth_position_unique.append(p)
+    truth_timestamp = truth_timestamp_unique
+    truth_position = truth_position_unique
+
+    # resample truth to event time (position already sampled correct)
+    for i in reversed(range(len(timestamp))):
+        if timestamp[i] < min(truth_timestamp) or timestamp[i] > max(truth_timestamp):
+            del timestamp[i]
+    truth_position_resampled = interpolate_positions(
+      timestamp, truth_timestamp, truth_position)
+
+    # convert truth to ENU
+    truth_position_resampled_enu = []
+    for pos in truth_position_resampled:
+        x, y, z = Geometry.lla2ecef(pos[0], pos[1], pos[2])
+        truth_position_resampled_enu.append(
+          Geometry.ecef2enu(x, y, z, 
+          radar4_lla[0], radar4_lla[1], radar4_lla[2]))
+
+    # plot x, y, z
+    plt.figure(figsize=(5,7))
+    for i in range(3):
+        yaxis_truth = [pos[i] for pos in truth_position_resampled_enu]
+        plt.subplot(3, 1, i+1)
+        plt.plot(timestamp, yaxis_truth, label="ADS-B Truth")
+    for method in position:
+        for i in range(3):
+            yaxis_target = [pos[i] for pos in position[method]["detections_enu"]]
+            plt.subplot(3, 1, i+1)
+            plt.plot(position[method]["timestamp"], yaxis_target, 'x', label=method)
+            plt.xlabel('Timestamp')
+            if i == 0:
+                plt.ylabel('ENU X (m)')
+            if i == 1:
+                plt.ylabel('ENU Y (m)')
+            if i == 2:
+                plt.ylabel('ENU Z (m)')
+    
+    plt.subplot(3, 1, 1)
+    plt.legend()
+    plt.tight_layout()
+    plt.savefig('save/plot_accuracy.png', bbox_inches='tight', pad_inches=0.01)
+
 if __name__ == "__main__":
     main()

script/requirements.txt

@@ -0,0 +1,2 @@
+numpy==1.26.4
+matplotlib==3.8.3