Skip to content

Commit

Permalink
Plot bundle adjust pointmap geodiff results (#27)
Browse files Browse the repository at this point in the history 
Typically, when running stereo processing, we also run `geodiff --csv-format "1:lon,2:lat,3:height_above_datum" pointmap.csv ref_dem.tif`. This PR resolves #7 by reading and plotting these files.
  • Loading branch information
@bpurinton
bpurinton committed Jul 10, 2024
1 parent 22f8a67 commit 4b95766
Show file tree
Hide file tree
Showing 9 changed files with 250 additions and 129 deletions.
69 changes: 51 additions & 18 deletions asp_plot/bundle_adjust.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 = [
Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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)",
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down
63 changes: 44 additions & 19 deletions asp_plot/cli/asp_plot.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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(
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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)"
Expand All @@ -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
Expand Down
144 changes: 92 additions & 52 deletions notebooks/bundle_adjust_plots.ipynb
View file

Large diffs are not rendered by default.

18 changes: 13 additions & 5 deletions notebooks/scene_plots.ipynb
View file

Large diffs are not rendered by default.

27 changes: 6 additions & 21 deletions notebooks/stereo_plots.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
]
},
{
Expand All @@ -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)"
]
},
{
Expand Down Expand Up @@ -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",
")"
]
},
Expand Down Expand Up @@ -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",
")"
]
},
Expand Down Expand Up @@ -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",
")"
]
},
Expand Down Expand Up @@ -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": {
Expand Down
Loading

0 comments on commit 4b95766

Please sign in to comment.