Case preparation from raw data

geoswe.data_prep turns the files an agency publishes into the arrays a Solver2D run takes: USGS 3DEP topography and NOAA CUDEM bathymetry into one projected bed, NLCD land cover into a Manning field, a NOAA CO-OPS CSV into a stage series, MRMS GRIB2 files into a rainfall series. Each function below hands back plain arrays; CaseData is a container to collect them in, which nothing in the package returns, so fill it yourself or pass the arrays to Solver2D directly.

Needs the io extra (rasterio, for the raster paths) and the forcings extra (pandas for the CSV and GRIB stamps, scipy for clean_dem):

pip install "geoswe[io,forcings]"

Every one of those imports is function-local, so import geoswe.data_prep works without them and the function you call names what it needs. mrms_to_uniform_timeseries needs more than either extra carries: xarray with the cfgrib engine and cfgrib’s own eccodes library, installed by hand. It raises rather than handing back a series a run cannot tell from a dry forecast.

Both of the DEM filters report what they moved, because a bed that was reshaped in silence is a run that succeeds on the wrong terrain. clip_range defaults to None and clips nothing; when you do pass a range, merge_dems_to_grid’s meta carries n_clipped and clipped_frac beside n_denormal and valid_frac, and clean_dem warns with the same numbers and the DEM’s own elevation range.

The grid

class geoswe.data_prep.CaseData(bed, manning, dx, dy, x0, y0, crs_wkt, rain_time_s=None, rain_rate_ms=None, stage_time_s=None, stage_m=None, stage_coastline_mask=None, gauges=<factory>, meta=<factory>)[source]

Bases: object

All inputs needed to instantiate a Solver2D simulation for a real case.

Parameters:
bed: ndarray
manning: ndarray
dx: float
dy: float
x0: float
y0: float
crs_wkt: str
rain_time_s: ndarray | None = None
rain_rate_ms: ndarray | None = None
stage_time_s: ndarray | None = None
stage_m: ndarray | None = None
stage_coastline_mask: ndarray | None = None
gauges: List[Tuple[str, float, float]]
meta: dict
geoswe.data_prep.merge_dems_to_grid(topo_paths, bathy_paths, bbox_latlon, target_dx, target_crs='EPSG:26917', nodata_fill=-9999.0, denormal_thr=0.001, clip_range=None)[source]

Merge USGS 3DEP topo and NOAA CUDEM bathy on a target UTM grid.

Bathymetry takes priority below MSL (z < 0); topo takes priority above. Returns (bed[nx,ny], metadata).

clip_range is an optional (low, high) pair in metres, handed on to clean_dem. It defaults to None, which clips nothing: the (-15, 50) it used to default to was the elevation range of the Pinellas County domain this pipeline was written for, and on terrain that rises above it the bed came back flattened with nothing in the run to say so. When a caller does pass a range, meta["n_clipped"] and meta["clipped_frac"] report how much of the bed it moved and clean_dem warns with the same numbers.

denormal_thr is the one knob for both near-zero filters, the per-layer one below and clean_dem’s; meta["n_denormal"] counts the cells the first one demoted. See clean_dem for why it is 1e-3 m and not the 1e-30 of a true float32 denormal, and for what that costs at a shoreline.

Parameters:
Return type:

Tuple[ndarray, dict]

geoswe.data_prep.clean_dem(bed, nodata=-9999.0, denormal_thr=0.001, abrupt_jump_m=8.0, clip_range=None, passes=3, stats=None)[source]

Postprocess a merged DEM in (rows, cols) orientation.

