Updating permafrost workflow (#98)

methods/secular/Secular_Requirement_Validation.ipynb

@@ -1984,7 +1984,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.13.3"
+   "version": "3.13.5"
   },
   "toc": {
    "base_numbering": 1,

my_sites.txt

@@ -611,6 +611,22 @@
             "SET_source" : "mintpy",
             "do_iono" : "True",
             "do_tropo" : "True",
-            "tropo_model" : "ERA5"}
+            "tropo_model" : "ERA5"},
+
+    "Info" : "######## PERMAFROST #########",
+
+    "NorthSlopeEastD102" : {
+                "calval_location" : "NorthSlopeEastD102",
+                "region_identifier" : "POINT(-149.37 69.09)",
+                "subset_region" : "[7620213:7686754, 641941:679925]",            
+                "download_start_date" : "20230525",
+                "download_end_date"  : "20230910",
+                "mintpy_ref_loc"      : "7651392, 666923",
+                "tempBaseMax" : "36",
+                "ifgExcludeList" : "auto",
+                "maskWater" : "True",
+                "sentinel_direction" : "DESCENDING",
+                "sentinel_path" : "102",
+                "sentinel_frame" : "362"}
 }
 }

prep/ARIA_prep.ipynb

@@ -671,7 +671,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.13.3"
+   "version": "3.13.5"
   }
  },
  "nbformat": 4,

prep/NISAR_prep.ipynb

@@ -305,9 +305,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "mintpy_dev",
+   "display_name": "Python (solid_earth_atbd_dev)",
    "language": "python",
-   "name": "mintpy_dev"
+   "name": "solid_earth_atbd_dev"
   },
   "language_info": {
    "codemirror_mode": {
@@ -319,7 +319,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.8"
+   "version": "3.13.5"
   }
  },
  "nbformat": 4,

solid_utils/plotting.py

