refactoring into package

This commit is a squash of about 30 messy commits I made while creating a package structure for asp_plot. In those commits, I significantly restructured the code, added a command line interface, and added tests. I also added a few new notebooks and scripts for plotting and processing stereo data. I also added a few new tests for th

.gitignore

@@ -109,6 +109,7 @@ venv/
 ENV/
 env.bak/
 venv.bak/
+.vscode/
 
 # Spyder project settings
 .spyderproject

README.md

@@ -6,3 +6,34 @@ Our objective is to release a set of standalone Python scripts that can be run a
 
 ## Status
 This is a work in progress, with initial notebooks compiled from recent projects using sample stereo images from the Maxar WorldView, Planet SkySat-C and BlackSky Global constellations. We plan to refine these tools in the coming year, and we welcome community contributions and input. 
+
+## Notes on tests
+
+The ASP uses this tarball for running tests (568 MB size):
+
+```
+wget https://github.com/NeoGeographyToolkit/StereoPipelineTest/releases/download/0.0.1/StereoPipelineTest.tar
+```
+
+[As specified here, it contains all tests scripts, **data**, and expected results.](https://github.com/NeoGeographyToolkit/StereoPipeline/blob/df18be8d32435f7f9db829b9fc951f59bda86e55/.github/workflows/build_test.sh#L200-L206)
+
+The `data/` directory in that tarball looks like:
+
+```
+8.9M - 13APR25_WV02_SEVZ_1030010021A8A100_10030010021A64500_DEM3-3m_10pct.tif
+ 55M - B17_016219_1978_XN_17N282W.8bit.cub
+ 37K - B17_016219_1978_XN_17N282W.8bit.json
+ 55M - B18_016575_1978_XN_17N282W.8bit.cub
+ 41K - B18_016575_1978_XN_17N282W.8bit.json
+3.1M - Copernicus_DSM.tif
+120M - M181058717LE.ce.cub
+171K - M181058717LE.ce.json
+120M - M181073012LE.ce.cub
+176K - M181073012LE.ce.json
+6.3M - Severnaya-Bedrock-UTM47-Ellipsoidal-Height.txt
+5.5M - basic_panchromatic.tif
+ 20M - ref-mars-PC.tif
+147K - ref-mars_yxz.csv
+```
+
+Can we use that data in the tests here as well? Or perhaps better to keep separate, some useful files maybe? (But probably not the `.cub` or `.json` files).

asp_plot/__init__.py

@@ -0,0 +1,4 @@
+import asp_plot.utils
+import asp_plot.processing_parameters
+import asp_plot.scenes
+import asp_plot.bundle_adjust

asp_plot/bundle_adjust.py

@@ -0,0 +1,130 @@
+import os
+import pandas as pd
+import geopandas as gpd
+import matplotlib.pyplot as plt
+import contextily as ctx
+from asp_plot.utils import ColorBar, Plotter
+
+
+class ReadResiduals:
+    def __init__(self, directory, ba_prefix):
+        self.directory = directory
+        self.ba_prefix = ba_prefix
+
+    def get_residual_csv_paths(self):
+        filenames = [
+            f"{self.ba_prefix}-initial_residuals_pointmap.csv",
+            f"{self.ba_prefix}-final_residuals_pointmap.csv",
+        ]
+
+        paths = [
+            os.path.join(os.path.expanduser(self.directory), filename)
+            for filename in filenames
+        ]
+
+        for path in paths:
+            if not os.path.isfile(path):
+                raise ValueError(f"Residuals CSV file not found: {path}")
+
+        init_resid_path, final_resid_path = paths
+
+        return init_resid_path, final_resid_path
+
+    def read_residuals_csv(self, csv_path):
+        resid_cols = [
+            "lon",
+            "lat",
+            "height_above_datum",
+            "mean_residual",
+            "num_observations",
+        ]
+
+        resid_df = pd.read_csv(csv_path, skiprows=2, names=resid_cols)
+
+        # Need the astype('str') to handle cases where column has dtype of int (without the # from DEM appended to some rows)
+        resid_df["from_DEM"] = (
+            resid_df["num_observations"].astype("str").str.contains("# from DEM")
+        )
+
+        resid_df["num_observations"] = (
+            resid_df["num_observations"]
+            .astype("str")
+            .str.split("#", expand=True)[0]
+            .astype(int)
+        )
+
+        resid_gdf = gpd.GeoDataFrame(
+            resid_df,
+            geometry=gpd.points_from_xy(
+                resid_df["lon"], resid_df["lat"], crs="EPSG:4326"
+            ),
+        )
+
+        resid_gdf.filename = os.path.basename(csv_path)
+
+        return resid_gdf
+
+    def get_residual_gdfs(self):
+        init_resid_path, final_resid_path = self.get_residual_csv_paths()
+        init_resid_gdf = self.read_residuals_csv(init_resid_path)
+        final_resid_gdf = self.read_residuals_csv(final_resid_path)
+        return init_resid_gdf, final_resid_gdf
+
+
+class PlotResiduals:
+    def __init__(self, geodataframes):
+        if not isinstance(geodataframes, list):
+            raise ValueError("Input must be a list of GeoDataFrames")
+        self.geodataframes = geodataframes
+
+    def plot_n_residuals(
+        self,
+        column_name="mean_residual",
+        clip_final=True,
+        lognorm=False,
+        clim=None,
+        cmap="inferno",
+        map_crs="EPSG:4326",
+        title_size=10,
+        **ctx_kwargs,
+    ):
+
+        # Get rows and columns and create subplots
+        n = len(self.geodataframes)
+        nrows = (n + 3) // 4
+        ncols = min(n, 4)
+        if n == 1:
+            f, axa = plt.subplots(1, 1, figsize=(8, 6))
+            axa = [axa]
+        else:
+            f, axa = plt.subplots(
+                nrows, ncols, figsize=(4 * ncols, 3 * nrows), sharex=True, sharey=True
+            )
+            axa = axa.flatten()
+
+        # Plot each GeoDataFrame
+        for i, gdf in enumerate(self.geodataframes):
+            if clim is None:
+                clim = ColorBar().get_clim(gdf[column_name])
+
+            gdf = gdf.sort_values(by=column_name).to_crs(map_crs)
+            plotter = Plotter(
+                ax=axa[i],
+                clim=clim,
+                label=column_name,
+            )
+
+            plotter.plot_geodataframe(gdf, column_name, lognorm)
+
+            ctx.add_basemap(ax=axa[i], **ctx_kwargs)
+
+            if clip_final and i == n - 1:
+                axa[i].autoscale(False)
+            axa[i].set_title(f"{column_name:} (n={gdf.shape[0]})", fontsize=title_size)
+
+        # Clean up axes and tighten layout
+        for i in range(n, nrows * ncols):
+            f.delaxes(axa[i])
+        plt.subplots_adjust(wspace=0.2, hspace=0.4)
+        plt.tight_layout()
+

asp_plot/cli/asp_plot.py

@@ -0,0 +1,5 @@
+# TODO: add command line interface; probably break into
+# separate scripts for each plot type
+
+if __name__ == "__main__":
+    main()

asp_plot/processing_parameters.py

@@ -0,0 +1,100 @@
+import matplotlib.pyplot as plt
+
+
+class ProcessingParameters:
+    def __init__(self, bundle_adjust_log, stereo_log, point2dem_log):
+        self.bundle_adjust_log = bundle_adjust_log
+        self.stereo_log = stereo_log
+        self.point2dem_log = point2dem_log
+        self.processing_parameters_dict = {}
+
+    def from_log_files(self):
+        with open(self.bundle_adjust_log, "r") as file:
+            content = file.readlines()
+
+        bundle_adjust_params = ""
+        processing_timestamp = ""
+
+        for line in content:
+            if "bundle_adjust" in line and not bundle_adjust_params:
+                bundle_adjust_params = line.strip()
+
+            if "[ console ]" in line and not processing_timestamp:
+                date, time = line.split()[0], line.split()[1]
+                processing_timestamp = f"{date}-{time[:5].replace(':', '')}"
+
+            if bundle_adjust_params and processing_timestamp:
+                break
+
+        with open(self.stereo_log, "r") as file:
+            content = file.readlines()
+
+        stereo_params = ""
+
+        for line in content:
+            if "stereo_tri" in line and not stereo_params:
+                stereo_params = line.strip()
+
+            if stereo_params:
+                break
+
+        with open(self.point2dem_log, "r") as file:
+            content = file.readlines()
+
+        point2dem_params = ""
+
+        for line in content:
+            if "point2dem" in line and not point2dem_params:
+                point2dem_params = line.strip()
+
+            if point2dem_params:
+                break
+
+        self.processing_parameters_dict = {
+            "bundle_adjust": bundle_adjust_params,
+            "stereo": stereo_params,
+            "point2dem": point2dem_params,
+            "processing_timestamp": processing_timestamp,
+        }
+
+        return self.processing_parameters_dict
+
+    def figure(self):
+        fig, axes = plt.subplots(3, 1, figsize=(10, 6))
+        ax1, ax2, ax3 = axes.flatten()
+
+        ax1.axis("off")
+        ax1.text(
+            0.5,
+            0.5,
+            f"Processed on: {self.processing_parameters_dict['processing_timestamp']:}\n\nBundle Adjust:\n{self.processing_parameters_dict['bundle_adjust']:}",
+            horizontalalignment="center",
+            verticalalignment="center",
+            fontsize=10,
+            wrap=True,
+        )
+
+        ax2.axis("off")
+        ax2.text(
+            0.5,
+            0.5,
+            f"Stereo:\n{self.processing_parameters_dict['stereo']:}",
+            horizontalalignment="center",
+            verticalalignment="center",
+            fontsize=10,
+            wrap=True,
+        )
+
+        ax3.axis("off")
+        ax3.text(
+            0.5,
+            0.5,
+            f"point2dem:\n{self.processing_parameters_dict['point2dem']:}",
+            horizontalalignment="center",
+            verticalalignment="center",
+            fontsize=10,
+            wrap=True,
+        )
+
+        plt.tight_layout()
+        plt.show()

asp_plot/scenes.py

@@ -0,0 +1,30 @@
+import matplotlib.pyplot as plt
+from asp_plot.utils import Raster, Plotter
+
+
+class ScenePlotter:
+    def __init__(self, left_ortho_fn, right_ortho_fn):
+        self.left_ortho_fn = left_ortho_fn
+        self.right_ortho_fn = right_ortho_fn
+
+    def plot_orthos(self, title=None):
+        f, axa = plt.subplots(1, 2, figsize=(10, 5), dpi=300)
+        f.suptitle(title, size=10)
+        axa = axa.ravel()
+
+        if self.left_ortho_fn:
+            ortho_ma = Raster(self.left_ortho_fn).read_array()
+            Plotter(ax=axa[0], cmap="gray", add_cbar=False).plot_array(ortho_ma)
+            axa[0].set_title("Left image")
+        else:
+            axa[0].axis("off")
+
+        if self.right_ortho_fn:
+            ortho_ma = Raster(self.right_ortho_fn).read_array()
+            Plotter(ax=axa[1], cmap="gray", add_cbar=False).plot_array(ortho_ma)
+            axa[1].set_title("Right image")
+        else:
+            axa[1].axis("off")
+
+        f.tight_layout()
+        plt.show()