Removes three known artifact classes:

  1. Non-zero elevations with |z| < denormal_thr, the near-zero values rasterio’s bilinear reproject leaves at tile edges. The default 1e-3 m is the value this pipeline was tuned with on the Pinellas tiles, not the 1e-30 of the float32 denormals the artifact is named after, so it also demotes real elevations within a millimetre of the datum: at a shoreline those are the beach, and they come back infilled from a neighbour. Pass 0.0 to keep every non-zero value, and do not tighten it to 1e-20 without the tiles in front of you, because the loose threshold is what keeps the 28-m tile-edge cliffs out. Exactly 0.0 m survives either way (the test is |z| > 0), a discontinuity on the datum.

  2. Single-pixel outliers that differ from their 8-neighbour median by more than abrupt_jump_m (tile-boundary cliffs from CUDEM).

  3. Out-of-range elevations, when clip_range is given: a (low, high) pair in metres, clipped to those bounds. It defaults to None, which clips nothing. The (-15, 50) it used to default to was the Pinellas County elevation range: on a 0 to 120 m hillslope it pinned 57.8% of the cells at exactly 50.0 m. Give it only when you know the elevation range of your own terrain: a bound below the ground flattens the landscape, and the run then succeeds on a wrong DEM. A clip that moves any cell warns with the count, the fraction and the DEM’s own range, rather than refusing, because a genuinely bounded study domain is a legitimate thing to ask for. A pair given high-first raises instead: it has no legitimate reading and would return one elevation for every cell.

Invalid cells are then infilled by nearest-neighbour from valid cells via a distance transform. Pass count controls how many outlier-detect/infill rounds we run (3 is enough for the corrupt-cells fraction we see).

stats, when a dict is given, is filled with n_denormal, n_clipped and clipped_frac: merge_dems_to_grid passes one and forwards the numbers into its meta.

Parameters:
Return type:

ndarray

geoswe.data_prep.landcover_to_manning_on_grid(nlcd_path, meta, nodata_manning=0.035)[source]

Read NLCD raster, reproject to the target grid, map class codes → Manning n.

Parameters:
Return type:

ndarray

The forcings and the coast

geoswe.data_prep.load_noaa_tide_csv(csv_path, t0_iso)[source]

Read a NOAA CO-OPS CSV; return (t_s_from_t0, stage_m), sorted by time.

The sort is not cosmetic: StageBoundary bisects the array it is handed, so one row out of order makes every lookup between it and its neighbour read the wrong sample with no error anywhere. A reorder that changes anything warns, because a CSV whose rows are not in time order is usually a download that went wrong rather than one to quietly repair.

Parameters:
Return type:

Tuple[ndarray, ndarray]

geoswe.data_prep.detect_coastline_cells(bed, mask_x_edge=True, mask_y_edge=False, shoreline_band_m=100.0, dx=30.0, min_depth_m=0.0)[source]

Detect cells along the open-water boundary where the stage BC applies.

Default min_depth_m=0 catches any cell with bed below MSL (the near-shore band). Setting min_depth_m=2 restricts to cells in water at least 2 m deep; this is the physically correct way to apply a stage Dirichlet BC, because those cells already contain enough water that overwriting (h, hu, hv) at the surge stage value doesn’t violate continuity. Coastal “wet=True at low tide, dry=True at high tide” cells with bed near 0 violate the BC and create spurious shocks.

shoreline_band_m is the width of the boundary strip (e.g. 200 m = 7 cells at dx=30); within this band, cells satisfying bed < -min_depth_m get the BC.

Parameters:
Return type:

ndarray

geoswe.data_prep.mrms_to_uniform_timeseries(grib_paths, t0_iso, bbox_latlon)[source]

Reduce a list of MRMS 24h-QPE GRIB files to a uniform-in-space time series of rainfall rate (m/s vs t_s_from_t0).

Each file gives the accumulated 24 h precipitation [mm] valid at the file’s stamp, and each file’s value is converted to a mean rate over its own 24 h window (accum / 24 h, NOT a difference of consecutive accumulations). This is only correct when the files are DAILY: for sub-daily files the overlapping 24 h windows would over-count rain by up to 24x, so this function raises unless the file spacing is ~24 h.

Needs xarray with the cfgrib GRIB2 engine and its eccodes library, which no extra of this package installs. A missing one raises, and so does a file list where every file failed to read: the all-zero pair this used to return is indistinguishable from a dry forecast, so the run went on to report “Rainfall (uniform): peak 0.00 mm/h” and flood nothing. Files that fail while others succeed warn with the count and are left out.

The returned pair is sorted by the stamps inside the files, not by filename, with a warning when that changes the order: every consumer of the pair bisects t_s.

Parameters:
Return type:

Tuple[ndarray, ndarray]