From 33d0d6b170861838021cc7143d6ba8b03bedfb43 Mon Sep 17 00:00:00 2001 From: mmacferrin Date: Tue, 1 Sep 2026 14:24:07 -0600 Subject: [PATCH] Import numpy, pandas and matplotlib under their conventional aliases Clears ICN001 from the ignore list (issue #40), taking it from 59 entries to 58. Nineteen imports and 262 usage sites across 12 files. The package was already inconsistent here: export_vector.py and plot_photon_clouds_v2.py used np and pd, and the latter aliased numpy and pandas while leaving matplotlib bare. This settles it the way the rule and the wider ecosystem expect. The rewrite targets ast.Name node positions, which by construction cover only real identifier references, so comments, docstrings and the import statements themselves were left alone. That matters here: twenty mentions sit in commented-out code, and five sit in docstrings where the full name is the correct prose form -- "returns a pandas.DataFrame" should not become "pd.DataFrame". `import matplotlib.pyplot as plt` and `from matplotlib import ticker` are untouched, since only bare `import matplotlib` was flagged. Equivalence was proven rather than inspected: both the old and new ASTs were normalised to a canonical form, with aliases expanded back to full module names and alias bindings stripped, and compared. All twelve files match. Since docstrings are AST constants, that also proves none were altered, and a separate check confirmed no comment line changed. Two knock-on effects came from the names simply getting shorter. An if/else in validate_dem.py that was previously too long to express as a ternary now fits, so SIM108 -- already enabled -- flagged it; reverting just that ternary makes the file's AST match too, confirming it is the only non-alias change in the diff. The formatter also repacked one file whose wrapped expressions now fit on single lines. --- pyproject.toml | 1 - src/ivert/cli.py | 4 +- src/ivert/convert_error_tif_to_vector.py | 6 +- src/ivert/icesat2_database_v2.py | 104 +++++++-------- src/ivert/icesat2_requests.py | 20 +-- src/ivert/plot_photon_clouds_v2.py | 4 +- src/ivert/plot_results_slope_centrality.py | 38 +++--- src/ivert/plot_validation_results.py | 76 +++++------ src/ivert/transform_points.py | 40 +++--- src/ivert/utils/cuboid_funcs.py | 14 +- src/ivert/utils/parallel_funcs.py | 42 +++--- src/ivert/validate_dem.py | 141 ++++++++++----------- src/ivert/validate_dem_collection.py | 36 +++--- 13 files changed, 260 insertions(+), 266 deletions(-) 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(