diff --git a/postprocessing/Geostacking.py b/postprocessing/Geostacking.py new file mode 100644 index 0000000..216c3c5 --- /dev/null +++ b/postprocessing/Geostacking.py @@ -0,0 +1,134 @@ +import xarray as xr +import numpy as np +from scipy.interpolate import interp1d + + +def add_distance_to_xarray(data_array, slow_time_values, gps_times, gps_distances): + """ + Add distance-along-track coordinate to xarray by interpolating GPS data. + + Parameters + ---------- + data_array : xr.DataArray + 2D xarray with dimensions (sample_idx, ...) and slow_time coordinate + slow_time_values : array-like + The slow_time values corresponding to each sample_idx + gps_times : array-like + Time values from GPS data + gps_distances : array-like + Distance along track values from GPS data + + Returns + ------- + xr.DataArray + Data array with added 'distance_along_track' coordinate + """ + # Create interpolation function for distance along track + # Use bounds_error=False to extrapolate for times outside GPS range + distance_interp = interp1d( + gps_times, + gps_distances, + kind='linear', + bounds_error=False, + fill_value='extrapolate' + ) + + # Interpolate distances for each slow_time + distances = distance_interp(slow_time_values) + + # Add as a coordinate to the data array + data_with_distance = data_array.assign_coords( + distance_along_track=('pulse_idx', distances) + ) + + return data_with_distance + + +def group_and_average_by_distance(data_array, bin_size=5.0, distance_coord='distance_along_track'): + """ + Group xarray columns by distance along track and compute mean. + + Parameters + ---------- + data_array : xr.DataArray + Data array with distance_along_track coordinate + bin_size : float + Size of distance bins in meters (default: 5.0) + distance_coord : str + Name of the distance coordinate (default: 'distance_along_track') + + Returns + ------- + xr.DataArray + Averaged data grouped by distance bins + """ + # Get the distance values + distances = data_array[distance_coord].values + + # Create bins + min_dist = np.floor(distances.min()) + max_dist = np.ceil(distances.max()) + bins = np.arange(min_dist, max_dist + bin_size, bin_size) + + # Assign each sample to a bin (using bin centers as labels) + bin_indices = np.digitize(distances, bins) - 1 + bin_centers = bins[:-1] + bin_size / 2 + + # Clip to valid range + bin_indices = np.clip(bin_indices, 0, len(bin_centers) - 1) + + # Add bin assignment as a coordinate + data_with_bins = data_array.assign_coords( + distance_bin=('pulse_idx', bin_centers[bin_indices]) + ) + + # Group by distance bin and compute mean + grouped = data_with_bins.groupby('distance_bin').mean(dim='pulse_idx') + + return grouped + +def stack(data: xr.Dataset, dist_stack: float): + #""" + #Stack (average) data along the slow time axis in chunks of n_stack chirps + #All relevant coordinates (i.e. slow_time, pulse_idx) are taken as their + #minimum value in the stack. + #""" + #return data.coarsen({'pulse_idx': n_stack}, + # boundary='trim', + # coord_func='min').mean() + return group_and_average_by_distance(data, + bin_size=dist_stack, + distance_coord='distance_along_track') + + +def process_data(data_array, slow_time_values, gps_times, gps_distances, bin_size=5.0): + """ + Complete pipeline: add distance coordinate and group by distance bins. + + Parameters + ---------- + data_array : xr.DataArray + 2D xarray with dimensions (sample_idx, ...) and slow_time coordinate + slow_time_values : array-like + The slow_time values corresponding to each sample_idx + gps_times : array-like + Time values from GPS data + gps_distances : array-like + Distance along track values from GPS data + bin_size : float + Size of distance bins in meters (default: 5.0) + + Returns + ------- + xr.DataArray + Averaged data grouped by distance bins + """ + # Step 1: Add distance coordinate + data_with_distance = add_distance_to_xarray( + data_array, slow_time_values, gps_times, gps_distances + ) + + # Step 2: Group and average by distance + result = group_and_average_by_distance(data_with_distance, bin_size=bin_size) + + return result \ No newline at end of file diff --git a/postprocessing/gps_processing.py b/postprocessing/gps_processing.py new file mode 100644 index 0000000..c43b8bd --- /dev/null +++ b/postprocessing/gps_processing.py @@ -0,0 +1,59 @@ +import sys +import json +import re +import pandas as pd + +def parse_line(line): + # Regex to extract UNIX time (e.g., 1764000892.487560) + unix_time_match = re.search(r'(\d+\.\d+)', line) + if not unix_time_match: + return None + unix_time = float(unix_time_match.group(1)) + + # Regex to extract the entire JSON object starting with "class":"TPV" + json_match = re.search(r'\{.*?"class":"TPV".*?\}', line, re.DOTALL) + if not json_match: + return None + json_str = json_match.group(0) + + try: + data = json.loads(json_str) + except json.JSONDecodeError: + return None + + # Extract the required fields + lat = data.get('lat') + lon = data.get('lon') + alt = data.get('alt') + altMSL = data.get('altMSL') + ecefx = data.get('ecefx') + ecefy = data.get('ecefy') + ecefz = data.get('ecefz') + + return { + 'unix_time': unix_time, + 'lat': lat, + 'lon': lon, + 'alt': alt, + 'altMSL': altMSL, + 'ecefx': ecefx, + 'ecefy': ecefy, + 'ecefz': ecefz, + } + +def extract_locations_from_gpspipe(lines): + results = [] + + for line in lines.split("\n"): + if 'TPV' in line: + result = parse_line(line) + if result: + results.append(result) + + # Convert list of dicts to DataFrame + df = pd.DataFrame(results) + + # Uncomment below to save as CSV + # df.to_csv('output.csv', index=False) + return df + \ No newline at end of file diff --git a/postprocessing/notebooks/Glass Field Processing.ipynb b/postprocessing/notebooks/Glass Field Processing.ipynb new file mode 100644 index 0000000..c8717e6 --- /dev/null +++ b/postprocessing/notebooks/Glass Field Processing.ipynb @@ -0,0 +1,841 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "9b3ac471-a786-4ccd-bdab-b0a7fdde8467", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6273bd01-f8a9-4b09-932a-52372d7b6678", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from dask.distributed import Client, LocalCluster\n", + "\n", + "client = Client() # Note that `memory_limit` is the limit **per worker**.\n", + "# n_workers=4,\n", + "# threads_per_worker=1,\n", + "# memory_limit='3GB'\n", + "client # If you click the dashboard link in the output, you can monitor real-time progress and get other cool visualizations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "770f50ea-f019-46de-a38c-aafcd483e9c3", + "metadata": {}, + "outputs": [], + "source": [ + "import copy\n", + "import sys\n", + "import xarray as xr\n", + "import numpy as np\n", + "import dask.array as da\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import hvplot.xarray\n", + "import hvplot.pandas\n", + "import holoviews as hv\n", + "from holoviews import opts\n", + "import scipy.constants\n", + "\n", + "sys.path.append(\"..\")\n", + "import processing_dask as pr\n", + "import gps_processing as gps\n", + "import Geostacking as geostack\n", + "import plot_dask\n", + "import processing as old_processing\n", + "\n", + "# for parsing GPS and stdout log files\n", + "import re\n", + "import pandas as pd\n", + "\n", + "sys.path.append(\"../../preprocessing/\")\n", + "from generate_chirp import generate_chirp" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a0f40923", + "metadata": {}, + "outputs": [], + "source": [ + "plt.rcParams.update({'font.size': 16})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9262b908-c466-47b1-b56f-5e1647540715", + "metadata": {}, + "outputs": [], + "source": [ + "def parse_stdout_times(line):\n", + " # Extract the timestamp from lines that contain specific phrases\n", + " if '[START] Beginning main loop' in line:\n", + " pattern = r'(\\d+\\.\\d+)'\n", + " match = re.search(pattern, line)\n", + " if match:\n", + " return float(match.group(1))\n", + " \n", + " elif '[TX] Closing file' in line:\n", + " pattern = r'(\\d+\\.\\d+)'\n", + " match = re.search(pattern, line)\n", + " if match:\n", + " return float(match.group(1))\n", + "\n", + " return None" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b1749d35-a1f7-49d4-bf07-d02dd6d5b76f", + "metadata": {}, + "outputs": [], + "source": [ + "def start_and_stop_from_log(logstring):\n", + " lines = logstring.split(\"\\n\")\n", + " results = []\n", + " for line in lines:\n", + " parsed = parse_stdout_times(line.strip())\n", + " if parsed:\n", + " results.append(parsed)\n", + " \n", + " return {\"start\": results[0], \"stop\": results[1]}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fa2332a9-1807-4eb5-bedb-d423590e6078", + "metadata": {}, + "outputs": [], + "source": [ + "def get_distance_along_track(df):\n", + " displacements = np.linalg.norm(df.iloc[1:,:].loc[:, [\"ecefx\", \"ecefy\", \"ecefz\"]].to_numpy() - \\\n", + " df.iloc[:-1,:].loc[:, [\"ecefx\", \"ecefy\", \"ecefz\"]].to_numpy(), axis=1)\n", + " return np.cumsum(displacements)" + ] + }, + { + "cell_type": "markdown", + "id": "e70b6b7e-dcd8-4be1-b94c-a3cad899cd73", + "metadata": {}, + "source": [ + "### Open and resave file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7793fbfe-8a8c-47cc-b631-0534e3e1c19a", + "metadata": {}, + "outputs": [], + "source": [ + "# file path to data and configs\n", + "\n", + "prefix = \"/Volumes/Extreme SSD/orca_paper/20240226_105437\" # Anna loopback spectrogram, no phase dithering\n", + "prefix = \"/home/thomas/Documents/StanfordGrad/RadioGlaciology/drone/radar_data/20230723-summit-day3-2start/20230723_104248\"\n", + "prefix = \"/Users/thatch/Projects/SORA/before-and-after-test/20251002_175426\"\n", + "prefix = \"/Volumes/sora-sshfs/20251127_142119\" # first shot - only went to 1km\n", + "prefix = \"/Volumes/sora-sshfs/20251127_143302\" # second shot, longer rx window, 4km\n", + "prefix = \"/Volumes/sora-sshfs/20251127_143818\" # thrid shot - Mt Rossman\n", + "prefix = \"/Volumes/sora-sshfs/20251127_144118\" # fourth shot - straight down - 1 sec or so\n", + "prefix = \"/Volumes/sora-sshfs/20251127_144528\" # fifth shot - straight down - 10 sec or so\n", + "prefix = \"/Users/thatch/Projects/SORA/data/20251127_144528\"\n", + "\n", + "prefix = \"/Volumes/sora-sshfs/20251130_124406\" # Nov 30 - first LO - OOPS! just attenuators! no cabl\n", + "\n", + "#path = \"/Volumes/sora-sshfs/\"\n", + "#path = \"/tmp/SORA-1130/\"\n", + "test = \"20251130_125618\" # first real LO at CI with attns and cable!! looks pretty good!\n", + "test = \"20251130_155708\" # lpda shot outside UG after SSD issue\n", + "#test = \"20251130_160010\" # lpda dragged by ~10 meters\n", + "#test = \"20251130_175353\" # | | pol 10 s\n", + "#test = \"20251130_175636\" # x pol 10 s\n", + "#test = \"20251130_175854\" # -- -- pol 10 s\n", + "#test = \"20251130_181559\" #10 sec loop;back\n", + "#test = \"20251130_183246\" # 10 sec -- -- pol\n", + "#test = \"20251130_183459\" # 10 sec | | pol\n", + "#test = \"20251130_234658\" # ref cable lo, 20db attn, rx only\n", + "#test = \"20251201_000404\" # suspected bad cable, white tape\n", + "#test = \"20251201_000952\" # suspected good cable, no tape\n", + "\n", + "#test = \"20251201_004608\" # suspected good\n", + "#test = \"20251201_005419\" # suspected bad\n", + "#test = \"20251201_005945\" # reference cable\n", + "\n", + "#test = \"20251201_011014\" # all up test in camp\n", + "#test = \"20251201_011420\" # and again\n", + "#test = \"20251201_011703\" # and agaibn\n", + "#test = \"20251201_011829\" # lower tx and rx gain\n", + "#test = \"20251201_012005\" # even less tx\n", + "#test = \"20251201_012130\" # even less tx gain\n", + "\n", + "#test = \"20251201_125531\" # low gain xpol\n", + "#test = \"20251201_125646\" # reg gain xpol\n", + "#test = \"20251201_125729\" # reg gain same pol\n", + "\n", + "#test = \"20251201_144352\" # rfspace antennas\n", + "\n", + "# 4 rfspace antenna tests\n", + "#test = \"20251201_155939\" # 450, rtk on\n", + "#test = \"20251201_160032\" # 390, rtk on\n", + "#test = \"20251201_160451\" # 390 rtk off\n", + "#test = \"20251201_160533\" # 450, rtk off\n", + "\n", + "#test = \"20251201_172432\" # 450, rtk off, lpdas\n", + "#test = \"20251201_172921\" # \"\" +5 tx gain\n", + "\n", + "path = \"/Volumes/CI1/\"\n", + "test = \"20251201_190207\" # static shot at 1km upstream\n", + "#test = \"20251201_191042\" # first traveerse, 200m\n", + "#test = \"20251201_192052\" # second traverse, 200m \n", + "#test = \"20251201_202215\" # 1km traverse upstream, ending at camp\n", + "\n", + "#dec 2\n", + "# parallel plane pol, tx gain50\n", + "# inplane pol, tx gain59\n", + "test = \"20251202_124842\" # inplane pol, tx gain50 - start of test\n", + "#test = \"20251202_125159\" # inplane, txgain 45\n", + "#test= \"20251202_125429\" # inplane, txgain 40\n", + "#test = \"20251202_125613\" # \", txgain 35\n", + "#test = \"20251202_125843\" # \", txgain 32\n", + "#test = \"20251202_130023\" # \", 27 db\n", + "#test = \"20251202_131542\" # same, farther from lehman\n", + "#test = \"20251202_131734\" # 22 db tx\n", + "test = \"20251202_131932\" # \"\", lower rx gain\n", + "test = \"20251202_132228\" # lower rx gain - still clipping\n", + "test = \"20251202_132409\" # 10 db rx gain\n", + "#test = \"20251202_132542\" # same ^\n", + "#test = \"20251202_134647\" # same but with rfspace antennas\n", + "#test = \"20251202_134809\" # higher tx rx gain\n", + "\n", + "test = \"20251202_180847\" # shot 3k afternoon dec 2\n", + "#test = \"20251202_181055\" # same +3db tx\n", + "#path = \"/tmp/SORA-1202/\"\n", + "#test = \"20251202_183128\" # ~1km traverse from shot 3k\n", + "#test = \"20251202_183653\" # station at line upstream of camp\n", + "#test = \"20251202_183833\" # same but +3db \n", + "test = \"20251202_190714\" # long transect -- sea water reflection!!!!\n", + "#test = \"20251202_190820\" # statijnary at end of tsect sea water!!!!!\n", + "#test = \"20251202_999999\"\n", + "\n", + "#path = \"/tmp/SORA-1203/\"\n", + "#test = \"20251203_182424\" # good conf from yesterday\n", + "#test =\"20251203_182518\" # same but 5us chrip\n", + "\n", + "#test = \"20251203_190524\"\n", + "#test = \"20251203_194038\" # long downstream transect dec 3, 7 kph\n", + "#test = \"20251203_195317\" # upstream transect to camp - 12 kph\n", + "\n", + "# rf space antenna test - hopefully less clipping\n", + "#test = \"20251203_210905\" # reference with lpdas\n", + "#test = \"20251203_211322\" # same conf with rf spaces tx gain 56\n", + "#test = \"20251203_211458\" # rf space, txgain 46\n", + "#test = \"20251203_211645\" # \", txgan 36\n", + "#test = \"20251203_211834\" # \"\" rxgain 26\n", + "#test = \"20251203_212053\" # rxgain 16\n", + "#test = \"20251203_212350\" # \"\", parallel pol\n", + "#test= \"20251203_212514\" # \"\", xpol\n", + "#test = \"20251203_212658\" # still xpol, big gains\n", + "#test = \"20251203_212821\" # parallel plane pol, big gains\n", + "\n", + "# dec 4\n", + "#test = \"20251204_212545\"\n", + "test = \"20251204_215913\" # 5km down das line\n", + "#test = \"20251204_221135\" # 2km side leg\n", + "#test = \"20251204_224207\" # 4 km upstream traverse on smart solo like (near smart solo 4016)\n", + "\n", + "\n", + "# Dec 5\n", + "#test = \"20251205_154838\" # 4 km down stream, 12 kph\n", + "#test = \"20251205_161625\" # 2km up starboard wing, 12 kph then 7 kph\n", + "#test = \"20251205_162708\" # diagonal down across GL, 12? kph\n", + "#test = \"20251205_164609\" # 3 km up stryde line\n", + "\n", + "# rx gain tests at camp\n", + "#test =\"20251205_170039\" # reference = tx gain 56 db, rx gain 36 db\n", + "#test = \"20251205_170508\" # 26 db\n", + "#test = \"20251205_170900\" # rx gain 166 !!!!!\n", + "#test = \"20251205_171306\" # rx gain 16 # 10 db SNR\n", + "#test = \"20251205_171525\" # rx gain 6\n", + "#test = \"20251205_171810\" # rx gain 8\n", + "\n", + "# Dec 6 - low rx gain!\n", + "# back 9\n", + "#test = \"20251206_134321\"\n", + "#test = \"20251206_140706\" # back 4\n", + "#test = \"20251206_142331\"\n", + "\n", + "# left wing\n", + "#test = \"20251206_161447\" # western boundary\n", + "#test = \"20251206_162356\" # ocean segment\n", + "#test = \"20251206_163720\" # diagonal upstream western GL\n", + "\n", + "# dec 7 left wing\n", + "#test = \"20251207_153056\" # down das line\n", + "#test = \"20251207_160617\" # up western boundary\n", + "#test = \"20251207_162558\" # down the dogleg\n", + "#test = \"20251207_164411\" # up the stride line\n", + "\n", + "# dec 7 opres test - no rfpa., 350 mhz\n", + "#test = \"20251207_170138\"\n", + "#test = \"20251207_170444\" # 70 txgain\n", + "#test = \"20251207_170748\" # + 4db rx gain\n", + "#test = \"20251207_171023\" # + 3 db rx\n", + "\n", + "#test = \"20251207_214057\" # test at start of run, rx gain36\n", + "#test = \"20251207_214310\" # rx gain26\n", + "\n", + "# Dec 7-8, opres overnight\n", + "#test = \"20251208_133603\" # great bed return! no real brightness variation... phase could be interesting!\n", + "\n", + "# dec 8\n", + "#test = \"20251208_155425\" # down das line, clipped tucker\n", + "#test = \"20251208_162754\" # starboard star offline segment\n", + "\n", + "# afternoon figure 8 track\n", + "#test = \"20251208_214018\" # igloo to ug04 (?)\n", + "#test = \"20251208_215131\" # gnss to right angle upstream\n", + "#test = \"20251208_220414\" # crossing line\n", + "#test = \"20251208_221734\" # upstream on west side\n", + "#test = \"20251208_222354\" #igloo to ug03\n", + "\n", + "# dec 9 !\n", + "#test = \"20251209_160157\" # ug03 to ug01, 350\n", + "#test = \"20251209_162232\" # ug01 to ug08 (out of box), 350\n", + "#test = \"20251209_162544\" # check at 450 at ug08, txgain 57\n", + "\n", + "# gl laps\n", + "\n", + "\n", + "#test = \"20251209_200750\" # leg 3 - down at 350, 58 db\n", + "\n", + "prefix = path + test\n", + "\n", + "# resave data as zarr for dask processing\n", + "zarr_base_location=\"/Users/thatch/Projects/SORA/before-and-after-test/test_tmp_zarr_cache/\"\n", + "zarr_path = pr.save_radar_data_to_zarr(prefix, zarr_base_location=zarr_base_location)\n", + "\n", + "# open zarr file, adjust chunk size to be 10 MB - 1 GB based on sample rate/bit depth\n", + "raw = xr.open_zarr(zarr_path, chunks={\"pulse_idx\": 1000})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c60de6cf-0145-436a-a612-b83357ddf1f8", + "metadata": {}, + "outputs": [], + "source": [ + "gpsdf = gps.extract_locations_from_gpspipe( raw.gpspipe_log )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "75130312-94cb-4876-8f42-7e1771f63490", + "metadata": {}, + "outputs": [], + "source": [ + "res = start_and_stop_from_log(raw.stdout_log)\n", + "res[\"stop\"] - res[\"start\"]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a55be0ee-f01d-4acb-8f29-7fe85482719a", + "metadata": {}, + "outputs": [], + "source": [ + "(raw.slow_time[-1] - raw.slow_time[1]).compute().values" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0dcaeffe-afca-43ee-916e-f833e8e2cd2d", + "metadata": {}, + "outputs": [], + "source": [ + "gpsdf[\"slow_time\"] = gpsdf[\"unix_time\"] - res[\"start\"]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "00cb998a-68ae-4aab-bd9c-e4d6a9a6b240", + "metadata": {}, + "outputs": [], + "source": [ + "georaw = geostack.add_distance_to_xarray(raw, \n", + " raw.slow_time.compute(), \n", + " gpsdf[\"slow_time\"].iloc[1:],\n", + " get_distance_along_track(gpsdf))" + ] + }, + { + "cell_type": "markdown", + "id": "6fb45450-15bd-4cf2-8311-709658955830", + "metadata": {}, + "source": [ + "### Enter processing parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "38423415-95ee-450f-a1d6-7a0fe84a3e58", + "metadata": {}, + "outputs": [], + "source": [ + "#zero_sample_idx = 36 # X310, fs = 20 MHz\n", + "#zero_sample_idx = 63 # X310, fs = 50 MHz\n", + "zero_sample_idx = 159 # B205mini, fs = 56 MHz\n", + "#zero_sample_idx = 166 # B205mini, fs = 20 MHz\n", + "\n", + "nstack = 1000 # number of pulses to stack\n", + "lenstack = 5 # number of meters to stack over\n", + "\n", + "modify_rx_window = False # set to true if you want to window the reference chirp only on receive, false uses ref chirp as transmitted in config file\n", + "rx_window = \"blackman\" # what you want to change the rx window to if modify_rx_window is true\n", + "\n", + "dielectric_constant = 3.17# ice (air = 1, 66% velocity coax = 2.2957)\n", + "#dielectric_constant = 2.2957 # COAX (air = 1, 66% velocity coax = 2.2957)\n", + "#dielectric_constant = 1.0\n", + "sig_speed = scipy.constants.c / np.sqrt(dielectric_constant)" + ] + }, + { + "cell_type": "markdown", + "id": "b6ed8d87-89fa-4368-b5ca-ac5b6e6519f3", + "metadata": {}, + "source": [ + "### Generate reference chirp" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "707ddadf-7434-4cb4-906a-112b2c2ac3f0", + "metadata": {}, + "outputs": [], + "source": [ + "if modify_rx_window:\n", + " config = copy.deepcopy(raw.config)\n", + " config['GENERATE']['window'] = rx_window\n", + "else:\n", + " config = raw.config\n", + "\n", + "chirp_ts, ref_chirp = generate_chirp(config)" + ] + }, + { + "cell_type": "markdown", + "id": "4316fce6-9d32-4f90-99ee-c9aef2cc3a51", + "metadata": {}, + "source": [ + "### View raw pulse in time domain to check for clipping" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19228216-c7f0-4f7c-a5c2-d63dc4973942", + "metadata": {}, + "outputs": [], + "source": [ + "single_pulse_raw = raw.radar_data[{'pulse_idx': 100}].compute()\n", + "plot1 = np.real(single_pulse_raw).hvplot.line(x='fast_time', color='red') * np.imag(single_pulse_raw).hvplot.line(x='fast_time')\n", + "\n", + "plot1 = plot1.opts(xlabel='Fast Time (s)', ylabel='Raw Amplitude', title=f\"Pulse 100 from {test}\")\n", + "plot1" + ] + }, + { + "cell_type": "markdown", + "id": "9b85eb39-f0f8-4329-83ce-456af3345685", + "metadata": {}, + "source": [ + "### Clean and stack data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe97ec17-ca91-4532-bfb7-c66f6384cfa9", + "metadata": {}, + "outputs": [], + "source": [ + "geostacked = pr.fill_errors(georaw, error_fill_value=0.0) # fill receiver errors with 0s\n", + "geostacked = geostack.stack(geostacked, lenstack)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fbe3f084-f365-4ff0-b9cf-8b0c070087c2", + "metadata": {}, + "outputs": [], + "source": [ + "stacked = pr.fill_errors(raw, error_fill_value=0.0) # fill receiver errors with 0s\n", + "stacked = pr.stack(stacked, nstack) # stack " + ] + }, + { + "cell_type": "markdown", + "id": "536a830d-9022-46ef-86c2-38578ce27a38", + "metadata": {}, + "source": [ + "### Pulse compress data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3affa606-7ff1-4a0f-b919-acb0cbe591a3", + "metadata": {}, + "outputs": [], + "source": [ + "compressed = pr.pulse_compress(stacked, ref_chirp,\n", + " fs=stacked.config['GENERATE']['sample_rate'],\n", + " zero_sample_idx=zero_sample_idx,\n", + " signal_speed=sig_speed)\n", + "\n", + "geocompressed = pr.pulse_compress(geostacked, ref_chirp,\n", + " fs=stacked.config['GENERATE']['sample_rate'],\n", + " zero_sample_idx=zero_sample_idx,\n", + " signal_speed=sig_speed)\n", + "\n", + "compressed_power = xr.apply_ufunc(\n", + " lambda x: 20*np.log10(np.abs(x)),\n", + " compressed,\n", + " dask=\"parallelized\"\n", + ")\n", + "\n", + "geocompressed_power = xr.apply_ufunc(\n", + " lambda x: 20*np.log10(np.abs(x)),\n", + " geocompressed,\n", + " dask=\"parallelized\"\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a669b2b1-b0e9-4d50-8a69-9b9c4fbc8a9c", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "0aa6043e-c392-437f-89bf-7a104240ec36", + "metadata": {}, + "source": [ + "### View 1D pulse compressed data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5347ab15-9815-4564-be5c-e4331947a102", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "plot1D = np.mean(compressed_power.radar_data, axis=0).hvplot.line(label=\"Mean of All Pulses\")\n", + "plot1D = plot1D * compressed_power.radar_data[0,:].hvplot.line(label=\"First Pulse\")\n", + "plot1D = plot1D * compressed_power.radar_data[-1,:].hvplot.line(label=\"Last Pulse\")\n", + "\n", + "# relevant options: xlim(-80,1000)\n", + "\n", + "plot1D = plot1D.opts(xlabel='Reflection Distance (m)', ylabel='Return Power (dB)')\n", + "plot1D = plot1D.opts(title=test+\" (ϵ=\"+str(dielectric_constant)+\")\")\n", + "plot1D.opts(xlim=(-50,2300), ylim=(-120, -25), show_grid=True)" + ] + }, + { + "cell_type": "markdown", + "id": "1981eed7-0a04-4013-ae25-503009b1d7f2", + "metadata": {}, + "source": [ + "### View 2D pulse compressed data (radargram)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "52502998-2657-489a-b093-31ca3263dda6", + "metadata": {}, + "outputs": [], + "source": [ + "# USING HOLOVIEWS (sometimes breaks)\n", + "plot2D = compressed_power.swap_dims({'pulse_idx': 'slow_time', 'travel_time': 'reflection_distance'}).hvplot.quadmesh(x='slow_time', cmap='inferno', flip_yaxis=True)\n", + "# relevant options: ylim=(100,-50), clim=(-90,-40)\n", + "\n", + "plot2D.opts(xlabel='Slow Time (s)', ylabel='Depth (m)', clabel='Return Power (dB)')\n", + "plot2D.opts(ylim=(-10, 2000), clim=(-115, -70))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fd5a7769-d01a-442d-88dd-697e01088ba4", + "metadata": {}, + "outputs": [], + "source": [ + "# Select data at reflection_distance ≈ 1616\n", + "depth_m = 1590\n", + "\n", + "line_data = compressed_power.swap_dims({\n", + " 'pulse_idx': 'slow_time', \n", + " 'travel_time': 'reflection_distance'\n", + "}).sel(reflection_distance=depth_m, method='nearest')\n", + "\n", + "# Plot as a line\n", + "plot1D = line_data.hvplot.line(x='slow_time')\n", + "plot1D.opts(xlabel='Slow Time (s)', ylabel='Return Power (dB)')\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2fd52a40-ca8d-4ff1-9533-3d436cc1d5f0", + "metadata": {}, + "outputs": [], + "source": [ + "# USING HOLOVIEWS (sometimes breaks)\n", + "geoplot2D = geocompressed_power.swap_dims({'travel_time': 'reflection_distance'}).hvplot\\\n", + " .quadmesh(x='distance_bin', cmap='inferno', flip_yaxis=True, width=1000)\n", + "# relevant options: ylim=(100,-50), clim=(-90,-40)\n", + "\n", + "geoplot2D.opts(xlabel='Distance along track (m)', ylabel='Depth (m)', clabel='Return Power (dB)')\n", + "geoplot2D.opts(ylim=(-10, 2000), clim=(-100, -65))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30f9a4f0-38ec-4dd1-abad-be7a6128f2ae", + "metadata": {}, + "outputs": [], + "source": [ + "camp_latitude = -79.58937\n", + "(gpsdf.loc[:, [\"slow_time\", \"lat\"]].hvplot.line(x='slow_time') * hv.HLine(camp_latitude) + \\\n", + " gpsdf.loc[:, [\"slow_time\", \"lon\"]].hvplot.line(x='slow_time', color='red')).cols(1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "306dc22e-7ae4-4d59-a84a-969f98e5f305", + "metadata": {}, + "outputs": [], + "source": [ + "# USING MATPLOTLIB (sometimes takes a while)\n", + "fig, ax = plt.subplots(1,1, figsize=(10,6), facecolor='white')\n", + "\n", + "p = ax.pcolormesh(compressed_power.slow_time, \\\n", + " compressed_power.reflection_distance, \\\n", + " compressed_power.radar_data.transpose(), \\\n", + " shading='auto', cmap='inferno',\\\n", + " # vmin = -140, vmax=-60\n", + " )\n", + "ax.invert_yaxis()\n", + "clb = fig.colorbar(p, ax=ax)\n", + "clb.set_label('Return Power (dB)')\n", + "ax.set_xlabel('Slow Time (s)')\n", + "ax.set_ylabel('Distance to Reflector (m)')\n", + "# relevant options: ax.set_ylim=(100,-50), ax.set_xlim=(0, 1), vmin=-90, vmax=40\n", + "ax.set_ylim(2000, -50)" + ] + }, + { + "cell_type": "markdown", + "id": "bbf58a46-2217-48d2-9c19-60021159567d", + "metadata": {}, + "source": [ + "### View spectrogram of stacked data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4e91cd22-855c-444f-842e-cacae16279e9", + "metadata": {}, + "outputs": [], + "source": [ + "inpt = raw\n", + "inpt[\"radar_data\"].shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c761f570-8984-4028-9d4e-02778deb3a28", + "metadata": {}, + "outputs": [], + "source": [ + "num_presums = raw.attrs[\"config\"][\"CHIRP\"][\"num_presums\"]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "94106b3f-7473-4697-aa00-2b44772f4703", + "metadata": {}, + "outputs": [], + "source": [ + "# data = stacked[\"radar_data\"].to_numpy()\n", + "n = 10\n", + "normalize = True\n", + "\n", + "pulse = pr.stack(inpt, n)[{'pulse_idx':0}][\"radar_data\"].to_numpy()\n", + "\n", + "f, t, S = scipy.signal.spectrogram(\n", + " pulse,\n", + " fs=raw.attrs[\"config\"][\"GENERATE\"][\"sample_rate\"],\n", + " window='flattop',\n", + " nperseg=128,\n", + " noverlap=64,\n", + " scaling='density', mode='psd',\n", + " return_onesided=False\n", + ")\n", + "\n", + "if normalize:\n", + " S /= np.max(S)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1629b3d5-c7b3-46eb-987b-fee76b1c62c5", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(facecolor='white', figsize=(10,6))\n", + "freq_mhz = (np.fft.fftshift(f) + raw.attrs['config']['RF0']['freq']) / 1e6\n", + "pcm = ax.pcolormesh(t, freq_mhz, 10*np.log10(np.abs(np.fft.fftshift(S, axes=0))), shading='nearest') # vmin=-420, vmax=-200\n", + "clb = fig.colorbar(pcm, ax=ax)\n", + "clb.set_label('Power [dB]')\n", + "ax.set_xlabel('Time [s]')\n", + "ax.set_ylabel('Frequency [MHz]')\n", + "#ax.set_title(f\"Spectrogram of received data with n_stack={n}\")\n", + "ax.text(0, 1.05, prefix.split(\"/\")[-1] + \"\\n\" + f\"n_stack * num_presums = {n * num_presums}\", horizontalalignment='left', verticalalignment='center', transform=ax.transAxes, fontdict={'size': 12})\n", + "fig.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "3954e602-4fb8-41cc-a627-6d3a116bf0f1", + "metadata": {}, + "source": [ + "### View Power Spectrum of All Received Data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "48418de1-e64f-4597-b40a-1ea8ee85e2f6", + "metadata": {}, + "outputs": [], + "source": [ + "single_stack = pr.stack(raw, raw.radar_data.shape[1])\n", + "\n", + "data_rx_fft = da.fft.fft(raw.radar_data, axis=0) / raw.radar_data.shape[0]\n", + "stacked_fft = da.fft.fft(stacked.radar_data, axis=0) / stacked.radar_data.shape[0]\n", + "full_fft = da.fft.fft(single_stack.radar_data, axis=0) / single_stack.radar_data.shape[0]\n", + "\n", + "data_rx_fft_pwr = 20*da.log10(da.abs(data_rx_fft))\n", + "stacked_fft_pwr = 20*da.log10(da.abs(stacked_fft))\n", + "full_fft_pwr = 20*da.log10(da.abs(full_fft))\n", + "\n", + "#data_rx_fft_pwr.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "263ebb96-9e43-4b86-99a5-533e24aa7d10", + "metadata": {}, + "outputs": [], + "source": [ + "# fig, axs = plt.subplots(2,1)\n", + "fig, axs = plt.subplots(facecolor='white', figsize=(10,6))\n", + "freqs = np.fft.fftshift(np.fft.fftfreq(data_rx_fft_pwr.shape[0], d=1/raw.config['GENERATE']['sample_rate']))\n", + "axs.plot(freqs/1e6, np.fft.fftshift(data_rx_fft_pwr[:,0]), label='Single Pulse')\n", + "axs.plot(freqs/1e6, np.fft.fftshift(stacked_fft_pwr[:,0]), label='Single Stack')\n", + "axs.plot(freqs/1e6, np.fft.fftshift(full_fft_pwr[:,0]), label='Full File')\n", + "axs.set_xlabel('Frequency [MHz]')\n", + "axs.set_ylabel('Power [dB]')\n", + "axs.set_title('Spectrum -- Power')\n", + "axs.grid()\n", + "axs.legend()\n", + "\n", + "# axs[1].plot(freqs/1e6, np.fft.fftshift(np.angle(data_rx_fft[:,0])))\n", + "# axs[1].plot(freqs/1e6, np.fft.fftshift(np.angle(stacked_fft[:,0])))\n", + "# axs[1].plot(freqs/1e6, np.fft.fftshift(np.angle(full_fft[:,0])))\n", + "# axs[1].set_xlabel('Frequency [MHz]')\n", + "# axs[1].set_ylabel('Phase [rad]')\n", + "# axs[1].set_title('Spectrum -- Phase')\n", + "# axs[1].grid()\n", + "# fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "93362d0a-e023-4b39-bc48-00d68319a5b6", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fdf7d33d-c8b0-403f-9fc5-4d5ecc5db65d", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/postprocessing/processing_dask.py b/postprocessing/processing_dask.py index 33c6456..606b634 100644 --- a/postprocessing/processing_dask.py +++ b/postprocessing/processing_dask.py @@ -91,6 +91,7 @@ def save_radar_data_to_zarr(prefix, skip_if_cached=True, zarr_base_location=None # Build filenames from prefix rx_samps_file = prefix + "_rx_samps.bin" log_file = prefix + "_uhd_stdout.log" + gpsd_file = prefix + "_gpspipe_stdout.log" # # Data loading @@ -137,6 +138,18 @@ def save_radar_data_to_zarr(prefix, skip_if_cached=True, zarr_base_location=None raise FileNotFoundError( f"Log file not found: {log_file}. If a log file is not required, set log_required=False") + gpsdlog = None + if os.path.exists(gpsd_file): + with open(gpsd_file, 'r') as gpsd_f: + gpsdlog = gpsd_f.read() + else: + print("No GPS log file found!!") + gpsdlog = None + #if log_required: + # raise FileNotFoundError( + # f"Log file not found: {log_file}. If a log file is not required, set log_required=False") + + # Save radar_data, slow_time, and fs to an xarray datarray data = xr.Dataset( data_vars={ @@ -151,6 +164,7 @@ def save_radar_data_to_zarr(prefix, skip_if_cached=True, zarr_base_location=None attrs={ "config": config, "stdout_log": log, + "gpspipe_log": gpsdlog, "prefix": prefix, "basename": basename } diff --git a/run.py b/run.py index a0e64b2..23f92a2 100644 --- a/run.py +++ b/run.py @@ -125,6 +125,9 @@ def run(self): if not self.setup_complete: raise Exception("Must call setup() before calling run(). If setup() does not complete successfully, you cannot call run().") + self.gpspipe_process = subprocess.Popen(["/usr/bin/gpspipe", "--json", "-uu", "--output", "gpspipe_stdout.log"]) + self.gpxlogger_process = subprocess.Popen(["/usr/bin/gpxlogger", "--reconnect", "-m2", "-f", "track.gpx", "-e", "sockets"]) + time.sleep(5) self.uhd_process = subprocess.Popen(["./radar", self.yaml_filename], stdout=subprocess.PIPE, bufsize=1, close_fds=True, text=True, cwd="sdr/build") self.uhd_output_reader_thread = threading.Thread(target=self.process_usrp_output, args=(self.uhd_process.stdout, open('uhd_stdout.log', 'w'), self.output_to_stdout)) self.uhd_output_reader_thread.daemon = True # thread dies with the program @@ -163,6 +166,9 @@ def stop(self, timeout = 10): was_force_killed = True self.is_running = False + self.gpspipe_process.kill() + self.gpxlogger_process.kill() + self.uhd_output_reader_thread.join() # If necessary, concatenate data files into a single file @@ -174,7 +180,11 @@ def stop(self, timeout = 10): # Save output print("Copying data files...") - file_prefix = save_data(self.yaml_filename, alternative_rx_samps_loc=alternative_rx_samps_loc, num_files=self.file_queue_size, extra_files={"uhd_stdout.log": "uhd_stdout.log"}) + file_prefix = save_data(self.yaml_filename, alternative_rx_samps_loc=alternative_rx_samps_loc, num_files=self.file_queue_size, + extra_files={"uhd_stdout.log": "uhd_stdout.log", + "gpspipe_stdout.log": "gpspipe_stdout.log", + "track.gpx": "track.gpx" + }) print("Finished copying data.") self.output_file = None