Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file not shown.
197 changes: 170 additions & 27 deletions methods/coseismic/Coseismic_Requirement_Validation.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand All @@ -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)"
]
},
{
Expand Down Expand Up @@ -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)"
]
},
{
Expand Down Expand Up @@ -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)"
]
},
Expand Down Expand Up @@ -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)"
]
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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)"
]
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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)"
]
},
{
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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)"
]
},
{
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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)"
]
},
Expand All @@ -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)"
]
},
{
Expand Down Expand Up @@ -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",
Expand Down
Loading