diff --git a/tutorials/mvp-1-geodesy/NB3-inspect-prepare-visualize-observations.ipynb b/tutorials/mvp-1-geodesy/NB3-inspect-prepare-visualize-observations.ipynb new file mode 100644 index 0000000..2c851ce --- /dev/null +++ b/tutorials/mvp-1-geodesy/NB3-inspect-prepare-visualize-observations.ipynb @@ -0,0 +1,1091 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "2abe5d1d-4a2f-4c96-a78c-6c3e5f9ac750", + "metadata": {}, + "source": [ + "# Notebook 3 - Inspecting, Preparing and Visualizing Observations" + ] + }, + { + "cell_type": "markdown", + "id": "35320cb4-aed2-4420-8f9c-b3deb365e792", + "metadata": {}, + "source": [ + "**Version:** 1.0 | **Last updated:** 2026-07-24 | \n", + "\n", + "**Author:** Eshanta Mishra | **Author institution:** EarthScope Consortium\n", + "\n", + "**Estimated Time:** ~ 30 minutes | **Pathway:** MVP1\n", + "\n", + "**License:** CC-BY-4.0" + ] + }, + { + "cell_type": "markdown", + "id": "b9b1d1c3-569f-459a-8595-e1ba758e4bc7", + "metadata": {}, + "source": [ + "## Introduction\n", + "\n", + "**What this notebook does:** It retrieves GNSS observations for several stations using the EarthScope SDK, then works through the checks you should run before using the data such as how complete it is, how it is sampled, where it is missing, and which values look wrong.\n", + "\n", + "**Why it is useful:** Raw GNSS observations are never uniformly complete. Satellites rise and set, receivers go offline, and some signals are never tracked. Positions, displacements, and reflector heights are all computed from these observations, so gaps and bad values carry through into the result. Checking coverage and data quality first tells you what the later numbers can support.\n", + "\n", + "**What you will accomplish:** By the end you will have mesured the time coverage and sampling interval of a real multi-station dataset, located its gaps, quantified per-satellite availability, aggregated the observations into time bins, flagged questionable signal-to-noise values without discarding them, and produced first-look plots that let you tell normal data from unusual data.\n", + "\n", + "---\n", + "\n", + "### Prerequisites\n", + "\n", + "* [ ] Have completed [Notebook 1 - Accessing GNSS Observations with the EarthScope SDK](NB1-access-gnss-via-SDK-checkpoint.ipynb])\n", + "* [ ] Be familiar with basic python and Polar dataframes.\n", + "\n", + "---\n", + "\n", + "## GeoLab Compute Resources\n", + "\n", + "| Setting | Recommended |\n", + "|---|---|\n", + "| Image | GeoLab (default image) |\n", + "| Server size | 4 GB RAM, ~0.5 CPUs (default server) |" + ] + }, + { + "cell_type": "markdown", + "id": "811fc68e-1e81-482c-82dd-5efa2a89d715", + "metadata": {}, + "source": [ + "## Learning Objectives\n", + "\n", + "By the end of this notebook, you will be able to:\n", + "\n", + "1. Inspect a retrieved observation dataframe to establish its time coverage, sampling interval, gaps, and per-satellite availability.\n", + "2. Prepare observations by accounting for missing values and aggregating them into time bins.\n", + "3. Flag implausible and low-quality SNR values without discarding them.\n", + "4. Visualize SNR over time and data availability, and compare across stations and satellites." + ] + }, + { + "cell_type": "markdown", + "id": "1ec3d79a-ad70-4c58-b8fc-287501adb7e9", + "metadata": {}, + "source": [ + "## Relevant Documentation & Resources\n", + "\n", + "* [EarthScope SDK - GNSS Observations tutorial](https://docs.earthscope.org/sdk/gnss-obs-tutorial)\n", + "* [EarthScope SDK - GNSS Satellite Ephemeris Positions tutorial](https://docs.earthscope.org/sdk/gnss-eph-pos-tutorial)\n", + "* [GeoLab Documentation](https://docs.earthscope.org/geolab)\n", + "* [Polars User Guide](https://docs.pola.rs/)\n", + "* [Altair](https://altair-viz.github.io/)" + ] + }, + { + "cell_type": "markdown", + "id": "076d400d-8c88-4113-a28d-91dac8648679", + "metadata": {}, + "source": [ + "## Contents\n", + "\n", + "1. [Setup & Imports](#id-1-setup-imports)\n", + "2. [Retrieve a Working Dataset](#id-2-retrieve-a-working-dataset)\n", + "3. [Inspecting the Data](#id-3-inspecting-the-data)\n", + "4. [Preparing the Data](#id-4-preparing-the-data)\n", + "5. [Flagging Outliers](#id-5-flagging-outliers)\n", + "6. [Visualizing Observations](#id-6-visualizing-observations)\n", + "7. [Comparing Across Stations and Satellites](#id-7-comparing-across-stations-and-satellites)\n", + "8. [Exploration Exercises](#id-8-exploration-exercises)\n", + "9. [Troubleshooting & Support](#id-9-troubleshooting-support)" + ] + }, + { + "cell_type": "markdown", + "id": "cca79b91-0348-4ee2-842e-740595288820", + "metadata": {}, + "source": [ + "## 1. Setup & Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8cd2903c-1b18-4877-a421-f19374f8fcf3", + "metadata": {}, + "outputs": [], + "source": [ + "# Standard library imports\n", + "import datetime as dt\n", + "\n", + "# Third-party imports\n", + "import altair as alt\n", + "import polars as pl\n", + "from earthscope_sdk import AsyncEarthScopeClient\n", + "\n", + "# Enable the Rust (vegafusion) backend so Altair can handle larger datasets efficiently\n", + "alt.data_transformers.enable(\"vegafusion\")\n", + "\n", + "es = AsyncEarthScopeClient()" + ] + }, + { + "cell_type": "markdown", + "id": "aa3e8c77-ceaf-4f6e-be2e-b99b0b7f0125", + "metadata": {}, + "source": [ + "### Configuration\n", + "\n", + "Set your parameters here before running the rest of the notebook. Every subsequent cell reads from these variables, so this is the only place you need to edit to point the notebook at\n", + "different data.\n", + "\n", + "The three stations below all belong to the permanent Alaska network. Two of them record continuously across the whole window whereas the third does not. Section 3 uses that difference to show what incomplete coverage looks like." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15ecb388-d718-4949-a072-5db6cd398d65", + "metadata": {}, + "outputs": [], + "source": [ + "# Modify these values before running the notebook.\n", + "STATIONS = [\"CLGO\", \"SELD\", \"PCHL\"] # three GNSS stations (4-character IDs)\n", + "SESSION = \"A\" # session name\n", + "SYSTEM = \"G\" # constellation: GPS\n", + "OBS_CODES = [\"1C\", \"2W\"] # L1 C/A and L2 P(Y) via Z-tracking\n", + "FIELDS = [\"snr\", \"phase\"] # measurement columns to return\n", + "\n", + "START = dt.datetime(2025, 7, 20) # query start (UTC)\n", + "END = dt.datetime(2025, 7, 24) # query end (UTC)\n", + "\n", + "NOMINAL_DT = dt.timedelta(seconds=15) # expected spacing between epochs\n", + "TIME_BIN = \"5m\" # aggregation window used in Section 4" + ] + }, + { + "cell_type": "markdown", + "id": "c1c52777-0410-4a6f-b3e6-198cfeebb318", + "metadata": {}, + "source": [ + "## 2. Retrieve a Working Dataset\n", + "\n", + "**What:** Four days of GPS observations from three Alaska stations, restricted to two signals and two measurement fields, returned as an Apache Arrow table and converted to a Polars dataframe. This is the same `gnss_observations()` call you met in Notebook 1, with the filters doing the work of keeping the request small. Each row is one satellite, one signal, one epoch.\n", + "\n", + "**Why:** Every filter we pass is pushed to the server, so only the data we asked for crosses the network. Dropping the constellation and observation-code filters would return several times as many rows for the same window, and requesting every field would widen each one. Filtering keeps the fetch to a minute or two.\n", + "\n", + "**Expected result:** On the order of one million rows across seven columns. The fetch typically takes a minute or two." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "475dd2e2-9db8-453d-b495-d6afa1417d8b", + "metadata": {}, + "outputs": [], + "source": [ + "# Describe the request, then .fetch() to run the query.\n", + "table = await es.data.gnss_observations(\n", + " start_datetime=START,\n", + " end_datetime=END,\n", + " station_name=STATIONS,\n", + " session_name=SESSION,\n", + " system=SYSTEM,\n", + " obs_code=OBS_CODES,\n", + " field=FIELDS,\n", + ").fetch()\n", + "\n", + "df = pl.from_arrow(table).sort(\"timestamp\")\n", + "print(f\"{len(df):,} rows x {df.width} columns\")\n", + "df.head()" + ] + }, + { + "cell_type": "markdown", + "id": "685f59bd-492c-4abe-8068-dfc684502cc1", + "metadata": {}, + "source": [ + "> **Check:** You should see a dataframe with `timestamp`, `satellite`, `obs_code`, `snr`, `phase`, `system`, and `igs` columns." + ] + }, + { + "cell_type": "markdown", + "id": "6069a914-e056-472a-a968-516aff9acf58", + "metadata": {}, + "source": [ + "## 3. Inspecting the data" + ] + }, + { + "cell_type": "markdown", + "id": "a0a9b34a-bb38-4172-92ab-07c3ec744f63", + "metadata": {}, + "source": [ + "These checks compare what the request returned against what it asked for. The five steps below go from the coarsest view of the dataset to the finest." + ] + }, + { + "cell_type": "markdown", + "id": "914d22b0-761e-4547-bfe8-d386c8c1d256", + "metadata": {}, + "source": [ + "### 3.1 Structure and contents" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "161efa31-b733-4d24-a5de-595b688bf9ae", + "metadata": {}, + "outputs": [], + "source": [ + "print(df.schema)\n", + "print()\n", + "print(\"Stations: \", df[\"igs\"].unique().sort().to_list())\n", + "print(\"Systems: \", df[\"system\"].unique().sort().to_list())\n", + "print(\"Obs codes: \", df[\"obs_code\"].unique().sort().to_list())\n", + "print(\"Satellites:\", sorted(df[\"satellite\"].unique().to_list()))" + ] + }, + { + "cell_type": "markdown", + "id": "90cf408c-1423-4c46-9f25-8bbe202a99ad", + "metadata": {}, + "source": [ + "**What to look for:** The stations, systems, and observation codes should match exactly what you requested in the Configuration cell. If a station you asked for is missing from this list, it returned no data at all for your window, which is a different problem from returning incomplete data." + ] + }, + { + "cell_type": "markdown", + "id": "e843c78a-cd11-46e7-a73e-cc1dc2e4cca8", + "metadata": {}, + "source": [ + "## 3.2 Time coverage per station" + ] + }, + { + "cell_type": "markdown", + "id": "1d41790f-5c81-46c2-8aec-f19b462425a5", + "metadata": {}, + "source": [ + "A dataframe can look healthy in aggregate while one station contributes almost nothing. Splitting the coverage per station is the first place that shows up.\n", + "\n", + "We compare the number of distinct epochs each station reported against the number we would expect if it had recorded continuously at the nominal interval." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c717985d-37d6-4103-815d-f6193c31bbe3", + "metadata": {}, + "outputs": [], + "source": [ + "expected_epochs = int((END - START) / NOMINAL_DT)\n", + "print(f\"Expected epochs per station over this window: {expected_epochs:,}\")\n", + "\n", + "coverage = (\n", + " df.group_by(\"igs\")\n", + " .agg(\n", + " pl.col(\"timestamp\").min().alias(\"first_epoch\"),\n", + " pl.col(\"timestamp\").max().alias(\"last_epoch\"),\n", + " pl.col(\"timestamp\").n_unique().alias(\"n_epochs\"),\n", + " pl.len().alias(\"n_rows\"),\n", + " )\n", + " .with_columns(\n", + " (pl.col(\"n_epochs\") / expected_epochs * 100).round(1).alias(\"completeness_pct\")\n", + " )\n", + " .sort(\"igs\")\n", + ")\n", + "coverage" + ] + }, + { + "cell_type": "markdown", + "id": "11d9cecd-adb3-4c5f-b219-249e2e604fa2", + "metadata": {}, + "source": [ + "> **Check:** Two of the three stations should sit close to 100% completeness and start at `2025-07-20 00:00:00`. The third should show a noticeably later `first_epoch` and a completeness figure well under 50%. That station was not reporting for the earlier part of the window. This is real, and it is exactly the kind of thing that would quietly bias a result computed across all three stations without checking." + ] + }, + { + "cell_type": "markdown", + "id": "cd100854-57d3-49f8-9464-b1c2f66e078e", + "metadata": {}, + "source": [ + "### 3.3 Sampling Interval" + ] + }, + { + "cell_type": "markdown", + "id": "499a7c45-2844-40f6-a6db-e34f5a5e82cd", + "metadata": {}, + "source": [ + "Notebook 2 found the sampling interval of the instantaneous position stream by taking the difference between consecutive timestamps. That worked because each position stream has exactly\n", + "one row per epoch.\n", + "\n", + "Observations are different. A single epoch contains one row per satellite, per signal, so a single timestamp is repeated many times over. Differencing the raw `timestamp` column\n", + "therefore returns mostly zeros and tells you nothing about the sampling rate. The cell below shows this." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c2e60dbd-ecdd-4099-8b18-39ce6c450d0f", + "metadata": {}, + "outputs": [], + "source": [ + "# This does NOT give the sampling interval, because many rows share the same epoch.\n", + "(\n", + " df.select(pl.col(\"timestamp\").diff().alias(\"dt\"))[\"dt\"]\n", + " .drop_nulls()\n", + " .value_counts()\n", + " .sort(\"count\", descending=True)\n", + " .head()\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "88f19b54-8b6f-4107-b243-d4a6a0bea140", + "metadata": {}, + "source": [ + "The fix is to reduce to the set of distinct epochs first, for one station at a time. Stations are independent recorders, so mixing them together would interleave their epochs and\n", + "give meaningless differences." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "71483f7f-414c-4ab2-9337-d7c688944a7d", + "metadata": {}, + "outputs": [], + "source": [ + "def epoch_spacing(frame: pl.DataFrame, station_igs: str) -> pl.DataFrame:\n", + " \"\"\"Distribution of time gaps between consecutive distinct epochs at one station.\"\"\"\n", + " return (\n", + " frame.filter(pl.col(\"igs\") == station_igs)\n", + " .select(\"timestamp\")\n", + " .unique()\n", + " .sort(\"timestamp\")\n", + " .select(pl.col(\"timestamp\").diff().alias(\"dt\"))[\"dt\"]\n", + " .drop_nulls()\n", + " .value_counts()\n", + " .sort(\"count\", descending=True)\n", + " )\n", + "\n", + "\n", + "for station_igs in df[\"igs\"].unique().sort():\n", + " print(f\"--- {station_igs} ---\")\n", + " print(epoch_spacing(df, station_igs).head())" + ] + }, + { + "cell_type": "markdown", + "id": "e462d737-d8db-4ef9-90fd-1d81b57a038e", + "metadata": {}, + "source": [ + "**What to look for:** For this window each station returns exactly one spacing, `15s`, which is the nominal sampling interval of these receivers. A single-row result is the healthy case.\n", + "\n", + "On other windows you may see more rows. A few large values scattered among the 15-second majority are ordinary interruptions, and Section 3.4 locates them." + ] + }, + { + "cell_type": "markdown", + "id": "cfb1409a-9df0-40d8-8c9a-32f29de6c303", + "metadata": {}, + "source": [ + "### 3.4 Locating the gaps" + ] + }, + { + "cell_type": "markdown", + "id": "e6f5357d-d589-45e0-8098-733d193baeac", + "metadata": {}, + "source": [ + "Whether gaps exist, and where, is a separate question from how complete a station is overall. The function below returns the start and end of every interruption longer than a threshold you choose.\n", + "\n", + "Expect an empty result on this window: all three stations record continuously once they start. An empty result is still a result. It tells you that the incomplete station's missing data is one contiguous block at the beginning, rather than the receiver dropping in and out repeatedly. Those are different faults and they call for different responses." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "007ed011-3a12-42fb-bf28-5969b047c625", + "metadata": {}, + "outputs": [], + "source": [ + "GAP_THRESHOLD = dt.timedelta(minutes=1)\n", + "\n", + "\n", + "def find_gaps(frame: pl.DataFrame, station_igs: str, threshold: dt.timedelta) -> pl.DataFrame:\n", + " \"\"\"Return every interior gap at one station that exceeds `threshold`.\"\"\"\n", + " epochs = (\n", + " frame.filter(pl.col(\"igs\") == station_igs)\n", + " .select(\"timestamp\")\n", + " .unique()\n", + " .sort(\"timestamp\")\n", + " )\n", + " return (\n", + " epochs.with_columns(pl.col(\"timestamp\").diff().alias(\"gap_length\"))\n", + " .filter(pl.col(\"gap_length\") > threshold)\n", + " .with_columns((pl.col(\"timestamp\") - pl.col(\"gap_length\")).alias(\"gap_start\"))\n", + " .rename({\"timestamp\": \"gap_end\"})\n", + " .select([\"gap_start\", \"gap_end\", \"gap_length\"])\n", + " )\n", + "\n", + "\n", + "for station_igs in df[\"igs\"].unique().sort():\n", + " gaps = find_gaps(df, station_igs, GAP_THRESHOLD)\n", + " print(f\"--- {station_igs}: {len(gaps)} gap(s) longer than {GAP_THRESHOLD} ---\")\n", + " if len(gaps):\n", + " print(gaps.head(10))" + ] + }, + { + "cell_type": "markdown", + "id": "e3e11250-04d1-4f8b-acd1-c24334984fd2", + "metadata": {}, + "source": [ + "> **Note:** This finds interior gaps only, and that limitation is the whole reason to run it alongside Section 3.2. A station that started recording late, or stopped early, produces no\n", + "> difference to detect, because there are no epochs on the missing side to difference against. PCHL is exactly that case: it is missing roughly three quarters of the window and this function reports nothing at all for it. Coverage tells you how much is missing. Gap detection tells you how the data is distributed in time. Both are needed" + ] + }, + { + "cell_type": "markdown", + "id": "998bb9dd-0543-4ca7-908b-6a40ed587af4", + "metadata": {}, + "source": [ + "### 3.5 Per-satellite availability" + ] + }, + { + "cell_type": "markdown", + "id": "f85f3404-29ed-47b8-840f-97f5986ef4d0", + "metadata": {}, + "source": [ + "Coverage can also be uneven within a station. A satellite that is tracked for far fewer epochs than its neighbours may have been obstructed from that site, or may have had its own problems.\n", + "\n", + "We restrict to a single observation code so that each row counts once per satellite per epoch." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cd509598-f2c0-487b-bede-796a65a0b7c4", + "metadata": {}, + "outputs": [], + "source": [ + "# Raw epoch counts are not comparable across stations here, because the three stations do not\n", + "# cover the same amount of wall-clock time. Normalising by each station's own total fixes that.\n", + "station_epochs = df.group_by(\"igs\").agg(pl.col(\"timestamp\").n_unique().alias(\"station_epochs\"))\n", + "\n", + "sat_availability = (\n", + " df.filter(pl.col(\"obs_code\") == OBS_CODES[0])\n", + " .group_by([\"igs\", \"satellite\"])\n", + " .agg(pl.col(\"timestamp\").n_unique().alias(\"n_epochs\"))\n", + " .join(station_epochs, on=\"igs\")\n", + " .with_columns(\n", + " (pl.col(\"n_epochs\") / pl.col(\"station_epochs\") * 100).round(1).alias(\"pct_of_window\")\n", + " )\n", + " .sort([\"igs\", \"satellite\"])\n", + ")\n", + "\n", + "# One row per satellite, one column per station: the share of that station's own recorded\n", + "# epochs during which the satellite was tracked.\n", + "sat_availability.pivot(on=\"igs\", index=\"satellite\", values=\"pct_of_window\").sort(\"satellite\")" + ] + }, + { + "cell_type": "markdown", + "id": "b3f044a4-6ee3-4aba-a4c9-74793839fd72", + "metadata": {}, + "source": [ + "**What to look for:** Each figure is the percentage of that station's own recorded epochs during which the satellite was tracked, so the three columns are directly comparable even though\n", + "one station covers far less wall-clock time than the others.\n", + "\n", + "Read the table two ways: *down a column* to spot a satellite tracked unusually little at one station, and *across a row* to see whether that satellite is under-represented everywhere.\n", + "Section 7 returns to this distinction, because the two cases have very different causes.\n", + "\n", + "> **Check:** Without the normalisation, PCHL's raw counts come out at roughly a third of CLGO's for every satellite, reflecting its shorter recording window rather than anything about the\n", + "> satellites themselves. Raw counts are not comparable across stations with unequal coverage." + ] + }, + { + "cell_type": "markdown", + "id": "c253b8a1-564c-4e66-be51-3bb99ec638ab", + "metadata": {}, + "source": [ + "## 4. Preparing the Data" + ] + }, + { + "cell_type": "markdown", + "id": "efb3f035-2105-4410-abd1-1bfca344667d", + "metadata": {}, + "source": [ + "### 4.1 Accounting for missing values" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ac3a55c5-8908-4258-988b-6c932e3faa67", + "metadata": {}, + "outputs": [], + "source": [ + "df.null_count()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f6c1bb13-556b-45e6-a488-9555f5589b4a", + "metadata": {}, + "outputs": [], + "source": [ + "(\n", + " df.group_by(\"obs_code\")\n", + " .agg(\n", + " pl.len().alias(\"n_rows\"),\n", + " pl.col(\"snr\").null_count().alias(\"snr_nulls\"),\n", + " pl.col(\"phase\").null_count().alias(\"phase_nulls\"),\n", + " )\n", + " .sort(\"obs_code\")\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "42de7c31-9908-4286-92fb-4afa03f5f297", + "metadata": {}, + "source": [ + "> **Note: not every null means the same thing.** Three distinct situations are easy to conflate, and should not be treated alike.\n", + ">\n", + "> * **A null in a measurement column**, such as `phase` on a row where `snr` is present, means the receiver observed the signal but did not produce that particular observable at that epoch. This is a per-signal tracking outcome, not a fault in the data.\n", + "> * **A column that is null throughout** means *nothing to report*.\n", + "> * **A missing row** is different in kind. If a satellite was below the horizon there is no row at all, and therefore no null to count. `null_count()` is blind to this entirely.\n", + "> \n", + "> Counting nulls tells you about the first two. Only comparing against expected epochs, as in Section 3, tells you about the third.\n", + "\n", + "**What to look for:** Compare the two signals rather than the totals. In this window `1C` carries tens of thousands of `phase` nulls while `2W` carries fewer than a hundred, despite the two having similar row counts. If you plan to use carrier phase, this tells you which signal actually delivers it here.\n", + "\n", + "The cause is not established in this notebook. It would require looking at receiver and firmware behavior, which is outside its scope." + ] + }, + { + "cell_type": "markdown", + "id": "68b7f540-9eba-489e-98b7-3ed435e50bf7", + "metadata": {}, + "source": [ + "### 4.2 Aggregating into time bins" + ] + }, + { + "cell_type": "markdown", + "id": "7cd29c4d-c457-40a7-b639-19d3ccb8394b", + "metadata": {}, + "source": [ + "**What it does:** Collapses the raw 15-second observations into summary statistics over fixed time windows, one per station, satellite, and signal.\n", + "\n", + "**Why it matters:** A million-row dataframe is awkward to plot and slow to reason about, and at 15-second resolution most of what you see is noise rather than structure. Binning to a few\n", + "minutes preserves the shape of each satellite pass while cutting the volume by an order of magnitude. Keeping `min`, `max`, and `n_obs` alongside the mean means the aggregation does not hide the variability it is smoothing over.\n", + "\n", + "**Expected output:** The same data reduced from roughly a million rows to a few tens of thousands." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "437f9c32-3ced-4d07-bc3f-fb505de554b4", + "metadata": {}, + "outputs": [], + "source": [ + "binned = (\n", + " df.with_columns(\n", + " # Truncate each timestamp down to the start of its bin\n", + " pl.col(\"timestamp\").dt.truncate(TIME_BIN).alias(\"time_bin\"),\n", + " # Composite label so each satellite/signal combination plots as its own series\n", + " pl.concat_str([\"system\", \"satellite\", \"obs_code\"], separator=\"-\").alias(\"system_sat_obs\"),\n", + " )\n", + " .group_by([\"igs\", \"time_bin\", \"system_sat_obs\", \"satellite\", \"obs_code\"])\n", + " .agg(\n", + " pl.col(\"snr\").mean().alias(\"mean_snr\"),\n", + " pl.col(\"snr\").min().alias(\"min_snr\"),\n", + " pl.col(\"snr\").max().alias(\"max_snr\"),\n", + " pl.len().alias(\"n_obs\"),\n", + " )\n", + " .sort([\"igs\", \"time_bin\"])\n", + ")\n", + "\n", + "print(f\"{len(df):,} raw rows -> {len(binned):,} binned rows\")\n", + "binned.head()" + ] + }, + { + "cell_type": "markdown", + "id": "50e345b1-4f88-43d9-ac32-225676315224", + "metadata": {}, + "source": [ + "> **Check:** `n_obs` should be close to 20 for a 5-minute bin at 15-second sampling. Bins with markedly fewer observations sit at the start or end of a satellite pass, or straddle a gap. Polars skips nulls when computing `mean`, so a bin's `mean_snr` is the average of the values that exist rather than null." + ] + }, + { + "cell_type": "markdown", + "id": "7a66e9a7-c8c7-4b1d-91da-0e4471397815", + "metadata": {}, + "source": [ + "## 5. Flagging Outliers\n", + "\n", + "Two quite different things get called outliers in SNR data, and it is worth separating them.\n", + "\n", + "* **Implausible values.** A signal-to-noise density at or below zero, or implausibly high, is not a weak signal. It is a broken number. These are bad data.\n", + "* **Low but real values.** Weak readings are correct measurements of an unfavourable situation. They are noisy, but they are not wrong. The next section looks at what drives them here.\n", + "\n", + "> **The thresholds below are author-chosen defaults, not standards.** `SNR_LOW = 30` and `SNR_MAX_PLAUSIBLE = 65` are starting points for this notebook, not published limits. Treat them as parameters to set for your own analysis and check them against your project's conventions before relying on them. \n", + "\n", + "### Why we flag rather than delete\n", + "\n", + "Low SNR observations are often filtered out, but they are not useless. GNSS reflectometry uses exactly these low-elevation data, because that is where the signal reflected off the ground\n", + "interferes with the direct signal strongly enough to measure.\n", + "\n", + "A `snr_quality` column marks them without removing them, so a later workflow can decide for itself which rows it needs." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "538f9587-fa9f-4e50-a746-b8d2b9c28f25", + "metadata": {}, + "outputs": [], + "source": [ + "SNR_MIN_PLAUSIBLE = 0.0 # dB-Hz, at or below this is not a real measurement\n", + "SNR_MAX_PLAUSIBLE = 65.0 # dB-Hz, above this is not physically expected\n", + "SNR_LOW = 30.0 # dB-Hz, below this is weak but real\n", + "\n", + "flagged = df.with_columns(\n", + " pl.when(pl.col(\"snr\").is_null())\n", + " .then(pl.lit(\"missing\"))\n", + " .when(\n", + " (pl.col(\"snr\") <= SNR_MIN_PLAUSIBLE) | (pl.col(\"snr\") > SNR_MAX_PLAUSIBLE)\n", + " )\n", + " .then(pl.lit(\"implausible\"))\n", + " .when(pl.col(\"snr\") < SNR_LOW)\n", + " .then(pl.lit(\"low\"))\n", + " .otherwise(pl.lit(\"ok\"))\n", + " .alias(\"snr_quality\")\n", + ")\n", + "\n", + "(\n", + " flagged.group_by(\"snr_quality\")\n", + " .agg(pl.len().alias(\"n\"))\n", + " .with_columns((pl.col(\"n\") / pl.col(\"n\").sum() * 100).round(2).alias(\"pct\"))\n", + " .sort(\"n\", descending=True)\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c8e9dba-2e99-4b91-a2cf-0ac7862778f1", + "metadata": {}, + "outputs": [], + "source": [ + "# The same breakdown per station, since signal environment is a property of the site.\n", + "(\n", + " flagged.group_by([\"igs\", \"snr_quality\"])\n", + " .agg(pl.len().alias(\"n\"))\n", + " .sort([\"igs\", \"snr_quality\"])\n", + " .pivot(on=\"snr_quality\", index=\"igs\", values=\"n\")\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "9eab0ba9-c339-4d03-8445-fd78bea98191", + "metadata": {}, + "source": [ + "> **Check:** `ok` should account for the large majority of rows and `low` for a minority. On this window the implausible and missing categories never occur, because every snr value is present and inside the plausible range. Note what that does to the pivoted table: pivot creates one column per value it finds, so the result has low and ok and nothing else. A window containing implausible readings would produce a third column.\n", + "\n", + "**What to look for:** The per-station shares come out close to one another, near a fifth of observations in each case. Whatever drives the `low` flag is affecting all three sites to much the same degree, which does not fit a site-specific cause such as local obstruction. The next cell splits the flag a different way." + ] + }, + { + "cell_type": "markdown", + "id": "848d030b-0e20-480a-b9bf-127679d00291", + "metadata": {}, + "source": [ + "### 5.1 Interrogating the flag" + ] + }, + { + "cell_type": "markdown", + "id": "85e163f4-e99c-4bdf-8329-a6637236491a", + "metadata": {}, + "source": [ + "Before using `snr_quality`, it is worth checking what it is actually responding to. Splitting it by station showed little variation, so the cell below splits it by signal instead." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9af6608f-2bd3-4ca0-bc5a-e57a8ca8d069", + "metadata": {}, + "outputs": [], + "source": [ + "(\n", + " flagged.group_by([\"obs_code\", \"snr_quality\"])\n", + " .agg(pl.len().alias(\"n\"))\n", + " .sort([\"obs_code\", \"snr_quality\"])\n", + " .pivot(on=\"snr_quality\", index=\"obs_code\", values=\"n\")\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "6dd5a8e8-83e9-4556-9320-014ad9726a6c", + "metadata": {}, + "source": [ + "**What to look for:** The `low` flag is not distributed evenly between the two signals. `1C` and `2W` are recorded from the same satellites, at the same epochs, by the same receivers, so a large imbalance here cannot be explained by satellite elevation or by sky obstruction. Those act on both signals at once.\n", + "\n", + "Notebook 1 describes what separates them: `1C` is the open L1 C/A signal, while `2W` is L2 P(Y) recovered by semi-codeless Z-tracking. They are acquired by different methods, and the table above indicates the two do not deliver comparable signal strength.\n", + "\n", + "A single threshold applied to both signals mostly separates `1C` from `2W`. It says little about sky conditions or station health. Two options are reasonable:\n", + "\n", + "* Set a threshold per observation code, so `low` means weak for that signal rather than weak compared with L1 C/A.\n", + "* Keep one threshold and always report results split by signal, so the difference stays visible." + ] + }, + { + "cell_type": "markdown", + "id": "c5008800-7811-4fbc-9721-9869da927fc2", + "metadata": {}, + "source": [ + "## 6. Visualizing Observations\n", + "\n", + "The plots below are drawn from data already aggregated in Polars, which keeps them responsive even though the underlying dataframe has a million rows." + ] + }, + { + "cell_type": "markdown", + "id": "70dcaf42-636f-4a65-9c4c-61dfaebe7fe5", + "metadata": {}, + "source": [ + "### 6.1 Satellite arcs" + ] + }, + { + "cell_type": "markdown", + "id": "927e80c0-c93d-484b-b620-9abeb2733df5", + "metadata": {}, + "source": [ + "The plot below shows mean SNR over time for a handful of satellites at a single station, on a single day.\n", + "\n", + "**What to look for:** Each satellite traces an arc. SNR climbs as the satellite rises above the horizon, peaks near its closest approach to the receiver, and falls away as it sets. Arcs begin and end abruptly because the satellite crosses the horizon. Departures from that smooth shape are the interesting part: a sudden drop in the middle of an arc suggests a temporary obstruction, while a persistently ragged arc suggests multipath from a reflective surface near the antenna.\n", + "\n", + "> **Note:** `DEMO_STATION` takes the alphabetically first station, which here is one with complete coverage. That is the right choice for seeing clean arcs, but it means the incomplete station never appears in this plot. Set `DEMO_STATION = \"PCHL00USA\"` and re-run the cell to see what a partial record looks like in the same view." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "edf9b94a-b045-4fcf-be34-f2e2f958f124", + "metadata": {}, + "outputs": [], + "source": [ + "DEMO_STATION = df[\"igs\"].unique().sort().to_list()[0] # first station alphabetically\n", + "DEMO_SATS = [5, 7, 11, 13]\n", + "\n", + "DEMO_START = dt.datetime(2025, 7, 23, tzinfo=dt.timezone.utc)\n", + "DEMO_END = dt.datetime(2025, 7, 24, tzinfo=dt.timezone.utc)\n", + "\n", + "arcs = binned.filter(\n", + " (pl.col(\"igs\") == DEMO_STATION)\n", + " & (pl.col(\"obs_code\") == OBS_CODES[0])\n", + " & (pl.col(\"satellite\").is_in(DEMO_SATS))\n", + " & (pl.col(\"time_bin\") >= DEMO_START)\n", + " & (pl.col(\"time_bin\") < DEMO_END)\n", + ")\n", + "\n", + "arcs.plot.point(x=\"time_bin\", y=\"mean_snr\", color=\"system_sat_obs\").properties(\n", + " width=800,\n", + " height=300,\n", + " title=f\"{DEMO_STATION}: mean SNR by satellite, {DEMO_START:%Y-%m-%d}\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "03fefc1c-4991-45d7-aa0f-7011a6a86c55", + "metadata": {}, + "source": [ + "### 6.2 Data availability heatmap" + ] + }, + { + "cell_type": "markdown", + "id": "c7abb785-35f4-433c-b526-de4d17809f5e", + "metadata": {}, + "source": [ + "The heatmap below puts satellites on the vertical axis and time on the horizontal, with colour showing how many observations arrived in each hour. It is the visual counterpart to the coverage table from Section 3.\n", + "\n", + "**What to look for:** Colored stripes slanting across the plot are individual satellite passes, which is the normal pattern. Blank regions show where data is absent. A blank vertical band across all satellites means the receiver was down. A blank horizontal band means one satellite was never tracked. A large blank block at one station and not the others is a station outage." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "954f0ec1-ddee-4651-ad31-0aa2710bf64f", + "metadata": {}, + "outputs": [], + "source": [ + "availability = (\n", + " df.filter(pl.col(\"obs_code\") == OBS_CODES[0])\n", + " .with_columns(pl.col(\"timestamp\").dt.truncate(\"1h\").alias(\"hour\"))\n", + " .group_by([\"igs\", \"hour\", \"satellite\"])\n", + " .agg(pl.len().alias(\"n_obs\"))\n", + ")\n", + "\n", + "\n", + "def availability_chart(frame: pl.DataFrame, station_igs: str) -> alt.Chart:\n", + " \"\"\"One availability heatmap for a single station.\"\"\"\n", + " return (\n", + " alt.Chart(frame.filter(pl.col(\"igs\") == station_igs))\n", + " .mark_rect()\n", + " .encode(\n", + " alt.X(\"hour:T\", title=\"Time (UTC)\"),\n", + " alt.Y(\"satellite:O\", title=\"Satellite\"),\n", + " alt.Color(\"n_obs:Q\", title=\"Obs / hour\", scale=alt.Scale(scheme=\"viridis\")),\n", + " )\n", + " .properties(width=700, height=260, title=f\"{station_igs}: observation availability\")\n", + " )\n", + "\n", + "\n", + "alt.vconcat(\n", + " *[availability_chart(availability, s) for s in df[\"igs\"].unique().sort()]\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "4d5fe3ce-4c6e-4ea9-80cb-90de10976051", + "metadata": {}, + "source": [ + "## 7. Comparing Across Stations and Satellites" + ] + }, + { + "cell_type": "markdown", + "id": "7ce22fd0-d673-4c11-9d5f-7c4eecdaf3b2", + "metadata": {}, + "source": [ + "A single SNR value is hard to judge on its own. Comparing stations against each other, and satellites against each other, gives you a reference for what is normal in this dataset." + ] + }, + { + "cell_type": "markdown", + "id": "78bfe008-e9fb-424c-8137-0d5faaf85e0a", + "metadata": {}, + "source": [ + "### 7.1 Signal strength distribution by station" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1bb7cc46-b4c9-40d7-aba3-68a34a013d73", + "metadata": {}, + "outputs": [], + "source": [ + "snr_stats = (\n", + " flagged.filter(\n", + " (pl.col(\"obs_code\") == OBS_CODES[0]) & (pl.col(\"snr_quality\") != \"implausible\")\n", + " )\n", + " .group_by(\"igs\")\n", + " .agg(\n", + " pl.col(\"snr\").mean().round(2).alias(\"mean\"),\n", + " pl.col(\"snr\").median().alias(\"median\"),\n", + " pl.col(\"snr\").std().round(2).alias(\"std\"),\n", + " pl.col(\"snr\").quantile(0.05).alias(\"p05\"),\n", + " pl.col(\"snr\").quantile(0.95).alias(\"p95\"),\n", + " )\n", + " .sort(\"igs\")\n", + ")\n", + "snr_stats" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15f19f9b-650e-4a45-9b59-e9c868cf9fb3", + "metadata": {}, + "outputs": [], + "source": [ + "# Pre-bin into 1 dB-Hz buckets in Polars so the chart stays light.\n", + "histogram = (\n", + " flagged.filter(\n", + " (pl.col(\"obs_code\") == OBS_CODES[0])\n", + " & (pl.col(\"snr_quality\") != \"implausible\")\n", + " & pl.col(\"snr\").is_not_null()\n", + " )\n", + " .with_columns(pl.col(\"snr\").floor().alias(\"snr_bin\"))\n", + " .group_by([\"igs\", \"snr_bin\"])\n", + " .agg(pl.len().alias(\"n\"))\n", + " .sort([\"igs\", \"snr_bin\"])\n", + ")\n", + "\n", + "histogram.plot.line(x=\"snr_bin\", y=\"n\", color=\"igs\").properties(\n", + " width=800, height=300, title=f\"SNR distribution by station ({OBS_CODES[0]})\"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "44328cad-50f8-43a2-be4f-b3db8dbb65db", + "metadata": {}, + "source": [ + "**What to look for:** On this window the three stations agree closely on center: the means land within a few tenths of a dB-Hz of one another. The differences are in spread instead, so compare `std` and `p05` in the table above rather than `mean`.\n", + "\n", + "Center and spread indicate different things:\n", + "\n", + "* A curve shifted left or right points at the station as a whole: a different antenna, a cable problem, a generally obstructed sky view.\n", + "* A curve with the same peak but a fatter low tail points at part of the sky only, which is what partial obstruction looks like.\n", + "* Similar center and similar spread means the sites are behaving comparably.\n", + "\n", + "Be careful comparing a station with much shorter coverage against the others. A record spanning a quarter of the window samples a different set of satellite geometries, and some of its difference in spread may come from that rather than from the site." + ] + }, + { + "cell_type": "markdown", + "id": "5eb70f80-9c47-4721-9d3c-abf591e09bef", + "metadata": {}, + "source": [ + "### 7.2 Per-satellite mean SNR across stations" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16226d52-e730-460e-991c-e2f4aa2212f7", + "metadata": {}, + "outputs": [], + "source": [ + "sat_comparison = (\n", + " flagged.filter(\n", + " (pl.col(\"obs_code\") == OBS_CODES[0]) & (pl.col(\"snr_quality\") != \"implausible\")\n", + " )\n", + " .group_by([\"igs\", \"satellite\"])\n", + " .agg(\n", + " pl.col(\"snr\").mean().round(2).alias(\"mean_snr\"),\n", + " pl.len().alias(\"n_obs\"),\n", + " )\n", + " .sort([\"satellite\", \"igs\"])\n", + ")\n", + "\n", + "sat_comparison.plot.point(x=\"satellite\", y=\"mean_snr\", color=\"igs\").properties(\n", + " width=800, height=300, title=\"Mean SNR per satellite, compared across stations\"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "999da740-58e5-435c-94dc-263e3831f0ea", + "metadata": {}, + "source": [ + "**What to look for:** This plot separates the two explanations for a weak satellite.\n", + "\n", + "* A satellite that is low at one station but normal at the others points at that station. Something about the site's sky view or hardware disadvantages that part of the sky.\n", + "* A satellite that is low at every station points at the satellite itself. Aging hardware and transmitter problems affect all receivers alike.\n", + "\n", + "Satellites with very few observations will have noisy means, so read `mean_snr` alongside `n_obs` rather than on its own." + ] + }, + { + "cell_type": "markdown", + "id": "1e147f5a-2352-4a16-b28b-3f5c1dbea902", + "metadata": {}, + "source": [ + "## 8. Exploration Exercises" + ] + }, + { + "cell_type": "markdown", + "id": "45f1beb7-097f-40f6-86a3-64b8d3a0d14e", + "metadata": {}, + "source": [ + "Now that you have completed the core workflow, try modifying the parameters below to explore how the results change.\n", + "\n", + "**Try these modifications:**\n", + "\n", + "1. **Change the bin width.** Set `TIME_BIN` in the Configuration section to `\"1m\"`, then `\"30m\"`, and re-run Sections 4 and 6.1. At what point does binning stop reducing noise and start erasing the shape of the satellite arcs?\n", + "\n", + "2. **Set thresholds per signal.** Section 5.1 showed that one global `SNR_LOW` mostly separates the two signals rather than separating good conditions from bad. Rewrite the `snr_quality` expression so the threshold depends on `obs_code`, picking a value for each signal from its own distribution in Section 7.1. Does the per-station picture change once the signal effect is taken out?\n", + "\n", + "3. **Move the quality threshold.** Set `SNR_LOW` to 25, then to 35, and re-run Section 5. What fraction of observations moves between `low` and `ok`? What does the sensitivity of that fraction tell you about publishing a fixed threshold as though it were a physical constant?\n", + "\n", + "4. **Find a worse station.** Replace one entry in `STATIONS` with another station, or swap the list for `network_name=\"PERM:Alaska\"` in the Section 2 request to pull the whole network at once. Which station has the largest share of `low` observations, and does its availability heatmap suggest why?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e24b46e4-799a-4754-81cd-9e719b2a2b9b", + "metadata": {}, + "outputs": [], + "source": [ + "# Exploration cell - use this space to experiment" + ] + }, + { + "cell_type": "markdown", + "id": "a0efc1c3-d60b-4683-9e8f-9772037131ea", + "metadata": {}, + "source": [ + "## 9. Troubleshooting & Support" + ] + }, + { + "cell_type": "markdown", + "id": "df7a004f-68e4-4966-9b14-8782df56518e", + "metadata": {}, + "source": [ + "### Common Issues\n", + "\n", + "| Error | Likely cause | Fix |\n", + "|---|---|---|\n", + "| Sampling interval comes back as `0s` | `timestamp` was differenced across all rows, but many rows share one epoch | Reduce to distinct epochs for a single station first: `.select(\"timestamp\").unique().sort(\"timestamp\")` |\n", + "| Dataframe has 0 rows | No data for that station, session, or window | Widen `START` / `END`, confirm the 4-character station IDs, or request `network_name` instead |\n", + "| `ColumnNotFoundError: phase` | The column was not requested | Add it to `FIELDS` and re-fetch. Server-side field selection means unrequested columns are absent, not null |\n", + "| Kernel dies or restarts during the fetch | Request too large for a 4 GB server | Narrow `OBS_CODES`, `FIELDS`, or the time window, or use a query plan as in Notebook 1 Section 6 |\n", + "| `mean_snr` is null for some bins | Every `snr` value in that bin was null | Expected at pass edges. Filter on `n_obs` if you need well-populated bins only |\n", + "\n", + "### Further Resources\n", + "\n", + "* [EarthScope SDK - GNSS Observations tutorial](https://docs.earthscope.org/sdk/gnss-obs-tutorial)\n", + "* [EarthScope SDK - GNSS Satellite Ephemeris Positions tutorial](https://docs.earthscope.org/sdk/gnss-eph-pos-tutorial)\n", + "* [GeoLab Documentation](https://docs.earthscope.org/geolab)\n", + "* [Polars User Guide](https://docs.pola.rs/)\n", + "* [Altair](https://altair-viz.github.io/)" + ] + } + ], + "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.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}