From 5d2be42bfd1f90a011d3d9cdaceeeb598dfb3caa Mon Sep 17 00:00:00 2001 From: Ben Purinton Date: Mon, 25 Nov 2024 19:41:36 -0800 Subject: [PATCH] Intelligent pc_align spawing, reporting, and writing out --- asp_plot/alignment.py | 91 +++++++++-------- asp_plot/altimetry.py | 49 +++++++++- notebooks/icesat2_plots.ipynb | 179 +++++++++++++++++++++++++--------- 3 files changed, 222 insertions(+), 97 deletions(-) diff --git a/asp_plot/alignment.py b/asp_plot/alignment.py index ffcb226..96011a2 100644 --- a/asp_plot/alignment.py +++ b/asp_plot/alignment.py @@ -1,5 +1,6 @@ import logging import os +import re import numpy as np from osgeo import gdal, osr @@ -15,7 +16,6 @@ def __init__( self, directory, dem_fn, - # aligned_dem_fn=None, **kwargs, ): self.directory = directory @@ -24,16 +24,6 @@ def __init__( raise ValueError(f"DEM file not found: {dem_fn}") self.dem_fn = dem_fn - # if aligned_dem_fn is not None and not os.path.exists(aligned_dem_fn): - # raise ValueError(f"Aligned DEM file not found: {aligned_dem_fn}") - # self.aligned_dem_fn = aligned_dem_fn - - # TODO: move all pc_align steps and dem translations to separate class - # call this alignment for different processing levels, with some minimum number - # of required points (a ~500 point parameter). - # Report all translation results to user via a dictionary and printing. - # If translation agrees within some tolerance (a 5% or 10% parameter) in XYZ components, - # apply the translation of the points found closest in time to the DEM. def pc_align_dem_to_atl06sr( self, max_displacement=20, @@ -78,50 +68,41 @@ def pc_align_report(self, output_prefix="pc_align/pc_align"): with open(pc_align_log, "r") as file: content = file.readlines() - report = "" + report = {} for line in content: - if "Input: error percentile of smallest errors (meters):" in line: - report_line = line.split("[ console ] : ")[1] - report += report_line - if "Input: mean of smallest errors (meters):" in line: - report_line = line.split("[ console ] : ")[1] - report += report_line - if "Output: error percentile of smallest errors (meters):" in line: - report_line = line.split("[ console ] : ")[1] - report += report_line - if "Output: mean of smallest errors (meters):" in line: - report_line = line.split("[ console ] : ")[1] - report += report_line + if "Input: error percentile" in line: + values = re.findall(r"(?:\d+%: )(\d+\.\d+)", line) + percentile_dict = { + "p16_beg": float(values[0]), + "p50_beg": float(values[1]), + "p84_beg": float(values[2]), + } + report = report | percentile_dict + if "Output: error percentile" in line: + values = re.findall(r"(?:\d+%: )(\d+\.\d+)", line) + percentile_dict = { + "p16_end": float(values[0]), + "p50_end": float(values[1]), + "p84_end": float(values[2]), + } + report = report | percentile_dict if "Translation vector (Cartesian, meters):" in line: - report_line = line.split("[ console ] : ")[1] - report += report_line ecef_shift = np.genfromtxt( [line.split("Vector3")[1][1:-2]], delimiter="," ) - report += f"\nECEF shift: {ecef_shift}\n" + xyz_shift_dict = { + "x_shift": ecef_shift[0], + "y_shift": ecef_shift[1], + "z_shift": ecef_shift[2], + } + report = report | xyz_shift_dict if "Translation vector magnitude (meters):" in line: - report_line = line.split("[ console ] : ")[1] - report += report_line + magnitude = re.findall(r"magnitude \(meters\): (\d+\.\d+)", line)[0] + report["translation_magnitude"] = float(magnitude) return report - def apply_dem_translation(self, output_prefix="pc_align/pc_align", inv_trans=True): - def get_proj_shift(src_c, src_shift, s_srs, t_srs, inv_trans=True): - if s_srs.IsSame(t_srs): - proj_shift = src_shift - else: - src_c_shift = src_c + src_shift - src2proj = osr.CoordinateTransformation(s_srs, t_srs) - proj_c = np.array(src2proj.TransformPoint(*src_c)) - proj_c_shift = np.array(src2proj.TransformPoint(*src_c_shift)) - if inv_trans: - proj_shift = proj_c - proj_c_shift - else: - proj_shift = proj_c_shift - proj_c - # Reduce unnecessary precision - proj_shift = np.around(proj_shift, 3) - return proj_shift - + def apply_dem_translation(self, output_prefix="pc_align/pc_align"): pc_align_log = glob_file(self.directory, f"{output_prefix}-log-pc_align*.txt") src = Raster(self.dem_fn) @@ -171,7 +152,7 @@ def get_proj_shift(src_c, src_shift, s_srs, t_srs, inv_trans=True): # Determine shift in original dataset coords t_srs = osr.SpatialReference() t_srs.ImportFromWkt(src.ds.crs.to_wkt()) - proj_shift = get_proj_shift(src_c, src_shift, s_srs, t_srs, inv_trans) + proj_shift = self.get_proj_shift(src_c, src_shift, s_srs, t_srs, inv_trans=True) aligned_dem_fn = self.dem_fn.replace(".tif", "_pc_align_translated.tif") print(f"\nWriting out: {aligned_dem_fn}\n") @@ -194,6 +175,22 @@ def get_proj_shift(src_c, src_shift, s_srs, t_srs, inv_trans=True): return aligned_dem_fn + def get_proj_shift(self, src_c, src_shift, s_srs, t_srs, inv_trans=True): + if s_srs.IsSame(t_srs): + proj_shift = src_shift + else: + src_c_shift = src_c + src_shift + src2proj = osr.CoordinateTransformation(s_srs, t_srs) + proj_c = np.array(src2proj.TransformPoint(*src_c)) + proj_c_shift = np.array(src2proj.TransformPoint(*src_c_shift)) + if inv_trans: + proj_shift = proj_c - proj_c_shift + else: + proj_shift = proj_c_shift - proj_c + # Reduce unnecessary precision + proj_shift = np.around(proj_shift, 3) + return proj_shift + # DEPRECATED: Very slow, better to use apply_dem_translation(). # def generate_translated_dem(self, pc_align_output, dem_out_fn): # """ diff --git a/asp_plot/altimetry.py b/asp_plot/altimetry.py index 9d4ea5d..d3c4239 100644 --- a/asp_plot/altimetry.py +++ b/asp_plot/altimetry.py @@ -208,7 +208,6 @@ def filter_esa_worldcover(self, filter_out="water"): mask = ~atl06sr["esa_worldcover.value"].isin(values) self.atl06sr_processing_levels_filtered[key] = atl06sr[mask] - # TODO: for pc_align: Spawn all four, see if they agree and provide similar translations) def predefined_temporal_filter_atl06sr(self, date=None): if date is None: date = StereopairMetadataParser(self.directory).get_pair_dict()["cdate"] @@ -285,7 +284,9 @@ def alignment_report( self, processing_level="high_confidence", minimum_points=500, + agreement_threshold=0.25, write_out_aligned_dem=False, + min_translation_threshold=0.1, key_for_aligned_dem="high_confidence", ): filtered_keys = [ @@ -318,11 +319,53 @@ def alignment_report( output_prefix=f"pc_align/pc_align_{key}", ) + report_data = [] for key in filtered_keys: report = alignment.pc_align_report(output_prefix=f"pc_align/pc_align_{key}") - print(f"\n{key} alignment report:\n{report}\n") + report_data.append({"key": key} | report) + alignment_report_df = pd.DataFrame(report_data) + + gsd = Raster(self.dem_fn).get_gsd() + if ( + alignment_report_df["translation_magnitude"].mean() + < min_translation_threshold * gsd + ): + write_out_aligned_dem = False + print( + f"\nTranslation magnitude is less than {min_translation_threshold*100}% of the DEM GSD. Skipping writing out aligned DEM.\n" + ) if write_out_aligned_dem: + # Calculate ranges and mean for each shift component + x_range = ( + alignment_report_df["x_shift"].max() + - alignment_report_df["x_shift"].min() + ) + y_range = ( + alignment_report_df["y_shift"].max() + - alignment_report_df["y_shift"].min() + ) + z_range = ( + alignment_report_df["z_shift"].max() + - alignment_report_df["z_shift"].min() + ) + x_mean = alignment_report_df["x_shift"].mean() + y_mean = alignment_report_df["y_shift"].mean() + z_mean = alignment_report_df["z_shift"].mean() + + # Check if range is more than X% of mean for any component + if ( + x_range > abs(x_mean * agreement_threshold) + or y_range > abs(y_mean * agreement_threshold) + or z_range > abs(z_mean * agreement_threshold) + ): + print( + f"\nWarning: Translation components vary by more than {agreement_threshold*100}% across temporal filters. The translation applied to the aligned DEM may be inaccurate.\n" + ) + print(f"X shift range: {x_range:.3f} m (mean: {x_mean:.3f} m)") + print(f"Y shift range: {y_range:.3f} m (mean: {y_mean:.3f} m)") + print(f"Z shift range: {z_range:.3f} m (mean: {z_mean:.3f} m)") + self.aligned_dem_fn = alignment.apply_dem_translation( output_prefix=f"pc_align/pc_align_{key_for_aligned_dem}", ) @@ -330,6 +373,8 @@ def alignment_report( f"\nWrote out {key_for_aligned_dem} aligned DEM to {self.aligned_dem_fn}\n" ) + self.alignment_report_df = alignment_report_df + def plot_atl06sr_time_stamps( self, key="high_confidence", diff --git a/notebooks/icesat2_plots.ipynb b/notebooks/icesat2_plots.ipynb index 14f9495..9d15176 100644 --- a/notebooks/icesat2_plots.ipynb +++ b/notebooks/icesat2_plots.ipynb @@ -448,7 +448,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 12, "metadata": {}, "outputs": [ { @@ -456,53 +456,11 @@ "output_type": "stream", "text": [ "\n", - "high_confidence alignment report:\n", - "Input: error percentile of smallest errors (meters): 16%: 0.379974, 50%: 0.613481, 84%: 0.823963\n", - "Input: mean of smallest errors (meters): 25%: 0.305042, 50%: 0.422038, 75%: 0.509095, 100%: 0.664736\n", - "Output: error percentile of smallest errors (meters): 16%: 0.0449552, 50%: 0.151377, 84%: 0.332819\n", - "Output: mean of smallest errors (meters): 25%: 0.0355416, 50%: 0.0730141, 75%: 0.116036, 100%: 0.267471\n", - "Translation vector (Cartesian, meters): Vector3(0.092759045,0.68025135,-0.5164948)\n", - "\n", - "ECEF shift: [ 0.09275904 0.68025135 -0.5164948 ]\n", - "Translation vector magnitude (meters): 0.85913504\n", - "\n", - "\n", - "\n", - "high_confidence_15_day_pad alignment report:\n", - "Input: error percentile of smallest errors (meters): 16%: 0.633969, 50%: 0.789984, 84%: 0.918273\n", - "Input: mean of smallest errors (meters): 25%: 0.577839, 50%: 0.65724, 75%: 0.71529, 100%: 0.854417\n", - "Output: error percentile of smallest errors (meters): 16%: 0.030055, 50%: 0.094727, 84%: 0.207692\n", - "Output: mean of smallest errors (meters): 25%: 0.0230105, 50%: 0.0470604, 75%: 0.073551, 100%: 0.228624\n", - "Translation vector (Cartesian, meters): Vector3(0.48710351,-0.0095034486,-0.67688551)\n", - "\n", - "ECEF shift: [ 0.48710351 -0.00950345 -0.67688551]\n", - "Translation vector magnitude (meters): 0.83398689\n", - "\n", - "\n", - "\n", - "high_confidence_45_day_pad alignment report:\n", - "Input: error percentile of smallest errors (meters): 16%: 0.641462, 50%: 0.79161, 84%: 0.923064\n", - "Input: mean of smallest errors (meters): 25%: 0.57232, 50%: 0.657621, 75%: 0.716122, 100%: 0.860368\n", - "Output: error percentile of smallest errors (meters): 16%: 0.028541, 50%: 0.0919282, 84%: 0.218268\n", - "Output: mean of smallest errors (meters): 25%: 0.0219495, 50%: 0.0453428, 75%: 0.0723805, 100%: 0.233776\n", - "Translation vector (Cartesian, meters): Vector3(0.51237937,0.16143899,-0.65319019)\n", - "\n", - "ECEF shift: [ 0.51237937 0.16143899 -0.65319019]\n", - "Translation vector magnitude (meters): 0.84572607\n", - "\n", - "\n", - "\n", - "high_confidence_seasonal alignment report:\n", - "Input: error percentile of smallest errors (meters): 16%: 0.566842, 50%: 0.732428, 84%: 0.888689\n", - "Input: mean of smallest errors (meters): 25%: 0.507818, 50%: 0.594075, 75%: 0.657678, 100%: 0.777453\n", - "Output: error percentile of smallest errors (meters): 16%: 0.03194, 50%: 0.109037, 84%: 0.240184\n", - "Output: mean of smallest errors (meters): 25%: 0.0251211, 50%: 0.0523517, 75%: 0.083187, 100%: 0.197556\n", - "Translation vector (Cartesian, meters): Vector3(0.22125086,-0.24869599,-0.73254096)\n", - "\n", - "ECEF shift: [ 0.22125086 -0.24869599 -0.73254096]\n", - "Translation vector magnitude (meters): 0.80462283\n", - "\n", + "Warning: Translation components vary by more than 25.0% across temporal filters. The translation applied to the aligned DEM may be inaccurate.\n", "\n", + "X shift range: 0.420 m (mean: 0.328 m)\n", + "Y shift range: 0.929 m (mean: 0.146 m)\n", + "Z shift range: 0.216 m (mean: -0.645 m)\n", "\n", "Writing out: /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/20220417_2252_1040010074793300_1040010075633C00-DEM_1m_pc_align_translated.tif\n", "\n", @@ -516,14 +474,139 @@ "icesat.alignment_report(\n", " processing_level=\"high_confidence\",\n", " minimum_points=500,\n", + " agreement_threshold=0.25,\n", " write_out_aligned_dem=True,\n", + " min_translation_threshold=0.1,\n", " key_for_aligned_dem=\"high_confidence_15_day_pad\",\n", ")" ] }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
keyp16_begp50_begp84_begp16_endp50_endp84_endx_shifty_shiftz_shifttranslation_magnitude
0high_confidence0.3799740.6134810.8239630.0449550.1513770.3328190.0927590.680251-0.5164950.859135
1high_confidence_15_day_pad0.6339690.7899840.9182730.0300550.0947270.2076920.487104-0.009503-0.6768860.833987
2high_confidence_45_day_pad0.6414620.7916100.9230640.0285410.0919280.2182680.5123790.161439-0.6531900.845726
3high_confidence_seasonal0.5668420.7324280.8886890.0319400.1090370.2401840.221251-0.248696-0.7325410.804623
\n", + "
" + ], + "text/plain": [ + " key p16_beg p50_beg p84_beg p16_end \\\n", + "0 high_confidence 0.379974 0.613481 0.823963 0.044955 \n", + "1 high_confidence_15_day_pad 0.633969 0.789984 0.918273 0.030055 \n", + "2 high_confidence_45_day_pad 0.641462 0.791610 0.923064 0.028541 \n", + "3 high_confidence_seasonal 0.566842 0.732428 0.888689 0.031940 \n", + "\n", + " p50_end p84_end x_shift y_shift z_shift translation_magnitude \n", + "0 0.151377 0.332819 0.092759 0.680251 -0.516495 0.859135 \n", + "1 0.094727 0.207692 0.487104 -0.009503 -0.676886 0.833987 \n", + "2 0.091928 0.218268 0.512379 0.161439 -0.653190 0.845726 \n", + "3 0.109037 0.240184 0.221251 -0.248696 -0.732541 0.804623 " + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "icesat.alignment_report_df" + ] + }, + { + "cell_type": "code", + "execution_count": 14, "metadata": {}, "outputs": [ {