diff --git a/methods/transient/Transient_Requirement_Validation.ipynb b/methods/transient/Transient_Requirement_Validation.ipynb index 1f94a9b..4aef3ae 100644 --- a/methods/transient/Transient_Requirement_Validation.ipynb +++ b/methods/transient/Transient_Requirement_Validation.ipynb @@ -120,7 +120,16 @@ "from solid_utils.plotting import display_coseismic_validation as display_transient_validation\n", "from solid_utils.configs import update_reference_point\n", "from solid_utils.corrections import run_cmd, pairwise_stack_from_timeseries\n", - "from solid_utils.saving import save_results" + "from solid_utils.saving import save_results\n", + "\n", + "from IPython.display import display, HTML\n", + "\n", + "def show_result(site, method, percent_pass):\n", + " display(HTML(f\"\"\"\n", + "
\n", + " Method {method} final result for {site}: {percent_pass}% of interferograms pass the requirement.\n", + "
\n", + " \"\"\"))" ] }, { @@ -139,12 +148,12 @@ "outputs": [], "source": [ "# === Basic Configuration ===\n", - "site = \"CalVal_NISAR_NewZealandD15\" # Cal/Val location ID\n", + "site = \"CalVal_NISAR_LosAngelesA34\" # Cal/Val location ID\n", "requirement = \"Transient\" # Options: 'Secular', 'Coseismic', 'Transient'\n", "dataset = \"NISAR\" # Dataset type: 'NISAR', 'ARIA_S1', 'ARIA_S1_new'\n", "\n", "\n", - "rundate = \"20260602\" # Date of this Cal/Val run\n", + "rundate = \"20260608\" # Date of this Cal/Val run\n", "version = \"1\" # Version of this Cal/Val run\n", "custom_sites = \"/home/jovyan/my_nisar_sites.txt\" # Path to custom site metadata\n", "\n", @@ -168,7 +177,31 @@ "except KeyError:\n", " raise ValueError(f\"Site ID '{site}' not found in {custom_sites}\")\n", "\n", - "print(f\"Loaded site: {site}\")" + "print(f\"Loaded site: {site}\")\n", + "\n", + "\n", + "def build_report_site_info(method, runtime_params=None):\n", + " \"\"\"Return site info with runtime-resolved values for validation report.\n", + "\n", + " Parameters from my_sites are treated as initial configuration.\n", + " Runtime-resolved values overwrite standard parameter names to reflect\n", + " what was actually used during notebook execution.\n", + " \"\"\"\n", + " runtime_params = runtime_params or {}\n", + " report_site_info = copy.deepcopy(site_info)\n", + "\n", + " # Add top-level run parameters (not in my_sites JSON)\n", + " report_site_info['site'] = site\n", + " report_site_info['requirement'] = requirement\n", + " report_site_info['dataset'] = dataset\n", + " report_site_info['custom_sites_file'] = custom_sites\n", + " report_site_info['validation_method'] = method\n", + "\n", + " # Update with runtime-resolved values under standard parameter names\n", + " for key, value in runtime_params.items():\n", + " report_site_info[key] = value\n", + "\n", + " return report_site_info" ] }, { @@ -292,20 +325,20 @@ "ifgramStack_file = os.path.join(mintpy_dir, 'inputs/ifgramStack.h5')\n", "\n", "# Modify network - base command\n", - "command = f\"modify_network.py {ifgramStack_file} -t {config_file} \"\n", + "command = [\"modify_network.py\", ifgramStack_file, \"-t\", config_file]\n", "\n", "# Check whether exclusions specified in my_sites file\n", "if site_info.get('ifgExcludePair') not in [None, 'auto', 'no']:\n", - " command += f\" --exclude-ifg {site_info.get('ifgExcludePair')}\"\n", + " command.extend([\"--exclude-ifg\", site_info.get('ifgExcludePair')])\n", "\n", "if site_info.get('ifgExcludeDate') not in [None, 'auto', 'no']:\n", - " command += f\" --exclude-date {site_info.get('ifgExcludeDate')}\"\n", + " command.extend([\"--exclude-date\", site_info.get('ifgExcludeDate')])\n", "\n", "if site_info.get('ifgExcludeIndex') not in [None, 'auto', 'no']:\n", - " command += f\" --exclude-ifg-index {site_info.get('ifgExcludeIndex')} \"\n", + " command.extend([\"--exclude-ifg-index\", str(site_info.get('ifgExcludeIndex'))])\n", "\n", "# Run command\n", - "process = subprocess.run(command, shell=True)" + "process = subprocess.run(command, check=True)" ] }, { @@ -315,34 +348,6 @@ "### Retrieve available date pairs" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "# Formulate ifgramStack file\n", - "ifgramStack_file = os.path.join(mintpy_dir, 'inputs/ifgramStack.h5')\n", - "\n", - "# Modify network - base command\n", - "command = f\"modify_network.py {ifgramStack_file} -t {config_file} \"\n", - "\n", - "# Check whether exclusions specified in my_sites file\n", - "if site_info.get('ifgExcludePair') not in [None, 'auto', 'no']:\n", - " command += f\" --exclude-ifg {site_info.get('ifgExcludePair')}\"\n", - "\n", - "if site_info.get('ifgExcludeDate') not in [None, 'auto', 'no']:\n", - " command += f\" --exclude-date {site_info.get('ifgExcludeDate')}\"\n", - "\n", - "if site_info.get('ifgExcludeIndex') not in [None, 'auto', 'no']:\n", - " command += f\" --exclude-ifg-index {site_info.get('ifgExcludeIndex')} \"\n", - "\n", - "# Run command\n", - "process = subprocess.run(command, shell=True)" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -420,12 +425,6 @@ " time_interval = (ifgs_date[i][1]-ifgs_date[i][0]).days\n", " if time_interval != 12:\n", " del_row_index.append(i)\n", - "while i" + "### 4.2. Set Validation parameters " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# === Validation Parameters ===\n", + "pixel_radius = 3\n", + "max_workers = 8\n", + "max_reference_site_retries = 100\n", + "\n", + "coherence_threshold = 0.3\n", + "pass_threshold = 0.683\n", + "\n", + "n_bins = 10\n", + "distance_min_km = 0.1\n", + "distance_max_km_method1 = 50.0\n", + "distance_max_km_method2 = 50.0\n", + "\n", + "num_pixel_pair_samples = 1000000" ] }, { @@ -1121,7 +1157,6 @@ "displacement = {}\n", "gnss_time_series = {}\n", "gnss_time_series_std = {}\n", - "pixel_radius = 3\n", "\n", "# coord once\n", "atr = readfile.read_attribute(ifgs_file)\n", @@ -1139,7 +1174,6 @@ "global_start_str = min(ifg_start).strftime(\"%Y%m%d\")\n", "global_end_str = max(ifg_end).strftime(\"%Y%m%d\")\n", "\n", - "# r = pixel_radius\n", "\n", "def process_station(site_name):\n", " \"\"\"Return per-station results so the main thread can merge dicts safely.\"\"\"\n", @@ -1208,11 +1242,18 @@ " return local_bad, local_disp, local_ts, local_ts_std\n", "\n", "# Run threads\n", - "max_workers = 8 # start here; try 8, 12, 16\n", "with ThreadPoolExecutor(max_workers=max_workers) as ex:\n", - " futures = [ex.submit(process_station, s) for s in site_names]\n", + " futures = {ex.submit(process_station, s): s for s in site_names}\n", " for fut in as_completed(futures):\n", - " local_bad, local_disp, local_ts, local_ts_std = fut.result()\n", + " site_name = futures[fut]\n", + "\n", + " try:\n", + " local_bad, local_disp, local_ts, local_ts_std = fut.result()\n", + " except Exception as e:\n", + " print(f\"Skipping GNSS station {site_name}: {e}\")\n", + " for ifg_ndx in range(n_ifg):\n", + " bad_stn[ifg_ndx].append(site_name)\n", + " continue\n", "\n", " for k, v in local_disp.items():\n", " displacement[k] = v\n", @@ -1243,7 +1284,10 @@ "gnss_time_series = dict(sorted(gnss_time_series.items()))\n", "gnss_time_series_std = dict(sorted(gnss_time_series_std.items()))\n", "displacement = dict(sorted(displacement.items()))\n", - "bad_stn = dict(sorted(bad_stn.items()))" + "bad_stn = dict(sorted(bad_stn.items()))\n", + "\n", + "# Method 1 prerequisites are consolidated after displacement_df is built.\n", + "method1_prereq_failures = []" ] }, { @@ -1266,11 +1310,16 @@ "outputs": [], "source": [ "# Convert displacement dictionaries to pandas dataframes\n", - "displacement = pd.DataFrame.from_dict(displacement, orient='index',\n", - " columns=['lat','lon','gnss_disp','insar_disp'])\n", + "if len(displacement) == 0:\n", + " method1_prereq_failures.append(\"No valid GNSS/InSAR displacements were collected.\")\n", + " displacement_df = pd.DataFrame(columns=['lat', 'lon', 'gnss_disp', 'insar_disp'])\n", + " displacement_df.index = pd.MultiIndex.from_tuples([], names=['ifg index', 'station'])\n", + "else:\n", + " displacement_df = pd.DataFrame.from_dict(displacement, orient='index',\n", + " columns=['lat','lon','gnss_disp','insar_disp'])\n", "\n", - "# Organize by IFG index and site name\n", - "displacement.index = pd.MultiIndex.from_tuples(displacement.index, names=['ifg index','station'])" + " # Organize by IFG index and site name\n", + " displacement_df.index = pd.MultiIndex.from_tuples(displacement_df.index, names=['ifg index','station'])" ] }, { @@ -1290,13 +1339,25 @@ "source": [ "# Drop IFGs with fewer than three stations\n", "drop_index = []\n", - "for i in displacement.index.get_level_values(0).unique():\n", - " if len(displacement.loc[i]) < 3:\n", + "for i in displacement_df.index.get_level_values(0).unique():\n", + " if len(displacement_df.loc[i]) < 3:\n", " drop_index.append(i)\n", - "displacement=displacement.drop(drop_index)\n", + "if drop_index:\n", + " displacement_df = displacement_df.drop(drop_index)\n", "\n", "# ifgs_date after drop for approach 1\n", - "ifgs_date_ap1=np.delete(ifgs_date,drop_index,axis=0)\n" + "method1_ifg_indices = list(displacement_df.index.get_level_values(0).unique())\n", + "ifgs_date_ap1 = [ifgs_date[i] for i in method1_ifg_indices]\n", + "if len(method1_ifg_indices) == 0:\n", + " method1_prereq_failures.append(\"No IFGs remain with at least 3 GNSS stations.\")\n", + "\n", + "method1_prereq_ok = len(method1_prereq_failures) == 0\n", + "method1_summary = \"Method 1 failed to complete.\"\n", + "\n", + "if not method1_prereq_ok:\n", + " print(\"\\033[1mMethod 1 cannot run:\\033[0m\")\n", + " for failure in method1_prereq_failures:\n", + " print(f\" - {failure}\")" ] }, { @@ -1335,41 +1396,87 @@ "metadata": {}, "outputs": [], "source": [ - "# Read reference site\n", - "gnss_ref_site_name = sitedata['sites'][site]['gps_ref_site_name']\n", - "print(f\"Using reference site: {gnss_ref_site_name:s}\")\n", + "if method1_prereq_ok:\n", + " try:\n", + " # Read reference site\n", + " configured_ref_site_name = sitedata['sites'][site]['gps_ref_site_name']\n", + " print(f\"Using reference site: {configured_ref_site_name:s}\")\n", "\n", - "# Loop through interferograms to re-reference\n", - "for ifg_ndx in displacement.index.get_level_values(0).unique():\n", - " # Determine reference site\n", - " if gnss_ref_site_name in ['auto', 'random']:\n", - " # Choose random GNSS site\n", - " gnss_ref_site_name = random.choice(displacement.loc[ifg_ndx].index.unique()) \n", - " print(f\"Ifg {ifg_ndx} reference site: {gnss_ref_site_name}\")\n", + " ifg_indices = list(displacement_df.index.get_level_values(0).unique())\n", "\n", - " \n", - " # Remove reference site values from GNSS and InSAR displacements\n", - " displacement.loc[ifg_ndx, 'gnss_disp'] = displacement.loc[ifg_ndx, 'gnss_disp'].values \\\n", - " - displacement.loc[(ifg_ndx, gnss_ref_site_name), 'gnss_disp']\n", - " displacement.loc[ifg_ndx, 'insar_disp'] = displacement.loc[ifg_ndx, 'insar_disp'].values \\\n", - " - displacement.loc[(ifg_ndx, gnss_ref_site_name), 'insar_disp']\n", + " def reference_site_has_all_data(site_name):\n", + " for ifg_ndx in ifg_indices:\n", + " if site_name not in displacement_df.loc[ifg_ndx].index:\n", + " return False\n", + " values = displacement_df.loc[(ifg_ndx, site_name), ['gnss_disp', 'insar_disp']]\n", + " if not np.isfinite(values.to_numpy(dtype=float)).all():\n", + " return False\n", + " return True\n", + "\n", + " if configured_ref_site_name in ['auto', 'random']:\n", + " candidate_sites = list(displacement_df.index.get_level_values('station').unique())\n", + " tried_sites = set()\n", + " ref_site_name = None\n", + " retry_count = 0\n", + "\n", + " while retry_count < max_reference_site_retries:\n", + " remaining_sites = [site_name for site_name in candidate_sites if site_name not in tried_sites]\n", + " if not remaining_sites:\n", + " break\n", "\n", - " # Reference point pixel coordinates\n", - " ref_y_value = int(atr['REF_Y'])\n", - " ref_x_value = int(atr['REF_X'])\n", - " # ref_x_value = round((displacement.loc[(ifg_ndx, gnss_ref_site_name),'lon'] - W)/lon_step)\n", - " # ref_y_value = round((displacement.loc[(ifg_ndx, gnss_ref_site_name),'lat'] - N)/lat_step)\n", + " candidate_ref_site = random.choice(remaining_sites)\n", + " tried_sites.add(candidate_ref_site)\n", + " retry_count += 1\n", "\n", - " # InSAR displacement values at site location for re-referencing\n", - " ref_disp_insar = insar_displacement[ifg_ndx,\n", - " ref_y_value-pixel_radius:ref_y_value+1+pixel_radius, \n", - " ref_x_value-pixel_radius:ref_x_value+1+pixel_radius]\n", + " if reference_site_has_all_data(candidate_ref_site):\n", + " ref_site_name = candidate_ref_site\n", + " break\n", "\n", - " # Re-referenced mean value at site location\n", - " ref_disp_insar = np.nanmean(ref_disp_insar)\n", + " if ref_site_name is None:\n", + " raise ValueError(\n", + " f\"No GNSS reference site found with finite GNSS/InSAR data for all \"\n", + " f\"{len(ifg_indices)} Method 1 interferograms after {retry_count} retries.\"\n", + " )\n", + " else:\n", + " ref_site_name = configured_ref_site_name\n", + " if not reference_site_has_all_data(ref_site_name):\n", + " raise ValueError(\n", + " f\"Reference site '{ref_site_name}' does not have finite GNSS/InSAR data \"\n", + " f\"for all {len(ifg_indices)} Method 1 interferograms.\"\n", + " )\n", "\n", - " # Subtract reference value from InSAR displacement\n", - " insar_displacement[ifg_ndx] -= ref_disp_insar" + " print(f\"Method 1 reference site: {ref_site_name}\")\n", + "\n", + " # Loop through interferograms to re-reference\n", + " for ifg_ndx in ifg_indices:\n", + " # Remove reference site values from GNSS and InSAR displacements\n", + " displacement_df.loc[ifg_ndx, 'gnss_disp'] = displacement_df.loc[ifg_ndx, 'gnss_disp'].values \\\n", + " - displacement_df.loc[(ifg_ndx, ref_site_name), 'gnss_disp']\n", + " displacement_df.loc[ifg_ndx, 'insar_disp'] = displacement_df.loc[ifg_ndx, 'insar_disp'].values \\\n", + " - displacement_df.loc[(ifg_ndx, ref_site_name), 'insar_disp']\n", + "\n", + " # Reference point pixel coordinates\n", + " ref_y_value = int(atr['REF_Y'])\n", + " ref_x_value = int(atr['REF_X'])\n", + " # ref_x_value = round((displacement.loc[(ifg_ndx, ref_site_name),'lon'] - W)/lon_step)\n", + " # ref_y_value = round((displacement.loc[(ifg_ndx, ref_site_name),'lat'] - N)/lat_step)\n", + "\n", + " # InSAR displacement values at site location for re-referencing\n", + " ref_disp_insar = insar_displacement[ifg_ndx,\n", + " ref_y_value-pixel_radius:ref_y_value+1+pixel_radius, \n", + " ref_x_value-pixel_radius:ref_x_value+1+pixel_radius]\n", + "\n", + " # Re-referenced mean value at site location\n", + " ref_disp_insar = np.nanmean(ref_disp_insar)\n", + "\n", + " # Subtract reference value from InSAR displacement\n", + " insar_displacement[ifg_ndx] -= ref_disp_insar\n", + " except Exception as e:\n", + " method1_prereq_ok = False\n", + " method1_prereq_failures.append(f\"Failed to re-reference GNSS/InSAR: {e}\")\n", + " print(f\"\\033[1mMethod 1 cannot run: {method1_prereq_failures[-1]}\\033[0m\")\n", + "else:\n", + " print(\"\\033[1mMethod 1 prerequisite failed. Skipping re-reference.\\033[0m\")" ] }, { @@ -1387,30 +1494,38 @@ }, "outputs": [], "source": [ - "# Set color values\n", - "cmap_obj = copy.copy(plt.get_cmap(cmap))\n", - "vmin_mm = vmin * 0.24 / 4 / np.pi * 1000\n", - "vmax_mm = vmax * 0.24 / 4 / np.pi * 1000\n", - "\n", - "# Loop through interferograms\n", "gnss_insar_figs = []\n", - "for ifg_ndx in displacement.index.get_level_values(0).unique():\n", - " fig, ax = plt.subplots(figsize = (8,8))\n", - " img1 = ax.imshow(insar_displacement[ifg_ndx],\n", - " cmap=cmap_str, vmin=vmin_mm, vmax=vmax_mm, interpolation='nearest',\n", - " extent=(W, E, S, N))\n", - " ax.set_title(f\"{ifgs_date[ifg_ndx][0].strftime('%Y%m%d')}-{ifgs_date[ifg_ndx][1].strftime('%Y%m%d')}\")\n", - " cbar1 = fig.colorbar(img1, ax=ax, orientation='horizontal')\n", - " cbar1.set_label('LOS displacement [mm]')\n", - "\n", - " for site_name in displacement.loc[ifg_ndx].index:\n", - " lon, lat = displacement.loc[(ifg_ndx, site_name), 'lon'], displacement.loc[(ifg_ndx, site_name), 'lat']\n", - " color = cmap_obj((displacement.loc[(ifg_ndx, site_name), 'gnss_disp']-vmin_mm)/(vmax_mm-vmin_mm))\n", - " ax.scatter(lon, lat, s=8**2, color=color, edgecolors='k')\n", - " ax.annotate(site_name, (lon,lat), color='black')\n", "\n", - " # Append figure to list\n", - " gnss_insar_figs.append(fig)" + "if method1_prereq_ok:\n", + " try:\n", + " # Set color values\n", + " cmap_obj = copy.copy(plt.get_cmap(cmap))\n", + " wavelength = float(insar_metadata['WAVELENGTH'])\n", + " vmin_mm = vmin * wavelength / 4 / np.pi * 1000\n", + " vmax_mm = vmax * wavelength / 4 / np.pi * 1000\n", + "\n", + " # Loop through interferograms\n", + " for ifg_ndx in displacement_df.index.get_level_values(0).unique():\n", + " fig, ax = plt.subplots(figsize = (8,8))\n", + " img1 = ax.imshow(insar_displacement[ifg_ndx],\n", + " cmap=cmap_str, vmin=vmin_mm, vmax=vmax_mm, interpolation='nearest',\n", + " extent=(W, E, S, N))\n", + " ax.set_title(f\"{ifgs_date[ifg_ndx][0].strftime('%Y%m%d')}-{ifgs_date[ifg_ndx][1].strftime('%Y%m%d')}\")\n", + " cbar1 = fig.colorbar(img1, ax=ax, orientation='horizontal')\n", + " cbar1.set_label('LOS displacement [mm]')\n", + "\n", + " for site_name in displacement_df.loc[ifg_ndx].index:\n", + " lon, lat = displacement_df.loc[(ifg_ndx, site_name), 'lon'], displacement_df.loc[(ifg_ndx, site_name), 'lat']\n", + " color = cmap_obj((displacement_df.loc[(ifg_ndx, site_name), 'gnss_disp']-vmin_mm)/(vmax_mm-vmin_mm))\n", + " ax.scatter(lon, lat, s=8**2, color=color, edgecolors='k')\n", + " ax.annotate(site_name, (lon,lat), color='black')\n", + "\n", + " # Append figure to list\n", + " gnss_insar_figs.append(fig)\n", + " except Exception as e:\n", + " print(f\"\\033[1mMethod 1 GNSS/InSAR plotting failed: {e}\\033[0m\")\n", + "else:\n", + " print(\"\\033[1mMethod 1 prerequisite failed. Skipping GNSS/InSAR station plots.\\033[0m\")" ] }, { @@ -1438,7 +1553,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We first pair up all GNSS stations and compare the relative measurement from both GNSS and InSAR. " + "Method 1 is wrapped in a guarded code cell so notebook execution can continue to Method 2 if it fails." ] }, { @@ -1447,53 +1562,97 @@ "metadata": {}, "outputs": [], "source": [ - "# Empty dictionaries for GNSS and InSAR measurements, etc.\n", - "insar_disp = {}\n", - "gnss_disp = {}\n", - "ddiff_dist = {}\n", - "ddiff_disp = {}\n", - "abs_ddiff_disp = {}\n", - "\n", - "# Define ellipsoid for distance calculation\n", - "geod = pyproj.Geod(ellps=\"WGS84\")\n", - "\n", - "# Loop through interferograms\n", - "for ifg_ndx in displacement.index.get_level_values(0).unique():\n", - " displacement_i = displacement.loc[ifg_ndx]\n", - " insar_disp_i = []\n", - " gnss_disp_i = []\n", - " ddiff_dist_i = []\n", - " ddiff_disp_i = []\n", - "\n", - " # Loop through site pairs\n", - " for sta1 in displacement_i.index:\n", - " for sta2 in displacement_i.index:\n", - " if sta2 == sta1:\n", - " break\n", - "\n", - " # Compute InSAR and GNSS displacement residuals\n", - " insar_disp_i.append(displacement_i.loc[sta1, 'insar_disp'] \\\n", - " - displacement_i.loc[sta2, 'insar_disp'])\n", - " gnss_disp_i.append(displacement_i.loc[sta1, 'gnss_disp'] \\\n", - " - displacement_i.loc[sta2, 'gnss_disp'])\n", - "\n", - " # Compute double-difference residual\n", - " ddiff_disp_i.append(gnss_disp_i[-1] - insar_disp_i[-1])\n", - "\n", - " # Compute distance between sites\n", - " _, _, distance = geod.inv(displacement_i.loc[sta1,'lon'], displacement_i.loc[sta1,'lat'],\n", - " displacement_i.loc[sta2,'lon'], displacement_i.loc[sta2,'lat'])\n", - " distance = distance / 1000 # convert unit from m to km\n", + "# Shared label for Methods 1 and 2\n", + "site_loc = sitedata['sites'][site]['calval_location']\n", "\n", - " # Record double difference\n", - " ddiff_dist_i.append(distance)\n", + "# Initialize Method 1 success flag\n", + "method1_ok = False\n", + " \n", "\n", - " # Record all site-to-site values within an IFG\n", - " insar_disp[ifg_ndx] = np.array(insar_disp_i)\n", - " gnss_disp[ifg_ndx] = np.array(gnss_disp_i)\n", - " ddiff_dist[ifg_ndx] = np.array(ddiff_dist_i)\n", - " ddiff_disp[ifg_ndx] = np.array(ddiff_disp_i)\n", - " abs_ddiff_disp[ifg_ndx] = abs(np.array(ddiff_disp_i))" + "if not method1_prereq_ok:\n", + " print(\"\\033[1mMethod 1 prerequisite failed. Continuing to Method 2...\\033[0m\")\n", + " for failure in method1_prereq_failures:\n", + " print(f\" - {failure}\")\n", + "else:\n", + " try:\n", + " # Empty dictionaries for GNSS and InSAR measurements, etc.\n", + " insar_disp = {}\n", + " gnss_disp = {}\n", + " ddiff_dist = {}\n", + " ddiff_disp = {}\n", + " abs_ddiff_disp = {}\n", + "\n", + " # Define ellipsoid for distance calculation\n", + " geod = pyproj.Geod(ellps=\"WGS84\")\n", + "\n", + " # Loop through interferograms\n", + " for ifg_ndx in displacement_df.index.get_level_values(0).unique():\n", + " displacement_i = displacement_df.loc[ifg_ndx]\n", + " insar_disp_i = []\n", + " gnss_disp_i = []\n", + " ddiff_dist_i = []\n", + " ddiff_disp_i = []\n", + "\n", + " # Loop through site pairs\n", + " for sta1 in displacement_i.index:\n", + " for sta2 in displacement_i.index:\n", + " if sta2 == sta1:\n", + " break\n", + "\n", + " # Compute InSAR and GNSS displacement residuals\n", + " insar_disp_i.append(displacement_i.loc[sta1, 'insar_disp'] \\\n", + " - displacement_i.loc[sta2, 'insar_disp'])\n", + " gnss_disp_i.append(displacement_i.loc[sta1, 'gnss_disp'] \\\n", + " - displacement_i.loc[sta2, 'gnss_disp'])\n", + "\n", + " # Compute double-difference residual\n", + " ddiff_disp_i.append(gnss_disp_i[-1] - insar_disp_i[-1])\n", + "\n", + " # Compute distance between sites\n", + " _, _, distance = geod.inv(displacement_i.loc[sta1,'lon'], displacement_i.loc[sta1,'lat'],\n", + " displacement_i.loc[sta2,'lon'], displacement_i.loc[sta2,'lat'])\n", + " distance = distance / 1000 # convert unit from m to km\n", + "\n", + " # Record double difference\n", + " ddiff_dist_i.append(distance)\n", + "\n", + " # Record all site-to-site values within an IFG\n", + " insar_disp[ifg_ndx] = np.array(insar_disp_i)\n", + " gnss_disp[ifg_ndx] = np.array(gnss_disp_i)\n", + " ddiff_dist[ifg_ndx] = np.array(ddiff_dist_i)\n", + " ddiff_disp[ifg_ndx] = np.array(ddiff_disp_i)\n", + " abs_ddiff_disp[ifg_ndx] = abs(np.array(ddiff_disp_i))\n", + "\n", + " # Set requirement thresholds\n", + " transient_distance_rqmt = (distance_min_km, distance_max_km_method1) # distances for evaluation\n", + " transient_threshold_rqmt = lambda L: 3 * (1 + np.sqrt(L)) # coseismic threshold in mm\n", + "\n", + " threshold = pass_threshold # fraction of Gaussian normal distribution for pass/fail\n", + "\n", + " # Loop through interferograms\n", + " method1_validation_figs = []\n", + " for ifg_ndx in displacement_df.index.get_level_values(0).unique():\n", + " # Start and end dates as strings\n", + " start_date = ifgs_date[ifg_ndx][0].strftime('%Y%m%d')\n", + " end_date = ifgs_date[ifg_ndx][1].strftime('%Y%m%d')\n", + "\n", + " # Validation figure and assessment\n", + " _, validation_fig_method1 = display_transient_validation(ddiff_dist[ifg_ndx], abs_ddiff_disp[ifg_ndx],\n", + " site_loc, start_date, end_date,\n", + " requirement=transient_threshold_rqmt,\n", + " distance_rqmt=transient_distance_rqmt,\n", + " n_bins=n_bins,\n", + " threshold=threshold,\n", + " sensor='NISAR',\n", + " validation_type=requirement.lower(),\n", + " validation_data='GNSS')\n", + " method1_validation_figs.append(validation_fig_method1)\n", + "\n", + " method1_ok = True\n", + "\n", + " except Exception as e:\n", + " print(f\"Method 1 failed: {e}\")\n", + " print(\"Continuing to Method 2...\")\n" ] }, { @@ -1507,7 +1666,6 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "\n", "## 5.3 Validate Requirement Based on Binned Measurement Residuals" ] }, @@ -1517,43 +1675,15 @@ "metadata": {}, "outputs": [], "source": [ - "# Set requirement thresholds\n", - "transient_distance_rqmt = (0.1, 50) # distances for evaluation\n", - "transient_threshold_rqmt = lambda L: 3 * (1 + np.sqrt(L)) # coseismic threshold in mm\n", - "\n", - "n_bins = 10 # number of distance bins for analysis\n", - "threshold = 0.683 # fraction of Gaussian normal distribution for pass/fail\n", - "\n", - "# Loop through interferograms\n", - "method1_validation_figs = []\n", - "site_loc = sitedata['sites'][site]['calval_location']\n", - "for ifg_ndx in displacement.index.get_level_values(0).unique():\n", - " # Start and end dates as strings\n", - " start_date = ifgs_date[ifg_ndx][0].strftime('%Y%m%d')\n", - " end_date = ifgs_date[ifg_ndx][1].strftime('%Y%m%d')\n", - "\n", - " # Validation figure and assessment\n", - " _, validation_fig_method1 = display_transient_validation(ddiff_dist[ifg_ndx], abs_ddiff_disp[ifg_ndx],\n", - " site_loc, start_date, end_date,\n", - " requirement=transient_threshold_rqmt,\n", - " distance_rqmt=transient_distance_rqmt,\n", - " n_bins=n_bins,\n", - " threshold=threshold,\n", - " sensor='NISAR',\n", - " validation_type=requirement.lower(),\n", - " validation_data='GNSS')\n", - " method1_validation_figs.append(validation_fig_method1)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Reformat double differences as list\n", - "ddiff_dist_ap1 = list(ddiff_dist.values())\n", - "abs_ddiff_disp_ap1 = list(abs_ddiff_disp.values())" + "if method1_ok:\n", + " # Reformat Method 1 outputs using the same IFG keys for dates and residuals\n", + " ifg_indices_ap1 = list(ddiff_dist.keys())\n", + " ddiff_dist_ap1 = [ddiff_dist[i] for i in ifg_indices_ap1]\n", + " abs_ddiff_disp_ap1 = [abs_ddiff_disp[i] for i in ifg_indices_ap1]\n", + " ifgs_date_ap1 = [ifgs_date[i] for i in ifg_indices_ap1]\n", + "else:\n", + " message = \"Method 1 did not complete successfully. Skipping Method 1 data reformatting.\"\n", + " print(f\"\\033[1m{message}\\033[0m\")" ] }, { @@ -1572,9 +1702,13 @@ "metadata": {}, "outputs": [], "source": [ - "# Define number of interferograms\n", - "n_ifgs = len(ddiff_dist_ap1)\n", - "print(f\"Analyzing {n_ifgs} interferograms\")" + "if method1_ok:\n", + " # Define number of interferograms\n", + " n_ifgs = len(ddiff_dist_ap1)\n", + " print(f\"Analyzing {n_ifgs} interferograms\")\n", + "else:\n", + " message = \"Method 1 failed. Skipping interferogram count.\"\n", + " print(f\"\\033[1m{message}\\033[0m\")" ] }, { @@ -1590,9 +1724,12 @@ "metadata": {}, "outputs": [], "source": [ - "# Define bins over distance requirement\n", - "n_bins = 10\n", - "bins = np.linspace(0.1, 50.0, num=n_bins+1)" + "if method1_ok:\n", + " # Define bins over distance requirement\n", + " bins = np.linspace(distance_min_km, distance_max_km_method1, num=n_bins+1)\n", + "else:\n", + " message = \"Method 1 failed. Skipping bin definition.\"\n", + " print(f\"\\033[1m{message}\\033[0m\")" ] }, { @@ -1601,33 +1738,37 @@ "metadata": {}, "outputs": [], "source": [ - "# Pre-allocate array for number of points for each IFG and bins\n", - "n_all = np.empty([n_ifgs, n_bins+1], dtype=int)\n", + "if method1_ok:\n", + " # Pre-allocate array for number of points for each IFG and bins\n", + " n_all = np.empty([n_ifgs, n_bins+1], dtype=int)\n", "\n", - "# Pre-allocate array for number of points that pass based on requirement\n", - "n_pass = np.empty([n_ifgs,n_bins+1], dtype=int)\n", + " # Pre-allocate array for number of points that pass based on requirement\n", + " n_pass = np.empty([n_ifgs,n_bins+1], dtype=int)\n", "\n", - "# Loop through interferograms\n", - "for i in range(n_ifgs):\n", - " # Determine bin indices\n", - " inds = np.digitize(ddiff_dist_ap1[i], bins)\n", + " # Loop through interferograms\n", + " for i in range(n_ifgs):\n", + " # Determine bin indices\n", + " inds = np.digitize(ddiff_dist_ap1[i], bins)\n", "\n", - " # Loop through bins\n", - " for j in range(1,n_bins+1):\n", - " # Evaluate requirement for the i-th IFG and j-th distance bin\n", - " rqmt = 3*(1+np.sqrt(ddiff_dist_ap1[i][inds==j]))\n", + " # Loop through bins\n", + " for j in range(1,n_bins+1):\n", + " # Evaluate requirement for the i-th IFG and j-th distance bin\n", + " rqmt = 3*(1+np.sqrt(ddiff_dist_ap1[i][inds==j]))\n", "\n", - " # Relative measurement of i-th IFG and j-th distance bin\n", - " rem = abs_ddiff_disp_ap1[i][inds==j]\n", - " assert len(rqmt) == len(rem)\n", - " n_all[i,j-1] = len(rem)\n", - " n_pass[i,j-1] = np.count_nonzero(rem 0, n_all, 1)\n", "\n", - "# Define threshold of data points in a bin that must pass\n", - "threshold = 0.683\n", + " # Define threshold of data points in a bin that must pass\n", + " threshold = pass_threshold\n", "\n", - "# The assumed nature of Gaussian distribution gives a probability of 0.683 of being within one standard deviation.\n", - "success_or_fail = ratio > threshold" + " # The pass threshold represents the one-standard-deviation probability for a Gaussian distribution.\n", + " success_or_fail = ratio > threshold\n", + "else:\n", + " message = \"Method 1 failed. Skipping ratio calculation.\"\n", + " print(f\"\\033[1m{message}\\033[0m\")" ] }, { @@ -1667,27 +1812,31 @@ "metadata": {}, "outputs": [], "source": [ - "def to_str(x:bool):\n", - " if x==True:\n", - " return 'true '\n", - " elif x==False:\n", - " return 'false '\n", + "if method1_ok:\n", + " def to_str(x:bool):\n", + " if x==True:\n", + " return 'true '\n", + " elif x==False:\n", + " return 'false '\n", "\n", - "success_or_fail_str = [list(map(to_str, x)) for x in success_or_fail]\n", + " success_or_fail_str = [list(map(to_str, x)) for x in success_or_fail]\n", "\n", - "columns = []\n", - "for i in range(n_bins):\n", - " columns.append(f'{bins[i]:.2f}-{bins[i+1]:.2f}')\n", - "columns.append('total')\n", + " columns = []\n", + " for i in range(n_bins):\n", + " columns.append(f'{bins[i]:.2f}-{bins[i+1]:.2f}')\n", + " columns.append('total')\n", "\n", - "index = []\n", - "for i in range(len(ifgs_date_ap1)):\n", - " index.append(ifgs_date_ap1[i][0].strftime('%Y%m%d')+'-'+ifgs_date_ap1[i][1].strftime('%Y%m%d'))\n", + " index = []\n", + " for i in range(len(ifgs_date_ap1)):\n", + " index.append(ifgs_date_ap1[i][0].strftime('%Y%m%d')+'-'+ifgs_date_ap1[i][1].strftime('%Y%m%d'))\n", "\n", - "n_all_pd = pd.DataFrame(n_all,columns=columns,index=index)\n", - "n_pass_pd = pd.DataFrame(n_pass,columns=columns,index=index)\n", - "ratio_pd = pd.DataFrame(ratio,columns=columns,index=index)\n", - "success_or_fail_pd = pd.DataFrame(success_or_fail_str,columns=columns,index=index)" + " n_all_pd = pd.DataFrame(n_all,columns=columns,index=index)\n", + " n_pass_pd = pd.DataFrame(n_pass,columns=columns,index=index)\n", + " ratio_pd = pd.DataFrame(ratio,columns=columns,index=index)\n", + " success_or_fail_pd = pd.DataFrame(success_or_fail_str,columns=columns,index=index)\n", + "else:\n", + " message = \"Method 1 failed. Skipping DataFrame creation.\"\n", + " print(f\"\\033[1m{message}\\033[0m\")" ] }, { @@ -1705,7 +1854,11 @@ }, "outputs": [], "source": [ - "n_all_pd" + "if method1_ok:\n", + " display(n_all_pd)\n", + "else:\n", + " message = \"Method 1 failed. No data to display.\"\n", + " print(f\"\\033[1m{message}\\033[0m\")" ] }, { @@ -1723,7 +1876,11 @@ }, "outputs": [], "source": [ - "n_pass_pd" + "if method1_ok:\n", + " display(n_pass_pd)\n", + "else:\n", + " message = \"Method 1 failed. No data to display.\"\n", + " print(f\"\\033[1m{message}\\033[0m\")" ] }, { @@ -1741,13 +1898,17 @@ }, "outputs": [], "source": [ - "# Stylized pandas table\n", - "validation_table_method1 = ratio_pd.style\n", - "validation_table_method1.set_table_styles([ # create internal CSS classes\n", - " {'selector': '.true', 'props': 'background-color: #e6ffe6;'},\n", - " {'selector': '.false', 'props': 'background-color: #ffe6e6;'},\n", - "], overwrite=False)\n", - "validation_table_method1.set_td_classes(success_or_fail_pd)" + "if method1_ok:\n", + " # Stylized pandas table\n", + " validation_table_method1 = ratio_pd.style\n", + " validation_table_method1.set_table_styles([ # create internal CSS classes\n", + " {'selector': '.true', 'props': 'background-color: #e6ffe6;'},\n", + " {'selector': '.false', 'props': 'background-color: #ffe6e6;'},\n", + " ], overwrite=False)\n", + " validation_table_method1.set_td_classes(success_or_fail_pd)\n", + "else:\n", + " message = \"Method 1 failed. No table to display.\"\n", + " print(f\"\\033[1m{message}\\033[0m\")" ] }, { @@ -1764,16 +1925,20 @@ "metadata": {}, "outputs": [], "source": [ - "percentage = np.count_nonzero(ratio_pd['total'] > threshold) / n_ifgs\n", + "if method1_ok:\n", + " percentage_m1 = (np.count_nonzero(ratio_pd['total'] > threshold) / n_ifgs) * 100\n", "\n", - "method_summary = f\"Percentage of interferograms passes the requirement: {percentage}\"\n", + " method1_summary = f\"Percentage of interferograms that pass the requirement: {percentage_m1}%\"\n", "\n", - "if percentage >= 0.70:\n", - " method_summary += \"\\nThe interferogram stack passes the requirement.\"\n", - "else:\n", - " method_summary += \"\\nThe interferogram stack fails the requirement.\"\n", + " if percentage_m1 >= 70:\n", + " method1_summary += \"\\nThe interferogram stack passes the requirement.\"\n", + " else:\n", + " method1_summary += \"\\nThe interferogram stack fails the requirement.\"\n", "\n", - "print(method_summary)" + " print(method1_summary)\n", + "else:\n", + " method1_summary = \"Method 1 failed to complete.\"\n", + " print(f\"\\033[1m{method1_summary}\\033[0m\")" ] }, { @@ -1782,37 +1947,50 @@ "metadata": {}, "outputs": [], "source": [ - "# Save Method 1 results to file\n", - "run_date = dt.now().strftime('%Y%m%dT%H%M%S')\n", - "save_fldr = f\"{run_date}-Transient-Method1\"\n", - "save_dir = os.path.join(mintpy_dir, save_fldr)\n", + "if method1_ok:\n", + " # Save Method 1 results to file\n", + " run_date = dt.now().strftime('%Y%m%dT%H%M%S')\n", + " save_fldr = f\"{run_date}-Transient-Method1\"\n", + " save_dir = os.path.join(mintpy_dir, save_fldr)\n", + "\n", + " # Build runtime parameters for report\n", + " runtime_params = {\n", + " 'gps_ref_site_name': ref_site_name,\n", + " }\n", "\n", - "save_params = {\n", - " 'save_dir': save_dir,\n", - " 'run_date': run_date,\n", - " 'requirement': requirement,\n", - " 'site': site,\n", - " 'method': '1',\n", - " 'sitedata': sitedata['sites'][site],\n", - " 'gnss_insar_figs': gnss_insar_figs,\n", - " 'validation_figs': method1_validation_figs,\n", - " 'validation_table': validation_table_method1,\n", - " 'summary': method_summary\n", - "}\n", - "save_results(**save_params)\n", + " save_params = {\n", + " 'save_dir': save_dir,\n", + " 'run_date': run_date,\n", + " 'requirement': requirement,\n", + " 'site': site,\n", + " 'method': '1',\n", + " 'sitedata': build_report_site_info('1', runtime_params),\n", + " 'gnss_insar_figs': gnss_insar_figs,\n", + " 'validation_figs': method1_validation_figs,\n", + " 'validation_table': validation_table_method1,\n", + " 'summary': method1_summary\n", + " }\n", + " save_results(**save_params)\n", "\n", - "# Save the report in the home directory as well\n", - "save_params['save_dir'] = home_dir\n", - "save_results(**save_params)" + " # Save the report in the home directory as well\n", + " save_params['save_dir'] = home_dir\n", + " save_results(**save_params)\n", + "else:\n", + " message = \"Method 1 failed. Skipping save.\"\n", + " print(f\"\\033[1m{method1_summary}\\033[0m\")" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "
\n", - "Approach 1 final result for CentralValleyA144: around 79% of interferograms passes the requirement.\n", - "
" + "if method1_ok:\n", + " show_result(site, \"1\", percentage_m1)\n", + "else:\n", + " message = \"Method 1 failed\"\n", + " print(f\"\\033[1m{message}\\033[0m\")" ] }, { @@ -1857,44 +2035,58 @@ "metadata": {}, "outputs": [], "source": [ - "# Construct dataset-layer names as lists\n", - "unwrapPhaseName = [f\"unwrapPhase-{date[0].strftime('%Y%m%d')}_{date[1].strftime('%Y%m%d')}\"\n", - " for date in ifgs_date]\n", - "coherenceName = [f\"coherence-{date[0].strftime('%Y%m%d')}_{date[1].strftime('%Y%m%d')}\"\n", - " for date in ifgs_date]\n", - "# Read unwrapped phase from selected interferograms\n", - "ifgs_unw, insar_metadata = readfile.read(ifgs_file, datasetName=unwrapPhaseName)\n", - "# Coherence is coming from the original ifgramStack.h5\n", - "insar_coherence, _ = readfile.read(os.path.join(mintpy_dir, \"inputs\", \"ifgramStack.h5\"), datasetName=coherenceName)\n", - "\n", - "# Convert phase to displacement in m and switch convention to positive range decrease\n", - "insar_displacement = -ifgs_unw*float(insar_metadata['WAVELENGTH']) / (4*np.pi)\n", - "\n", - "# Convert displacement units from m to mm\n", - "insar_displacement = insar_displacement * 1000.\n", - "if singleStack:\n", - " insar_displacement = insar_displacement[np.newaxis, :, :]\n", - " insar_coherence = insar_coherence[np.newaxis, :, :]\n", + "method2_ok = True\n", + "method2_summary = \"Method 2 failed to complete.\"\n", + "hist_insar_figs = []\n", + "method2_validation_figs = []\n", "\n", - "# Read 2D mask array\n", - "msk, _ = readfile.read(msk_file, datasetName=\"waterMask\")\n", + "def mark_method2_failed(stage, error):\n", + " global method2_ok, method2_summary\n", + " method2_ok = False\n", + " method2_summary = f\"Method 2 failed during {stage}: {error}\"\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")\n", "\n", - "# Repeat mask array for each interferogram\n", - "msk = np.stack([msk] * insar_displacement.shape[0], axis=0)\n", - " \n", - "# Set masked pixels to NaN\n", - "insar_displacement[msk == 0] = np.nan\n", - "insar_displacement[insar_displacement==0.0] = np.nan\n", + "try:\n", + " # Construct dataset-layer names as lists\n", + " unwrapPhaseName = [f\"unwrapPhase-{date[0].strftime('%Y%m%d')}_{date[1].strftime('%Y%m%d')}\"\n", + " for date in ifgs_date]\n", + " coherenceName = [f\"coherence-{date[0].strftime('%Y%m%d')}_{date[1].strftime('%Y%m%d')}\"\n", + " for date in ifgs_date]\n", + " # Read unwrapped phase from selected interferograms\n", + " ifgs_unw, insar_metadata = readfile.read(ifgs_file, datasetName=unwrapPhaseName)\n", + " # Coherence is coming from the original ifgramStack.h5\n", + " insar_coherence, _ = readfile.read(os.path.join(mintpy_dir, \"inputs\", \"ifgramStack.h5\"), datasetName=coherenceName)\n", + "\n", + " # Convert phase to displacement in m and switch convention to positive range decrease\n", + " insar_displacement = -ifgs_unw*float(insar_metadata['WAVELENGTH']) / (4*np.pi)\n", + "\n", + " # Convert displacement units from m to mm\n", + " insar_displacement = insar_displacement * 1000.\n", + " if singleStack:\n", + " insar_displacement = insar_displacement[np.newaxis, :, :]\n", + " insar_coherence = insar_coherence[np.newaxis, :, :]\n", + "\n", + " # Read 2D mask array\n", + " msk, _ = readfile.read(msk_file, datasetName=\"waterMask\")\n", + "\n", + " # Repeat mask array for each interferogram\n", + " msk = np.stack([msk] * insar_displacement.shape[0], axis=0)\n", + " \n", + " # Set masked pixels to NaN\n", + " insar_displacement[msk == 0] = np.nan\n", + " insar_displacement[insar_displacement==0.0] = np.nan\n", "\n", - "# Clean up phase-only IFGs to avoid future confusion\n", - "del ifgs_unw\n", + " # Clean up phase-only IFGs to avoid future confusion\n", + " del ifgs_unw\n", "\n", - "# Define number of interferograms\n", - "n_ifgs = len(ddiff_dist_ap1)\n", - "print(f\"Analyzing {n_ifgs} interferograms\")\n", + " # Define number of interferograms for Method 2\n", + " n_ifgs_method2 = len(ifgs_date)\n", + " print(f\"Analyzing {n_ifgs_method2} interferograms\")\n", "\n", - "# Masking pixels with low coherence\n", - "insar_displacement[insar_coherence <0.3] = np.nan" + " # Masking pixels with low coherence\n", + " insar_displacement[insar_coherence < coherence_threshold] = np.nan\n", + "except Exception as e:\n", + " mark_method2_failed(\"InSAR array loading\", e)" ] }, { @@ -1912,15 +2104,21 @@ }, "outputs": [], "source": [ - "cmap_obj = copy.copy(plt.get_cmap('gray'))\n", + "if method2_ok:\n", + " try:\n", + " cmap_obj = copy.copy(plt.get_cmap('gray'))\n", "\n", - "for ifg_ndx in range(n_ifgs):\n", - " fig, ax = plt.subplots(figsize=[18, 5.5])\n", - " img1 = ax.imshow(insar_coherence[ifg_ndx], cmap=cmap_obj, interpolation='nearest', extent=(W, E, S, N))\n", - " ax.set_title(f\"Coherence \"\n", - " f\"\\n Date range {ifgs_date[ifg_ndx][0].strftime('%Y%m%d')}-{ifgs_date[ifg_ndx][1].strftime('%Y%m%d')}\")\n", - " cbar1 = fig.colorbar(img1, ax=ax)\n", - " cbar1.set_label('coherence')" + " for ifg_ndx in range(n_ifgs_method2):\n", + " fig, ax = plt.subplots(figsize=[18, 5.5])\n", + " img1 = ax.imshow(insar_coherence[ifg_ndx], cmap=cmap_obj, interpolation='nearest', extent=(W, E, S, N))\n", + " ax.set_title(f\"Coherence \"\n", + " f\"\\n Date range {ifgs_date[ifg_ndx][0].strftime('%Y%m%d')}-{ifgs_date[ifg_ndx][1].strftime('%Y%m%d')}\")\n", + " cbar1 = fig.colorbar(img1, ax=ax)\n", + " cbar1.set_label('coherence')\n", + " except Exception as e:\n", + " print(f\"\\033[1mMethod 2 warning: coherence plotting failed: {e}\\033[0m\")\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -1931,15 +2129,21 @@ }, "outputs": [], "source": [ - "cmap_obj = copy.copy(plt.get_cmap(cmap))\n", + "if method2_ok:\n", + " try:\n", + " cmap_obj = copy.copy(plt.get_cmap(cmap))\n", "\n", - "for ifg_ndx in range(n_ifgs):\n", - " fig, ax = plt.subplots(figsize=[18, 5.5])\n", - " img1 = ax.imshow(insar_displacement[ifg_ndx], cmap=cmap_obj, interpolation='nearest', extent=(W, E, S, N))\n", - " ax.set_title(f\"Interferogram \"\n", - " f\"\\n Date range {ifgs_date[i][0].strftime('%Y%m%d')}-{ifgs_date[i][1].strftime('%Y%m%d')}\")\n", - " cbar1 = fig.colorbar(img1, ax=ax)\n", - " cbar1.set_label('LOS displacement [mm]')" + " for ifg_ndx in range(n_ifgs_method2):\n", + " fig, ax = plt.subplots(figsize=[18, 5.5])\n", + " img1 = ax.imshow(insar_displacement[ifg_ndx], cmap=cmap_obj, interpolation='nearest', extent=(W, E, S, N))\n", + " ax.set_title(f\"Interferogram \"\n", + " f\"\\n Date range {ifgs_date[ifg_ndx][0].strftime('%Y%m%d')}-{ifgs_date[ifg_ndx][1].strftime('%Y%m%d')}\")\n", + " cbar1 = fig.colorbar(img1, ax=ax)\n", + " cbar1.set_label('LOS displacement [mm]')\n", + " except Exception as e:\n", + " print(f\"\\033[1mMethod 2 warning: interferogram plotting failed: {e}\\033[0m\")\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -1965,11 +2169,17 @@ }, "outputs": [], "source": [ - "if 'NISAR' in dataset:\n", - " X0,Y0 = load_geo_utm(insar_metadata)\n", + "if method2_ok:\n", + " try:\n", + " if 'NISAR' in dataset:\n", + " X0,Y0 = load_geo_utm(insar_metadata)\n", + " else:\n", + " X0,Y0 = load_geo(insar_metadata)\n", + " X0_2d, Y0_2d = np.meshgrid(X0, Y0)\n", + " except Exception as e:\n", + " mark_method2_failed(\"coordinate grid creation\", e)\n", "else:\n", - " X0,Y0 = load_geo(insar_metadata)\n", - "X0_2d, Y0_2d = np.meshgrid(X0, Y0)" + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -1985,12 +2195,18 @@ "metadata": {}, "outputs": [], "source": [ - "# Determine the distance and phase difference between site pairs\n", - "dist = []; rel_measure = []\n", - "for ifg_ndx in range(n_ifgs):\n", - " dist_i, rel_measure_i = samp_pair(X0_2d, Y0_2d, insar_displacement[ifg_ndx], num_samples=1000000)\n", - " dist.append(dist_i/1000)\n", - " rel_measure.append(rel_measure_i)" + "if method2_ok:\n", + " try:\n", + " # Determine the distance and phase difference between site pairs\n", + " dist = []; rel_measure = []\n", + " for ifg_ndx in range(n_ifgs_method2):\n", + " dist_i, rel_measure_i = samp_pair(X0_2d, Y0_2d, insar_displacement[ifg_ndx], num_samples=num_pixel_pair_samples)\n", + " dist.append(dist_i/1000)\n", + " rel_measure.append(rel_measure_i)\n", + " except Exception as e:\n", + " mark_method2_failed(\"pixel-pair sampling\", e)\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -2008,15 +2224,23 @@ }, "outputs": [], "source": [ - "# Plot histogram of distances between pixel pairs\n", - "for ifg_ndx in range(n_ifgs):\n", - " fig, ax = plt.subplots(figsize=[18, 5.5])\n", - " img1 = ax.hist(dist[ifg_ndx], bins=100)\n", - " ax.set_title(f\"Histogram of distance \"\n", - " f\"\\n Date range {ifgs_date[ifg_ndx][0].strftime('%Y%m%d')}-{ifgs_date[ifg_ndx][1].strftime('%Y%m%d')}\")\n", - " ax.set_xlabel(r'Distance ($km$)')\n", - " ax.set_ylabel('Frequency')\n", - " ax.set_xlim(0, 50)" + "if method2_ok:\n", + " try:\n", + " # Plot histogram of distances between pixel pairs\n", + " hist_insar_figs = []\n", + " for ifg_ndx in range(n_ifgs_method2):\n", + " fig, ax = plt.subplots(figsize=[18, 5.5])\n", + " img1 = ax.hist(dist[ifg_ndx], bins=100)\n", + " ax.set_title(f\"Histogram of distance \"\n", + " f\"\\n Date range {ifgs_date[ifg_ndx][0].strftime('%Y%m%d')}-{ifgs_date[ifg_ndx][1].strftime('%Y%m%d')}\")\n", + " ax.set_xlabel(r'Distance ($km$)')\n", + " ax.set_ylabel('Frequency')\n", + " ax.set_xlim(0, 50)\n", + " hist_insar_figs.append(fig)\n", + " except Exception as e:\n", + " print(f\"\\033[1mMethod 2 warning: distance histogram plotting failed: {e}\\033[0m\")\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -2027,14 +2251,21 @@ }, "outputs": [], "source": [ - "# Plot histogram of relative measurements\n", - "for ifg_ndx in range(n_ifgs):\n", - " fig, ax = plt.subplots(figsize=[18, 5.5])\n", - " img1 = ax.hist(rel_measure[ifg_ndx], bins=100)\n", - " ax.set_title(f\"Histogram of Relative Measurement \"\n", - " f\"\\n Date range {ifgs_date[ifg_ndx][0].strftime('%Y%m%d')}-{ifgs_date[ifg_ndx][1].strftime('%Y%m%d')}\")\n", - " ax.set_xlabel(r'Relative Measurement ($mm$)')\n", - " ax.set_ylabel('Frequency')" + "if method2_ok:\n", + " try:\n", + " # Plot histogram of relative measurements\n", + " for ifg_ndx in range(n_ifgs_method2):\n", + " fig, ax = plt.subplots(figsize=[18, 5.5])\n", + " img1 = ax.hist(rel_measure[ifg_ndx], bins=100)\n", + " ax.set_title(f\"Histogram of Relative Measurement \"\n", + " f\"\\n Date range {ifgs_date[ifg_ndx][0].strftime('%Y%m%d')}-{ifgs_date[ifg_ndx][1].strftime('%Y%m%d')}\")\n", + " ax.set_xlabel(r'Relative Measurement ($mm$)')\n", + " ax.set_ylabel('Frequency')\n", + " hist_insar_figs.append(fig)\n", + " except Exception as e:\n", + " print(f\"\\033[1mMethod 2 warning: relative-measurement histogram plotting failed: {e}\\033[0m\")\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -2061,9 +2292,15 @@ "metadata": {}, "outputs": [], "source": [ - "# Define number of interferograms\n", - "n_ifgs = len(ddiff_dist_ap1)\n", - "print(f\"Analyzing {n_ifgs} interferograms\")" + "if method2_ok:\n", + " try:\n", + " # Define number of interferograms for Method 2\n", + " n_ifgs_method2 = len(ifgs_date)\n", + " print(f\"Analyzing {n_ifgs_method2} interferograms\")\n", + " except Exception as e:\n", + " mark_method2_failed(\"interferogram counting\", e)\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -2079,9 +2316,14 @@ "metadata": {}, "outputs": [], "source": [ - "# Define bins over distance requirement\n", - "n_bins = 10\n", - "bins = np.linspace(0.1, 50.0, num=n_bins+1)" + "if method2_ok:\n", + " try:\n", + " # Define bins over distance requirement\n", + " bins = np.linspace(distance_min_km, distance_max_km_method2, num=n_bins+1)\n", + " except Exception as e:\n", + " mark_method2_failed(\"bin definition\", e)\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -2090,33 +2332,39 @@ "metadata": {}, "outputs": [], "source": [ - "# Number of points for each ifgs and bins\n", - "n_all = np.empty([n_ifgs, n_bins+1], dtype=int)\n", + "if method2_ok:\n", + " try:\n", + " # Number of points for each ifgs and bins\n", + " n_all = np.empty([n_ifgs_method2, n_bins+1], dtype=int)\n", "\n", - "# Number of points pass\n", - "n_pass = np.empty([n_ifgs,n_bins+1], dtype=int)\n", + " # Number of points pass\n", + " n_pass = np.empty([n_ifgs_method2,n_bins+1], dtype=int)\n", "\n", - "# Loop through interferograms\n", - "for i in range(n_ifgs):\n", - " # Determine bin indices\n", - " inds = np.digitize(dist[i], bins)\n", + " # Loop through interferograms\n", + " for i in range(n_ifgs_method2):\n", + " # Determine bin indices\n", + " inds = np.digitize(dist[i], bins)\n", "\n", - " # Loop through bins\n", - " for j in range(1, n_bins+1):\n", - " # Evaluate requirement for the i-th IFG and j-th distance bin\n", - " rqmt = 3*(1+np.sqrt(dist[i][inds==j])) # mission requirement for i-th ifgs and j-th bins\n", + " # Loop through bins\n", + " for j in range(1, n_bins+1):\n", + " # Evaluate requirement for the i-th IFG and j-th distance bin\n", + " rqmt = 3*(1+np.sqrt(dist[i][inds==j])) # mission requirement for i-th ifgs and j-th bins\n", "\n", - " # Relative measurement of i-th IFG and j-th distance bin\n", - " rem = rel_measure[i][inds==j] # relative measurement\n", - " assert len(rqmt) == len(rem)\n", - " n_all[i,j-1] = len(rem)\n", - " n_pass[i,j-1] = np.count_nonzero(rem 0, n_all, 1)\n", + " mean_ratio = np.array([np.mean(ratio[:,:-1],axis=1)])\n", + " ratio = np.hstack((ratio,mean_ratio.T))\n", "\n", - "# Define threshold of data points in a bin that must pass\n", - "thresthod = 0.683\n", + " # Define threshold of data points in a bin that must pass\n", + " threshold = pass_threshold\n", "\n", - "#The assumed nature of Gaussian distribution gives a probability of 0.683 of being within one standard deviation.\n", - "success_or_fail = ratio > thresthod" + " # The pass threshold represents the one-standard-deviation probability for a Gaussian distribution.\n", + " success_or_fail = ratio > threshold\n", + " except Exception as e:\n", + " mark_method2_failed(\"ratio calculation\", e)\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -2143,33 +2397,37 @@ "metadata": {}, "outputs": [], "source": [ - "# Set requirement thresholds\n", - "transient_distance_rqmt = (0.1, 100) # distances for evaluation\n", - "transient_threshold_rqmt = lambda L: 3 * (1 + np.sqrt(L)) # coseismic threshold in mm\n", + "if method2_ok:\n", + " try:\n", + " # Set requirement thresholds\n", + " transient_distance_rqmt = (distance_min_km, distance_max_km_method2) # distances for evaluation\n", + " transient_threshold_rqmt = lambda L: 3 * (1 + np.sqrt(L)) # coseismic threshold in mm\n", "\n", - "n_bins = 10 # number of distance bins for analysis\n", - "threshold = 0.683 # fraction of Gaussian normal distribution for pass/fail\n", + " threshold = pass_threshold # fraction of Gaussian normal distribution for pass/fail\n", "\n", "\n", - "# Loop through interferograms\n", - "method2_validation_figs = []\n", - "for ifg_ndx in displacement.index.get_level_values(0).unique():\n", - " # Start and end dates as strings\n", - " start_date = ifgs_date[ifg_ndx][0].strftime('%Y%m%d')\n", - " end_date = ifgs_date[ifg_ndx][1].strftime('%Y%m%d')\n", + " # Loop through interferograms\n", + " for ifg_ndx in range(n_ifgs_method2):\n", + " # Start and end dates as strings\n", + " start_date = ifgs_date[ifg_ndx][0].strftime('%Y%m%d')\n", + " end_date = ifgs_date[ifg_ndx][1].strftime('%Y%m%d')\n", "\n", - " # Validation figure and assessment\n", - " _, validation_fig_method2 = display_transient_validation(dist[ifg_ndx], rel_measure[ifg_ndx],\n", - " site_loc, start_date, end_date,\n", - " requirement=transient_threshold_rqmt,\n", - " distance_rqmt=transient_distance_rqmt,\n", - " n_bins=n_bins,\n", - " threshold=threshold,\n", - " sensor='NISAR',\n", - " validation_type=requirement.lower(),\n", - " validation_data='INSAR')\n", + " # Validation figure and assessment\n", + " _, validation_fig_method2 = display_transient_validation(dist[ifg_ndx], rel_measure[ifg_ndx],\n", + " site_loc, start_date, end_date,\n", + " requirement=transient_threshold_rqmt,\n", + " distance_rqmt=transient_distance_rqmt,\n", + " n_bins=n_bins,\n", + " threshold=threshold,\n", + " sensor='NISAR',\n", + " validation_type=requirement.lower(),\n", + " validation_data='INSAR')\n", "\n", - " method2_validation_figs.append(validation_fig_method2)" + " method2_validation_figs.append(validation_fig_method2)\n", + " except Exception as e:\n", + " print(f\"\\033[1mMethod 2 warning: validation plotting failed: {e}\\033[0m\")\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -2193,28 +2451,34 @@ "metadata": {}, "outputs": [], "source": [ - "# Format evaluation success/failure as string\n", - "def to_str(x:bool):\n", - " if x==True:\n", - " return 'true '\n", - " elif x==False:\n", - " return 'false '\n", + "if method2_ok:\n", + " try:\n", + " # Format evaluation success/failure as string\n", + " def to_str(x:bool):\n", + " if x==True:\n", + " return 'true '\n", + " elif x==False:\n", + " return 'false '\n", "\n", - "success_or_fail_str = [list(map(to_str, x)) for x in success_or_fail]\n", + " success_or_fail_str = [list(map(to_str, x)) for x in success_or_fail]\n", "\n", - "columns = []\n", - "for i in range(n_bins):\n", - " columns.append(f\"{bins[i]:.2f}-{bins[i+1]:.2f}\")\n", - "columns.append('total')\n", + " columns = []\n", + " for i in range(n_bins):\n", + " columns.append(f\"{bins[i]:.2f}-{bins[i+1]:.2f}\")\n", + " columns.append('total')\n", "\n", - "index = []\n", - "for i in range(len(ifgs_date)):\n", - " index.append(f\"{ifgs_date[i][0].strftime('%Y%m%d')}-{ifgs_date[i][1].strftime('%Y%m%d')}\")\n", + " index = []\n", + " for i in range(len(ifgs_date)):\n", + " index.append(f\"{ifgs_date[i][0].strftime('%Y%m%d')}-{ifgs_date[i][1].strftime('%Y%m%d')}\")\n", "\n", - "n_all_pd = pd.DataFrame(n_all, columns=columns, index=index)\n", - "n_pass_pd = pd.DataFrame(n_pass, columns=columns, index=index)\n", - "ratio_pd = pd.DataFrame(ratio, columns=columns+['mean'], index=index)\n", - "success_or_fail_pd = pd.DataFrame(success_or_fail_str, columns=columns+['mean'], index=index)" + " n_all_pd = pd.DataFrame(n_all, columns=columns, index=index)\n", + " n_pass_pd = pd.DataFrame(n_pass, columns=columns, index=index)\n", + " ratio_pd = pd.DataFrame(ratio, columns=columns+['mean'], index=index)\n", + " success_or_fail_pd = pd.DataFrame(success_or_fail_str, columns=columns+['mean'], index=index)\n", + " except Exception as e:\n", + " mark_method2_failed(\"DataFrame creation\", e)\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -2232,7 +2496,10 @@ }, "outputs": [], "source": [ - "n_all_pd" + "if method2_ok:\n", + " display(n_all_pd)\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -2250,7 +2517,10 @@ }, "outputs": [], "source": [ - "n_pass_pd" + "if method2_ok:\n", + " display(n_pass_pd)\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -2268,13 +2538,20 @@ }, "outputs": [], "source": [ - "# Stylized pandas table\n", - "validation_table_method2 = ratio_pd.style\n", - "validation_table_method2.set_table_styles([ # create internal CSS classes\n", - " {'selector': '.true', 'props': 'background-color: #e6ffe6;'},\n", - " {'selector': '.false', 'props': 'background-color: #ffe6e6;'},\n", - "], overwrite=False)\n", - "validation_table_method2.set_td_classes(success_or_fail_pd)" + "if method2_ok:\n", + " try:\n", + " # Stylized pandas table\n", + " validation_table_method2 = ratio_pd.style\n", + " validation_table_method2.set_table_styles([ # create internal CSS classes\n", + " {'selector': '.true', 'props': 'background-color: #e6ffe6;'},\n", + " {'selector': '.false', 'props': 'background-color: #ffe6e6;'},\n", + " ], overwrite=False)\n", + " validation_table_method2.set_td_classes(success_or_fail_pd)\n", + " display(validation_table_method2)\n", + " except Exception as e:\n", + " mark_method2_failed(\"table styling\", e)\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -2298,7 +2575,13 @@ "metadata": {}, "outputs": [], "source": [ - "percentage = np.count_nonzero(ratio_pd['mean'] > thresthod) / n_ifgs" + "if method2_ok:\n", + " try:\n", + " percentage_m2 = (np.count_nonzero(ratio_pd['mean'] > threshold) / n_ifgs_method2) * 100\n", + " except Exception as e:\n", + " mark_method2_failed(\"percentage calculation\", e)\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -2307,11 +2590,19 @@ "metadata": {}, "outputs": [], "source": [ - "print(f\"Percentage of interferograms passes the requirement (70%): {percentage}.\")\n", - "if percentage >= 0.70:\n", - " print('The interferogram stack passes the requirement.')\n", + "if method2_ok:\n", + " try:\n", + " method2_summary = f\"Percentage of interferograms that pass the requirement: {percentage_m2}%.\"\n", + " if percentage_m2 >= 70:\n", + " method2_summary += \"\\nThe interferogram stack passes the requirement.\"\n", + " else:\n", + " method2_summary += \"\\nThe interferogram stack fails the requirement.\"\n", + "\n", + " print(method2_summary)\n", + " except Exception as e:\n", + " mark_method2_failed(\"summary creation\", e)\n", "else:\n", - " print('The interferogram stack fails the requirement.')" + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -2320,36 +2611,52 @@ "metadata": {}, "outputs": [], "source": [ - "# Save Method 2 results to file\n", - "save_fldr = f\"{dt.now().strftime('%Y%m%dT%H%M%S')}-Transient-Method2\"\n", - "save_dir = os.path.join(mintpy_dir, save_fldr)\n", + "if method2_ok:\n", + " try:\n", + " # Save Method 2 results to file\n", + " run_date_method2 = dt.now().strftime('%Y%m%dT%H%M%S')\n", + " save_fldr = f\"{run_date_method2}-Transient-Method2\"\n", + " save_dir = os.path.join(mintpy_dir, save_fldr)\n", "\n", - "save_params = {\n", - " 'save_dir': save_dir,\n", - " 'run_date': run_date,\n", - " 'requirement': requirement,\n", - " 'site': site,\n", - " 'method': '2',\n", - " 'sitedata': sitedata['sites'][site],\n", - " 'gnss_insar_figs': gnss_insar_figs,\n", - " 'validation_figs': method2_validation_figs,\n", - " 'validation_table': validation_table_method2,\n", - " 'summary': method_summary\n", - "}\n", - "save_results(**save_params)\n", + " # Build runtime parameters for report\n", + " # Method 2 does not use GPS reference site\n", + " runtime_params = {\n", + " 'gps_ref_site_name': 'not used by Method 2',\n", + " }\n", "\n", - "# Save the report in the home directory as well\n", - "save_params['save_dir'] = home_dir\n", - "save_results(**save_params)" + " save_params = {\n", + " 'save_dir': save_dir,\n", + " 'run_date': run_date_method2,\n", + " 'requirement': requirement,\n", + " 'site': site,\n", + " 'method': '2',\n", + " 'sitedata': build_report_site_info('2', runtime_params),\n", + " 'gnss_insar_figs': hist_insar_figs,\n", + " 'validation_figs': method2_validation_figs,\n", + " 'validation_table': validation_table_method2,\n", + " 'summary': method2_summary\n", + " }\n", + " save_results(**save_params)\n", + "\n", + " # Save the report in the home directory as well\n", + " save_params['save_dir'] = home_dir\n", + " save_results(**save_params)\n", + " except Exception as e:\n", + " mark_method2_failed(\"result saving\", e)\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "
\n", - "Approach 2 final result for CentralValleyA144: 100% of interferograms passes the requirement.\n", - "
" + "if method2_ok:\n", + " show_result(site, \"2\", percentage_m2)\n", + "else:\n", + " print(f\"\\033[1m{method2_summary}\\033[0m\")" ] }, { @@ -2373,26 +2680,30 @@ "metadata": {}, "outputs": [], "source": [ - "# Loop through interferograms\n", - "for ifg_ndx in displacement.index.get_level_values(0).unique():\n", - " # Initialize figure\n", - " plt.figure(figsize=(11,7))\n", - "\n", - " # Define range of displacement values\n", - " disp_range = (min([*insar_disp[ifg_ndx],*gnss_disp[ifg_ndx]]), max([*insar_disp[ifg_ndx],*gnss_disp[ifg_ndx]]))\n", - "\n", - " # Plot histograms of InSAR and GNSS displacement values\n", - " plt.hist(insar_disp[ifg_ndx], bins=100, range=disp_range, color = \"green\", label='D_InSAR')\n", - " plt.hist(gnss_disp[ifg_ndx], bins=100, range=disp_range, color=\"orange\", label='D_GNSS', alpha=0.5)\n", - "\n", - " # Format figure\n", - " plt.legend(loc='upper right')\n", - " plt.title(f\"Displacements \"\n", - " f\"\\n Date range {ifgs_date[ifg_ndx][0].strftime('%Y%m%d')}-{ifgs_date[ifg_ndx][1].strftime('%Y%m%d')} \"\n", - " f\"\\n Number of station pairs used: {len(insar_disp[ifg_ndx])}\")\n", - " plt.xlabel('LOS Displacement (mm)')\n", - " plt.ylabel('Number of Station Pairs')\n", - " plt.show()" + "if method1_ok:\n", + " # Loop through interferograms\n", + " for ifg_ndx in displacement_df.index.get_level_values(0).unique():\n", + " # Initialize figure\n", + " plt.figure(figsize=(11,7))\n", + "\n", + " # Define range of displacement values\n", + " disp_range = (min([*insar_disp[ifg_ndx],*gnss_disp[ifg_ndx]]), max([*insar_disp[ifg_ndx],*gnss_disp[ifg_ndx]]))\n", + "\n", + " # Plot histograms of InSAR and GNSS displacement values\n", + " plt.hist(insar_disp[ifg_ndx], bins=100, range=disp_range, color = \"green\", label='D_InSAR')\n", + " plt.hist(gnss_disp[ifg_ndx], bins=100, range=disp_range, color=\"orange\", label='D_GNSS', alpha=0.5)\n", + "\n", + " # Format figure\n", + " plt.legend(loc='upper right')\n", + " plt.title(f\"Displacements \"\n", + " f\"\\n Date range {ifgs_date[ifg_ndx][0].strftime('%Y%m%d')}-{ifgs_date[ifg_ndx][1].strftime('%Y%m%d')} \"\n", + " f\"\\n Number of station pairs used: {len(insar_disp[ifg_ndx])}\")\n", + " plt.xlabel('LOS Displacement (mm)')\n", + " plt.ylabel('Number of Station Pairs')\n", + " plt.show()\n", + "else:\n", + " message = \"Method 1 failed. Skipping displacement histogram plots.\"\n", + " print(f\"\\033[1m{message}\\033[0m\")" ] }, { @@ -2401,22 +2712,26 @@ "metadata": {}, "outputs": [], "source": [ - "# Loop through interferograms\n", - "for ifg_ndx in displacement.index.get_level_values(0).unique():\n", - " # Initialize figure\n", - " plt.figure(figsize=(11,7))\n", - "\n", - " # Plot histogram\n", - " plt.hist(ddiff_disp[ifg_ndx], bins=100, color='darkblue', linewidth=1, label='D_gnss - D_InSAR')\n", - "\n", - " # Format figure\n", - " plt.legend(loc='upper right')\n", - " plt.title(f\"Residuals\"\n", - " f\"\\n Date range {ifgs_date[ifg_ndx][0].strftime('%Y%m%d')}-{ifgs_date[ifg_ndx][1].strftime('%Y%m%d')}\"\n", - " f\"\\n Number of stations pairs used: {len(ddiff_disp[i])}\")\n", - " plt.xlabel('Displacement Residual (mm)')\n", - " plt.ylabel('N Stations')\n", - " plt.show()" + "if method1_ok:\n", + " # Loop through interferograms\n", + " for ifg_ndx in displacement_df.index.get_level_values(0).unique():\n", + " # Initialize figure\n", + " plt.figure(figsize=(11,7))\n", + "\n", + " # Plot histogram\n", + " plt.hist(ddiff_disp[ifg_ndx], bins=100, color='darkblue', linewidth=1, label='D_gnss - D_InSAR')\n", + "\n", + " # Format figure\n", + " plt.legend(loc='upper right')\n", + " plt.title(f\"Residuals\"\n", + " f\"\\n Date range {ifgs_date[ifg_ndx][0].strftime('%Y%m%d')}-{ifgs_date[ifg_ndx][1].strftime('%Y%m%d')}\"\n", + " f\"\\n Number of stations pairs used: {len(ddiff_disp[ifg_ndx])}\")\n", + " plt.xlabel('Displacement Residual (mm)')\n", + " plt.ylabel('N Stations')\n", + " plt.show()\n", + "else:\n", + " message = \"Method 1 failed. Skipping residual histogram plots.\"\n", + " print(f\"\\033[1m{message}\\033[0m\")" ] }, { @@ -2425,64 +2740,73 @@ "metadata": {}, "outputs": [], "source": [ - "# Loop through interferograms\n", - "for ifg_ndx in displacement.index.get_level_values(0).unique():\n", - " # Initialize figure\n", - " plt.figure(figsize=(11,7))\n", - "\n", - " # Draw distance threshold\n", - " dist_th = np.linspace(min(ddiff_dist[ifg_ndx]), max(ddiff_dist[ifg_ndx]),100)\n", - " acpt_error = 3*(1+np.sqrt(dist_th))\n", - "\n", - " # Plot residuals\n", - " plt.scatter(ddiff_dist[ifg_ndx], abs_ddiff_disp[ifg_ndx], s=1)\n", - " plt.plot(dist_th, acpt_error, 'r')\n", - "\n", - " # Format plot\n", - " plt.xlabel(\"Distance (km)\")\n", - " plt.ylabel(\"Amplitude of Displacement Residuals (mm)\")\n", - " plt.title(f\"Residuals \"\n", - " f\"\\n Date range {ifgs_date[ifg_ndx][0].strftime('%Y%m%d')}-{ifgs_date[ifg_ndx][1].strftime('%Y%m%d')} \"\n", - " f\"\\n Number of stations pairs used: {len(ddiff_dist[ifg_ndx])}\")\n", - " plt.legend([\"Measurement\", \"Mission Reqiurement\"])\n", - " plt.show()" + "if method1_ok:\n", + " # Loop through interferograms\n", + " for ifg_ndx in displacement_df.index.get_level_values(0).unique():\n", + " # Initialize figure\n", + " plt.figure(figsize=(11,7))\n", + "\n", + " # Draw distance threshold\n", + " dist_th = np.linspace(min(ddiff_dist[ifg_ndx]), max(ddiff_dist[ifg_ndx]),100)\n", + " acpt_error = 3*(1+np.sqrt(dist_th))\n", + "\n", + " # Plot residuals\n", + " plt.scatter(ddiff_dist[ifg_ndx], abs_ddiff_disp[ifg_ndx], s=1)\n", + " plt.plot(dist_th, acpt_error, 'r')\n", + "\n", + " # Format plot\n", + " plt.xlabel(\"Distance (km)\")\n", + " plt.ylabel(\"Amplitude of Displacement Residuals (mm)\")\n", + " plt.title(f\"Residuals \"\n", + " f\"\\n Date range {ifgs_date[ifg_ndx][0].strftime('%Y%m%d')}-{ifgs_date[ifg_ndx][1].strftime('%Y%m%d')} \"\n", + " f\"\\n Number of stations pairs used: {len(ddiff_dist[ifg_ndx])}\")\n", + " plt.legend([\"Measurement\", \"Mission Reqiurement\"])\n", + " plt.show()\n", + "else:\n", + " message = \"Method 1 failed. Skipping residual scatter plots.\"\n", + " print(f\"\\033[1m{message}\\033[0m\") " ] }, { "cell_type": "code", "execution_count": null, "metadata": { + "scrolled": true, "tags": [] }, "outputs": [], "source": [ - "# Loop through interferograms\n", - "for ifg_ndx in range(insar_displacement.shape[0]):\n", - " # Define start and end dates\n", - " start_time_str = ifgs_date[ifg_ndx][0].strftime('%Y%m%d')\n", - " end_time_str = ifgs_date[ifg_ndx][1].strftime('%Y%m%d')\n", + "if method1_ok:\n", + " # Loop through Method 1 interferograms with GNSS time-series data\n", + " for ifg_ndx in gnss_time_series.columns.get_level_values(0).unique():\n", + " # Define start and end dates\n", + " start_time_str = ifgs_date[ifg_ndx][0].strftime('%Y%m%d')\n", + " end_time_str = ifgs_date[ifg_ndx][1].strftime('%Y%m%d')\n", "\n", - " # Loop through stations in IFG\n", - " for site_name in gnss_time_series[ifg_ndx].columns:\n", - " print(f\"Plotting GPS postion from {start_time_str} to {end_time_str} at station: {site_name}\")\n", + " # Loop through stations in IFG\n", + " for site_name in gnss_time_series[ifg_ndx].columns:\n", + " print(f\"Plotting GPS postion from {start_time_str} to {end_time_str} at station: {site_name}\")\n", "\n", - " # Retrieve GNSS time-series\n", - " series = gnss_time_series[ifg_ndx, site_name]\n", + " # Retrieve GNSS time-series\n", + " series = gnss_time_series[(ifg_ndx, site_name)]\n", "\n", - " # Initialize figure\n", - " plt.figure(figsize=(15,5))\n", + " # Initialize figure\n", + " plt.figure(figsize=(15,5))\n", "\n", - " # Plot time-series data\n", - " plt.scatter(pd.date_range(start=ifgs_date[ifg_ndx][0], end=ifgs_date[ifg_ndx][1]),series)\n", + " # Plot time-series data\n", + " plt.scatter(pd.date_range(start=ifgs_date[ifg_ndx][0], end=ifgs_date[ifg_ndx][1]),series)\n", "\n", - " # Foramt figure\n", - " plt.title(f\"station name: {site_name}\")\n", - " plt.xlabel('Time')\n", - " plt.ylabel('Relative position in LOS direction (mm)')\n", + " # Foramt figure\n", + " plt.title(f\"station name: {site_name}\")\n", + " plt.xlabel('Time')\n", + " plt.ylabel('Relative position in LOS direction (mm)')\n", "\n", - " # Save figure\n", - " plt.savefig(os.path.join(work_dir, f\"{start_time_str}_{end_time_str}_{site_name}.jpg\"))\n", - " plt.close()" + " # Save figure\n", + " plt.savefig(os.path.join(work_dir, f\"{start_time_str}_{end_time_str}_{site_name}.jpg\"))\n", + " plt.close()\n", + "else:\n", + " message = \"Method 1 failed. Skipping GNSS time series plots.\"\n", + " print(f\"\\033[1m{message}\\033[0m\") " ] }, {