@@ -15,7 +15,7 @@ def display_validation(pair_distance: NDArray, pair_difference: NDArray,
                        requirement: float = 2, distance_rqmt: list = [0.1, 50],
                        n_bins: int = 10, threshold: float = 0.683, 
                        sensor:str ='Sentinel-1', validation_type:str='secular',
-                       validation_data:str='GNSS'):
+                       validation_data:str='Field_meas'):
 
     """Display double-difference validation results. Bin the double differences
     as a function of point separation, compute the stastics for each bin, and
@@ -514,5 +514,151 @@ def style_specific_cells(val):
     # plt.savefig(f'transient_validation_{index[i]}.png', bbox_inches='tight', transparent=True)
 
     plt.close()
-    return styled_df, fig  
+    return styled_df, fig 
+
+def display_permafrost_validation(pair_distance: NDArray, pair_difference: NDArray,
+                       site_name: str, start_date: str, end_date: str,
+                       requirement: float = 2, req_dist_fcn: bool = False,
+                       distance_rqmt: list = [0.1, 50], n_bins: int = 10, threshold: float = 0.8, 
+                       sensor:str ='Sentinel-1', validation_type:str='permafrost',
+                       validation_data:str='field'):
+   '''
+    Parameters:
+      pair_distance : array      - 1d array of pair distances used in validation
+      pair_difference : array    - 1d array 0f pair double differenced velocity residuals
+      site_name : str            - name of the cal/val site
+      start_date  : str          - data record start date, eg. 20190101
+      end_date : str             - data record end date, eg. 20200101
+      requirement : float        - value required for test to pass
+                                    e.g, 2 mm/yr for 3 years of data over distance requiremeent
+      req_dist_fcn: bool         - flag to scale requirement with distance. If True, requirement = value * (1+sqrt(L))
+      distance_rqmt : list       - distance over requirement is tested, eg. length scales of 0.1-50 km
+      n_bins : int               - number of bins
+      threshold : float          - threshold represents percentile of Gaussian normal distribution
+                                    within residuals are expected to be to pass the test
+                                    e.g. 0.8 for 80%
+      sensor : str               - sensor used in validation, e.g Sentinel-1 or NISAR
+      validation_type : str      - type of validation: permafrost
+      validation_data : str      - data used to validate against; field or INSAR
+
+   Return
+      validation_table
+      validation_figure
+   '''
+   # init dataframe
+   pair_req_met = np.array(pair_difference < requirement*(1+float(req_dist_fcn)*np.sqrt(pair_distance)))
+   pair_req_met[np.isnan(pair_difference)]=np.nan
+   df = pd.DataFrame(np.vstack([pair_distance,
+                                pair_difference,
+                                pair_req_met]).T,
+                                columns=['distance', 'double_diff','req_met'])
+
+   # remove nans
+   df_nonan = df.dropna(subset=['double_diff'])
+   bins = np.linspace(*distance_rqmt, num=n_bins+1)
+   bin_centers = (bins[:-1] + bins[1:]) / 2
+   binned_df = df_nonan.groupby(pd.cut(df_nonan['distance'], bins),
+                                observed=False)[['req_met']]
+   binned_df_diff = df_nonan.groupby(pd.cut(df_nonan['distance'], bins),
+                                observed=False)[['double_diff']]
+
+   # get binned validation table 
+   bin_req = requirement*(1+float(req_dist_fcn)*np.sqrt(bin_centers))
+   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))
+
+   # Add total at the end
+   validation = pd.concat([validation, pd.DataFrame(validation.sum(axis=0)).T])
+   validation['passed_pc'] = validation['passed_req.[#]'] / validation['total_count[#]']
+   validation['success_fail'] = validation['passed_pc'] > threshold
+   validation.index.name = 'distance[km]'
+   # Rename last row
+   validation.rename({validation.iloc[-1].name:'Total'}, inplace=True)
+
+   # Figure
+   fig, ax = plt.subplots(1, figsize=(9, 3), layout="none", dpi=200)
+   ymax = 20
+   if validation_type=='permafrost':
+       ymax=50
+
+
+   # Plot residuals
+   ms = 8 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')
+
+   for i,r in enumerate(bin_req):
+       ibin = [bins[i],bins[i+1]]
+       ax.fill_between(ibin, 0, r, color='#e6ffe6', zorder=0, alpha=0.6)
+       ax.fill_between(ibin, r, 51, color='#ffe6e6', zorder=0, alpha=0.6)
+   ax.vlines(bins, 0, ymax+1, linewidth=0.3, color='gray', zorder=1)
+
+   req_line_x = np.linspace(*distance_rqmt, num=200)
+   req_line_y = requirement*(1+req_dist_fcn*np.sqrt(req_line_x))
+   ax.plot(req_line_x,req_line_y, color='k', linestyle='--', zorder=3)
+
+   # Bar plot for each bin
+   quantile_th = binned_df_diff.quantile(q=threshold)['double_diff'].values
+   for bin_center, quantile, flag in zip(bin_centers,
+                                         quantile_th,
+                                         validation['success_fail']):
+      if flag:
+         color = '#227522'
+      else:
+         color = '#7c1b1b'
+      ax.bar(bin_center, quantile, align='center', width=np.diff(bins)[0],
+            color='None', edgecolor=color, linewidth=2, zorder=3)
+
+   # Add legend with data info
+   legend_kwargs = dict(transform=ax.transAxes, verticalalignment='top')
+   props = dict(boxstyle='square', facecolor='white', alpha=1, linewidth=0.4)
+   textstr = f'Sensor: {sensor} \n{validation_data}-InSAR point pairs\n'
+   textstr += f'Record: {start_date}-{end_date}'
+
+   # place a text box in upper left in axes coords
+   ax.text(0.02, 0.95, textstr, fontsize=8, bbox=props, **legend_kwargs)
+
+   # Add legend with validation info 
+   textstr = f'{validation_type.capitalize()} requirement\n'
+   textstr += f'Site: {site_name}\n'
+   if validation.loc['Total']['success_fail']:
+      validation_flag = 'PASSED'
+      validation_color = '#239d23'
+   else: 
+      validation_flag ='FAILED'
+      validation_color = '#bc2e2e'
+
+   props = {**props, **{'facecolor':'none', 'edgecolor':'none'}}
+   ax.text(0.818, 0.93, textstr, fontsize=8, bbox=props, **legend_kwargs)
+   ax.text(0.852, 0.82,  f"{validation_flag}",
+           fontsize=10, weight='bold',
+           bbox=props, **legend_kwargs)
+
+   rect = patches.Rectangle((0.8, 0.75), 0.19, 0.2,
+                           linewidth=1, edgecolor='black',
+                           facecolor=validation_color,
+                           transform=ax.transAxes)
+   ax.add_patch(rect)
+
+   # Title & labels
+   fig.suptitle(f"{validation_type.capitalize()} requirement: {site_name}", fontsize=10)
+   ax.set_xlabel("Distance (km)", fontsize=8)
+   if validation_data == 'GNSS':
+       txt = "Double-Differenced \nVelocity Residual (mm/yr)"
+   else:
+       txt = "Relative displacement measurement (mm)"    
+   ax.set_ylabel(txt, fontsize=8)
+   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,ymax)
+   ax.set_xlim(*distance_rqmt)
+
+   validation = validation.rename(columns={'success_fail': f'passed_req [>{threshold*100:.1f}%]'})
+
+   return validation, fig
 

solid_utils/saving.py

@@ -127,3 +127,129 @@ def save_results(save_dir:str, run_date:str, requirement:str, site:str, method:s
 
     # Save to zip file
     shutil.make_archive(save_dir, 'zip', save_dir)
+
+def save_results_permafrost(save_dir:str, run_date:str, requirement:str, site:str, method:str,
+                 sitedata:dict, gnss_insar_figs, validation_figs:list, validation_table:pd.DataFrame,
+                 ts_functions:dict=None, summary:str=""):
+    """Save the input parameters and results to an output folder.
+    The results folder will include image files, and an HTML file
+    recording the inputs, essential figures, and validation table.
+    """
+    # Create new directory for outputs
+    if os.path.exists(save_dir):
+        print(f"Directory {save_dir} exists")
+    else:
+        print(f"Creating {os.path.basename(save_dir)}")
+        os.makedirs(save_dir)
+
+    # Begin HTML string
+    html_str = """<!DOCTYPE html>
+<html lang="en">
+<head>
+    <meta charset="UTF-8">
+    <title>{requirement:s} Method {method}</title>
+</head>
+<body>
+    <h1>{requirement:s} requirement validation 
+{site:s} site, Method {method}</h1>
+</body>
+""".format(requirement=requirement.title(),
+           site=site,
+           method=method)
+
+    # Run date time
+    html_str += """<body>
+    <h2>Run time:</h2>
+    {run_date:s}
+    </body>
+""".format(run_date=run_date)
+
+    # Site parameters
+    html_str += "<body>"
+    html_str += """<h2>Setup parameters:</h2>
+<ul>
+"""
+    for key, value in sitedata.items():
+        html_str += "<li>{}: {}</li>".format(key, value)
+    html_str += """</ul>
+</body>"""
+
+    # Timeseries basis functions
+    if ts_functions is not None:
+        html_str += "<body>"
+        html_str += """<h3>Timeseries basis functions</h3>
+    <ul>
+    """
+        for key, value in ts_functions.items():
+            html_str += "<li>{}: {}</li>".format(key, value)
+        html_str += """</ul>
+</body>"""
+
+    # Processing results
+    html_str += """<body>
+<h2>InSAR and GNSS LOS Velocities</h2>
+Visual comparison of GNSS and InSAR LOS velocities
+<br>
+</body>
+""".format(method=method)
+
+    # Save GNSS InSAR figures
+    for i, gnss_insar_fig in enumerate(gnss_insar_figs):
+        fig_name = f"{requirement:s}_Method{method}_gnss_insar_figure{i + 1:d}.png"
+        fig_path = os.path.join(save_dir, fig_name)
+        print(f"Saving GNSS-InSAR figure to: {fig_name}")
+        gnss_insar_fig.savefig(fig_path, bbox_inches='tight', transparent=True,
+                               dpi=300)
+
+        # Embed GNSS-InSAR figures
+        html_str += """<body>
+<img src="{fig_name:s}" alt="Method {method} GNSS InSAR Image" width="800">
+</body>
+""".format(fig_name=fig_name,
+           method=method)
+
+    # Validation results
+    html_str += """<body>
+<h2>Method {method} Results</h2>
+</body>
+""".format(method=method)
+
+    # Save validation figures
+    for i, validation_fig in enumerate(validation_figs):
+        fig_name = f"{requirement:s}_Method{method}_validation_figure{i + 1:d}.png"
+        fig_path = os.path.join(save_dir, fig_name)
+        print(f"Saving validation figure to: {fig_name}")
+        validation_fig.savefig(fig_path, bbox_inches='tight', transparent=True,
+                               dpi=300)
+
+        # Embed validation figure
+        html_str += """<body>
+    <img src="{fig_name:s}" alt="Method {method} Validation Image" width="800">
+    </body>
+    """.format(fig_name=fig_name,
+               method=method)
+
+    # Validation table
+    html_str += """<body>
+{}
+</body>""".format(validation_table.to_html())
+
+    # Save summary
+    html_str += """<body>
+<h2>Summary</h2>
+{:s}
+</body>""".format(summary)
+
+    # Write to HTML file
+    html_name = f"{requirement:s}_Method{method}_validation_report.html"
+    html_path = os.path.join(save_dir, html_name)
+    with open(html_path, 'w') as html_file:
+        html_file.write(html_str)
+
+    # pdf_file = html_path.split(".")[0]+".pdf"
+    # HTML(html_path).write_pdf(pdf_file)
+    # print(f"Saved PDF version of report to: {pdf_file}")
+    print(f"Saved parameters and results to: {save_dir:s}")
+
+    # Save to zip file
+    shutil.make_archive(save_dir, 'zip', save_dir)