Plot bundle adjust pointmap geodiff results

asp_plot/bundle_adjust.py

@@ -8,30 +8,31 @@
 from asp_plot.utils import ColorBar, Plotter, save_figure
 
 
-class ReadResiduals:
+class ReadBundleAdjustFiles:
     def __init__(self, directory, bundle_adjust_directory):
         self.directory = directory
         self.bundle_adjust_directory = bundle_adjust_directory
 
-    def get_init_final_residuals_csvs(self):
+    def get_csv_paths(self, geodiff_files=False):
         filenames = [
             "*-initial_residuals_pointmap.csv",
             "*-final_residuals_pointmap.csv",
         ]
 
+        if geodiff_files:
+            filenames = [f.replace(".csv", "-diff.csv") for f in filenames]
+
         paths = [
-            glob.glob(
-                os.path.join(self.directory, self.bundle_adjust_directory, filename)
-            )[0]
-            for filename in filenames
+            glob.glob(os.path.join(self.directory, self.bundle_adjust_directory, f))[0]
+            for f in filenames
         ]
 
         for path in paths:
             if not os.path.isfile(path):
-                raise ValueError(f"Residuals CSV file not found: {path}")
+                raise ValueError(f"CSV file not found: {path}")
 
-        resid_init_path, resid_final_path = paths
-        return resid_init_path, resid_final_path
+        initial, final = paths
+        return initial, final
 
     def get_residuals_gdf(self, csv_path):
         cols = [
@@ -66,11 +67,38 @@ def get_residuals_gdf(self, csv_path):
         resid_gdf.filename = os.path.basename(csv_path)
         return resid_gdf
 
-    def get_init_final_residuals_gdfs(self):
-        resid_init_path, resid_final_path = self.get_init_final_residuals_csvs()
-        resid_init_gdf = self.get_residuals_gdf(resid_init_path)
+    def get_geodiff_gdf(self, csv_path):
+        cols = [
+            "lon",
+            "lat",
+            "height_diff_meters",
+        ]
+
+        geodiff_df = pd.read_csv(csv_path, skiprows=7, names=cols)
+
+        geodiff_gdf = gpd.GeoDataFrame(
+            geodiff_df,
+            geometry=gpd.points_from_xy(
+                geodiff_df["lon"], geodiff_df["lat"], crs="EPSG:4326"
+            ),
+        )
+
+        geodiff_gdf.filename = os.path.basename(csv_path)
+        return geodiff_gdf
+
+    def get_initial_final_residuals_gdfs(self):
+        resid_initial_path, resid_final_path = self.get_csv_paths()
+        resid_initial_gdf = self.get_residuals_gdf(resid_initial_path)
         resid_final_gdf = self.get_residuals_gdf(resid_final_path)
-        return resid_init_gdf, resid_final_gdf
+        return resid_initial_gdf, resid_final_gdf
+
+    def get_initial_final_geodiff_gdfs(self):
+        geodiff_initial_path, geodiff_final_path = self.get_csv_paths(
+            geodiff_files=True
+        )
+        geodiff_initial_gdf = self.get_geodiff_gdf(geodiff_initial_path)
+        geodiff_final_gdf = self.get_geodiff_gdf(geodiff_final_path)
+        return geodiff_initial_gdf, geodiff_final_gdf
 
     def get_mapproj_residuals_gdf(self):
         path = glob.glob(
@@ -123,18 +151,18 @@ def get_propagated_triangulation_uncert_df(self):
         return resid_triangulation_uncert_df
 
 
-class PlotResiduals(Plotter):
+class PlotBundleAdjustFiles(Plotter):
     def __init__(self, geodataframes, **kwargs):
         super().__init__(**kwargs)
         if not isinstance(geodataframes, list):
             raise ValueError("Input must be a list of GeoDataFrames")
         self.geodataframes = geodataframes
 
-    def get_residual_stats(self, gdf, column_name="mean_residual"):
+    def gdf_percentile_stats(self, gdf, column_name="mean_residual"):
         stats = gdf[column_name].quantile([0.25, 0.50, 0.84, 0.95]).round(2).tolist()
         return stats
 
-    def plot_n_residuals(
+    def plot_n_gdfs(
         self,
         column_name="mean_residual",
         cbar_label="Mean Residual (m)",
@@ -175,10 +203,15 @@ def plot_n_residuals(
                     clim=clim,
                     column_name=column_name,
                     cbar_label=cbar_label,
+                    cmap=cmap,
                 )
             else:
                 self.plot_geodataframe(
-                    ax=axa[i], gdf=gdf, column_name=column_name, cbar_label=cbar_label
+                    ax=axa[i],
+                    gdf=gdf,
+                    column_name=column_name,
+                    cbar_label=cbar_label,
+                    cmap=cmap,
                 )
 
             ctx.add_basemap(ax=axa[i], **ctx_kwargs)
@@ -187,7 +220,7 @@ def plot_n_residuals(
                 axa[i].autoscale(False)
 
             # Show some statistics and information
-            stats = self.get_residual_stats(gdf, column_name)
+            stats = self.gdf_percentile_stats(gdf, column_name)
             stats_text = f"(n={gdf.shape[0]})\n" + "\n".join(
                 f"{quantile*100:.0f}th: {stat}"
                 for quantile, stat in zip([0.25, 0.50, 0.84, 0.95], stats)

asp_plot/cli/asp_plot.py

@@ -3,9 +3,10 @@
 import subprocess
 import click
 import contextily as ctx
+from itertools import count
 from asp_plot.processing_parameters import ProcessingParameters
 from asp_plot.scenes import ScenePlotter, SceneGeometryPlotter
-from asp_plot.bundle_adjust import ReadResiduals, PlotResiduals
+from asp_plot.bundle_adjust import ReadBundleAdjustFiles, PlotBundleAdjustFiles
 from asp_plot.stereo import StereoPlotter
 from asp_plot.utils import compile_report
 
@@ -44,7 +45,7 @@
 @click.option(
     "--plots_directory",
     prompt=True,
-    default="asp_plots",
+    default="asp_plots_for_report",
     help="Directory to put output plots. Default: asp_plots",
 )
 @click.option(
@@ -68,20 +69,23 @@ def main(
     os.makedirs(plots_directory, exist_ok=True)
     report_pdf_path = os.path.join(directory, report_filename)
 
+    figure_counter = count(0)
+
     # Geometry plot
     plotter = SceneGeometryPlotter(directory)
 
-    plotter.dg_geom_plot(save_dir=plots_directory, fig_fn="00_geometry.png")
+    plotter.dg_geom_plot(save_dir=plots_directory, fig_fn=f"{next(figure_counter):02}.png")
 
     # Scene plot
     plotter = ScenePlotter(directory, stereo_directory, title="Mapprojected Scenes")
 
-    plotter.plot_orthos(save_dir=plots_directory, fig_fn="01_orthos.png")
+    plotter.plot_orthos(save_dir=plots_directory, fig_fn=f"{next(figure_counter):02}.png")
 
     # Bundle adjustment plots
-    residuals = ReadResiduals(directory, bundle_adjust_directory)
-    resid_init_gdf, resid_final_gdf = residuals.get_init_final_residuals_gdfs()
-    resid_mapprojected_gdf = residuals.get_mapproj_residuals_gdf()
+    ba_files = ReadBundleAdjustFiles(directory, bundle_adjust_directory)
+    resid_initial_gdf, resid_final_gdf = ba_files.get_initial_final_residuals_gdfs()
+    geodiff_initial_gdf, geodiff_final_gdf = ba_files.get_initial_final_geodiff_gdfs()
+    resid_mapprojected_gdf = ba_files.get_mapproj_residuals_gdf()
 
     ctx_kwargs = {
         "crs": map_crs,
@@ -90,48 +94,69 @@ def main(
         "alpha": 0.5,
     }
 
-    plotter = PlotResiduals(
-        [resid_init_gdf, resid_final_gdf],
+    plotter = PlotBundleAdjustFiles(
+        [resid_initial_gdf, resid_final_gdf],
         lognorm=True,
         title="Bundle Adjust Initial and Final Residuals (Log Scale)",
     )
 
-    plotter.plot_n_residuals(
+    plotter.plot_n_gdfs(
         column_name="mean_residual",
         cbar_label="Mean Residual (m)",
         map_crs=map_crs,
         save_dir=plots_directory,
-        fig_fn="02_ba_residuals_log.png",
+        fig_fn=f"{next(figure_counter):02}.png",
         **ctx_kwargs,
     )
 
     plotter.lognorm = False
     plotter.title = "Bundle Adjust Initial and Final Residuals (Linear Scale)"
 
-    plotter.plot_n_residuals(
+    plotter.plot_n_gdfs(
         column_name="mean_residual",
         cbar_label="Mean Residual (m)",
         common_clim=False,
         map_crs=map_crs,
         save_dir=plots_directory,
-        fig_fn="03_ba_residuals_linear.png",
+        fig_fn=f"{next(figure_counter):02}.png",
         **ctx_kwargs,
     )
 
-    plotter = PlotResiduals(
+    plotter = PlotBundleAdjustFiles(
         [resid_mapprojected_gdf],
         title="Bundle Adjust Midpoint distance between\nfinal interest points projected onto reference DEM",
     )
 
-    plotter.plot_n_residuals(
+    plotter.plot_n_gdfs(
         column_name="mapproj_ip_dist_meters",
         cbar_label="Interest point distance (m)",
         map_crs=map_crs,
         save_dir=plots_directory,
-        fig_fn="04_ba_residuals_mapproj_dist.png",
+        fig_fn=f"{next(figure_counter):02}.png",
         **ctx_kwargs,
     )
 
+    plotter = PlotBundleAdjustFiles(
+        [geodiff_initial_gdf, geodiff_final_gdf],
+        lognorm=False,
+        title="Bundle Adjust Initial and Final Geodiff vs. Reference DEM"
+    )
+
+    clim = (float(geodiff_initial_gdf["height_diff_meters"].quantile(0.05)), float(geodiff_initial_gdf["height_diff_meters"].quantile(0.95)))
+    abs_max = max(abs(clim[0]), abs(clim[1]))
+    clim = (-abs_max, abs_max)
+
+    plotter.plot_n_gdfs(
+        column_name="height_diff_meters",
+        cbar_label="Height difference (m)",
+        map_crs=map_crs,
+        cmap="RdBu",
+        clim=clim,
+        save_dir=plots_directory,
+        fig_fn=f"{next(figure_counter):02}.png",
+        **ctx_kwargs
+    )
+
     # Stereo plots
     plotter = StereoPlotter(
         directory,
@@ -143,7 +168,7 @@ def main(
 
     plotter.plot_match_points(
         save_dir=plots_directory,
-        fig_fn="05_stereo_match_points.png",
+        fig_fn=f"{next(figure_counter):02}.png",
     )
 
     plotter.title = "Disparity (pixels)"
@@ -152,14 +177,14 @@ def main(
         unit="pixels",
         quiver=True,
         save_dir=plots_directory,
-        fig_fn="06_disparity_pixels.png",
+        fig_fn=f"{next(figure_counter):02}.png",
     )
 
     plotter.title = "Stereo DEM Results"
 
     plotter.plot_dem_results(
         save_dir=plots_directory,
-        fig_fn="07_stereo_dem_results.png",
+        fig_fn=f"{next(figure_counter):02}.png",
     )
 
     # Compile report

notebooks/stereo_plots.ipynb

@@ -8,8 +8,7 @@
    "source": [
     "import os\n",
     "from asp_plot.stereo import StereoPlotter\n",
-    "from asp_plot.processing_parameters import ProcessingParameters\n",
-    "from asp_plot.utils import compile_report"
+    "from asp_plot.processing_parameters import ProcessingParameters"
    ]
   },
   {
@@ -23,10 +22,8 @@
     "ba_directory = \"ba/\"\n",
     "reference_dem = \"/Users/ben/Dropbox/UW_Shean/COP/COP30_utqiagvik_lzw-adj_proj.tif\"\n",
     "\n",
-    "plots_directory = os.path.join(directory, \"stereo_report_plots\")\n",
-    "os.makedirs(plots_directory, exist_ok=True)\n",
-    "\n",
-    "report_pdf_path = os.path.join(directory, \"asp_report_WV03_20220417.pdf\")\n"
+    "plots_directory = os.path.join(directory, \"asp_plots\")\n",
+    "os.makedirs(plots_directory, exist_ok=True)"
    ]
   },
   {
@@ -63,7 +60,7 @@
     "\n",
     "plotter.plot_match_points(\n",
     "  save_dir=plots_directory,\n",
-    "  fig_fn=\"05_stereo_match_points.png\",\n",
+    "  fig_fn=\"stereo_match_points.png\",\n",
     ")"
    ]
   },
@@ -122,7 +119,7 @@
     "  unit=\"pixels\",\n",
     "  quiver=True,\n",
     "  save_dir=plots_directory,\n",
-    "  fig_fn=\"06_disparity_pixels.png\",\n",
+    "  fig_fn=\"disparity_pixels.png\",\n",
     ")"
    ]
   },
@@ -155,7 +152,7 @@
     "\n",
     "plotter.plot_dem_results(\n",
     "  save_dir=plots_directory,\n",
-    "  fig_fn=\"07_stereo_dem_results.png\",\n",
+    "  fig_fn=\"stereo_dem_results.png\",\n",
     ")"
    ]
   },
@@ -191,18 +188,6 @@
     "  diff_clim=(-5, 5),\n",
     ")"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 3,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "processing_parameters = ProcessingParameters(directory, ba_directory, stereo_directory)\n",
-    "processing_parameters_dict = processing_parameters.from_log_files()\n",
-    "\n",
-    "compile_report(plots_directory, processing_parameters_dict, report_pdf_path)"
-   ]
   }
  ],
  "metadata": {