diff --git a/pyproject.toml b/pyproject.toml
index 07d91db..ef7c5a6 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -94,7 +94,6 @@ ignore = [
     "FBT001",   # boolean-type-hint-positional-argument
     "FBT002",   # boolean-default-value-positional-argument
     "FIX002",   # line-contains-todo
-    "ICN001",   # unconventional-import-alias
     "PLC0415",  # import-outside-top-level
     "PLR0912",  # too-many-branches
     "PLR0913",  # too-many-arguments
diff --git a/src/ivert/cli.py b/src/ivert/cli.py
index b2ac7a1..7778f21 100644
--- a/src/ivert/cli.py
+++ b/src/ivert/cli.py
@@ -1673,7 +1673,7 @@ def database_export(
 
     # --- Read, subset, and merge granules. ---
     import geopandas
-    import pandas
+    import pandas as pd
 
     dt_min = _yyyymmdd_to_delta_time(tmin) if date_filtering else None
     dt_max = _yyyymmdd_to_delta_time(tmax) if date_filtering else None
@@ -1706,7 +1706,7 @@ def database_export(
         raise click.ClickException("No photons found to export after filtering.")
 
     merged = geopandas.GeoDataFrame(
-        pandas.concat(gdfs, ignore_index=True),
+        pd.concat(gdfs, ignore_index=True),
         crs=ev.WGS84_EPSG,
     )
 
diff --git a/src/ivert/convert_error_tif_to_vector.py b/src/ivert/convert_error_tif_to_vector.py
index 4bc1afb..2ee4b5a 100644
--- a/src/ivert/convert_error_tif_to_vector.py
+++ b/src/ivert/convert_error_tif_to_vector.py
@@ -2,7 +2,7 @@
 import sys
 
 import geopandas
-import numpy
+import numpy as np
 import rasterio
 import shapely.geometry
 
@@ -29,8 +29,8 @@ def convert_ivert_error_map_to_vector(
 
         ndv = src.nodatavals[0]
 
-        good_rows, good_cols = numpy.where(
-            (~array.isnan()) if numpy.isnan(ndv) else (array != ndv),
+        good_rows, good_cols = np.where(
+            (~array.isnan()) if np.isnan(ndv) else (array != ndv),
         )
         good_xs, good_ys = rasterio.transform.xy(
             src.transform,
diff --git a/src/ivert/icesat2_database_v2.py b/src/ivert/icesat2_database_v2.py
index 367745a..7cfc188 100644
--- a/src/ivert/icesat2_database_v2.py
+++ b/src/ivert/icesat2_database_v2.py
@@ -20,8 +20,8 @@
 import fetchez.core
 import fetchez.spatial
 import globato
-import numpy
-import pandas
+import numpy as np
+import pandas as pd
 import xarray
 from fetchez.modules.earthdata import IceSat2 as _FetchezIceSat2
 
@@ -193,7 +193,7 @@ def create_new_database(
         self,
         populate: bool = True,
         overwrite: bool = False,
-    ) -> pandas.DataFrame:
+    ) -> pd.DataFrame:
         """Create a new database from scratch.
 
         Parameters
@@ -232,15 +232,15 @@ def create_new_database(
                     records.append(meta)
 
             if records:
-                gdf = pandas.DataFrame(records)[list(self._empty_db_dict().keys())]
+                gdf = pd.DataFrame(records)[list(self._empty_db_dict().keys())]
             else:
-                gdf = pandas.DataFrame(self._empty_db_dict()).drop(
+                gdf = pd.DataFrame(self._empty_db_dict()).drop(
                     labels=0,
                     axis="rows",
                 )
 
         else:
-            gdf = pandas.DataFrame(self._empty_db_dict()).drop(labels=0, axis="rows")
+            gdf = pd.DataFrame(self._empty_db_dict()).drop(labels=0, axis="rows")
 
         self._write_index(gdf)
         if os.path.exists(self.db_fname):
@@ -443,7 +443,7 @@ def _source_granule_from_filename(cls, filename: str) -> str:
         return cls._NC_SUFFIX_RE.sub("", filename)
 
     @staticmethod
-    def _h5_along_track_m(h5_fn: str, beams) -> pandas.DataFrame:
+    def _h5_along_track_m(h5_fn: str, beams) -> pd.DataFrame:
         """Return per-photon cumulative along-track distance (m) for the given beams.
 
         Cumulative distance = sum of segment_length values up to each photon's segment
@@ -468,20 +468,20 @@ def _h5_along_track_m(h5_fn: str, beams) -> pandas.DataFrame:
                     continue
 
                 # Cumulative distance at the start of each segment
-                seg_cumul_start = numpy.concatenate(
-                    [[0.0], numpy.cumsum(seg_length[:-1])],
+                seg_cumul_start = np.concatenate(
+                    [[0.0], np.cumsum(seg_length[:-1])],
                 )
 
                 # Map each photon to its segment via ph_index_beg
                 n = len(delta_time)
-                seg_of_ph = numpy.clip(
-                    numpy.searchsorted(ph_index_beg, numpy.arange(n), side="right") - 1,
+                seg_of_ph = np.clip(
+                    np.searchsorted(ph_index_beg, np.arange(n), side="right") - 1,
                     0,
                     len(ph_index_beg) - 1,
                 )
 
                 dfs.append(
-                    pandas.DataFrame(
+                    pd.DataFrame(
                         {
                             "laser": beam,
                             "delta_time": delta_time,
@@ -493,9 +493,9 @@ def _h5_along_track_m(h5_fn: str, beams) -> pandas.DataFrame:
                 )
 
         return (
-            pandas.concat(dfs, ignore_index=True)
+            pd.concat(dfs, ignore_index=True)
             if dfs
-            else pandas.DataFrame(
+            else pd.DataFrame(
                 columns=["laser", "delta_time", "x", "y", "along_track_m"],
             )
         )
@@ -602,12 +602,12 @@ def _process_h5_to_nc(
             use_external_masks=use_external_masks,
         )
 
-        chunks = [pandas.DataFrame(chunk) for chunk in stream]
+        chunks = [pd.DataFrame(chunk) for chunk in stream]
 
         if not chunks:
             return None
 
-        df = pandas.concat(chunks, ignore_index=True)
+        df = pd.concat(chunks, ignore_index=True)
         df = df.rename(columns={"ph_h_classed": "class_code"})
 
         # Temporal filter
@@ -657,16 +657,16 @@ def _process_h5_to_nc(
             "data_bbox": [xmin, xmax, ymin, ymax, tmin, tmax],
             "zbounds": [zmin, zmax],
             "numphotons": len(df),
-            "numphotons_unclassified": int(numpy.count_nonzero(cc == -1)),
-            "numphotons_noise": int(numpy.count_nonzero(cc == 0)),
-            "numphotons_ground": int(numpy.count_nonzero(cc == 1)),
-            "numphotons_canopy": int(numpy.count_nonzero(cc == 2)),
-            "numphotons_canopy_top": int(numpy.count_nonzero(cc == 3)),
-            "numphotons_ice_surface": int(numpy.count_nonzero(cc == 6)),
-            "numphotons_buildings": int(numpy.count_nonzero(cc == 7)),
-            "numphotons_bathy_floor": int(numpy.count_nonzero(cc == 40)),
-            "numphotons_bathy_surface": int(numpy.count_nonzero(cc == 41)),
-            "numphotons_inland_water_surface": int(numpy.count_nonzero(cc == 42)),
+            "numphotons_unclassified": int(np.count_nonzero(cc == -1)),
+            "numphotons_noise": int(np.count_nonzero(cc == 0)),
+            "numphotons_ground": int(np.count_nonzero(cc == 1)),
+            "numphotons_canopy": int(np.count_nonzero(cc == 2)),
+            "numphotons_canopy_top": int(np.count_nonzero(cc == 3)),
+            "numphotons_ice_surface": int(np.count_nonzero(cc == 6)),
+            "numphotons_buildings": int(np.count_nonzero(cc == 7)),
+            "numphotons_bathy_floor": int(np.count_nonzero(cc == 40)),
+            "numphotons_bathy_surface": int(np.count_nonzero(cc == 41)),
+            "numphotons_inland_water_surface": int(np.count_nonzero(cc == 42)),
             "downloaded_on_utc": int(
                 datetime.datetime.now(datetime.UTC).strftime("%Y%m%d"),
             ),
@@ -725,30 +725,30 @@ def _write_index(self, df) -> None:
 
         for col in self._INDEX_STR_COLS:
             values = (
-                numpy.array(
+                np.array(
                     ["" if v is None else str(v) for v in df[col]],
                     dtype=object,
                 )
                 if n
-                else numpy.array([], dtype=object)
+                else np.array([], dtype=object)
             )
             data_vars[col] = ("record", values)
 
         int_cols = (*self._INDEX_INT_COLS, *self._INDEX_BBOX_INT_COLS)
         for col in int_cols:
             values = (
-                numpy.asarray(df[col].to_list(), dtype="int64")
+                np.asarray(df[col].to_list(), dtype="int64")
                 if n
-                else numpy.array([], dtype="int64")
+                else np.array([], dtype="int64")
             )
             data_vars[col] = ("record", values)
 
         float_cols = (*self._INDEX_BBOX_FLOAT_COLS, *self._INDEX_ZBOUNDS_COLS)
         for col in float_cols:
             values = (
-                numpy.asarray(df[col].to_list(), dtype="float64")
+                np.asarray(df[col].to_list(), dtype="float64")
                 if n
-                else numpy.array([], dtype="float64")
+                else np.array([], dtype="float64")
             )
             data_vars[col] = ("record", values)
 
@@ -770,7 +770,7 @@ def _write_index(self, df) -> None:
         ds.to_netcdf(self.db_fname, encoding=encoding)
 
     @classmethod
-    def read_index_file(cls, index_fname: str) -> pandas.DataFrame:
+    def read_index_file(cls, index_fname: str) -> pd.DataFrame:
         """Read any IVERT NetCDF index file into a plain pandas DataFrame.
 
         Every column comes back as a numpy array with no per-row Python loops and
@@ -793,7 +793,7 @@ def read_index_file(cls, index_fname: str) -> pandas.DataFrame:
         ):
             data[col] = ds[cls._stored_col_name(col, ds)].to_numpy()
 
-        df = pandas.DataFrame(data)
+        df = pd.DataFrame(data)
         # Restore the canonical column order used elsewhere in the codebase.
         return df[list(cls._empty_db_dict().keys())]
 
@@ -886,7 +886,7 @@ def read_database_file(
     def omit_photons_from_exclusion_bbox(
         dataframe,
         bbox_to_exclude,
-    ) -> pandas.DataFrame:
+    ) -> pd.DataFrame:
         """Exclude any photons that fall within the particular bounding box."""
         x = dataframe["x"]
         y = dataframe["y"]
@@ -924,7 +924,7 @@ def read_granule(
         granule_fn: str,
         subset_bbox: list | tuple | None = None,
         photon_classes: list | tuple | None = None,
-    ) -> pandas.DataFrame:
+    ) -> pd.DataFrame:
         """Read classified photons from a NetCDF granule file.
 
         Parameters
@@ -987,7 +987,7 @@ def query_photons(
         min_confidence_level: int = 1,
         omit_bboxes=None,
         # download_new_data: bool = False,
-    ) -> pandas.DataFrame | None:
+    ) -> pd.DataFrame | None:
         """Query the database for photons in a given bounding box and date range.
 
         Parameters
@@ -1025,7 +1025,7 @@ def query_photons(
         )
         logger.info(
             "%d granules exist with %s ground photons and %s bathy_floor photons.",
-            numpy.count_nonzero(fnames.apply(os.path.exists)),
+            np.count_nonzero(fnames.apply(os.path.exists)),
             f"{gdf_subset['numphotons_ground'].sum():,}",
             f"{gdf_subset['numphotons_bathy_floor'].sum():,}",
         )
@@ -1046,7 +1046,7 @@ def query_photons(
         if len(granule_dfs) == 0:
             return None
 
-        photons_df = pandas.concat(granule_dfs, ignore_index=True)
+        photons_df = pd.concat(granule_dfs, ignore_index=True)
 
         if min_bathy_confidence > 0.0:
             photons_df = photons_df[
@@ -1064,7 +1064,7 @@ def query_photons(
             omit_bboxes = []
 
         # If we're given a single bounding box of exclusions as a 4- or 6-tuple of numbers (not iterables), put it in a 1-length list.
-        if len(omit_bboxes) in (4, 6) and not numpy.any(
+        if len(omit_bboxes) in (4, 6) and not np.any(
             [self.is_iterable(num) for num in omit_bboxes],
         ):
             omit_bboxes = [omit_bboxes]
@@ -1078,8 +1078,8 @@ def query_photons(
                 "Trimmed granules from %s to %s photons (%s ground, %s bathy).",
                 f"{gdf_subset['numphotons'].sum():,}",
                 f"{len(photons_df):,}",
-                f"{numpy.count_nonzero(photons_df['class_code'] == 1):,}",
-                f"{numpy.count_nonzero(photons_df['class_code'] == 40):,}",
+                f"{np.count_nonzero(photons_df['class_code'] == 1):,}",
+                f"{np.count_nonzero(photons_df['class_code'] == 40):,}",
             )
         else:
             logger.info("No photons in bbox.")
@@ -1126,7 +1126,7 @@ def convert_date_to_yyyymmdd(
                 "Date must be an int, str, datetime.datetime, or datetime.date.",
             )
 
-    def query_granules(self, bbox: list | tuple) -> pandas.DataFrame | None:
+    def query_granules(self, bbox: list | tuple) -> pd.DataFrame | None:
         """Return a sub-dataframe of granules in the database that possibly intersect the bounding box, using data bounding boxes."""
         gdf = self.open_gdf()
         if gdf is None or len(gdf) == 0:
@@ -1462,7 +1462,7 @@ def download_new_granules(
             if not new_records:
                 continue
 
-            new_gdf = pandas.DataFrame(new_records)[list(self._empty_db_dict().keys())]
+            new_gdf = pd.DataFrame(new_records)[list(self._empty_db_dict().keys())]
             if existing_gdf is None or len(existing_gdf) == 0:
                 self.gdf = new_gdf
             else:
@@ -1490,7 +1490,7 @@ def download_new_granules(
                     # granule but a non-overlapping query bbox are legitimately distinct data
                     # (e.g. a different date range or region) and must not be dropped, since
                     # doing so would discard original data rather than de-duplicate it.
-                    bbox_overlaps = pandas.Series(
+                    bbox_overlaps = pd.Series(
                         ivert.utils.cuboid_funcs.cuboids_intersect_vectorized(
                             existing_gdf["query_bbox_xmin"].to_numpy(),
                             existing_gdf["query_bbox_xmax"].to_numpy(),
@@ -1515,7 +1515,7 @@ def download_new_granules(
                             if os.path.exists(old_fpath):
                                 os.remove(old_fpath)
                         existing_gdf = existing_gdf[~is_replaced]
-                self.gdf = pandas.concat(
+                self.gdf = pd.concat(
                     [existing_gdf, new_gdf],
                     ignore_index=True,
                 )
@@ -1586,7 +1586,7 @@ def bounds(self, axis: str, data_or_query: str = "data") -> tuple | None:
 
     def unique_bboxes(
         self,
-        gdf: pandas.DataFrame | None = None,
+        gdf: pd.DataFrame | None = None,
         data_or_query: str = "query",
     ) -> list | None:
         """Return a numpy array of unique query bounding boxes in the database.
@@ -1796,22 +1796,22 @@ def split_bbox_into_parts(
     xmin, xmax, ymin, ymax = bbox
     max_deg_size = tile_size_deg * max_tile_scale_factor
 
-    xbins = numpy.arange(xmin, xmax, tile_size_deg)
-    ybins = numpy.arange(ymin, ymax, tile_size_deg)
+    xbins = np.arange(xmin, xmax, tile_size_deg)
+    ybins = np.arange(ymin, ymax, tile_size_deg)
 
     if xbins[-1] < xmax:
         if len(xbins) == 1 or ((xmax - xbins[-2]) > max_deg_size):
-            xbins = numpy.append(xbins, xmax)
+            xbins = np.append(xbins, xmax)
         else:
             xbins[-1] = xmax
 
     if ybins[-1] < ymax:
         if len(ybins) == 1 or ((ymax - ybins[-2]) > max_deg_size):
-            ybins = numpy.append(ybins, ymax)
+            ybins = np.append(ybins, ymax)
         else:
             ybins[-1] = ymax
 
-    binxs, binys = numpy.meshgrid(xbins, ybins)
+    binxs, binys = np.meshgrid(xbins, ybins)
     bin_xmins = binxs[:-1, :-1].flatten()
     bin_xmaxs = binxs[1:, 1:].flatten()
     bin_ymins = binys[:-1, :-1].flatten()
diff --git a/src/ivert/icesat2_requests.py b/src/ivert/icesat2_requests.py
index 1d5d8f1..7547404 100644
--- a/src/ivert/icesat2_requests.py
+++ b/src/ivert/icesat2_requests.py
@@ -11,8 +11,8 @@
 import time
 
 import dateparser
-import numpy
-import pandas
+import numpy as np
+import pandas as pd
 
 import ivert.utils.configfile
 
@@ -49,7 +49,7 @@ def find_matching_request(
         only_unexpired: bool = True,
         tolerance: float = 1e-9,
         return_rows: bool = False,
-    ) -> dict | pandas.DataFrame | None:
+    ) -> dict | pd.DataFrame | None:
         """Return the cached Harmony JSON for a matching request, or None.
 
         Parameters
@@ -79,7 +79,7 @@ def find_matching_request(
                 self._is_expired,
             )
 
-        if numpy.any(matching_mask):
+        if np.any(matching_mask):
             if return_rows:
                 return self.df[matching_mask]
             return self._read_json(self.df[matching_mask].iloc[0]["json"])
@@ -104,7 +104,7 @@ def add_record(
         if isinstance(json_dict, str):
             json_dict = self._read_json(json_dict)
 
-        new_row = pandas.DataFrame(
+        new_row = pd.DataFrame(
             [
                 {
                     "atl_dataset": atl_dataset,
@@ -117,7 +117,7 @@ def add_record(
             ],
         )
 
-        self.df = pandas.concat([self.df, new_row], ignore_index=True)
+        self.df = pd.concat([self.df, new_row], ignore_index=True)
 
         if write_file:
             self.export()
@@ -172,9 +172,9 @@ def open(self, read_again: bool = False, create_if_nonexistent: bool = True):
             num_tries = 0
             while num_tries < 20:
                 try:
-                    self.df = pandas.read_csv(self.csv_file, index_col=False)
+                    self.df = pd.read_csv(self.csv_file, index_col=False)
                     break
-                except (TypeError, pandas.errors.ParserError):
+                except (TypeError, pd.errors.ParserError):
                     num_tries += 1
                     if num_tries >= 20:
                         raise
@@ -198,7 +198,7 @@ def clean_csv(self, verbose: bool = True):
             self.open()
 
         expired = self.df["expiration_date"].apply(self._is_expired)
-        if numpy.any(expired):
+        if np.any(expired):
             if verbose:
                 print(f"Removing {expired.sum()} expired Harmony request record(s).")
             self.df = self.df[~expired]
@@ -211,7 +211,7 @@ def clean_csv(self, verbose: bool = True):
 
     def _create_empty(self):
         """Create an empty CSV with the correct columns."""
-        self.df = pandas.DataFrame(
+        self.df = pd.DataFrame(
             columns=[
                 "atl_dataset",
                 "bbox",
diff --git a/src/ivert/plot_photon_clouds_v2.py b/src/ivert/plot_photon_clouds_v2.py
index 42695e9..b2b12d3 100644
--- a/src/ivert/plot_photon_clouds_v2.py
+++ b/src/ivert/plot_photon_clouds_v2.py
@@ -18,11 +18,11 @@
 import sys
 
 import click
-import matplotlib
+import matplotlib as mpl
 import numpy as np
 import pandas as pd
 
-matplotlib.use("Agg")
+mpl.use("Agg")
 import matplotlib.pyplot as plt
 import netCDF4
 
diff --git a/src/ivert/plot_results_slope_centrality.py b/src/ivert/plot_results_slope_centrality.py
index 8f154a5..ea88f24 100644
--- a/src/ivert/plot_results_slope_centrality.py
+++ b/src/ivert/plot_results_slope_centrality.py
@@ -6,8 +6,8 @@
 # Include the base /src/ directory of thie project, to add all the other modules.
 import import_parent_dir
 import matplotlib.pyplot as plt
-import numpy
-import pandas
+import numpy as np
+import pandas as pd
 import rasterio
 from matplotlib import ticker
 
@@ -21,8 +21,8 @@
 
 def add_lat_lons(df):
     """From the filename, get the lat/lon of each grid-cell point, from the filename and the i,j position of the pixel."""
-    lat = numpy.empty((len(df),), dtype=float)
-    lon = numpy.empty((len(df),), dtype=float)
+    lat = np.empty((len(df),), dtype=float)
+    lon = np.empty((len(df),), dtype=float)
 
     for idx, ((i, j), row) in enumerate(df.iterrows()):
         fname = row["filename"]
@@ -50,7 +50,7 @@ def get_slopes(df, files_dirname):
     }
 
     fn_array_dict = {}
-    slopes = numpy.empty((len(df),), dtype=float)
+    slopes = np.empty((len(df),), dtype=float)
     for idx, ((i, j), row) in enumerate(df.iterrows()):
         fname = row.filename
         if fname in fn_array_dict:
@@ -85,7 +85,7 @@ def plot_errors_against_slope_centrality(
         "total_results.h5",
     )
     if os.path.exists(total_results_h5):
-        data = pandas.read_hdf(total_results_h5)
+        data = pd.read_hdf(total_results_h5)
         if verbose:
             print(os.path.basename(total_results_h5), "read.")
     else:
@@ -131,15 +131,15 @@ def plot_errors_against_slope_centrality(
     # slope           = data["slope"]
 
     fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, dpi=dpi, figsize=(8, 5))
-    low, hi = numpy.percentile(meandiff, [2.5, 97.5])
+    low, hi = np.percentile(meandiff, [2.5, 97.5])
     ax1.hist(meandiff, bins=100, range=(low, hi))
     ax1.set_ylabel("% of data cells")
     ax1.set_xlabel("Elevation difference (m)")
     ax1.yaxis.set_major_formatter(ticker.PercentFormatter(len(meandiff), decimals=0))
 
     # Add the lines for mean +- std
-    center = numpy.mean(meandiff)
-    std = numpy.std(meandiff)
+    center = np.mean(meandiff)
+    std = np.std(meandiff)
     ax1.axvline(x=center, color="darkgreen", linewidth=0.75)
     ax1.axvline(x=center + std, color="darkgreen", linestyle="--", linewidth=0.5)
     ax1.axvline(x=center - std, color="darkgreen", linestyle="--", linewidth=0.5)
@@ -171,15 +171,15 @@ def plot_errors_against_slope_centrality(
     ax2.set_ylabel("Error (m)")
 
     # There definitely seems to be a relationship between errors and % coverage of each pixel. Plot it over this.
-    unique_coverages = numpy.unique(coverage_pct)
-    coverage_lo = numpy.zeros((len(unique_coverages),))
-    coverage_hi = numpy.zeros((len(unique_coverages),))
-    coverage_rmse = numpy.zeros(
+    unique_coverages = np.unique(coverage_pct)
+    coverage_lo = np.zeros((len(unique_coverages),))
+    coverage_hi = np.zeros((len(unique_coverages),))
+    coverage_rmse = np.zeros(
         len(
             unique_coverages,
         ),
     )
-    coverage_counts = numpy.zeros(
+    coverage_counts = np.zeros(
         len(
             unique_coverages,
         ),
@@ -187,8 +187,8 @@ def plot_errors_against_slope_centrality(
     for i, c in enumerate(unique_coverages):
         cmask = coverage_pct == c
         diff_subset = meandiff[cmask]
-        coverage_lo[i], coverage_hi[i] = numpy.percentile(diff_subset, (2, 98))
-        coverage_rmse[i] = numpy.sqrt(numpy.mean(diff_subset**2))
+        coverage_lo[i], coverage_hi[i] = np.percentile(diff_subset, (2, 98))
+        coverage_rmse[i] = np.sqrt(np.mean(diff_subset**2))
         coverage_counts[i] = diff_subset.size
 
     ax2.fill_between(
@@ -221,7 +221,7 @@ def plot_errors_against_slope_centrality(
     meandiff_good = meandiff[good_coverage_mask]
     # coverage_good = coverage_pct[good_coverage_mask]
 
-    low, hi = numpy.percentile(meandiff, [5, 90])
+    low, hi = np.percentile(meandiff, [5, 90])
     ax4.hist(meandiff_good, bins=100, range=(low, hi), color="darkred", alpha=0.6)
     ax4.set_ylabel("% of data cells")
     ax4.set_xlabel("Elevation difference (m)")
@@ -230,8 +230,8 @@ def plot_errors_against_slope_centrality(
     )
 
     # Add the lines for mean +- std
-    center = numpy.mean(meandiff_good)
-    std = numpy.std(meandiff_good)
+    center = np.mean(meandiff_good)
+    std = np.std(meandiff_good)
     ax4.axvline(x=center, color="darkred", linewidth=0.75)
     ax4.axvline(x=center + std, color="darkred", linestyle="--", linewidth=0.5)
     ax4.axvline(x=center - std, color="darkred", linestyle="--", linewidth=0.5)
diff --git a/src/ivert/plot_validation_results.py b/src/ivert/plot_validation_results.py
index b64d762..62ba950 100644
--- a/src/ivert/plot_validation_results.py
+++ b/src/ivert/plot_validation_results.py
@@ -2,10 +2,10 @@
 import math
 import os
 
-import matplotlib
+import matplotlib as mpl
 import matplotlib.pyplot as plt
-import numpy
-import pandas
+import numpy as np
+import pandas as pd
 import six
 import tqdm
 from matplotlib import ticker
@@ -29,10 +29,10 @@ def get_data_from_h5_or_list(
     empty_val: float = ivert_config.dem_default_ndv,
     include_filenames: bool = False,
     verbose: bool = True,
-) -> pandas.DataFrame:
+) -> pd.DataFrame:
     """Return the data either from a single hdf5 results file, or a list of them. Filter out empty (bad data) values."""
     if type(h5_name_or_list) is str:
-        data = pandas.read_hdf(h5_name_or_list)
+        data = pd.read_hdf(h5_name_or_list)
         if include_filenames:
             if orig_filenames is None:
                 data["filename"] = os.path.basename(h5_name_or_list)
@@ -55,7 +55,7 @@ def get_data_from_h5_or_list(
             ),
         ):
             if os.path.exists(h5_file):
-                temp_data = pandas.read_hdf(h5_file)
+                temp_data = pd.read_hdf(h5_file)
                 if include_filenames:
                     if orig_filenames is None:
                         temp_data["filename"] = os.path.basename(h5_file)
@@ -65,7 +65,7 @@ def get_data_from_h5_or_list(
 
                 data_list.append(temp_data)
 
-        data = pandas.concat(data_list)
+        data = pd.concat(data_list)
     else:
         raise TypeError(
             "Non-iterable value for parameter 'results_h5_name_or_list':",
@@ -115,9 +115,9 @@ def plot_histograms_and_line(
     # If we're writing a PNG file, use the "Agg" backend (no display).
     # This helps avoid errors.
     if os.path.splitext(output_figure_name)[1].lower() == ".png":
-        matplotlib.use("Agg")
+        mpl.use("Agg")
 
-    if type(results_h5_or_list_or_df) is pandas.DataFrame:
+    if type(results_h5_or_list_or_df) is pd.DataFrame:
         data = results_h5_or_list_or_df
     else:
         data = get_data_from_h5_or_list(results_h5_or_list_or_df, empty_val=empty_val)
@@ -147,7 +147,7 @@ def plot_histograms_and_line(
 
     # Generate figure. Scale width proportionally to panel count.
     if figsize is None:
-        figsize = matplotlib.rcParams["figure.figsize"]
+        figsize = mpl.rcParams["figure.figsize"]
     scaled_figsize = (figsize[0] * ncols / 3, figsize[1])
     fig, axes = plt.subplots(
         1,
@@ -186,15 +186,15 @@ def plot_histograms_and_line(
         )
 
         # Add the lines for mean +- std
-        center = numpy.mean(meandiff_land)
-        std = numpy.std(meandiff_land)
+        center = np.mean(meandiff_land)
+        std = np.std(meandiff_land)
         ax1.axvline(x=center, color="black", linewidth=0.75)
         ax1.axvline(x=center + std, color="black", linestyle="--", linewidth=0.5)
         ax1.axvline(x=center - std, color="black", linestyle="--", linewidth=0.5)
 
         # Crop the left & right (only if greater than 20 points)
         if len(meandiff_land) >= 20:
-            cutoffs = numpy.percentile(meandiff_land, [1, 99])
+            cutoffs = np.percentile(meandiff_land, [1, 99])
         else:
             cutoffs = [min(meandiff_land), max(meandiff_land)]
 
@@ -211,7 +211,7 @@ def plot_histograms_and_line(
             cutoffs[0] = center - (std * hist_cutoff_num_stddevs)
 
         # Just error checking, if any of the cutoffs come back with NaN or Inf, just clip it to -1, 1, debug later.
-        if numpy.any(numpy.isnan(cutoffs) | numpy.isinf(cutoffs)):
+        if np.any(np.isnan(cutoffs) | np.isinf(cutoffs)):
             cutoffs = [-1, 1]
 
         ax1.set_xlim(cutoffs)
@@ -240,7 +240,7 @@ def plot_histograms_and_line(
 
         # If requested, add the RMSE value to the figure.
         if also_add_rmse_to_hist:
-            rmse = numpy.sqrt(numpy.mean(meandiff_land**2))
+            rmse = np.sqrt(np.mean(meandiff_land**2))
             txt_std = ax1.text(
                 0.97,
                 0.95,
@@ -287,15 +287,15 @@ def plot_histograms_and_line(
         )
 
         # Add the lines for mean +- std
-        center = numpy.mean(meandiff_bathy)
-        std = numpy.std(meandiff_bathy)
+        center = np.mean(meandiff_bathy)
+        std = np.std(meandiff_bathy)
         ax2.axvline(x=center, color="black", linewidth=0.75)
         ax2.axvline(x=center + std, color="black", linestyle="--", linewidth=0.5)
         ax2.axvline(x=center - std, color="black", linestyle="--", linewidth=0.5)
 
         # Crop the left & right (only if greater than 20 points)
         if len(meandiff_bathy) >= 20:
-            cutoffs = numpy.percentile(meandiff_bathy, [1, 99])
+            cutoffs = np.percentile(meandiff_bathy, [1, 99])
         else:
             cutoffs = [min(meandiff_bathy), max(meandiff_bathy)]
 
@@ -312,7 +312,7 @@ def plot_histograms_and_line(
             cutoffs[0] = center - (std * hist_cutoff_num_stddevs)
 
         # Just error checking, if any of the cutoffs come back with NaN or Inf, just clip it to -1, 1, debug later.
-        if numpy.any(numpy.isnan(cutoffs) | numpy.isinf(cutoffs)):
+        if np.any(np.isnan(cutoffs) | np.isinf(cutoffs)):
             cutoffs = [-1, 1]
 
         ax2.set_xlim(cutoffs)
@@ -341,7 +341,7 @@ def plot_histograms_and_line(
 
         # If requested, add the RMSE value to the figure.
         if also_add_rmse_to_hist:
-            rmse = numpy.sqrt(numpy.mean(meandiff_bathy**2))
+            rmse = np.sqrt(np.mean(meandiff_bathy**2))
             txt_std = ax2.text(
                 0.97,
                 0.95,
@@ -427,7 +427,7 @@ def plot_histograms_and_line(
     if place_name is None:
         place_name = "DEM"
 
-    rmse = (numpy.sum(meandiff**2) / len(meandiff)) ** 0.5
+    rmse = (np.sum(meandiff**2) / len(meandiff)) ** 0.5
 
     fig.suptitle(
         f"{place_name}: Errors and Distributions\nRMSE = {rmse:0.3f} m,   N = {len(meandiff):,} cells",
@@ -476,9 +476,9 @@ def plot_histogram_and_error_stats_4_panels(
     # If we're writing a PNG file, use the "Agg" backend (no display).
     # This helps avoid errors.
     if os.path.splitext(output_figure_name)[1].lower() == ".png":
-        matplotlib.use("Agg")
+        mpl.use("Agg")
 
-    if type(results_h5_or_list_or_df) is pandas.DataFrame:
+    if type(results_h5_or_list_or_df) is pd.DataFrame:
         data = results_h5_or_list_or_df
     else:
         data = get_data_from_h5_or_list(results_h5_or_list_or_df, empty_val=empty_val)
@@ -517,7 +517,7 @@ def plot_histogram_and_error_stats_4_panels(
 
     # Generate figure. If a figure size isn't given, use the matplotlib.rcParams default.
     if figsize is None:
-        figsize = matplotlib.rcParams["figure.figsize"]
+        figsize = mpl.rcParams["figure.figsize"]
     fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(
         2,
         2,
@@ -545,14 +545,14 @@ def plot_histogram_and_error_stats_4_panels(
     )
 
     # Add the lines for mean +- std
-    center = numpy.mean(meandiff)
-    std = numpy.std(meandiff)
+    center = np.mean(meandiff)
+    std = np.std(meandiff)
     ax1.axvline(x=center, color="darkred", linewidth=0.75)
     ax1.axvline(x=center + std, color="darkred", linestyle="--", linewidth=0.5)
     ax1.axvline(x=center - std, color="darkred", linestyle="--", linewidth=0.5)
 
     # Crop the left & right
-    cutoffs = numpy.percentile(meandiff, [1, 99])
+    cutoffs = np.percentile(meandiff, [1, 99])
     # cutoffs = [-max(numpy.abs(cutoffs)), max(numpy.abs(cutoffs))]
     # Do not crop the photo to make the stddev lines fall outside the plot.
     # If they do, reset the min/max cutoff to be 2 stddev away from the mean on that side.
@@ -584,7 +584,7 @@ def plot_histogram_and_error_stats_4_panels(
 
     # If requested, as the RMSE value to the figure.
     if also_add_rmse_to_hist:
-        rmse = numpy.sqrt(numpy.mean(meandiff**2))
+        rmse = np.sqrt(np.mean(meandiff**2))
         txt_std = ax1.text(
             0.97,  # 0.12 if text_left else 0.97,
             0.95,
@@ -672,17 +672,17 @@ def plot_histogram_and_error_stats_4_panels(
     )
 
     # Crop the left & right, right at 98 percentile.
-    cutoff = numpy.percentile(numphotons_intd, 99)
+    cutoff = np.percentile(numphotons_intd, 99)
     xmin = ax3.get_xlim()[0]
     ax3.set_xlim(xmin * 0.5, cutoff)
 
-    center = numpy.mean(numphotons_intd)
-    std = numpy.std(numphotons_intd)
+    center = np.mean(numphotons_intd)
+    std = np.std(numphotons_intd)
 
     ax3.text(
         0.95,
         0.95,
-        f"{int(numpy.round(center)):d} $\\pm$ {int(numpy.round(std)):d}\nphotons per cell",
+        f"{int(np.round(center)):d} $\\pm$ {int(np.round(std)):d}\nphotons per cell",
         ha="right",
         va="top",
         fontsize="small",
@@ -713,7 +713,7 @@ def plot_histogram_and_error_stats_4_panels(
     )
 
     # Crop the right edge at the 99th percentile
-    cutoff = numpy.percentile(canopy_fraction, 99)
+    cutoff = np.percentile(canopy_fraction, 99)
     # Add just a bit of padding on the left to make more room for the "D" label.
     xmin = ax4.get_xlim()[0] - 0.025 * (cutoff - ax4.get_xlim()[0])
     ax4.set_xlim(xmin, cutoff)
@@ -723,14 +723,14 @@ def plot_histogram_and_error_stats_4_panels(
     # median = numpy.median(canopy_fraction)
     canopy_mask = (
         (canopy_fraction > 0.0)
-        & numpy.isfinite(canopy_fraction)
-        & ~numpy.isnan(canopy_fraction)
+        & np.isfinite(canopy_fraction)
+        & ~np.isnan(canopy_fraction)
     )
 
     ax4.text(
         0.95,
         0.95,
-        f"{numpy.count_nonzero(canopy_mask) * 100 / canopy_fraction.size:0.1f} % of cells have >0 cover:\n{numpy.mean(canopy_fraction[canopy_mask]):0.1f} $\\pm$ {numpy.std(canopy_fraction[canopy_mask]):0.1f} % canopy cover\nin non-zero cells",
+        f"{np.count_nonzero(canopy_mask) * 100 / canopy_fraction.size:0.1f} % of cells have >0 cover:\n{np.mean(canopy_fraction[canopy_mask]):0.1f} $\\pm$ {np.std(canopy_fraction[canopy_mask]):0.1f} % canopy cover\nin non-zero cells",
         ha="right",
         va="top",
         transform=ax4.transAxes,
@@ -763,7 +763,7 @@ def plot_histogram_and_error_stats_4_panels(
         print(output_figure_name, "written.")
 
         # Compute the RMSE and spit that out too.
-        rmse = (numpy.sum(meandiff**2) / len(meandiff)) ** 0.5
+        rmse = (np.sum(meandiff**2) / len(meandiff)) ** 0.5
         print(f"\tRMSE: {rmse:0.3f} m")
 
     # Clear the figure and close the plot.
@@ -862,7 +862,7 @@ def plot_histogram_and_error_stats_4_panels(
 
 
 if __name__ == "__main__":
-    results_df = pandas.read_hdf(
+    results_df = pd.read_hdf(
         "/home/mmacferrin/Research/DEMs/CUDEMs_1_9_Oregon_2025/dems/2025.10.10_w_metadata/icesat2/ncei19_n45x50_w124x00_2025v1_results.h5",
     )
     plot_histograms_and_line(
diff --git a/src/ivert/transform_points.py b/src/ivert/transform_points.py
index 75a5e08..8acafb0 100644
--- a/src/ivert/transform_points.py
+++ b/src/ivert/transform_points.py
@@ -2,7 +2,7 @@
 import math
 import os
 
-import numpy
+import numpy as np
 import pyproj
 import rasterio
 import transformez
@@ -19,12 +19,12 @@
 
 
 def transform_points(
-    x: list | tuple | numpy.ndarray,
-    y: list | tuple | numpy.ndarray,
-    z: list | tuple | numpy.ndarray,
+    x: list | tuple | np.ndarray,
+    y: list | tuple | np.ndarray,
+    z: list | tuple | np.ndarray,
     src_epsg: str | int,
     dst_epsg: str | int,
-    src_region: list | tuple | numpy.ndarray | None = None,
+    src_region: list | tuple | np.ndarray | None = None,
     cache_dir: str | None = None,
 ) -> tuple:
     """Transform a set of 3D points from one coordinate reference system to another.
@@ -62,9 +62,9 @@ def transform_points(
     if src_crs.is_exact_same(dst_crs):
         return x, y, z
 
-    x = numpy.asarray(x, dtype=float)
-    y = numpy.asarray(y, dtype=float)
-    z = numpy.asarray(z, dtype=float)
+    x = np.asarray(x, dtype=float)
+    y = np.asarray(y, dtype=float)
+    z = np.asarray(z, dtype=float)
 
     src_horz, src_vert_epsg = _decompose_crs(src_crs)
     dst_horz, dst_vert_epsg = _decompose_crs(dst_crs)
@@ -165,14 +165,14 @@ def _grid_covers_region(grid_fn: str, region_bounds: list[float], tol=1e-9) -> b
 
 
 def _apply_vertical_transform(
-    x: numpy.ndarray,
-    y: numpy.ndarray,
-    z: numpy.ndarray,
+    x: np.ndarray,
+    y: np.ndarray,
+    z: np.ndarray,
     src_vert_epsg: str,
     dst_vert_epsg: str,
-    src_region: list | tuple | numpy.ndarray | None,
+    src_region: list | tuple | np.ndarray | None,
     cache_dir: str | None,
-) -> numpy.ndarray:
+) -> np.ndarray:
     """Compute and apply a vertical datum shift to z via a cached transformez grid."""
     from scipy.interpolate import RegularGridInterpolator
 
@@ -240,11 +240,11 @@ def _apply_vertical_transform(
         shift_data = src.read(1).astype(float)
         grid_bounds = src.bounds
         if src.nodata is not None:
-            shift_data[numpy.isclose(shift_data, src.nodata, atol=1e-4)] = numpy.nan
+            shift_data[np.isclose(shift_data, src.nodata, atol=1e-4)] = np.nan
 
     height, width = shift_data.shape
-    lons = numpy.linspace(grid_bounds.left, grid_bounds.right, width)
-    lats = numpy.linspace(grid_bounds.bottom, grid_bounds.top, height)
+    lons = np.linspace(grid_bounds.left, grid_bounds.right, width)
+    lats = np.linspace(grid_bounds.bottom, grid_bounds.top, height)
 
     # A point outside the grid must never come back as a 0.0 shift: an unconverted
     # height is indistinguishable from a legitimate zero separation, and silently
@@ -258,7 +258,7 @@ def _apply_vertical_transform(
         | (y > grid_bounds.top)
     )
     if outside.any():
-        n_outside = int(numpy.count_nonzero(outside))
+        n_outside = int(np.count_nonzero(outside))
         raise ValueError(
             f"Vertical transform failed: {n_outside:,} of {outside.size:,} points fall "
             f"outside the shift grid {grid_fn} "
@@ -274,15 +274,15 @@ def _apply_vertical_transform(
         shift_data[::-1, :],
         method="linear",
         bounds_error=False,
-        fill_value=numpy.nan,
+        fill_value=np.nan,
     )
 
-    shifts = interp(numpy.column_stack([y, x]))
+    shifts = interp(np.column_stack([y, x]))
 
     # In-bounds NaNs mean the grid itself has nodata there (a genuine gap in the
     # datum model), not a caching problem. Those heights cannot be converted, so
     # they stay NaN and are dropped downstream rather than silently passed through.
-    n_nan = int(numpy.count_nonzero(numpy.isnan(shifts)))
+    n_nan = int(np.count_nonzero(np.isnan(shifts)))
     if n_nan:
         logger.warning(
             "%d of %d points fall on nodata cells of shift grid %s; "
diff --git a/src/ivert/utils/cuboid_funcs.py b/src/ivert/utils/cuboid_funcs.py
index cd6360c..a593773 100644
--- a/src/ivert/utils/cuboid_funcs.py
+++ b/src/ivert/utils/cuboid_funcs.py
@@ -1,4 +1,4 @@
-import numpy
+import numpy as np
 
 # from itertools import product
 
@@ -477,12 +477,12 @@ def cuboids_intersect_vectorized(
             f"Invalid bbox_order: {bbox_order}. Must be 'point' or 'axis'.",
         )
 
-    xmin = numpy.asarray(xmin)
-    xmax = numpy.asarray(xmax)
-    ymin = numpy.asarray(ymin)
-    ymax = numpy.asarray(ymax)
-    zmin = numpy.asarray(zmin)
-    zmax = numpy.asarray(zmax)
+    xmin = np.asarray(xmin)
+    xmax = np.asarray(xmax)
+    ymin = np.asarray(ymin)
+    ymax = np.asarray(ymax)
+    zmin = np.asarray(zmin)
+    zmax = np.asarray(zmax)
 
     return (
         (xmin < qxmax - tol)
diff --git a/src/ivert/utils/parallel_funcs.py b/src/ivert/utils/parallel_funcs.py
index 7db6cf7..d967a1b 100644
--- a/src/ivert/utils/parallel_funcs.py
+++ b/src/ivert/utils/parallel_funcs.py
@@ -4,7 +4,7 @@
 import sys
 import time
 
-import numpy
+import numpy as np
 import psutil
 import tqdm
 
@@ -30,27 +30,27 @@ def physical_cpu_count():
 # For integers & floats... does not handle character/string arrays.
 # Reference: https://docs.python.org/3/library/array.html
 dtypes_dict = {
-    numpy.int8: "b",
-    numpy.uint8: "B",
-    numpy.int16: "h",
-    numpy.uint16: "H",
-    numpy.int32: "l",
-    numpy.uint32: "L",
-    numpy.int64: "q",
-    numpy.uint64: "Q",
-    numpy.float32: "f",
-    numpy.float64: "d",
+    np.int8: "b",
+    np.uint8: "B",
+    np.int16: "h",
+    np.uint16: "H",
+    np.int32: "l",
+    np.uint32: "L",
+    np.int64: "q",
+    np.uint64: "Q",
+    np.float32: "f",
+    np.float64: "d",
     # Repeat for these expressions of dtype as well.
-    numpy.dtype("int8"): "b",
-    numpy.dtype("uint8"): "B",
-    numpy.dtype("int16"): "h",
-    numpy.dtype("uint16"): "H",
-    numpy.dtype("int32"): "l",
-    numpy.dtype("uint32"): "L",
-    numpy.dtype("int64"): "q",
-    numpy.dtype("uint64"): "Q",
-    numpy.dtype("float32"): "f",
-    numpy.dtype("float64"): "d",
+    np.dtype("int8"): "b",
+    np.dtype("uint8"): "B",
+    np.dtype("int16"): "h",
+    np.dtype("uint16"): "H",
+    np.dtype("int32"): "l",
+    np.dtype("uint32"): "L",
+    np.dtype("int64"): "q",
+    np.dtype("uint64"): "Q",
+    np.dtype("float32"): "f",
+    np.dtype("float64"): "d",
 }
 
 
diff --git a/src/ivert/validate_dem.py b/src/ivert/validate_dem.py
index e85c551..56a2225 100644
--- a/src/ivert/validate_dem.py
+++ b/src/ivert/validate_dem.py
@@ -23,8 +23,8 @@
 import click
 import geopandas
 import numexpr
-import numpy
-import pandas
+import numpy as np
+import pandas as pd
 import pyproj
 import rasterio
 import shapely
@@ -50,7 +50,7 @@
 INTERDECILE_MIN_PHOTONS = 5
 
 
-def read_dataframe_file(df_filename: str) -> pandas.DataFrame:
+def read_dataframe_file(df_filename: str) -> pd.DataFrame:
     """Read a dataframe file, either from a picklefile, HDF, CSV, or feather.
 
     (Can handle other formats by adding more "elif ..." statements in the function.)
@@ -59,13 +59,13 @@ def read_dataframe_file(df_filename: str) -> pandas.DataFrame:
     ext = os.path.splitext(df_filename)[1]
     ext = ext.lower()
     if ext == ".pickle":
-        dataframe = pandas.read_pickle(df_filename)
+        dataframe = pd.read_pickle(df_filename)
     elif ext in (".h5", ".hdf"):
-        dataframe = pandas.read_hdf(df_filename, mode="r")
+        dataframe = pd.read_hdf(df_filename, mode="r")
     elif ext in (".csv", ".txt"):
-        dataframe = pandas.read_csv(df_filename)
+        dataframe = pd.read_csv(df_filename)
     elif ext == ".feather":
-        dataframe = pandas.read_feather(df_filename)
+        dataframe = pd.read_feather(df_filename)
     else:
         raise NotImplementedError(
             f"ERROR: Unknown dataframe file extension '{ext}'. (Currently supporting .pickle, .h5, .hdf, .csv, .txt, or .feather)",
@@ -195,31 +195,31 @@ def validate_dem_child_process(
 
     # Define shared memory arrays here.
     h_shm = shared_memory.SharedMemory(name=height_array_name)
-    heights = numpy.ndarray(array_shape, dtype=height_dtype, buffer=h_shm.buf)
+    heights = np.ndarray(array_shape, dtype=height_dtype, buffer=h_shm.buf)
 
     pi_shm = shared_memory.SharedMemory(name=i_array_name)
-    photon_i = numpy.ndarray(  # noqa: F841 (used by name inside numexpr.evaluate() below)
+    photon_i = np.ndarray(  # noqa: F841 (used by name inside numexpr.evaluate() below)
         array_shape,
         dtype=i_dtype,
         buffer=pi_shm.buf,
     )
 
     pj_shm = shared_memory.SharedMemory(name=j_array_name)
-    photon_j = numpy.ndarray(  # noqa: F841 (used by name inside numexpr.evaluate() below)
+    photon_j = np.ndarray(  # noqa: F841 (used by name inside numexpr.evaluate() below)
         array_shape,
         dtype=j_dtype,
         buffer=pj_shm.buf,
     )
 
     pc_shm = shared_memory.SharedMemory(name=code_array_name)
-    ph_codes = numpy.ndarray(array_shape, dtype=code_dtype, buffer=pc_shm.buf)
+    ph_codes = np.ndarray(array_shape, dtype=code_dtype, buffer=pc_shm.buf)
 
     if measure_coverage:
         x_shm = shared_memory.SharedMemory(name=x_array_name)
-        ph_x = numpy.ndarray(array_shape, dtype=x_dtype, buffer=x_shm.buf)
+        ph_x = np.ndarray(array_shape, dtype=x_dtype, buffer=x_shm.buf)
 
         y_shm = shared_memory.SharedMemory(name=y_array_name)
-        ph_y = numpy.ndarray(array_shape, dtype=y_dtype, buffer=y_shm.buf)
+        ph_y = np.ndarray(array_shape, dtype=y_dtype, buffer=y_shm.buf)
     else:
         x_shm = None
         y_shm = None
@@ -272,22 +272,19 @@ def validate_dem_child_process(
             # Do the work.
             # r_keep marks the cells that had enough photons to validate. Cells left
             # False are dropped from the results below, not reported as empty.
-            r_keep = numpy.zeros((n,), dtype=bool)
-            r_mean = numpy.zeros((n,), dtype=float)
-            r_numphotons = numpy.zeros((n,), dtype=numpy.uint32)
+            r_keep = np.zeros((n,), dtype=bool)
+            r_mean = np.zeros((n,), dtype=float)
+            r_numphotons = np.zeros((n,), dtype=np.uint32)
             r_numphotons_bathy = r_numphotons.copy()
             r_numphotons_intd = r_numphotons.copy()
-            r_std = numpy.zeros((n,), dtype=float)
-            r_interdecile = numpy.zeros((n,), float)
-            r_range = numpy.zeros((n,), heights.dtype)
-            r_10p = numpy.zeros((n,), float)
-            r_90p = numpy.zeros((n,), float)
-            r_dem_elev = numpy.zeros((n,), dtype=float)
-            r_mean_diff = numpy.zeros((n,), dtype=float)
-            if measure_coverage:
-                r_coverage_frac = numpy.zeros((n,), dtype=float)
-            else:
-                r_coverage_frac = None
+            r_std = np.zeros((n,), dtype=float)
+            r_interdecile = np.zeros((n,), float)
+            r_range = np.zeros((n,), heights.dtype)
+            r_10p = np.zeros((n,), float)
+            r_90p = np.zeros((n,), float)
+            r_dem_elev = np.zeros((n,), dtype=float)
+            r_mean_diff = np.zeros((n,), dtype=float)
+            r_coverage_frac = np.zeros((n,), dtype=float) if measure_coverage else None
 
             # 'i' and 'j' look unused, but numexpr.evaluate() resolves the names in
             # its expression string against this frame's locals, so they are read
@@ -298,7 +295,7 @@ def validate_dem_child_process(
                 # Using numexpr.evaluate here is far more memory-and-time efficient than just doing it with the numpy arrays.
                 ph_subset_mask = numexpr.evaluate("(photon_i == i) & (photon_j == j)")
                 # Generate a small pandas dataframe from the subset
-                subset_df = pandas.DataFrame(
+                subset_df = pd.DataFrame(
                     {
                         "height": heights[ph_subset_mask],
                         "ph_code": ph_codes[ph_subset_mask],
@@ -324,10 +321,10 @@ def validate_dem_child_process(
                     # Equal to the geotransform, the y-value starts at the top (max) and iterate downward (negative step.)
                     cell_ystep = (cell_ymin - cell_ymax) / num_subdivisions
 
-                    subset_df["subset_i"] = numpy.floor(
+                    subset_df["subset_i"] = np.floor(
                         (subset_df.ycoord - cell_ymax) / cell_ystep,
                     ).astype(int)
-                    subset_df["subset_j"] = numpy.floor(
+                    subset_df["subset_j"] = np.floor(
                         (subset_df.xcoord - cell_xmin) / cell_xstep,
                     ).astype(int)
 
@@ -359,7 +356,7 @@ def validate_dem_child_process(
                 r_keep[counter] = True
                 r_numphotons[counter] = n_photons
                 r_dem_elev[counter] = dem_elev_list[counter]
-                r_numphotons_bathy[counter] = numpy.count_nonzero(
+                r_numphotons_bathy[counter] = np.count_nonzero(
                     subset_df.ph_code == 40,
                 )
                 r_range[counter] = subset_df.height.max() - subset_df.height.min()
@@ -395,10 +392,10 @@ def validate_dem_child_process(
             # Generate a little dataframe of the outputs for the grid cells that had
             # enough photons to validate. Cells below 'min_photons' were skipped in
             # the loop above and are omitted here rather than reported as empty.
-            results_df = pandas.DataFrame(
+            results_df = pd.DataFrame(
                 {
-                    "i": numpy.asarray(dem_i_list)[r_keep],
-                    "j": numpy.asarray(dem_j_list)[r_keep],
+                    "i": np.asarray(dem_i_list)[r_keep],
+                    "j": np.asarray(dem_j_list)[r_keep],
                     "mean": r_mean[r_keep],
                     "stddev": r_std[r_keep],
                     "numphotons": r_numphotons[r_keep],
@@ -518,10 +515,10 @@ def subdivide_dem(
 
 
 def reset_results_indexes_after_merge(
-    sub_results_df: pandas.DataFrame,
+    sub_results_df: pd.DataFrame,
     sub_dem_fname: str,
     parent_dem_fname: str,
-) -> pandas.DataFrame:
+) -> pd.DataFrame:
     """DEM results dataframes are indexed by (i, j).  Reset the index after merging."""
     if (
         "i" not in sub_results_df.columns and "i" not in sub_results_df.index.names
@@ -806,7 +803,7 @@ class are dropped before any statistics are computed. Default: [1, 6, 40], meani
             # Concatenate the results dataframes.
             output_dfs = []
             for fname in all_fnames:
-                dem_results_df = pandas.read_hdf(fname)
+                dem_results_df = pd.read_hdf(fname)
                 # Now I gotta reset the i,j indexes.
                 sub_dem_name = fname.replace("_results.h5", ".tif")
                 parent_dem_name = dem_name
@@ -817,7 +814,7 @@ class are dropped before any statistics are computed. Default: [1, 6, 40], meani
                 )
                 output_dfs.append(dem_results_df)
 
-            shared_results_df = pandas.concat(output_dfs, ignore_index=False, axis=0)
+            shared_results_df = pd.concat(output_dfs, ignore_index=False, axis=0)
             # After we've combined all the resutls, *then* filter out outliers if they exist.
             if outliers_sd_threshold is not None:
                 assert type(outliers_sd_threshold) in (int, float)
@@ -908,8 +905,8 @@ class are dropped before any statistics are computed. Default: [1, 6, 40], meani
                 for i in range(len(sub_shared_ret_values))
                 if common_key in sub_shared_ret_values[i]
             ]
-            results_df = pandas.concat(
-                [pandas.read_hdf(fname) for fname in all_fnames],
+            results_df = pd.concat(
+                [pd.read_hdf(fname) for fname in all_fnames],
                 ignore_index=True,
                 axis=0,
             )
@@ -1300,15 +1297,15 @@ def _compute_photon_overlap(
             )
             if verbose and excluded_mask.any():
                 print(
-                    f"{numpy.count_nonzero(excluded_mask):,}",
+                    f"{np.count_nonzero(excluded_mask):,}",
                     "photons excluded by exclusion zone(s).",
                 )
             photon_df = photon_df[~excluded_mask]
 
     xstart, xstep, xrot, ystart, yrot, ystep = dem_ds.transform.to_gdal()
     _check_dem_geotransform(dem_ds.name, xstep, ystep, xrot, yrot)
-    photon_df["i"] = numpy.floor((photon_df["dem_y"] - ystart) / ystep).astype(int)
-    photon_df["j"] = numpy.floor((photon_df["dem_x"] - xstart) / xstep).astype(int)
+    photon_df["i"] = np.floor((photon_df["dem_y"] - ystart) / ystep).astype(int)
+    photon_df["j"] = np.floor((photon_df["dem_x"] - xstart) / xstep).astype(int)
 
     photon_df = photon_df[
         (photon_df["i"] >= 0)
@@ -1320,11 +1317,11 @@ def _compute_photon_overlap(
     # Keep only the photon classes requested for this validation. Doing it here,
     # rather than inside each child process, means the height/class arrays copied
     # into shared memory carry only the photons that will actually be used.
-    class_mask = numpy.isin(photon_df["class_code"], classes)
+    class_mask = np.isin(photon_df["class_code"], classes)
     if not class_mask.all():
         if verbose:
             print(
-                f"{numpy.count_nonzero(~class_mask):,} photons dropped as outside the",
+                f"{np.count_nonzero(~class_mask):,} photons dropped as outside the",
                 "requested photon classes",
                 f"({'/'.join(str(c) for c in classes)}).",
             )
@@ -1347,21 +1344,21 @@ def _compute_photon_overlap(
         if dem_ndv is None:
             dem_ndv = EMPTY_VAL
 
-    if numpy.isnan(dem_ndv):
-        dem_goodpixel_mask = ~numpy.isnan(dem_array)
+    if np.isnan(dem_ndv):
+        dem_goodpixel_mask = ~np.isnan(dem_array)
     else:
         dem_goodpixel_mask = dem_array != dem_ndv
 
     photon_df = photon_df.set_index(["i", "j"], drop=False)
-    dem_mask_w_photons = numpy.zeros(dem_array.shape, dtype=bool)
+    dem_mask_w_photons = np.zeros(dem_array.shape, dtype=bool)
     dem_mask_w_photons[photon_df.i, photon_df.j] = 1
 
     dem_overlap_mask = dem_goodpixel_mask & dem_mask_w_photons
-    dem_overlap_i, dem_overlap_j = numpy.where(dem_overlap_mask)
+    dem_overlap_i, dem_overlap_j = np.where(dem_overlap_mask)
     dem_overlap_elevs = dem_array[dem_overlap_mask]
 
     if verbose:
-        num_goodpixels = numpy.count_nonzero(dem_goodpixel_mask)
+        num_goodpixels = np.count_nonzero(dem_goodpixel_mask)
         print(f"{num_goodpixels:,}", "land cells exist in the DEM.")
         if num_goodpixels == 0:
             print(
@@ -1371,10 +1368,10 @@ def _compute_photon_overlap(
         print(
             f"{len(photon_df):,} ICESat-2 photons overlap",
             f"{len(dem_overlap_i):,}",
-            f"DEM cells ({numpy.count_nonzero(dem_overlap_mask) * 100 / num_goodpixels:0.2f}% of total DEM data).",
+            f"DEM cells ({np.count_nonzero(dem_overlap_mask) * 100 / num_goodpixels:0.2f}% of total DEM data).",
         )
 
-    if numpy.count_nonzero(dem_overlap_mask) == 0:
+    if np.count_nonzero(dem_overlap_mask) == 0:
         if verbose:
             print(
                 "No overlapping ICESat-2 data with valid land cells. Stopping and moving on.",
@@ -1418,9 +1415,9 @@ def _run_photon_level_validation(
         print("Performing photon-level validation...")
         print("\tGenerating DEM elevation dataframe... ", end="")
 
-    dem_elev_df = pandas.DataFrame(
+    dem_elev_df = pd.DataFrame(
         {"dem_elevation": dem_overlap_elevs},
-        index=pandas.MultiIndex.from_arrays(
+        index=pd.MultiIndex.from_arrays(
             (dem_overlap_i, dem_overlap_j),
             names=("i", "j"),
         ),
@@ -1431,7 +1428,7 @@ def _run_photon_level_validation(
 
     photon_df_with_dem_elevs = photon_df.join(dem_elev_df, how="left")
     photon_df_with_dem_elevs = photon_df_with_dem_elevs[
-        pandas.notna(photon_df_with_dem_elevs["dem_elevation"])
+        pd.notna(photon_df_with_dem_elevs["dem_elevation"])
     ]
     if verbose:
         print(f"Done with {len(photon_df_with_dem_elevs)} records.")
@@ -1828,13 +1825,11 @@ def _write_validation_outputs(
     if len(results_dataframes_list) == 0:
         return files_to_export
 
-    results_dataframe = pandas.concat(results_dataframes_list)
+    results_dataframe = pd.concat(results_dataframes_list)
     # Cells with too few photons were already dropped by the child processes, which
     # enforce 'min_photons_per_cell'. This only guards against a non-finite mean
     # arising from bad photon elevations.
-    results_dataframe = results_dataframe[
-        numpy.isfinite(results_dataframe["mean"])
-    ].copy()
+    results_dataframe = results_dataframe[np.isfinite(results_dataframe["mean"])].copy()
 
     # Drop cells below the requested minimum ICESat-2 coverage. Coverage is a
     # per-cell property, so filtering here (per subset, before any outlier removal)
@@ -2149,17 +2144,17 @@ def _format_stat(value) -> str:
 
     """
     x = float(value)
-    if not numpy.isfinite(x):
+    if not np.isfinite(x):
         return str(x)
     if x != 0 and abs(x) < 0.10:
         # Enough decimals to show 2 significant digits for small magnitudes.
-        decimals = 1 - int(numpy.floor(numpy.log10(abs(x))))
+        decimals = 1 - int(np.floor(np.log10(abs(x))))
         return f"{x:.{decimals}f}"
     return f"{x:.2f}"
 
 
 def write_summary_stats_file(
-    results_df: pandas.DataFrame,
+    results_df: pd.DataFrame,
     statsfile_name: str,
     verbose: bool = True,
 ) -> None:
@@ -2207,12 +2202,12 @@ def write_summary_stats_file(
         f"Mean bias error (DEM - ICESat-2) (m): {_format_stat(mean_diff.mean())}",
     )
     lines.append(
-        f"RMSE (m): {_format_stat(numpy.sqrt(numpy.mean(numpy.power(mean_diff, 2))))}",
+        f"RMSE (m): {_format_stat(np.sqrt(np.mean(np.power(mean_diff, 2))))}",
     )
 
     lines.append(
         "Number of cells with bathymetry photons: {:d}".format(
-            numpy.count_nonzero(results_df["numphotons_bathy"] > 0),
+            np.count_nonzero(results_df["numphotons_bathy"] > 0),
         ),
     )
 
@@ -2230,7 +2225,7 @@ def write_summary_stats_file(
     )
 
     percentile_levels = [0, 1, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 99, 100]
-    percentile_values = numpy.percentile(mean_diff, percentile_levels)
+    percentile_values = np.percentile(mean_diff, percentile_levels)
     for level, v in zip(percentile_levels, percentile_values, strict=True):
         lines.append(f"    {level:>3d} percentile error level (m): {_format_stat(v)}")
 
@@ -2246,10 +2241,10 @@ def write_summary_stats_file(
         )
         for pct_of_cells in range(100, 0, -10):
             # The coverage threshold that retains this fraction of the best-covered cells.
-            coverage_threshold = numpy.percentile(coverage_frac, 100 - pct_of_cells)
+            coverage_threshold = np.percentile(coverage_frac, 100 - pct_of_cells)
             mask = coverage_frac >= coverage_threshold
             subset_diff = mean_diff[mask]
-            rmse = numpy.sqrt(numpy.mean(numpy.power(subset_diff, 2)))
+            rmse = np.sqrt(np.mean(np.power(subset_diff, 2)))
             lines.append(
                 f"    RMSE for grid cells with >{coverage_threshold * 100:0.1f}% coverage ({pct_of_cells:d}% of cells) (m): {_format_stat(rmse)}",
             )
@@ -2282,7 +2277,7 @@ def generate_result_geotiff(
     """
     xsize, ysize = dem_ds.width, dem_ds.height
     emptyval = float(EMPTY_VAL)
-    result_array = numpy.zeros([ysize, xsize], dtype=numpy.float32) + emptyval
+    result_array = np.zeros([ysize, xsize], dtype=np.float32) + emptyval
 
     indices = results_dataframe.index.to_numpy()
     ivals = [idx[0] for idx in indices]
@@ -2336,8 +2331,8 @@ def _results_cell_centers(results_dataframe, dem_ds):
     """
     gt = dem_ds.transform.to_gdal()
     indices = results_dataframe.index.to_numpy()
-    ivals = numpy.array([idx[0] for idx in indices], dtype=float)
-    jvals = numpy.array([idx[1] for idx in indices], dtype=float)
+    ivals = np.array([idx[0] for idx in indices], dtype=float)
+    jvals = np.array([idx[1] for idx in indices], dtype=float)
     x = gt[0] + (jvals + 0.5) * gt[1] + (ivals + 0.5) * gt[2]
     y = gt[3] + (jvals + 0.5) * gt[4] + (ivals + 0.5) * gt[5]
     return x, y
@@ -2365,7 +2360,7 @@ def _export_errors_vector(results_dataframe, dem_ds, out_fname, fmt, verbose=Tru
     for name, col in fields:
         vals = results_dataframe[col].to_numpy()
         data[name] = (
-            vals.astype(numpy.int32)
+            vals.astype(np.int32)
             if col.startswith("numphotons")
             else vals.astype(float)
         )
@@ -2387,9 +2382,9 @@ def _export_errors_xyz(results_dataframe, dem_ds, out_fname, verbose=True):
     """Write a whitespace-delimited 'x y error' text file, one cell-center point per line."""
     x_centers, y_centers = _results_cell_centers(results_dataframe, dem_ds)
     errors = results_dataframe["diff_mean"].to_numpy()
-    numpy.savetxt(
+    np.savetxt(
         out_fname,
-        numpy.column_stack([x_centers, y_centers, errors]),
+        np.column_stack([x_centers, y_centers, errors]),
         fmt="%.8g",
     )
     if verbose:
diff --git a/src/ivert/validate_dem_collection.py b/src/ivert/validate_dem_collection.py
index 66fc249..a2c2581 100644
--- a/src/ivert/validate_dem_collection.py
+++ b/src/ivert/validate_dem_collection.py
@@ -7,8 +7,8 @@
 import traceback
 
 import click
-import numpy
-import pandas
+import numpy as np
+import pandas as pd
 
 import ivert.icesat2_database_v2
 import ivert.utils.query_yes_no as yes_no
@@ -16,16 +16,16 @@
 
 
 def write_summary_csv_file(
-    total_results_df_or_file: pandas.DataFrame | str,
+    total_results_df_or_file: pd.DataFrame | str,
     list_of_empty_files: list[str] | tuple[str],
     csv_name: str,
     verbose: bool = True,
-) -> pandas.DataFrame:
+) -> pd.DataFrame:
     """Write a summary csv of all the results in a collection, after they've been run."""
     if type(total_results_df_or_file) is str:
-        total_df = pandas.read_hdf(total_results_df_or_file)
+        total_df = pd.read_hdf(total_results_df_or_file)
     else:
-        assert isinstance(total_results_df_or_file, pandas.DataFrame)
+        assert isinstance(total_results_df_or_file, pd.DataFrame)
         total_df = total_results_df_or_file
 
     if "filename" not in total_df.columns:
@@ -35,11 +35,11 @@ def write_summary_csv_file(
     all_filenames = list(unique_files) + list(list_of_empty_files)
     n = len(all_filenames)
 
-    means = numpy.empty((n,), dtype=float)
-    stds = numpy.empty((n,), dtype=float)
-    rmses = numpy.empty((n,), dtype=float)
-    n_cells = numpy.empty((n,), dtype=int)
-    photons_per_cell = numpy.empty((n,), dtype=float)
+    means = np.empty((n,), dtype=float)
+    stds = np.empty((n,), dtype=float)
+    rmses = np.empty((n,), dtype=float)
+    n_cells = np.empty((n,), dtype=int)
+    photons_per_cell = np.empty((n,), dtype=float)
     # canopy_mean = numpy.empty((n,), dtype=float)
     # canopy_mean_gt0 = numpy.empty((n,), dtype=float)
 
@@ -58,15 +58,15 @@ def write_summary_csv_file(
         else:
             # For files with no results, just list n/a for this.
             assert fname in list_of_empty_files
-            means[i] = numpy.nan
-            stds[i] = numpy.nan
-            rmses[i] = numpy.nan
+            means[i] = np.nan
+            stds[i] = np.nan
+            rmses[i] = np.nan
             n_cells[i] = 0
-            photons_per_cell[i] = numpy.nan
+            photons_per_cell[i] = np.nan
             # canopy_mean[i] = numpy.nan
             # canopy_mean_gt0[i] = numpy.nan
 
-    output_df = pandas.DataFrame(
+    output_df = pd.DataFrame(
         data={
             "filename": all_filenames,
             "rmse": rmses,
@@ -231,7 +231,7 @@ def validate_list_of_dems(
         results_df = None
 
         if not os.path.exists(statsfile_name):
-            results_df = pandas.read_hdf(results_h5)
+            results_df = pd.read_hdf(results_h5)
             if verbose:
                 print(results_df, "read.")
             validate_dem.write_summary_stats_file(
@@ -242,7 +242,7 @@ def validate_list_of_dems(
 
         if not os.path.exists(plot_file_name):
             if results_df is None:
-                results_df = pandas.read_hdf(results_h5)
+                results_df = pd.read_hdf(results_h5)
                 if verbose:
                     print(results_df, "read.")
             plot_validation_results.plot_histograms_and_line(