diff --git a/docs/On-demand_SE_setup_instructions.pptx b/docs/On-demand_SE_setup_instructions_082025.pptx similarity index 72% rename from docs/On-demand_SE_setup_instructions.pptx rename to docs/On-demand_SE_setup_instructions_082025.pptx index 577c58e..f60faf5 100644 Binary files a/docs/On-demand_SE_setup_instructions.pptx and b/docs/On-demand_SE_setup_instructions_082025.pptx differ diff --git a/methods/coseismic/Coseismic_Requirement_Validation.ipynb b/methods/coseismic/Coseismic_Requirement_Validation.ipynb index d990640..cf91f1f 100644 --- a/methods/coseismic/Coseismic_Requirement_Validation.ipynb +++ b/methods/coseismic/Coseismic_Requirement_Validation.ipynb @@ -146,7 +146,7 @@ "dataset = \"ARIA_S1_new\" # Dataset type: 'ARIA_S1', 'ARIA_S1_new'\n", "aria_gunw_version = \"3_0_1\"\n", "\n", - "rundate = \"20250820\" # Date of this Cal/Val run\n", + "rundate = \"20250826\" # Date of this Cal/Val run\n", "version = \"1\" # Version of this Cal/Val run\n", "custom_sites = \"/home/jovyan/my_sites.txt\" # Path to custom site metadata\n", "\n", @@ -173,8 +173,9 @@ "print(f\"Loaded site: {site}\")\n", "\n", "# === Plot Parameters ===\n", - "vmin, vmax = -100, 100 # mm/yr\n", - "cmap = plt.get_cmap('RdBu_r')" + "vmin, vmax = -25, 25 # mm/yr\n", + "cmap_str = 'RdBu_r'\n", + "cmap = plt.get_cmap(cmap_str)" ] }, { @@ -416,8 +417,21 @@ "metadata": {}, "outputs": [], "source": [ - "scp_args = f\"velocity.h5 velocity -v {vmin} {vmax} --colormap RdBu_r --figtitle LOS_Velocity --unit mm/yr -m {msk_file}\"\n", - "view.main(scp_args.split());" + "# Visualize velocity before any corrections\n", + "title = 'LOS velocity, no corrections'\n", + "view_opts = [\n", + " velocity_filename, \"velocity\",\n", + " \"-m\", msk_file,\n", + " \"--dem\", f\"{mintpy_dir}/inputs/geometryGeo.h5\", \"--dem-blend\", \"--dem-nocontour\",\n", + " \"--lalo-label\",\n", + " \"--coastline\", \"10m\",\n", + " \"-c\", cmap_str,\n", + " \"-u\", \"mm\",\n", + " \"-v\", str(vmin), str(vmax),\n", + " \"--title\", title,\n", + " \"--noverbose\"\n", + "]\n", + "view.main(view_opts)" ] }, { @@ -494,10 +508,43 @@ "outputs": [], "source": [ "# Visualize the corrections\n", - "if 'do_SET' in site_info.keys() and site_info.get('do_SET') != \"False\":\n", - " view.main([set_cor_file, '-m', msk_file])\n", - "else:\n", - " site_info['do_SET'] = \"False\"\n", + "if site_info.get('do_SET') != \"False\":\n", + " view_opts = [\n", + " set_cor_file,\n", + " \"-m\", msk_file,\n", + " \"-c\", cmap_str,\n", + " \"--noverbose\"\n", + " ]\n", + " view.main(view_opts)\n", + "else: \n", + " print('#'*10, 'Solid Earth Tide Correction set to False', '#'*10)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Visualize velocity\n", + "print(velocity_filename)\n", + "if site_info.get('do_SET') != \"False\":\n", + " view_opts = [\n", + " velocity_filename, \"velocity\",\n", + " \"-m\", msk_file,\n", + " \"--dem\", f\"{mintpy_dir}/inputs/geometryGeo.h5\",\n", + " \"--dem-blend\", # blend the DEM shade with input image to have a GMT-like impression.\n", + " \"--dem-nocontour\", # do not show DEM contour lines\n", + " \"--coastline\", \"10m\", # Draw coastline with specified resolution (10m, 50m, 110m). This enables '--lalo-label'\n", + " \"--lalo-label\", # Show N, S, E, W tick label for plot in geo-coordinate.\n", + " \"-c\", cmap_str,\n", + " \"-u\", \"mm\",\n", + " \"-v\", str(vmin), str(vmax),\n", + " \"--title\", 'Velocity after correction for SET',\n", + " \"--noverbose\"\n", + " ]\n", + " view.main(view_opts)\n", + "else: \n", " print('#'*10, 'Solid Earth Tide Correction set to False', '#'*10)" ] }, @@ -550,7 +597,24 @@ "source": [ "# Visualize the corrections\n", "if site_info.get('do_iono') != \"False\":\n", - " view.main([iono_cor_file, '-m', msk_file])\n", + " view_opts = [\n", + " iono_stack_file,\n", + " \"--wrap\",\n", + " \"-m\", msk_file,\n", + " \"-c\", cmap_str,\n", + " \"--noverbose\"\n", + " ]\n", + " print('#'*10, 'Wrapped pair-wise ionosphere', '#'*10)\n", + " view.main(view_opts)\n", + "\n", + " view_opts = [\n", + " iono_stack_file,\n", + " \"-m\", msk_file,\n", + " \"-c\", cmap_str,\n", + " \"--noverbose\"\n", + " ]\n", + " print('#'*10, 'Unwrapped pair-wise ionosphere', '#'*10)\n", + " view.main(view_opts)\n", "else: \n", " print('#'*10, 'Ionosphere Correction set to False', '#'*10)" ] @@ -603,6 +667,7 @@ "\n", " tropo_cor_file = os.path.join(mintpy_dir, \"inputs\", \"ERA5.h5\")\n", " timeseries_filename = output_timeseries\n", + " velocity_filename = output_velocity\n", "\n", " else:\n", " # HRRR-based correction\n", @@ -634,7 +699,13 @@ "source": [ "# Visualize the corrections\n", "if site_info.get('do_tropo') != \"False\":\n", - " view.main([tropo_cor_file, '-m', msk_file])\n", + " view_opts = [\n", + " tropo_cor_file,\n", + " \"-m\", msk_file,\n", + " \"-c\", cmap_str,\n", + " \"--noverbose\"\n", + " ]\n", + " view.main(view_opts)\n", "else: \n", " print('#'*10, 'Troposhere Correction set to False', '#'*10)" ] @@ -793,7 +864,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ "# Parse the TS model into argument string for command line execution\n", @@ -823,7 +896,8 @@ " function_str += f' {val:f}'\n", "\n", "# Run command\n", - "command = 'timeseries2velocity.py ' + timeseries_filename + function_str\n", + "vel_file = 'velocity_nostep.h5'\n", + "command = 'timeseries2velocity.py ' + timeseries_filename + ' -o ' + vel_file + function_str\n", "process = subprocess.run(command, shell=True)\n", "\n", "# Load velocity file\n", @@ -851,8 +925,21 @@ }, "outputs": [], "source": [ - "scp_args = f\"velocity.h5 velocity -v {vmin/2} {vmax/2} --colormap RdBu --figtitle LOS_Velocity --unit mm/yr\"\n", - "view.main(scp_args.split())" + "# Visualize velocity\n", + "title = \"LOS velocity, linear fit\"\n", + "view_opts = [\n", + " vel_file, \"velocity\",\n", + " \"-m\", msk_file,\n", + " \"--dem\", f\"{mintpy_dir}/inputs/geometryGeo.h5\", \"--dem-blend\", \"--dem-nocontour\",\n", + " \"--lalo-label\",\n", + " \"--coastline\", \"10m\",\n", + " \"-c\", cmap_str,\n", + " \"-u\", \"mm\",\n", + " \"-v\", str(vmin), str(vmax),\n", + " \"--title\", title,\n", + " \"--noverbose\"\n", + "]\n", + "view.main(view_opts)" ] }, { @@ -888,7 +975,8 @@ }, "outputs": [], "source": [ - "command = 'timeseries2velocity.py ' + timeseries_filename + function_str\n", + "vel_file = \"velocity_coseismic.h5\"\n", + "command = 'timeseries2velocity.py ' + timeseries_filename + ' -o ' + vel_file + function_str\n", "process = subprocess.run(command, shell=True)\n", "\n", "EQdataset = 'step' + sitedata['sites'][site]['earthquakeDate']\n", @@ -918,10 +1006,36 @@ }, "outputs": [], "source": [ - "scp_args = f\"velocity.h5 step{sitedata['sites'][site]['earthquakeDate']} -v {vmin} {vmax} --colormap RdBu --unit mm --figtitle LOS_Coseismic\"\n", - "view.main(scp_args.split())\n", - "scp_args = f\"velocity.h5 velocity -v {vmin/2} {vmax/2} --colormap RdBu --unit mm/yr --figtitle LOS_Velocity\"\n", - "view.main(scp_args.split())" + "# Visualize velocity\n", + "title = \"LOS Coseismic, earthquake date\"\n", + "view_opts = [\n", + " vel_file, \"step\"+sitedata['sites'][site]['earthquakeDate'],\n", + " \"-m\", msk_file,\n", + " \"--dem\", f\"{mintpy_dir}/inputs/geometryGeo.h5\", \"--dem-blend\", \"--dem-nocontour\",\n", + " \"--lalo-label\",\n", + " \"--coastline\", \"10m\",\n", + " \"-c\", cmap_str,\n", + " \"-u\", \"mm\",\n", + " \"-v\", str(vmin), str(vmax),\n", + " \"--title\", title,\n", + " \"--noverbose\"\n", + "]\n", + "view.main(view_opts)\n", + "\n", + "title = \"LOS Linear Velocity Fit\"\n", + "view_opts = [\n", + " vel_file, \"velocity\",\n", + " \"-m\", msk_file,\n", + " \"--dem\", f\"{mintpy_dir}/inputs/geometryGeo.h5\", \"--dem-blend\", \"--dem-nocontour\",\n", + " \"--lalo-label\",\n", + " \"--coastline\", \"10m\",\n", + " \"-c\", cmap_str,\n", + " \"-u\", \"mm/yr\",\n", + " \"-v\", str(vmin), str(vmax),\n", + " \"--title\", title,\n", + " \"--noverbose\"\n", + "]\n", + "view.main(view_opts)" ] }, { @@ -1436,10 +1550,11 @@ "double_diff_rel_measure_err = np.array([diff_res.dis_err for diff_res in double_diffs.values()])\n", "\n", "# Validation figure and assessment\n", + "site_loc = sitedata['sites'][site]['calval_location']\n", "(validation_table_method1,\n", " validation_fig_method1) = display_coseismic_validation(site_dist, # binned distance for point\n", " double_diff_rel_measure, # binned double-difference velocities mm/yr\n", - " site, # cal/val site name\n", + " site_loc, # cal/val site name\n", " start_date, # start date of InSAR dataset\n", " end_date, # end date of InSAR dataset \n", " requirement=coseismic_threshold_rqmt, # measurement requirement to meet, e.g 2 mm/yr for 3 years of data over 0.1-50km\n", @@ -1543,7 +1658,8 @@ "print(f\"Present directory: {mintpy_dir:s}\")\n", "\n", "# Step function fit when there was no earthquake\n", - "command = 'timeseries2velocity.py ' + timeseries_filename + ' --step ' + sitedata['sites'][site]['noEarthquakeDate']\n", + "vel_file = \"velocity_noearthquake.h5\"\n", + "command = 'timeseries2velocity.py ' + timeseries_filename + ' -o ' + vel_file + ' --step ' + sitedata['sites'][site]['noEarthquakeDate']\n", "process = subprocess.run(command, shell=True)" ] }, @@ -1560,10 +1676,36 @@ "metadata": {}, "outputs": [], "source": [ - "scp_args = f\"velocity.h5 step{sitedata['sites'][site]['noEarthquakeDate']} -v {vmin} {vmax} --colormap RdBu --unit mm --figtitle LOS_Coseismic\"\n", - "view.main(scp_args.split())\n", - "scp_args = f\"velocity.h5 velocity -v {vmin/2} {vmax/2} --colormap RdBu --unit mm/yr --figtitle LOS_Velocity\"\n", - "view.main(scp_args.split())" + "# Visualize velocity\n", + "title = \"LOS Coseismic, no earthquake date\"\n", + "view_opts = [\n", + " vel_file, \"step\"+sitedata['sites'][site]['noEarthquakeDate'],\n", + " \"-m\", msk_file,\n", + " \"--dem\", f\"{mintpy_dir}/inputs/geometryGeo.h5\", \"--dem-blend\", \"--dem-nocontour\",\n", + " \"--lalo-label\",\n", + " \"--coastline\", \"10m\",\n", + " \"-c\", cmap_str,\n", + " \"-u\", \"mm\",\n", + " \"-v\", str(vmin), str(vmax),\n", + " \"--title\", title,\n", + " \"--noverbose\"\n", + "]\n", + "view.main(view_opts)\n", + "\n", + "title = \"LOS Linear Velocity Fit\"\n", + "view_opts = [\n", + " vel_file, \"velocity\",\n", + " \"-m\", msk_file,\n", + " \"--dem\", f\"{mintpy_dir}/inputs/geometryGeo.h5\", \"--dem-blend\", \"--dem-nocontour\",\n", + " \"--lalo-label\",\n", + " \"--coastline\", \"10m\",\n", + " \"-c\", cmap_str,\n", + " \"-u\", \"mm/yr\",\n", + " \"-v\", str(vmin), str(vmax),\n", + " \"--title\", title,\n", + " \"--noverbose\"\n", + "]\n", + "view.main(view_opts)" ] }, { @@ -1733,10 +1875,11 @@ "outputs": [], "source": [ "# Validation figure and assessment\n", + "site_loc = sitedata['sites'][site]['calval_location']\n", "(validation_table_method2,\n", " validation_fig_method2) = display_coseismic_validation(insar_sample_dist, # binned distance for point\n", " insar_rel_measure, # binned double-difference velocities mm/yr\n", - " site, # cal/val site name\n", + " site_loc, # cal/val site name\n", " start_date, # start date of InSAR dataset\n", " end_date, # end date of InSAR dataset \n", " requirement=coseismic_threshold_rqmt, # measurement requirement to meet, e.g 2 mm/yr for 3 years of data over 0.1-50km\n", diff --git a/methods/secular/Secular_Requirement_Validation.ipynb b/methods/secular/Secular_Requirement_Validation.ipynb index 6257a6d..fd82084 100644 --- a/methods/secular/Secular_Requirement_Validation.ipynb +++ b/methods/secular/Secular_Requirement_Validation.ipynb @@ -149,7 +149,7 @@ "dataset = \"ARIA_S1_new\" # Dataset type: 'ARIA_S1', 'ARIA_S1_new'\n", "aria_gunw_version = \"3_0_1\"\n", "\n", - "rundate = \"20250820\" # Date of this Cal/Val run\n", + "rundate = \"20250826\" # Date of this Cal/Val run\n", "version = \"1\" # Version of this Cal/Val run\n", "custom_sites = \"/home/jovyan/my_sites.txt\" # Path to custom site metadata\n", "\n", @@ -561,7 +561,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ "# Visualize the corrections\n", @@ -675,7 +677,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ "if 'do_tropo' in site_info.keys() and site_info.get(\"do_tropo\") != \"False\":\n", @@ -970,7 +974,9 @@ " function_str += f' {val:f}'\n", "\n", "# Run command\n", - "command = 'timeseries2velocity.py ' + timeseries_filename + function_str\n", + "vel_file = \"velocity_final_fit.h5\"\n", + "command = 'timeseries2velocity.py ' + timeseries_filename + ' -o ' + vel_file + function_str\n", + "print(command)\n", "process = subprocess.run(command, shell=True)\n", "\n", "# Load velocity file\n", @@ -996,8 +1002,21 @@ "metadata": {}, "outputs": [], "source": [ - "scp_args = f\"velocity.h5 velocity -v {vmin} {vmax} --colormap RdBu_r --figtitle LOS_Velocity --unit mm/yr -m {msk_file}\"\n", - "view.main(scp_args.split());" + "view_opts = [\n", + " vel_file, \"velocity\",\n", + " \"-m\", msk_file,\n", + " \"--dem\", f\"{mintpy_dir}/inputs/geometryGeo.h5\",\n", + " \"--dem-blend\", # blend the DEM shade with input image to have a GMT-like impression.\n", + " \"--dem-nocontour\", # do not show DEM contour lines\n", + " \"--coastline\", \"10m\", # Draw coastline with specified resolution (10m, 50m, 110m). This enables '--lalo-label'\n", + " \"--lalo-label\", # Show N, S, E, W tick label for plot in geo-coordinate.\n", + " \"-c\", cmap_str,\n", + " \"-u\", \"mm\",\n", + " \"-v\", str(vmin), str(vmax),\n", + " \"--title\", \"Secular velocity\",\n", + " \"--noverbose\"\n", + "]\n", + "view.main(view_opts)" ] }, { @@ -1165,14 +1184,14 @@ "# Plot valid stations\n", "for i, gnss_stn in enumerate(gnss_stns.values()):\n", " # Plot station\n", - " ax.scatter(gnss_stn.site_lon, gnss_stn.site_lat, s=8**2, c='k',\n", + " ax.scatter(gnss_stn.site_lon, gnss_stn.site_lat, s=20, c='k',\n", " label=('valid' if i == 0 else None))\n", " ax.annotate(gnss_stn.site, (gnss_stn.site_lon, gnss_stn.site_lat))\n", "\n", "# Plot rejected stations\n", "for i, gnss_stn in enumerate(bad_stns.values()):\n", " # Plot station\n", - " ax.scatter(gnss_stn.site_lon, gnss_stn.site_lat, s=8**2,\n", + " ax.scatter(gnss_stn.site_lon, gnss_stn.site_lat, s=40,\n", " marker='X', facecolor='none', edgecolor='k',\n", " label=('rejected' if i == 0 else None))\n", " ax.annotate(gnss_stn.site, (gnss_stn.site_lon, gnss_stn.site_lat))\n", @@ -1478,7 +1497,7 @@ "source": [ "# Define requirement\n", "secular_threshold_rqmt = 2 # mm/yr\n", - "secular_distance_rqmt = [0.1, 50.0] # km\n", + "secular_distance_rqmt = [0.1, 100.0] # km\n", "\n", "# Validation parameters\n", "n_bins = 10\n", @@ -1489,10 +1508,11 @@ "double_diff_rel_measure = np.abs(np.array([diff_res.vel for diff_res in double_diffs.values()]))\n", "\n", "# Validation figure and assessment\n", + "siteloc = sitedata['sites'][site]['calval_location']\n", "(validation_table_method1,\n", " validation_fig_method1) = display_validation(site_dist, # binned distance for point\n", " double_diff_rel_measure, # binned double-difference velocities mm/yr\n", - " site, # cal/val site name\n", + " siteloc, # cal/val site name\n", " start_date, # start date of InSAR dataset\n", " end_date, # end date of InSAR dataset \n", " requirement=secular_threshold_rqmt, # measurement requirement to meet, e.g 2 mm/yr for 3 years of data over 0.1-50km\n", @@ -1670,12 +1690,14 @@ "# Validation parameters\n", "n_bins = 10\n", "threshold = 0.683 \n", + "secular_distance_rqmt = [0.1, 100.0] # km\n", "\n", "# Validation figure and assessment\n", + "siteloc = sitedata['sites'][site]['calval_location']\n", "(validation_table_method2,\n", " validation_fig_method2) = display_validation(insar_sample_dist, # binned distance for point\n", " insar_rel_measure, # binned relative velocities mm/yr\n", - " site, # cal/val site name\n", + " siteloc, # cal/val site name\n", " start_date, # start date of InSAR dataset\n", " end_date, # end date of InSAR dataset \n", " requirement=secular_threshold_rqmt, # measurement requirement to meet, e.g 2 mm/yr for 3 years of data over 0.1-50km\n", @@ -1936,6 +1958,13 @@ "prog_bar.close()\n", "plt.show()" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/methods/transient/Transient_Requirement_Validation.ipynb b/methods/transient/Transient_Requirement_Validation.ipynb index 9352e97..66f79e4 100644 --- a/methods/transient/Transient_Requirement_Validation.ipynb +++ b/methods/transient/Transient_Requirement_Validation.ipynb @@ -138,7 +138,7 @@ "outputs": [], "source": [ "# === Basic Configuration ===\n", - "site = \"test\"\n", + "site = \"CalVal_S1_LosAngelesA64\"\n", "requirement = \"Transient\"\n", "dataset = 'ARIA_S1_new' # For Sentinel-1 testing with aria-tools\n", "aria_gunw_version = \"3_0_1\"\n", @@ -1330,6 +1330,7 @@ "\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", @@ -1337,7 +1338,7 @@ "\n", " # Validation figure and assessment\n", " _, validation_fig_method1 = display_transient_validation(ddiff_dist[ifg_ndx], abs_ddiff_disp[ifg_ndx],\n", - " site, start_date, end_date,\n", + " site_loc, start_date, end_date,\n", " requirement=transient_threshold_rqmt,\n", " distance_rqmt=transient_distance_rqmt,\n", " n_bins=n_bins,\n", diff --git a/prep/ARIA_prep.ipynb b/prep/ARIA_prep.ipynb index 87a09d0..a085736 100644 --- a/prep/ARIA_prep.ipynb +++ b/prep/ARIA_prep.ipynb @@ -95,11 +95,11 @@ "source": [ "# === Basic Configuration ===\n", "site = \"test\" # Cal/Val location ID\n", - "requirement = \"Transient\" # Options: 'Secular', 'Coseismic', 'Transient'\n", + "requirement = \"Secular\" # Options: 'Secular', 'Coseismic', 'Transient'\n", "dataset = \"ARIA_S1_new\" # Dataset type: 'ARIA_S1', 'ARIA_S1_new'\n", "aria_gunw_version = \"3_0_1\"\n", "\n", - "rundate = \"20250820\" # Date of this Cal/Val run\n", + "rundate = \"20250826\" # Date of this Cal/Val run\n", "version = \"1\" # Version of this Cal/Val run\n", "custom_sites = \"/home/jovyan/my_sites.txt\" # Path to custom site metadata\n", "\n", @@ -471,8 +471,8 @@ "- mintpy.load.corFile = ../stack/cohStack.vrt\n", "- mintpy.load.connCompFile = ../stack/connCompStack.vrt\n", "- mintpy.load.demFile = ../DEM/SRTM_3arcsec.dem\n", - "- mintpy.load.incAngleFile = ../incidenceAngle/{download_start_date}_{download_edn_date}.vrt\n", - "- mintpy.load.azAngleFile = ../azimuthAngle/{download_start_date}_{download_edn_date}.vrt\n", + "- mintpy.load.incAngleFile = ../incidenceAngle/{download_start_date}_{download_end_date}.vrt\n", + "- mintpy.load.azAngleFile = ../azimuthAngle/{download_start_date}_{download_end_date}.vrt\n", "- mintpy.load.waterMaskFile = ../mask/watermask.msk\n", "- mintpy.reference.lalo = auto, or somewhere in your bounding box\n", "- mintpy.topographicResidual.pixelwiseGeometry = no\n", @@ -653,22 +653,6 @@ "print('Mintpy input files:')\n", "[x for x in os.listdir('inputs') if x.endswith('.h5')]" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6c846c09-26ea-4ea9-901c-55fd77bfcf5f", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9e9678ca-dd7d-43fb-903e-a2ce39664238", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/solid_utils/plotting.py b/solid_utils/plotting.py index 4fabee8..edac857 100644 --- a/solid_utils/plotting.py +++ b/solid_utils/plotting.py @@ -49,7 +49,9 @@ def display_validation(pair_distance: NDArray, pair_difference: NDArray, # Remove nans df_nonan = df.dropna(subset=['double_diff']) - bins = np.linspace(*distance_rqmt, num=n_bins+1) + distance_rqmt_log = [0.011,np.log10(distance_rqmt[1])] + + bins = np.logspace(*distance_rqmt_log, num=n_bins+1) bin_centers = (bins[:-1] + bins[1:]) / 2 binned_df = df_nonan.groupby(pd.cut(df_nonan['distance'], bins), observed=False)[['double_diff']] @@ -58,6 +60,7 @@ def display_validation(pair_distance: NDArray, pair_difference: NDArray, validation = pd.DataFrame([]) validation['total_count[#]'] = binned_df.apply(lambda x: np.ma.masked_invalid(x).count()) validation['passed_req.[#]'] = binned_df.apply(lambda x: np.count_nonzero(x < requirement)) + bin_counts = validation['total_count[#]'].values # Add total at the end validation = pd.concat([validation, pd.DataFrame(validation.sum(axis=0)).T]) @@ -70,28 +73,32 @@ def display_validation(pair_distance: NDArray, pair_difference: NDArray, # Figure fig, ax = plt.subplots(1, figsize=(9, 3), layout="none", dpi=200) - + ax.set_xscale('log') + # Plot residuals - ms = 8 if pair_difference.shape[0] < 1e4 else 0.3 + distance = np.log(df_nonan.distance) + ms = 6 if pair_difference.shape[0] < 1e4 else 0.3 alpha = 0.6 if pair_difference.shape[0] < 1e4 else 0.2 ax.scatter(df_nonan.distance, df_nonan.double_diff, color='black', s=ms, zorder=1, alpha=alpha, edgecolor='None') ax.fill_between(distance_rqmt, 0, requirement, color='#e6ffe6', zorder=0, alpha=0.6) - ax.fill_between(distance_rqmt, requirement, 21, color='#ffe6e6', zorder=0, alpha=0.6) - ax.vlines(bins, 0, 21, linewidth=0.3, color='gray', zorder=1) + ax.fill_between(distance_rqmt, requirement, 31, color='#ffe6e6', zorder=0, alpha=0.6) + ax.vlines(bins, 0, 31, linewidth=0.3, color='gray', zorder=1) ax.axhline(requirement, color='k', linestyle='--', zorder=3) # Bar plot for each bin quantile_th = binned_df.quantile(q=threshold)['double_diff'].values - for bin_center, quantile, flag in zip(bin_centers, + for bin_center, quantile, flag, bindiff in zip(bin_centers, quantile_th, - validation['success_fail']): + validation['success_fail'],np.diff(bins)): + if flag: color = '#227522' else: color = '#7c1b1b' - ax.bar(bin_center, quantile, align='center', width=np.diff(bins)[0], + + ax.bar(bin_center, quantile, align='center', width=bindiff, color='None', edgecolor=color, linewidth=2, zorder=3) # Add legend with data info @@ -125,6 +132,12 @@ def display_validation(pair_distance: NDArray, pair_difference: NDArray, transform=ax.transAxes) ax.add_patch(rect) + # Plot number of points in bin + ax2 = ax.twinx() + ax2.plot(bin_centers, bin_counts, marker='o', linestyle='-') + ax2.set_yscale('log') + ax2.set_ylim(1+np.min(bin_counts), 10*np.max(bin_counts)) + # Title & labels fig.suptitle(f"{validation_type.capitalize()} requirement: {site_name}", fontsize=10) ax.set_xlabel("Distance (km)", fontsize=8) @@ -136,10 +149,9 @@ def display_validation(pair_distance: NDArray, pair_difference: NDArray, ax.minorticks_on() ax.tick_params(axis='x', which='minor', length=4, direction='in', top=False, width=1.5) ax.tick_params(axis='both', labelsize=8) - ax.set_xticks(bin_centers, minor=True) - ax.set_xticks(np.arange(0,55,5)) - ax.set_ylim(0,20) - ax.set_xlim(*distance_rqmt) + ax.set_xticks([1,10,100]) + ax.set_ylim(0,30) + ax.set_xlim([1,100]) validation = validation.rename(columns={'success_fail': f'passed_req [>{threshold*100:.1f}%]'})