From cff5ae04de7c4a3f2c623f6c6e5d87689464e7a9 Mon Sep 17 00:00:00 2001 From: Matthew Love Date: Fri, 27 Mar 2026 22:43:31 -0700 Subject: [PATCH 1/2] ruff formatting --- src/transformez/cli.py | 1 - src/transformez/transform.py | 3 ++- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/transformez/cli.py b/src/transformez/cli.py index 921f46c..ef9f730 100644 --- a/src/transformez/cli.py +++ b/src/transformez/cli.py @@ -15,7 +15,6 @@ import click import logging from transformez import api -from transformez import grid_engine logger = logging.getLogger(__name__) diff --git a/src/transformez/transform.py b/src/transformez/transform.py index 725ca41..5970b7f 100644 --- a/src/transformez/transform.py +++ b/src/transformez/transform.py @@ -413,7 +413,8 @@ def _step_from_hub(self, epsg, ref_type, geoid=None, epoch=None): if proxy_name: # Revert the erroneous HTDP shift to NAD83 (since global is WGS84) total_out -= htdp_shift - if desc_parts: desc_parts.pop() + if desc_parts: + desc_parts.pop() s, d = self._get_global_chain(proxy_name, model='fes2014') if s is not None: From c1144fee5ad511fd4cbef313e8834283f1a03420 Mon Sep 17 00:00:00 2001 From: Matthew Love Date: Fri, 27 Mar 2026 22:44:00 -0700 Subject: [PATCH 2/2] ruff formatting --- src/transformez/__init__.py | 4 + src/transformez/api.py | 57 ++-- src/transformez/cli.py | 71 +++-- src/transformez/definitions.py | 475 +++++++++++++++++++++++++-------- src/transformez/grid_engine.py | 91 +++++-- src/transformez/hooks.py | 30 ++- src/transformez/htdp.py | 56 ++-- src/transformez/modules.py | 45 ++-- src/transformez/spatial.py | 83 +++--- src/transformez/srs.py | 124 +++++---- src/transformez/transform.py | 210 ++++++++++----- src/transformez/utils.py | 13 +- src/transformez/vdatum.py | 40 ++- 13 files changed, 877 insertions(+), 422 deletions(-) diff --git a/src/transformez/__init__.py b/src/transformez/__init__.py index 0588399..e5fd382 100644 --- a/src/transformez/__init__.py +++ b/src/transformez/__init__.py @@ -33,11 +33,13 @@ # Expose the module for fetchez from .modules import TransformezMod + def _find_proj_lib(): """Locate the best available PROJ_LIB path.""" try: import rasterio + r_path = os.path.join(os.path.dirname(rasterio.__file__), "proj_data") if os.path.exists(os.path.join(r_path, "proj.db")): return r_path @@ -53,6 +55,7 @@ def _find_proj_lib(): try: import pyproj + p_path = pyproj.datadir.get_data_dir() if os.path.exists(os.path.join(p_path, "proj.db")): return p_path @@ -61,6 +64,7 @@ def _find_proj_lib(): return None + target_proj_lib = _find_proj_lib() if "PROJ_LIB" in os.environ: diff --git a/src/transformez/api.py b/src/transformez/api.py index fc55db6..17316eb 100644 --- a/src/transformez/api.py +++ b/src/transformez/api.py @@ -46,12 +46,10 @@ def _parse_datum(datum_arg: str) -> Tuple[Optional[int], Optional[str]]: if not datum_arg: return None, None s = str(datum_arg) - if ':' in s: - parts = s.split(':') + if ":" in s: + parts = s.split(":") geoid_part = parts[1] - geoid = ( - geoid_part.split('=')[1] if 'geoid=' in geoid_part else geoid_part - ) + geoid = geoid_part.split("=")[1] if "geoid=" in geoid_part else geoid_part return Datums.get_vdatum_by_name(parts[0]), geoid return Datums.get_vdatum_by_name(s), None @@ -75,8 +73,7 @@ def plot_grid(grid_array, region, title="Vertical Shift Preview"): region_obj = regions[0] masked_data = np.ma.masked_where( - (np.isnan(grid_array)) | (grid_array == -9999) | (grid_array == 0), - grid_array + (np.isnan(grid_array)) | (grid_array == -9999) | (grid_array == 0), grid_array ) if masked_data.count() == 0: @@ -93,13 +90,19 @@ def plot_grid(grid_array, region, title="Vertical Shift Preview"): plt.title(title) plt.xlabel("Longitude") plt.ylabel("Latitude") - plt.grid(True, linestyle=':', alpha=0.6) + plt.grid(True, linestyle=":", alpha=0.6) - stats = (f"Min: {masked_data.min():.3f} m\n" - f"Max: {masked_data.max():.3f} m\n" - f"Mean: {masked_data.mean():.3f} m") - plt.annotate(stats, xy=(0.02, 0.02), xycoords='axes fraction', - bbox=dict(boxstyle="round,pad=0.3", fc="white", ec="gray", alpha=0.8)) + stats = ( + f"Min: {masked_data.min():.3f} m\n" + f"Max: {masked_data.max():.3f} m\n" + f"Mean: {masked_data.mean():.3f} m" + ) + plt.annotate( + stats, + xy=(0.02, 0.02), + xycoords="axes fraction", + bbox=dict(boxstyle="round,pad=0.3", fc="white", ec="gray", alpha=0.8), + ) logger.info("Displaying preview... Close the plot window to continue.") plt.show() @@ -156,12 +159,15 @@ def generate_grid( vt = VerticalTransform( region=region_obj, - nx=nx, ny=ny, - epsg_in=epsg_in, epsg_out=epsg_out, - geoid_in=geoid_in, geoid_out=geoid_out, + nx=nx, + ny=ny, + epsg_in=epsg_in, + epsg_out=epsg_out, + geoid_in=geoid_in, + geoid_out=geoid_out, decay_pixels=decay_pixels, cache_dir=cache_dir, - verbose=verbose + verbose=verbose, ) shift_array, _ = vt._vertical_transform(vt.epsg_in, vt.epsg_out) @@ -184,7 +190,7 @@ def transform_raster( decay_pixels: int = 100, output_raster: Optional[str] = None, cache_dir: Optional[str] = None, - verbose: bool = False + verbose: bool = False, ) -> Optional[str]: """Apply a vertical datum transformation directly to an existing raster file. @@ -209,9 +215,7 @@ def transform_raster( with rasterio.open(input_raster) as src: bounds = src.bounds - region_obj = ( - Region(bounds.left, bounds.right, bounds.bottom, bounds.top) - ) + region_obj = Region(bounds.left, bounds.right, bounds.bottom, bounds.top) nx, ny = src.width, src.height epsg_in, geoid_in = _parse_datum(datum_in) @@ -227,12 +231,15 @@ def transform_raster( vt = VerticalTransform( region=region_obj, - nx=nx, ny=ny, - epsg_in=epsg_in, epsg_out=epsg_out, - geoid_in=geoid_in, geoid_out=geoid_out, + nx=nx, + ny=ny, + epsg_in=epsg_in, + epsg_out=epsg_out, + geoid_in=geoid_in, + geoid_out=geoid_out, decay_pixels=decay_pixels, cache_dir=cache_dir, - verbose=verbose + verbose=verbose, ) shift_array, _ = vt._vertical_transform(vt.epsg_in, vt.epsg_out) diff --git a/src/transformez/cli.py b/src/transformez/cli.py index ef9f730..12f4d99 100644 --- a/src/transformez/cli.py +++ b/src/transformez/cli.py @@ -20,7 +20,7 @@ @click.group(name="transform") -@click.version_option(package_name='transformez') +@click.version_option(package_name="transformez") def transformez_cli(): """Apply vertical datum transformations and generate shift grids.""" @@ -29,14 +29,32 @@ def transformez_cli(): @transformez_cli.command("run") @click.argument("input_file", required=False) -@click.option("-R", "--region", help="Bounding box or location string (if no input file).") -@click.option("-E", "--increment", help="Resolution (e.g., 1s, 30m) (if no input file).") -@click.option("-I", "--input-datum", required=True, help="Source Datum (e.g., 'mllw', '5703').") -@click.option("-O", "--output-datum", required=True, help="Target Datum (e.g., '4979', '5703:g2012b').") +@click.option( + "-R", "--region", help="Bounding box or location string (if no input file)." +) +@click.option( + "-E", "--increment", help="Resolution (e.g., 1s, 30m) (if no input file)." +) +@click.option( + "-I", "--input-datum", required=True, help="Source Datum (e.g., 'mllw', '5703')." +) +@click.option( + "-O", + "--output-datum", + required=True, + help="Target Datum (e.g., '4979', '5703:g2012b').", +) @click.option("--out", "-o", help="Output filename (default: auto-named).") -@click.option("--decay-pixels", type=int, default=100, help="Number of pixels to decay tidal shifts inland.") +@click.option( + "--decay-pixels", + type=int, + default=100, + help="Number of pixels to decay tidal shifts inland.", +) @click.option("--preview", is_flag=True, help="Preview the transformation output.") -def transform_run(input_file, region, increment, input_datum, output_datum, out, decay_pixels, preview): +def transform_run( + input_file, region, increment, input_datum, output_datum, out, decay_pixels, preview +): """Transform a raster's vertical datum or generate a standalone shift grid. If an INPUT_FILE is provided, that specific raster is transformed in place. @@ -57,17 +75,23 @@ def transform_run(input_file, region, increment, input_datum, output_datum, out, datum_out=output_datum, decay_pixels=decay_pixels, output_raster=out, - verbose=True + verbose=True, ) if result: - click.secho(f"Successfully transformed raster: {result}", fg="green", bold=True) + click.secho( + f"Successfully transformed raster: {result}", fg="green", bold=True + ) else: click.secho("Failed to transform raster.", fg="red") sys.exit(1) elif region and increment: - click.secho(f"Generating vertical shift grid for region: {region}...", fg="cyan", bold=True) + click.secho( + f"Generating vertical shift grid for region: {region}...", + fg="cyan", + bold=True, + ) click.echo(f" Shift: {input_datum} āž” {output_datum} @ {increment}") # Auto-generate an output name if one wasn't provided @@ -80,20 +104,25 @@ def transform_run(input_file, region, increment, input_datum, output_datum, out, datum_out=output_datum, decay_pixels=decay_pixels, out_fn=out_fn, - verbose=True + verbose=True, ) if preview: api.plot_grid(result, region) if result is not None: - click.secho(f"Successfully generated shift grid: {out_fn}", fg="green", bold=True) + click.secho( + f"Successfully generated shift grid: {out_fn}", fg="green", bold=True + ) else: click.secho("Failed to generate shift grid.", fg="red") sys.exit(1) else: - click.secho("Error: You must provide either an INPUT_FILE or both --region and --increment.", fg="red") + click.secho( + "Error: You must provide either an INPUT_FILE or both --region and --increment.", + fg="red", + ) sys.exit(1) @@ -106,7 +135,7 @@ def transform_list(): click.secho("\nšŸŒŠ Supported Tidal Surfaces:", fg="cyan", bold=True) # For tidal datums, the user types the dictionary key (e.g., 'mllw') for key, v in Datums.SURFACES.items(): - region_str = v.get('region', 'global').upper() + region_str = v.get("region", "global").upper() click.echo(f" {key:<12} : {v.get('name', key):<30} [{region_str}]") click.secho("\nšŸŒ Ellipsoidal / Frame Datums (EPSG):", fg="cyan", bold=True) @@ -118,19 +147,23 @@ def transform_list(): # For orthometric, the key in Datums.CDN is typically the EPSG code (e.g., '5703') for epsg_key, v in Datums.CDN.items(): # Fallback to the key if 'epsg' isn't explicitly defined in the dict - epsg_code = v.get('epsg', epsg_key) - geoid_str = v.get('default_geoid', 'None') - click.echo(f" {str(epsg_code):<12} : {v.get('name', 'Unknown'):<30} (Default Geoid: {geoid_str})") + epsg_code = v.get("epsg", epsg_key) + geoid_str = v.get("default_geoid", "None") + click.echo( + f" {str(epsg_code):<12} : {v.get('name', 'Unknown'):<30} (Default Geoid: {geoid_str})" + ) click.secho("\nšŸŒ Available Geoids:", fg="cyan", bold=True) click.echo(f" {', '.join(Datums.GEOIDS.keys())}") click.secho("\nšŸ’” Pro-Tip:", fg="yellow", bold=True, nl=False) - click.echo(" Combine an EPSG and a specific Geoid using a colon (e.g., -O 5703:g2012b)\n") + click.echo( + " Combine an EPSG and a specific Geoid using a colon (e.g., -O 5703:g2012b)\n" + ) except ImportError: click.secho("Error: Could not load Transformez datum definitions.", fg="red") -if __name__ == '__main__': +if __name__ == "__main__": transformez_cli() diff --git a/src/transformez/definitions.py b/src/transformez/definitions.py index c80b3a7..8ea72e0 100644 --- a/src/transformez/definitions.py +++ b/src/transformez/definitions.py @@ -25,105 +25,366 @@ class Datums: # ========================================================================= SURFACES = { # --- Tidal Datums --- - 1089: {'name': 'mllw', 'description': 'Mean Lower Low Water', 'region': 'usa'}, - 5866: {'name': 'mllw', 'description': 'Mean Lower Low Water', 'region': 'usa'}, - 1091: {'name': 'mlw', 'description': 'Mean Low Water', 'region': 'usa'}, - 5869: {'name': 'mhhw', 'description': 'Mean Higher High Water', 'region': 'usa'}, - 5868: {'name': 'mhw', 'description': 'Mean High Water', 'region': 'usa'}, - 5714: {'name': 'msl', 'description': 'Mean Sea Level', 'region': 'usa'}, - 5713: {'name': 'mtl', 'description': 'Mean Tide Level', 'region': 'usa'}, - + 1089: {"name": "mllw", "description": "Mean Lower Low Water", "region": "usa"}, + 5866: {"name": "mllw", "description": "Mean Lower Low Water", "region": "usa"}, + 1091: {"name": "mlw", "description": "Mean Low Water", "region": "usa"}, + 5869: { + "name": "mhhw", + "description": "Mean Higher High Water", + "region": "usa", + }, + 5868: {"name": "mhw", "description": "Mean High Water", "region": "usa"}, + 5714: {"name": "msl", "description": "Mean Sea Level", "region": "usa"}, + 5713: {"name": "mtl", "description": "Mean Tide Level", "region": "usa"}, # --- Hydraulic / River Datums --- # Columbia River Datum (No standard EPSG, using 0 placeholder or custom) - 0: {'name': 'crd', 'description': 'Columbia River Datum', 'uncertainty': 0, 'epsg': 0, 'region': 'usa'}, - + 0: { + "name": "crd", + "description": "Columbia River Datum", + "uncertainty": 0, + "epsg": 0, + "region": "usa", + }, # IGLD 1985 (Dynamic Height) - 5609: {'name': 'IGLD85', 'description': 'International Great Lakes Datum 1985', 'uncertainty': 0, 'epsg': 5609, 'region': 'usa'}, - + 5609: { + "name": "IGLD85", + "description": "International Great Lakes Datum 1985", + "uncertainty": 0, + "epsg": 5609, + "region": "usa", + }, # IGLD Low Water Datum (Chart Datum for Lakes) # VDatum uses 'LWD_IGLD85' string - 9000: {'name': 'LWD_IGLD85', 'description': 'IGLD85 Low Water Datum', 'uncertainty': 0, 'epsg': 5609, 'region': 'usa'}, - + 9000: { + "name": "LWD_IGLD85", + "description": "IGLD85 Low Water Datum", + "uncertainty": 0, + "epsg": 5609, + "region": "usa", + }, # --- Legacy Vertical --- # NGVD29 is often best handled via VDatum (VERTCON) if PROJ isn't configured - 5702: {'name': 'NGVD29', 'description': 'National Geodetic Vertical Datum 1929', 'uncertainty': 0.05, 'epsg': 5702}, - + 5702: { + "name": "NGVD29", + "description": "National Geodetic Vertical Datum 1929", + "uncertainty": 0.05, + "epsg": 5702, + }, # --- Global Tidal Datums (DTU/FES) --- # Using pseudo-EPSG codes or negative placeholders for custom logic - 9001: {'name': 'lat', 'description': 'Lowest Astronomical Tide', 'region': 'global'}, - 9002: {'name': 'hat', 'description': 'Highest Astronomical Tide', 'region': 'global'}, - 9003: {'name': 'mss', 'description': 'Mean Sea Surface', 'region': 'global'}, + 9001: { + "name": "lat", + "description": "Lowest Astronomical Tide", + "region": "global", + }, + 9002: { + "name": "hat", + "description": "Highest Astronomical Tide", + "region": "global", + }, + 9003: {"name": "mss", "description": "Mean Sea Surface", "region": "global"}, } GLOBAL_ALIASES = { - 'mllw': 'lat', # Mean Lower Low Water -> Lowest Astronomical Tide - 'mlw': 'lat', # Mean Low Water -> LAT (Conservative) - 'mhhw': 'hat', # Mean Higher High Water -> Highest Astronomical Tide - 'mhw': 'hat', # Mean High Water -> HAT (Conservative) - 'msl': 'mss', # Mean Sea Level -> Mean Sea Surface - 'mtl': 'mss', # Mean Tide Level -> MSS - 'dtl': 'mss', # Diurnal Tide Level -> MSS + "mllw": "lat", # Mean Lower Low Water -> Lowest Astronomical Tide + "mlw": "lat", # Mean Low Water -> LAT (Conservative) + "mhhw": "hat", # Mean Higher High Water -> Highest Astronomical Tide + "mhw": "hat", # Mean High Water -> HAT (Conservative) + "msl": "mss", # Mean Sea Level -> Mean Sea Surface + "mtl": "mss", # Mean Tide Level -> MSS + "dtl": "mss", # Diurnal Tide Level -> MSS } HTDP = { - 4269: {'name': 'NAD_83(2011/CORS96/2007)', 'description': '(North American plate fixed)', 'htdp_id': 1, 'uncertainty': .02, 'epoch': 1997.0}, - 6781: {'name': 'NAD_83(2011/CORS96/2007)', 'description': '(North American plate fixed)', 'htdp_id': 1, 'uncertainty': .02, 'epoch': 1997.0}, - 6319: {'name': 'NAD_83(2011/CORS96/2007)', 'description': '(North American plate fixed)', 'htdp_id': 1, 'uncertainty': .02, 'epoch': 1997.0}, - 6321: {'name': 'NAD_83(PA11/PACP00)', 'description': '(Pacific plate fixed)', 'htdp_id': 2, 'uncertainty': .02, 'epoch': 1997.0}, - 6324: {'name': 'NAD_83(MA11/MARP00)', 'description': '(Mariana plate fixed)', 'htdp_id': 3, 'uncertainty': .02, 'epoch': 1997.0}, - 4979: {'name': 'WGS_84(original)', 'description': '(NAD_83(2011) used)', 'htdp_id': 4, 'uncertainty': 0, 'epoch': 1997.0}, - 7815: {'name': 'WGS_84(original)', 'description': '(NAD_83(2011) used)', 'htdp_id': 4, 'uncertainty': 0, 'epoch': 1997.0}, - 7816: {'name': 'WGS_84(original)', 'description': '(NAD_83(2011) used)', 'htdp_id': 4, 'uncertainty': 0, 'epoch': 1997.0}, - 7656: {'name': 'WGS_84(G730)', 'description': '(ITRF91 used)', 'htdp_id': 5, 'uncertainty': 0, 'epoch': 1997.0}, - 7657: {'name': 'WGS_84(G730)', 'description': '(ITRF91 used)', 'htdp_id': 5, 'uncertainty': 0, 'epoch': 1997.0}, - 7658: {'name': 'WGS_84(G873)', 'description': '(ITRF94 used)', 'htdp_id': 6, 'uncertainty': 0, 'epoch': 1997.0}, - 7659: {'name': 'WGS_84(G873)', 'description': '(ITRF94 used)', 'htdp_id': 6, 'uncertainty': 0, 'epoch': 1997.0}, - 7660: {'name': 'WGS_84(G1150)', 'description': '(ITRF2000 used)', 'htdp_id': 7, 'uncertainty': 0, 'epoch': 1997.0}, - 7661: {'name': 'WGS_84(G1150)', 'description': '(ITRF2000 used)', 'htdp_id': 7, 'uncertainty': 0, 'epoch': 1997.0}, - 7662: {'name': 'WGS_84(G1674)', 'description': '(ITRF2008 used)', 'htdp_id': 8, 'uncertainty': 0, 'epoch': 2000.0}, - 7663: {'name': 'WGS_84(G1674)', 'description': '(ITRF2008 used)', 'htdp_id': 8, 'uncertainty': 0, 'epoch': 2000.0}, - 7664: {'name': 'WGS_84(G1762)', 'description': '(IGb08 used)', 'htdp_id': 9, 'uncertainty': 0, 'epoch': 2000.0}, - 7665: {'name': 'WGS_84(G1762)', 'description': '(IGb08 used)', 'htdp_id': 9, 'uncertainty': 0, 'epoch': 2000.0}, - 7666: {'name': 'WGS_84(G2139)', 'description': '(ITRF2014=IGS14=IGb14 used)', 'htdp_id': 10, 'uncertainty': 0, 'epoch': 1997.0}, - 7667: {'name': 'WGS_84(G2139)', 'description': '(ITRF2014=IGS14=IGb14 used)', 'htdp_id': 10, 'uncertainty': 0, 'epoch': 1997.0}, - 4910: {'name': 'ITRF88', 'description': '', 'htdp_id': 11, 'uncertainty': 0, 'epoch': 1988.0}, - 4911: {'name': 'ITRF89', 'description': '', 'htdp_id': 12, 'uncertainty': 0, 'epoch': 1988.0}, - 7901: {'name': 'ITRF89', 'description': '', 'htdp_id': 12, 'uncertainty': 0, 'epoch': 1988.0}, - 7902: {'name': 'ITRF90', 'description': '(PNEOS90/NEOS90)', 'htdp_id': 13, 'uncertainty': 0, 'epoch': 1988.0}, - 7903: {'name': 'ITRF91', 'description': '', 'htdp_id': 14, 'uncertainty': 0, 'epoch': 1988.0}, - 7904: {'name': 'ITRF92', 'description': '', 'htdp_id': 15, 'uncertainty': 0, 'epoch': 1988.0}, - 7905: {'name': 'ITRF93', 'description': '', 'htdp_id': 16, 'uncertainty': 0, 'epoch': 1988.0}, - 7906: {'name': 'ITRF94', 'description': '', 'htdp_id': 17, 'uncertainty': 0, 'epoch': 1988.0}, - 7907: {'name': 'ITRF96', 'description': '', 'htdp_id': 18, 'uncertainty': 0, 'epoch': 1996.0}, - 7908: {'name': 'ITRF97', 'description': 'IGS97', 'htdp_id': 19, 'uncertainty': 0, 'epoch': 1997.0}, - 7909: {'name': 'ITRF2000', 'description': 'IGS00/IGb00', 'htdp_id': 20, 'uncertainty': 0, 'epoch': 2000.0}, - 7910: {'name': 'ITRF2005', 'description': 'IGS05', 'htdp_id': 21, 'uncertainty': 0, 'epoch': 2000.0}, - 7911: {'name': 'ITRF2008', 'description': 'IGS08/IGb08', 'htdp_id': 22, 'uncertainty': 0, 'epoch': 2000.0}, - 7912: {'name': 'ELLIPSOID', 'description': 'IGS14/IGb14/WGS84/ITRF2014 Ellipsoid', 'htdp_id': 23, 'uncertainty': 0, 'epoch': 2000.0}, - 1322: {'name': 'ITRF2020', 'description': 'IGS20', 'htdp_id': 24, 'uncertainty': 0, 'epoch': 2000.0}, - + 4269: { + "name": "NAD_83(2011/CORS96/2007)", + "description": "(North American plate fixed)", + "htdp_id": 1, + "uncertainty": 0.02, + "epoch": 1997.0, + }, + 6781: { + "name": "NAD_83(2011/CORS96/2007)", + "description": "(North American plate fixed)", + "htdp_id": 1, + "uncertainty": 0.02, + "epoch": 1997.0, + }, + 6319: { + "name": "NAD_83(2011/CORS96/2007)", + "description": "(North American plate fixed)", + "htdp_id": 1, + "uncertainty": 0.02, + "epoch": 1997.0, + }, + 6321: { + "name": "NAD_83(PA11/PACP00)", + "description": "(Pacific plate fixed)", + "htdp_id": 2, + "uncertainty": 0.02, + "epoch": 1997.0, + }, + 6324: { + "name": "NAD_83(MA11/MARP00)", + "description": "(Mariana plate fixed)", + "htdp_id": 3, + "uncertainty": 0.02, + "epoch": 1997.0, + }, + 4979: { + "name": "WGS_84(original)", + "description": "(NAD_83(2011) used)", + "htdp_id": 4, + "uncertainty": 0, + "epoch": 1997.0, + }, + 7815: { + "name": "WGS_84(original)", + "description": "(NAD_83(2011) used)", + "htdp_id": 4, + "uncertainty": 0, + "epoch": 1997.0, + }, + 7816: { + "name": "WGS_84(original)", + "description": "(NAD_83(2011) used)", + "htdp_id": 4, + "uncertainty": 0, + "epoch": 1997.0, + }, + 7656: { + "name": "WGS_84(G730)", + "description": "(ITRF91 used)", + "htdp_id": 5, + "uncertainty": 0, + "epoch": 1997.0, + }, + 7657: { + "name": "WGS_84(G730)", + "description": "(ITRF91 used)", + "htdp_id": 5, + "uncertainty": 0, + "epoch": 1997.0, + }, + 7658: { + "name": "WGS_84(G873)", + "description": "(ITRF94 used)", + "htdp_id": 6, + "uncertainty": 0, + "epoch": 1997.0, + }, + 7659: { + "name": "WGS_84(G873)", + "description": "(ITRF94 used)", + "htdp_id": 6, + "uncertainty": 0, + "epoch": 1997.0, + }, + 7660: { + "name": "WGS_84(G1150)", + "description": "(ITRF2000 used)", + "htdp_id": 7, + "uncertainty": 0, + "epoch": 1997.0, + }, + 7661: { + "name": "WGS_84(G1150)", + "description": "(ITRF2000 used)", + "htdp_id": 7, + "uncertainty": 0, + "epoch": 1997.0, + }, + 7662: { + "name": "WGS_84(G1674)", + "description": "(ITRF2008 used)", + "htdp_id": 8, + "uncertainty": 0, + "epoch": 2000.0, + }, + 7663: { + "name": "WGS_84(G1674)", + "description": "(ITRF2008 used)", + "htdp_id": 8, + "uncertainty": 0, + "epoch": 2000.0, + }, + 7664: { + "name": "WGS_84(G1762)", + "description": "(IGb08 used)", + "htdp_id": 9, + "uncertainty": 0, + "epoch": 2000.0, + }, + 7665: { + "name": "WGS_84(G1762)", + "description": "(IGb08 used)", + "htdp_id": 9, + "uncertainty": 0, + "epoch": 2000.0, + }, + 7666: { + "name": "WGS_84(G2139)", + "description": "(ITRF2014=IGS14=IGb14 used)", + "htdp_id": 10, + "uncertainty": 0, + "epoch": 1997.0, + }, + 7667: { + "name": "WGS_84(G2139)", + "description": "(ITRF2014=IGS14=IGb14 used)", + "htdp_id": 10, + "uncertainty": 0, + "epoch": 1997.0, + }, + 4910: { + "name": "ITRF88", + "description": "", + "htdp_id": 11, + "uncertainty": 0, + "epoch": 1988.0, + }, + 4911: { + "name": "ITRF89", + "description": "", + "htdp_id": 12, + "uncertainty": 0, + "epoch": 1988.0, + }, + 7901: { + "name": "ITRF89", + "description": "", + "htdp_id": 12, + "uncertainty": 0, + "epoch": 1988.0, + }, + 7902: { + "name": "ITRF90", + "description": "(PNEOS90/NEOS90)", + "htdp_id": 13, + "uncertainty": 0, + "epoch": 1988.0, + }, + 7903: { + "name": "ITRF91", + "description": "", + "htdp_id": 14, + "uncertainty": 0, + "epoch": 1988.0, + }, + 7904: { + "name": "ITRF92", + "description": "", + "htdp_id": 15, + "uncertainty": 0, + "epoch": 1988.0, + }, + 7905: { + "name": "ITRF93", + "description": "", + "htdp_id": 16, + "uncertainty": 0, + "epoch": 1988.0, + }, + 7906: { + "name": "ITRF94", + "description": "", + "htdp_id": 17, + "uncertainty": 0, + "epoch": 1988.0, + }, + 7907: { + "name": "ITRF96", + "description": "", + "htdp_id": 18, + "uncertainty": 0, + "epoch": 1996.0, + }, + 7908: { + "name": "ITRF97", + "description": "IGS97", + "htdp_id": 19, + "uncertainty": 0, + "epoch": 1997.0, + }, + 7909: { + "name": "ITRF2000", + "description": "IGS00/IGb00", + "htdp_id": 20, + "uncertainty": 0, + "epoch": 2000.0, + }, + 7910: { + "name": "ITRF2005", + "description": "IGS05", + "htdp_id": 21, + "uncertainty": 0, + "epoch": 2000.0, + }, + 7911: { + "name": "ITRF2008", + "description": "IGS08/IGb08", + "htdp_id": 22, + "uncertainty": 0, + "epoch": 2000.0, + }, + 7912: { + "name": "ELLIPSOID", + "description": "IGS14/IGb14/WGS84/ITRF2014 Ellipsoid", + "htdp_id": 23, + "uncertainty": 0, + "epoch": 2000.0, + }, + 1322: { + "name": "ITRF2020", + "description": "IGS20", + "htdp_id": 24, + "uncertainty": 0, + "epoch": 2000.0, + }, } CDN = { # CONUS / Alaska / Hawaii / PR / VI - 5703: {'name': 'NAVD88 height', 'vdatum_id': 'navd88:m:height', 'default_geoid': 'g2018', 'ellipsoid': 6319}, - 6360: {'name': 'NAVD88 height (usFt)', 'default_geoid': 'g2018'}, - 8228: {'name': 'NAVD88 height (Ft)', 'default_geoid': 'g2018'}, - + 5703: { + "name": "NAVD88 height", + "vdatum_id": "navd88:m:height", + "default_geoid": "g2018", + "ellipsoid": 6319, + }, + 6360: {"name": "NAVD88 height (usFt)", "default_geoid": "g2018"}, + 8228: {"name": "NAVD88 height (Ft)", "default_geoid": "g2018"}, # Puerto Rico - 6641: {'name': 'PRVD02 height', 'vdatum_id': 'prvd02:m:height', 'default_geoid': 'g2018', 'ellipsoid': 6319}, - + 6641: { + "name": "PRVD02 height", + "vdatum_id": "prvd02:m:height", + "default_geoid": "g2018", + "ellipsoid": 6319, + }, # Virgin Islands - 6642: {'name': 'VIVD09 height', 'vdatum_id': 'vivd09:m:height', 'default_geoid': 'g2018', 'ellipsoid': 6319}, - + 6642: { + "name": "VIVD09 height", + "vdatum_id": "vivd09:m:height", + "default_geoid": "g2018", + "ellipsoid": 6319, + }, # Canada (CGVD2013 uses CGG2013 geoid) # Note: You need to ensure 'CGG2013' is fetchable via your fetcher or map it to a filename - 6647: {'name': 'CGVD2013(CGG2013)', 'vdatum_id': 'cgvd2013:m:height', 'default_geoid': 'CGG2013'}, - + 6647: { + "name": "CGVD2013(CGG2013)", + "vdatum_id": "cgvd2013:m:height", + "default_geoid": "CGG2013", + }, # Global EGM - 3855: {'name': 'EGM2008 height', 'vdatum_id': 'egm2008:m:height', 'default_geoid': 'egm2008'}, - 5773: {'name': 'EGM96 height', 'vdatum_id': 'egm96:m:height', 'default_geoid': 'egm96'}, - + 3855: { + "name": "EGM2008 height", + "vdatum_id": "egm2008:m:height", + "default_geoid": "egm2008", + }, + 5773: { + "name": "EGM96 height", + "vdatum_id": "egm96:m:height", + "default_geoid": "egm96", + }, # # Ellipsoidal (Hubs) - No Geoid needed # 6319: {'name': 'NAD83(2011)', 'vdatum_id': 'nad83_2011:m:height'}, # 4979: {'name': 'WGS84', 'vdatum_id': 'wgs84:m:height'}, @@ -131,39 +392,21 @@ class Datums: GEOIDS = { # Standard PROJ-CDN Geoids (Default provider is 'proj') - 'g2018': {'name': 'geoid 2018', 'uncertainty': .0127, 'provider': 'proj'}, - 'g2012b': {'name': 'geoid 2012b', 'uncertainty': .017, 'provider': 'proj'}, - 'geoid09': {'name': 'geoid 2009', 'uncertainty': .05, 'provider': 'proj'}, - + "g2018": {"name": "geoid 2018", "uncertainty": 0.0127, "provider": "proj"}, + "g2012b": {"name": "geoid 2012b", "uncertainty": 0.017, "provider": "proj"}, + "geoid09": {"name": "geoid 2009", "uncertainty": 0.05, "provider": "proj"}, # New XGEOIDs via VDatum (Provider is 'vdatum') - 'xgeoid20b': {'name': 'xgeoid20b', 'uncertainty': .02, 'provider': 'vdatum'}, - 'xgeoid19b': {'name': 'xgeoid19b', 'uncertainty': .02, 'provider': 'vdatum'}, - - 'egm2008': {'name': 'EGM2008', 'uncertainty': 0, 'provider': 'proj'}, - 'egm96': {'name': 'EGM96', 'uncertainty': 0, 'provider': 'proj'}, - - 'CGG2013': {'name': 'CGG2013', 'uncertainty': 0.01, 'provider': 'proj'}, + "xgeoid20b": {"name": "xgeoid20b", "uncertainty": 0.02, "provider": "vdatum"}, + "xgeoid19b": {"name": "xgeoid19b", "uncertainty": 0.02, "provider": "vdatum"}, + "egm2008": {"name": "EGM2008", "uncertainty": 0, "provider": "proj"}, + "egm96": {"name": "EGM96", "uncertainty": 0, "provider": "proj"}, + "CGG2013": {"name": "CGG2013", "uncertainty": 0.01, "provider": "proj"}, } MODELS = { - 'fes2014': { - 'provider': 'seanoe', - 'grids': { - 'lat': 'LAT', - 'mss': 'MSL' - } - }, - 'dtu10': { - 'provider': 'dtu', - 'grids': { - 'mss': 'mss', - 'mdt': 'mdt' - } - }, - 'egm2008': { - 'provider': 'proj', - 'grid': 'egm2008' - } + "fes2014": {"provider": "seanoe", "grids": {"lat": "LAT", "mss": "MSL"}}, + "dtu10": {"provider": "dtu", "grids": {"mss": "mss", "mdt": "mdt"}}, + "egm2008": {"provider": "proj", "grid": "egm2008"}, } @classmethod @@ -188,7 +431,7 @@ def get_vdatum_by_name(cls, datum_name, region_check=None): # Direct Match (e.g. 'mllw' -> 5866) for epsg, info in cls.SURFACES.items(): - if s_name == info['name']: + if s_name == info["name"]: return epsg # We should check for aliases here @@ -200,12 +443,12 @@ def get_global_proxy(cls, datum_name): s_name = str(datum_name).lower() - if s_name in ['lat', 'hat', 'mss']: + if s_name in ["lat", "hat", "mss"]: return s_name for epsg, info in cls.SURFACES.items(): if str(epsg) == s_name: - s_name = info['name'] + s_name = info["name"] break return cls.GLOBAL_ALIASES.get(s_name) @@ -216,14 +459,14 @@ def get_frame_type(cls, epsg): if epsg in cls.SURFACES: # Distinguish between NOAA VDatum and Global - if cls.SURFACES[epsg].get('region') == 'global': - return 'global_tidal' - return 'surface' + if cls.SURFACES[epsg].get("region") == "global": + return "global_tidal" + return "surface" if epsg in cls.HTDP: - return 'htdp' + return "htdp" if epsg in cls.CDN: - return 'cdn' + return "cdn" return None @@ -237,7 +480,7 @@ def get_default_geoid(cls, epsg): return None if e_int in cls.CDN: - return cls.CDN[e_int].get('default_geoid') + return cls.CDN[e_int].get("default_geoid") return None @classmethod @@ -245,12 +488,12 @@ def get_vdatum_id(cls, epsg): """Retrieve the NOAA VDatum CLI string for an EPSG.""" if epsg in cls.SURFACES: - return cls.SURFACES[epsg].get('vdatum_id') + return cls.SURFACES[epsg].get("vdatum_id") if epsg in cls.CDN: - return cls.CDN[epsg].get('vdatum_id') + return cls.CDN[epsg].get("vdatum_id") if epsg == 6319: - return 'nad83_2011:m:height' + return "nad83_2011:m:height" return None diff --git a/src/transformez/grid_engine.py b/src/transformez/grid_engine.py index edcaa5e..f837c85 100644 --- a/src/transformez/grid_engine.py +++ b/src/transformez/grid_engine.py @@ -34,8 +34,7 @@ def plot_grid(grid_array, region, title="Vertical Shift Preview"): return masked_data = np.ma.masked_where( - (np.isnan(grid_array)) | (grid_array == -9999) | (grid_array == 0), - grid_array + (np.isnan(grid_array)) | (grid_array == -9999) | (grid_array == 0), grid_array ) if masked_data.count() == 0: @@ -52,13 +51,19 @@ def plot_grid(grid_array, region, title="Vertical Shift Preview"): plt.title(title) plt.xlabel("Longitude") plt.ylabel("Latitude") - plt.grid(True, linestyle=':', alpha=0.6) + plt.grid(True, linestyle=":", alpha=0.6) - stats = (f"Min: {masked_data.min():.3f} m\n" - f"Max: {masked_data.max():.3f} m\n" - f"Mean: {masked_data.mean():.3f} m") - plt.annotate(stats, xy=(0.02, 0.02), xycoords='axes fraction', - bbox=dict(boxstyle="round,pad=0.3", fc="white", ec="gray", alpha=0.8)) + stats = ( + f"Min: {masked_data.min():.3f} m\n" + f"Max: {masked_data.max():.3f} m\n" + f"Mean: {masked_data.mean():.3f} m" + ) + plt.annotate( + stats, + xy=(0.02, 0.02), + xycoords="axes fraction", + bbox=dict(boxstyle="round,pad=0.3", fc="white", ec="gray", alpha=0.8), + ) logger.info("Displaying preview... Close the plot window to continue.") plt.show() @@ -68,9 +73,14 @@ class GridEngine: def load_and_interpolate(source_files, target_region, nx, ny, decay_pixels=100): """Composites grids using GDAL Warper.""" - xmin, xmax, ymin, ymax = target_region.xmin, target_region.xmax, target_region.ymin, target_region.ymax + xmin, xmax, ymin, ymax = ( + target_region.xmin, + target_region.xmax, + target_region.ymin, + target_region.ymax, + ) dst_transform = from_bounds(xmin, ymin, xmax, ymax, nx, ny) - dst_crs = 'EPSG:4326' + dst_crs = "EPSG:4326" mosaic = np.full((ny, nx), np.nan, dtype=np.float32) @@ -85,7 +95,7 @@ def load_and_interpolate(source_files, target_region, nx, ny, decay_pixels=100): if src_nodata is not None: src_data[np.isclose(src_data, src_nodata, atol=1e-4)] = np.nan - if fn.endswith('.gtx'): + if fn.endswith(".gtx"): src_data[np.isclose(src_data, -88.8888, atol=1e-2)] = np.nan temp_buffer = np.full((ny, nx), np.nan, dtype=np.float32) @@ -95,12 +105,12 @@ def load_and_interpolate(source_files, target_region, nx, ny, decay_pixels=100): source=src_data, destination=temp_buffer, src_transform=src.transform, - src_crs=src.crs or 'EPSG:4326', + src_crs=src.crs or "EPSG:4326", src_nodata=np.nan, dst_transform=dst_transform, dst_crs=dst_crs, dst_nodata=np.nan, - resampling=Resampling.bilinear + resampling=Resampling.bilinear, ) valid_mask = ~np.isnan(temp_buffer) @@ -131,7 +141,9 @@ def smart_blend(in_grid, background_grid, blend_pixels=50): dist = ndimage.distance_transform_edt(mask) alpha = np.clip(dist / blend_pixels, 0.0, 1.0) - nearest_indices = ndimage.distance_transform_edt(mask, return_distances=False, return_indices=True) + nearest_indices = ndimage.distance_transform_edt( + mask, return_distances=False, return_indices=True + ) extended_vdatum = in_grid.copy() extended_vdatum[mask] = in_grid[tuple(nearest_indices)][mask] @@ -140,7 +152,13 @@ def smart_blend(in_grid, background_grid, blend_pixels=50): return blended_data @staticmethod - def coastal_aware_composite(vdatum_grid, global_grid, decay_pixels=100, buffer_pixels=10, max_discontinuity=0.5): + def coastal_aware_composite( + vdatum_grid, + global_grid, + decay_pixels=100, + buffer_pixels=10, + max_discontinuity=0.5, + ): """Intelligently handles inland decay vs. offshore blending, while filtering out low-resolution global artifacts. """ @@ -148,10 +166,14 @@ def coastal_aware_composite(vdatum_grid, global_grid, decay_pixels=100, buffer_p final_grid = vdatum_grid.copy() vdatum_mask = np.isnan(vdatum_grid) if not vdatum_mask.all(): - nearest_idx = ndimage.distance_transform_edt(vdatum_mask, return_distances=False, return_indices=True) + nearest_idx = ndimage.distance_transform_edt( + vdatum_mask, return_distances=False, return_indices=True + ) nearest_vdatum_vals = vdatum_grid[tuple(nearest_idx)] - fes_anomaly_mask = np.abs(global_grid - nearest_vdatum_vals) > max_discontinuity + fes_anomaly_mask = ( + np.abs(global_grid - nearest_vdatum_vals) > max_discontinuity + ) global_grid[fes_anomaly_mask] = np.nan @@ -162,11 +184,15 @@ def coastal_aware_composite(vdatum_grid, global_grid, decay_pixels=100, buffer_p is_offshore = ~is_vdatum & is_ocean if is_offshore.any(): - blended_ocean = GridEngine.smart_blend(vdatum_grid, global_grid, blend_pixels=50) + blended_ocean = GridEngine.smart_blend( + vdatum_grid, global_grid, blend_pixels=50 + ) final_grid[is_offshore] = blended_ocean[is_offshore] if is_inland.any(): - decayed_inland = GridEngine.fill_nans(vdatum_grid, decay_pixels=decay_pixels, buffer_pixels=buffer_pixels) + decayed_inland = GridEngine.fill_nans( + vdatum_grid, decay_pixels=decay_pixels, buffer_pixels=buffer_pixels + ) final_grid[is_inland] = decayed_inland[is_inland] return final_grid @@ -205,7 +231,9 @@ def apply_vertical_shift(src_dem, shift_array, dst_dem): data = src.read(1) if data.shape != shift_array.shape: - raise ValueError(f"Dimension mismatch: DEM {data.shape} vs Shift {shift_array.shape}") + raise ValueError( + f"Dimension mismatch: DEM {data.shape} vs Shift {shift_array.shape}" + ) nodata = src.nodata if src.nodata is not None else -9999 profile.update(nodata=nodata) @@ -214,7 +242,7 @@ def apply_vertical_shift(src_dem, shift_array, dst_dem): data[valid_mask] += shift_array[valid_mask] data[~valid_mask] = nodata - with rasterio.open(dst_dem, 'w', **profile) as dst: + with rasterio.open(dst_dem, "w", **profile) as dst: dst.write(data, 1) logger.info(f"Successfully wrote transformed DEM to: {dst_dem}") return True @@ -222,6 +250,7 @@ def apply_vertical_shift(src_dem, shift_array, dst_dem): logger.error(f"Failed to apply shift to DEM: {e}") return False + class GridWriter: @staticmethod def write(filename, data, region): @@ -231,8 +260,8 @@ def write(filename, data, region): if dirname and not os.path.exists(dirname): os.makedirs(dirname) - if not filename.endswith('.tif'): - filename = os.path.splitext(filename)[0] + '.tif' + if not filename.endswith(".tif"): + filename = os.path.splitext(filename)[0] + ".tif" rows, cols = data.shape xmin, xmax, ymin, ymax = region.xmin, region.xmax, region.ymin, region.ymax @@ -242,9 +271,17 @@ def write(filename, data, region): transform = rasterio.transform.from_origin(xmin, ymax, res_x, res_y) with rasterio.open( - filename, 'w', driver='GTiff', height=rows, width=cols, count=1, - dtype='float32', crs='EPSG:4326', transform=transform, - compress='deflate', tiled=True + filename, + "w", + driver="GTiff", + height=rows, + width=cols, + count=1, + dtype="float32", + crs="EPSG:4326", + transform=transform, + compress="deflate", + tiled=True, ) as dst: - dst.write(data.astype('float32'), 1) + dst.write(data.astype("float32"), 1) return filename diff --git a/src/transformez/hooks.py b/src/transformez/hooks.py index 3520ba8..e5e82dd 100644 --- a/src/transformez/hooks.py +++ b/src/transformez/hooks.py @@ -39,9 +39,16 @@ class TransformezHook(FetchHook): stage = "pre" desc = "Generate a vertical transformation shift grid." - def __init__(self, datum_in="5703", datum_out="6319", increment="3s", - output_grid="shift.tif", keep_grid=True, apply=False, - **kwargs): + def __init__( + self, + datum_in="5703", + datum_out="6319", + increment="3s", + output_grid="shift.tif", + keep_grid=True, + apply=False, + **kwargs, + ): super().__init__(**kwargs) self.datum_in = datum_in self.datum_out = datum_out @@ -57,7 +64,9 @@ def run(self, entries): module = entries[0][0] region = getattr(module, "region", None) if not region: - logger.warning("Module has no region defined. Cannot generate shift grid in PRE stage.") + logger.warning( + "Module has no region defined. Cannot generate shift grid in PRE stage." + ) return entries logger.info(f"Generating vertical shift grid for region: {region}") @@ -75,7 +84,9 @@ def _run_file(self, entries): """Apply the shift grid to specific files.""" if not os.path.exists(self.output_grid): - logger.warning(f"Shift grid {self.output_grid} not found. Skipping transform.") + logger.warning( + f"Shift grid {self.output_grid} not found. Skipping transform." + ) return entries for mod, entry in entries: @@ -114,10 +125,11 @@ def _generate_grid(self, region): vt = VerticalTransform( extent=region, - nx=nx, ny=ny, + nx=nx, + ny=ny, epsg_in=self.datum_in, epsg_out=self.datum_out, - verbose=False + verbose=False, ) shift_array, _ = vt._vertical_transform(vt.epsg_in, vt.epsg_out) @@ -125,9 +137,7 @@ def _generate_grid(self, region): if shift_array is None: logger.error("Transformation failed to generate a grid.") else: - GridWriter.write( - self.output_grid, shift_array, region.to_list() - ) + GridWriter.write(self.output_grid, shift_array, region.to_list()) logger.info(f"Saved shift grid to {self.output_grid}") def _apply_raster(self, src, grid): diff --git a/src/transformez/htdp.py b/src/transformez/htdp.py index 76c2527..bec68f5 100644 --- a/src/transformez/htdp.py +++ b/src/transformez/htdp.py @@ -29,13 +29,18 @@ ae = ".exe" if sys.platform == "win32" else "" # Validating HTDP existence -htdp_cmd = "echo 0 | htdp 2>&1" if sys.platform == "win32" else "echo 0 | htdp 2>&1 | grep SOFTWARE | awk '{print $3}'" -HAS_HTDP = utils.cmd_check(f'htdp{ae}', htdp_cmd) +htdp_cmd = ( + "echo 0 | htdp 2>&1" + if sys.platform == "win32" + else "echo 0 | htdp 2>&1 | grep SOFTWARE | awk '{print $3}'" +) +HAS_HTDP = utils.cmd_check(f"htdp{ae}", htdp_cmd) + class HTDP: """Wrapper for the NGS HTDP software.""" - def __init__(self, htdp_bin: str = 'htdp', verbose: bool = True): + def __init__(self, htdp_bin: str = "htdp", verbose: bool = True): self.htdp_bin = htdp_bin self.verbose = verbose if not HAS_HTDP: @@ -54,16 +59,14 @@ def run_grid(self, region, nx, ny, epsg_in, epsg_out, epoch_in, epoch_out): # transform.py passes ints (EPSGs), we need HTDP internal IDs def get_id(epsg): if epsg in Datums.HTDP: - return Datums.HTDP[epsg]['htdp_id'] + return Datums.HTDP[epsg]["htdp_id"] # Fallback for common codes if not in dictionary if epsg == 6319: return 1 # NAD83(2011) if epsg == 4979: - return 10 # WGS84(G1762) - raise ValueError( - f"EPSG {epsg} not defined in HTDP dictionary." - ) + return 10 # WGS84(G1762) + raise ValueError(f"EPSG {epsg} not defined in HTDP dictionary.") try: id_in = get_id(epsg_in) @@ -83,7 +86,7 @@ def get_id(epsg): lats = np.linspace(region.ymin, region.ymax, ny) # Write input file - with open(in_fn, 'w') as f: + with open(in_fn, "w") as f: for y_idx, lat in enumerate(lats): for x_idx, lon in enumerate(lons): # "Lat Lon Height TextID" @@ -92,9 +95,7 @@ def get_id(epsg): # Write Control File self._write_control( - ctl_fn, out_fn, in_fn, - id_in, epoch_in, - id_out, epoch_out + ctl_fn, out_fn, in_fn, id_in, epoch_in, id_out, epoch_out ) # Run HTDP @@ -113,13 +114,13 @@ def _read_grid(self, filename: str, shape: Tuple[int, int]) -> np.ndarray: """Parses HTDP output, mapping PNT_x_y tags to grid indices.""" grid = np.zeros(shape) - with open(filename, 'r') as f: + with open(filename, "r") as f: for line in f: if "PNT_" not in line: continue try: - parts = line.replace('"', ' ').split() + parts = line.replace('"', " ").split() # HTDP Output Format: Lat, Lon, Height, Text idx_off = 1 if parts[0] == "*" else 0 @@ -139,8 +140,9 @@ def _read_grid(self, filename: str, shape: Tuple[int, int]) -> np.ndarray: return grid - def _write_control(self, control_fn, out_fn, in_fn, - id_in, epoch_in, id_out, epoch_out): + def _write_control( + self, control_fn, out_fn, in_fn, id_in, epoch_in, id_out, epoch_out + ): """Writes the batch control file.""" # 4 = Transform Positions # 2 = Input file @@ -169,19 +171,19 @@ def _write_control(self, control_fn, out_fn, in_fn, f"0\n" f"0\n" ) - with open(control_fn, 'w') as f: + with open(control_fn, "w") as f: f.write(content) def run_cmd(self, control_fn): """Executes the binary.""" try: - with open(control_fn, 'r') as stdin: + with open(control_fn, "r") as stdin: subprocess.run( [self.htdp_bin], stdin=stdin, check=True, stdout=subprocess.DEVNULL, - stderr=subprocess.PIPE + stderr=subprocess.PIPE, ) except Exception as e: logger.error(f"HTDP Runtime Error: {e}") @@ -199,13 +201,11 @@ def download_htdp(target_dir=None): os.makedirs(target_dir, exist_ok=True) - if sys.platform == 'win32': + if sys.platform == "win32": # Windows: Download pre-compiled EXE url = "https://geodesy.noaa.gov/TOOLS/Htdp/htdp.exe" out_path = os.path.join(target_dir, "htdp.exe") - logger.info( - "Downloading HTDP executable for Windows from NOAA..." - ) + logger.info("Downloading HTDP executable for Windows from NOAA...") try: urllib.request.urlretrieve(url, out_path) logger.info(f"Success! Downloaded to: {out_path}") @@ -221,13 +221,11 @@ def download_htdp(target_dir=None): zip_path = os.path.join(target_dir, "HTDP-download.zip") src_dir = os.path.join(target_dir, "htdp_source") - logger.info( - "Downloading HTDP source code for Unix from NOAA..." - ) + logger.info("Downloading HTDP source code for Unix from NOAA...") try: urllib.request.urlretrieve(url, zip_path) - with zipfile.ZipFile(zip_path, 'r') as zip_ref: + with zipfile.ZipFile(zip_path, "r") as zip_ref: zip_ref.extractall(src_dir) os.remove(zip_path) @@ -236,11 +234,11 @@ def download_htdp(target_dir=None): logger.info( "HTDP requires compilation on Linux/macOS. Run the following commands:" ) - print("\n" + "="*50) + print("\n" + "=" * 50) print(f"cd {src_dir}") print("gfortran -o htdp htdp.f") print("sudo mv htdp /usr/local/bin/") - print("="*50 + "\n") + print("=" * 50 + "\n") except Exception as e: logger.error(f"Failed to download or extract HTDP source: {e}") diff --git a/src/transformez/modules.py b/src/transformez/modules.py index 615bb9f..65f9613 100644 --- a/src/transformez/modules.py +++ b/src/transformez/modules.py @@ -27,7 +27,7 @@ src_datum="Source Datum (e.g. 'mllw', '5703', '4979').", dst_datum="Destination Datum (e.g. '5703:geoid=g2012b').", increment="Grid resolution (default: 3s).", - output_name="Optional output filename override." + output_name="Optional output filename override.", ) class TransformezMod(FetchModule): """A dynamic Fetchez module that generates vertical shift grids on demand. @@ -37,38 +37,50 @@ class TransformezMod(FetchModule): """ name = "transformez" - meta_desc = 'Generate vertical datum shift grids on-demand.' + meta_desc = "Generate vertical datum shift grids on-demand." meta_category = "Tools" meta_tags = ["vdatum", "transformation", "shift-grid"] - def __init__(self, src_datum='5703', dst_datum='4979', increment='3s', output_name=None, **kwargs): + def __init__( + self, + src_datum="5703", + dst_datum="4979", + increment="3s", + output_name=None, + **kwargs, + ): super().__init__(name="transformez", **kwargs) self.src_datum = src_datum self.dst_datum = dst_datum self.increment = increment self.output_name = output_name - s_name = str(self.src_datum).replace(':', '_') - d_name = str(self.dst_datum).replace(':', '_') + s_name = str(self.src_datum).replace(":", "_") + d_name = str(self.dst_datum).replace(":", "_") w, e, s, n = self.region - self.dst_fn = os.path.join(self._outdir, f"shift_{s_name}_to_{d_name}_{w}_{s}.tif") + self.dst_fn = os.path.join( + self._outdir, f"shift_{s_name}_to_{d_name}_{w}_{s}.tif" + ) def run(self): from fetchez import utils + try: inc_val = utils.str2inc(self.increment) nx = int(self.region.width / inc_val) ny = int(self.region.height / inc_val) except Exception: - logger.warning(f"Invalid increment '{self.increment}', defaulting to 3s (~0.000833).") + logger.warning( + f"Invalid increment '{self.increment}', defaulting to 3s (~0.000833)." + ) # Default roughly 3 arc-seconds (approx 90m) nx = int(self.region.width / 0.00083333333) ny = int(self.region.height / 0.00083333333) def parse_d(d_str): - if ':' in str(d_str): - parts = d_str.split(':') - geoid = parts[1].split('=')[1] if 'geoid=' in parts[1] else parts[1] + if ":" in str(d_str): + parts = d_str.split(":") + geoid = parts[1].split("=")[1] if "geoid=" in parts[1] else parts[1] return parts[0], geoid return d_str, None @@ -77,9 +89,12 @@ def parse_d(d_str): vt = VerticalTransform( region=self.region, - nx=nx, ny=ny, - epsg_in=epsg_in, epsg_out=epsg_out, - geoid_in=geoid_in, geoid_out=geoid_out, + nx=nx, + ny=ny, + epsg_in=epsg_in, + epsg_out=epsg_out, + geoid_in=geoid_in, + geoid_out=geoid_out, ) logger.info(f"Generating shift grid: {self.src_datum} -> {self.dst_datum}...") @@ -98,6 +113,6 @@ def parse_d(d_str): meta={ "src_datum": self.src_datum, "dst_datum": self.dst_datum, - "generator": "transformez" - } + "generator": "transformez", + }, ) diff --git a/src/transformez/spatial.py b/src/transformez/spatial.py index cef863a..47bc8d3 100644 --- a/src/transformez/spatial.py +++ b/src/transformez/spatial.py @@ -26,7 +26,6 @@ class TransRegion(Region): - # gtl is (geo-transform, xcount, ycount) @classmethod def from_geo_transform(cls, gtl): @@ -41,7 +40,7 @@ def from_geo_transform(cls, gtl): return cls(xmin, xmax, ymin, ymax) - def srcwin(self, geo_transform, x_count, y_count, node='grid'): + def srcwin(self, geo_transform, x_count, y_count, node="grid"): """Output the appropriate GDAL srcwin (xoff, yoff, xsize, ysize).""" ul = _geo2pixel(self.xmin, self.ymax, geo_transform, node=node) @@ -60,11 +59,13 @@ def srcwin(self, geo_transform, x_count, y_count, node='grid'): y_size = y_stop - y_start if x_size <= 0 or y_size <= 0: - return 0, 0, 0, 0 + return 0, 0, 0, 0 return int(x_start), int(y_start), int(x_size), int(y_size) - def geo_transform(self, x_inc: float = 0, y_inc: Optional[float] = None, node: str = 'grid'): + def geo_transform( + self, x_inc: float = 0, y_inc: Optional[float] = None, node: str = "grid" + ): """Return dimensions and a geotransform based on the region and a cellsize. Returns: @@ -80,7 +81,11 @@ def geo_transform(self, x_inc: float = 0, y_inc: Optional[float] = None, node: s this_origin = _geo2pixel(self.xmin, self.ymax, dst_gt, node=node) this_end = _geo2pixel(self.xmax, self.ymin, dst_gt, node=node) - return int(this_end[0] - this_origin[0]), int(this_end[1] - this_origin[1]), dst_gt + return ( + int(this_end[0] - this_origin[0]), + int(this_end[1] - this_origin[1]), + dst_gt, + ) def geo_transform_from_count(self, x_count: int = 0, y_count: int = 0): x_inc = (self.xmax - self.xmin) / x_count @@ -133,16 +138,18 @@ def densify_edges(self, density=20): xs.extend(np.linspace(self.xmax, self.xmin, density)) ys.extend([self.ymin] * density) - #return list(zip(xs, ys)) + # return list(zip(xs, ys)) return xs, ys def transform_densify(self, transformer=None, transform_direction="FORWARD"): if transformer is None or not self.valid_p(): - logger.error(f'Could not perform region transformation; {self}') + logger.error(f"Could not perform region transformation; {self}") return self points_x, points_y = self.densify_edges(20) - trans_points_x, trans_points_y = transformer.transform(points_x, points_y, direction=transform_direction) + trans_points_x, trans_points_y = transformer.transform( + points_x, points_y, direction=transform_direction + ) self.xmin = min(trans_points_x) self.xmax = max(trans_points_x) @@ -150,35 +157,39 @@ def transform_densify(self, transformer=None, transform_direction="FORWARD"): self.ymax = max(trans_points_y) # set the new SRS - #self.src_srs = d_srs.ExportToWkt() + # self.src_srs = d_srs.ExportToWkt() return self def transform(self, transformer=None, transform_direction="FORWARD"): if transformer is None or not self.valid_p(): - logger.error(f'Could not perform region transformation; {self}') + logger.error(f"Could not perform region transformation; {self}") return self - self.xmin, self.ymin = transformer.transform(self.xmin, self.ymin, direction=transform_direction) - self.xmax, self.ymax = transformer.transform(self.xmax, self.ymax, direction=transform_direction) + self.xmin, self.ymin = transformer.transform( + self.xmin, self.ymin, direction=transform_direction + ) + self.xmax, self.ymax = transformer.transform( + self.xmax, self.ymax, direction=transform_direction + ) return self - def warp(self, dst_srs='epsg:4326'): + def warp(self, dst_srs="epsg:4326"): """Transform region horizontally to a new CRS.""" if str_or(self.srs) is None: - logger.warning(f'Region has no valid associated srs: {self.srs}') + logger.warning(f"Region has no valid associated srs: {self.srs}") return self with warnings.catch_warnings(): - warnings.simplefilter('ignore') + warnings.simplefilter("ignore") transform = srs.SRSParser(src_srs=self.srs, dst_srs=dst_srs) - if transform.tc['src_horz_crs'] and transform.tc['dst_horz_crs']: - pipeline_str = '+proj=pipeline +step {} +inv +step {}'.format( - transform.tc['src_horz_crs'].to_proj4(), - transform.tc['dst_horz_crs'].to_proj4() + if transform.tc["src_horz_crs"] and transform.tc["dst_horz_crs"]: + pipeline_str = "+proj=pipeline +step {} +inv +step {}".format( + transform.tc["src_horz_crs"].to_proj4(), + transform.tc["dst_horz_crs"].to_proj4(), ) transformer = pyproj.Transformer.from_pipeline(pipeline_str) @@ -189,28 +200,36 @@ def warp(self, dst_srs='epsg:4326'): return self -def _geo2pixel(geo_x, geo_y, geo_transform, node='grid'): +def _geo2pixel(geo_x, geo_y, geo_transform, node="grid"): """Convert a geographic x,y value to a pixel location.""" if geo_transform[2] + geo_transform[4] == 0: pixel_x = (geo_x - geo_transform[0]) / geo_transform[1] pixel_y = (geo_y - geo_transform[3]) / geo_transform[5] - if node == 'grid': - pixel_x += .5 - pixel_y += .5 + if node == "grid": + pixel_x += 0.5 + pixel_y += 0.5 else: pixel_x, pixel_y = _apply_gt(geo_x, geo_y, _invert_gt(geo_transform)) return int(pixel_x), int(pixel_y) -def _apply_gt(in_x, in_y, geo_transform, node='pixel'): +def _apply_gt(in_x, in_y, geo_transform, node="pixel"): """Apply geotransform to in_x, in_y.""" - offset = 0.5 if node == 'pixel' else 0 - - out_x = geo_transform[0] + (in_x + offset) * geo_transform[1] + (in_y + offset) * geo_transform[2] - out_y = geo_transform[3] + (in_x + offset) * geo_transform[4] + (in_y + offset) * geo_transform[5] + offset = 0.5 if node == "pixel" else 0 + + out_x = ( + geo_transform[0] + + (in_x + offset) * geo_transform[1] + + (in_y + offset) * geo_transform[2] + ) + out_y = ( + geo_transform[3] + + (in_x + offset) * geo_transform[4] + + (in_y + offset) * geo_transform[5] + ) return out_x, out_y @@ -227,8 +246,12 @@ def _invert_gt(geo_transform): out_gt[4] = -geo_transform[4] * inv_det out_gt[2] = -geo_transform[2] * inv_det out_gt[5] = geo_transform[1] * inv_det - out_gt[0] = (geo_transform[2] * geo_transform[3] - geo_transform[0] * geo_transform[5]) * inv_det - out_gt[3] = (-geo_transform[1] * geo_transform[3] + geo_transform[0] * geo_transform[4]) * inv_det + out_gt[0] = ( + geo_transform[2] * geo_transform[3] - geo_transform[0] * geo_transform[5] + ) * inv_det + out_gt[3] = ( + -geo_transform[1] * geo_transform[3] + geo_transform[0] * geo_transform[4] + ) * inv_det return out_gt diff --git a/src/transformez/srs.py b/src/transformez/srs.py index c8ad73d..ca67a8d 100644 --- a/src/transformez/srs.py +++ b/src/transformez/srs.py @@ -31,7 +31,9 @@ class SRSParser: - Vertical: Z + Shift_Grid """ - def __init__(self, src_srs, dst_srs, region=None, vert_grid=None, cache_dir='.', **kwargs): + def __init__( + self, src_srs, dst_srs, region=None, vert_grid=None, cache_dir=".", **kwargs + ): self.src_srs_input = src_srs self.dst_srs_input = dst_srs self.region = region @@ -39,24 +41,24 @@ def __init__(self, src_srs, dst_srs, region=None, vert_grid=None, cache_dir='.', self.cache_dir = cache_dir self.tc = { - 'src_crs': None, - 'dst_crs': None, - 'src_vert_epsg': None, - 'dst_vert_epsg': None, - 'src_geoid': None, - 'dst_geoid': None, - 'want_vertical': False, - 'trans_fn': None, + "src_crs": None, + "dst_crs": None, + "src_vert_epsg": None, + "dst_vert_epsg": None, + "src_geoid": None, + "dst_geoid": None, + "want_vertical": False, + "trans_fn": None, } self._parse_srs() def _extract_geoid(self, srs_str): - parts = str(srs_str).split('+geoid:') + parts = str(srs_str).split("+geoid:") return parts[0], (parts[1] if len(parts) > 1 else None) def _extract_vertical(self, srs_str): - parts = str(srs_str).split('+') + parts = str(srs_str).split("+") return parts[0], (parts[1] if len(parts) > 1 else None) def _get_epsg_int(self, crs): @@ -68,12 +70,12 @@ def _get_epsg_int(self, crs): return None def _parse_srs(self): - clean_src, self.tc['src_geoid'] = self._extract_geoid(self.src_srs_input) - clean_dst, self.tc['dst_geoid'] = self._extract_geoid(self.dst_srs_input) + clean_src, self.tc["src_geoid"] = self._extract_geoid(self.src_srs_input) + clean_dst, self.tc["dst_geoid"] = self._extract_geoid(self.dst_srs_input) try: - self.tc['src_crs'] = CRS.from_user_input(clean_src) - self.tc['dst_crs'] = CRS.from_user_input(clean_dst) + self.tc["src_crs"] = CRS.from_user_input(clean_src) + self.tc["dst_crs"] = CRS.from_user_input(clean_dst) vert_epsg_src = None vert_epsg_dst = None except Exception: @@ -81,50 +83,57 @@ def _parse_srs(self): clean_dst, vert_epsg_dst = self._extract_vertical(self.dst_srs_input) try: - self.tc['src_crs'] = CRS.from_user_input(clean_src) - self.tc['dst_crs'] = CRS.from_user_input(clean_dst) + self.tc["src_crs"] = CRS.from_user_input(clean_src) + self.tc["dst_crs"] = CRS.from_user_input(clean_dst) except Exception as e: logger.error(f"Invalid SRS: {e}") return # Extract vertical components before flattening - if self.tc['src_crs'].is_compound: - self.tc['src_vert_epsg'] = self._get_epsg_int(self.tc['src_crs'].sub_crs_list[1]) + if self.tc["src_crs"].is_compound: + self.tc["src_vert_epsg"] = self._get_epsg_int( + self.tc["src_crs"].sub_crs_list[1] + ) # Strip to Horizontal for PROJ Transformer - self.tc['src_crs'] = self.tc['src_crs'].sub_crs_list[0] - elif self.tc['src_crs'].is_vertical: - self.tc['src_vert_epsg'] = self._get_epsg_int(self.tc['src_crs']) + self.tc["src_crs"] = self.tc["src_crs"].sub_crs_list[0] + elif self.tc["src_crs"].is_vertical: + self.tc["src_vert_epsg"] = self._get_epsg_int(self.tc["src_crs"]) - if self.tc['dst_crs'].is_compound: - self.tc['dst_vert_epsg'] = self._get_epsg_int(self.tc['dst_crs'].sub_crs_list[1]) - self.tc['dst_crs'] = self.tc['dst_crs'].sub_crs_list[0] - elif self.tc['dst_crs'].is_vertical: - self.tc['dst_vert_epsg'] = self._get_epsg_int(self.tc['dst_crs']) + if self.tc["dst_crs"].is_compound: + self.tc["dst_vert_epsg"] = self._get_epsg_int( + self.tc["dst_crs"].sub_crs_list[1] + ) + self.tc["dst_crs"] = self.tc["dst_crs"].sub_crs_list[0] + elif self.tc["dst_crs"].is_vertical: + self.tc["dst_vert_epsg"] = self._get_epsg_int(self.tc["dst_crs"]) - if self.tc['src_vert_epsg'] is None: - self.tc['src_vert_epsg'] = vert_epsg_src - if self.tc['dst_vert_epsg'] is None: - self.tc['dst_vert_epsg'] = vert_epsg_dst + if self.tc["src_vert_epsg"] is None: + self.tc["src_vert_epsg"] = vert_epsg_src + if self.tc["dst_vert_epsg"] is None: + self.tc["dst_vert_epsg"] = vert_epsg_dst # Lookup default geoids # If we have a vertical EPSG but no manual geoid, look it up in definitions.py - if self.tc['src_vert_epsg'] and not self.tc['src_geoid']: - self.tc['src_geoid'] = Datums.get_default_geoid(self.tc['src_vert_epsg']) + if self.tc["src_vert_epsg"] and not self.tc["src_geoid"]: + self.tc["src_geoid"] = Datums.get_default_geoid(self.tc["src_vert_epsg"]) - if self.tc['dst_vert_epsg'] and not self.tc['dst_geoid']: - self.tc['dst_geoid'] = Datums.get_default_geoid(self.tc['dst_vert_epsg']) + if self.tc["dst_vert_epsg"] and not self.tc["dst_geoid"]: + self.tc["dst_geoid"] = Datums.get_default_geoid(self.tc["dst_vert_epsg"]) # We want vertical if we have explicit vertical EPSGs OR manual geoids - has_src_vert = (self.tc['src_vert_epsg'] is not None) or (self.tc['src_geoid'] is not None) - has_dst_vert = (self.tc['dst_vert_epsg'] is not None) or (self.tc['dst_geoid'] is not None) - - self.tc['want_vertical'] = has_src_vert or has_dst_vert + has_src_vert = (self.tc["src_vert_epsg"] is not None) or ( + self.tc["src_geoid"] is not None + ) + has_dst_vert = (self.tc["dst_vert_epsg"] is not None) or ( + self.tc["dst_geoid"] is not None + ) + self.tc["want_vertical"] = has_src_vert or has_dst_vert def set_vertical_transform(self): """Generates the vertical shift grid using VerticalTransform.""" - if not self.region or not self.tc['want_vertical']: + if not self.region or not self.tc["want_vertical"]: return try: @@ -134,39 +143,40 @@ def set_vertical_transform(self): proc_region = Region.from_list(self.region) proc_region.buffer(pct=5) - s_ident = self.tc['src_vert_epsg'] - d_ident = self.tc['dst_vert_epsg'] + s_ident = self.tc["src_vert_epsg"] + d_ident = self.tc["dst_vert_epsg"] - if not s_ident and self.tc['src_geoid']: + if not s_ident and self.tc["src_geoid"]: s_ident = 6319 - if not d_ident and self.tc['dst_geoid']: + if not d_ident and self.tc["dst_geoid"]: d_ident = 6319 if not s_ident or not d_ident: return - s_name = str(s_ident).replace(':', '_').replace(' ', '_').replace('/', '_') - d_name = str(d_ident).replace(':', '_').replace(' ', '_').replace('/', '_') + s_name = str(s_ident).replace(":", "_").replace(" ", "_").replace("/", "_") + d_name = str(d_ident).replace(":", "_").replace(" ", "_").replace("/", "_") grid_name = f"_vdatum_trans_{s_name}_{d_name}_{proc_region.format('fn')}.tif" - self.tc['trans_fn'] = grid_name.replace('\\', '/') + self.tc["trans_fn"] = grid_name.replace("\\", "/") - if not os.path.exists(self.tc['trans_fn']): - logger.info(f"Generating vertical grid: {s_ident} -> {d_ident} : {self.tc['trans_fn']} :") + if not os.path.exists(self.tc["trans_fn"]): + logger.info( + f"Generating vertical grid: {s_ident} -> {d_ident} : {self.tc['trans_fn']} :" + ) vt = VerticalTransform( proc_region, nx=max(10, int(proc_region.width / 0.0008333)), ny=max(10, int(proc_region.height / 0.0008333)), epsg_in=s_ident, epsg_out=d_ident, - geoid_in=self.tc['src_geoid'], - geoid_out=self.tc['dst_geoid'] + geoid_in=self.tc["src_geoid"], + geoid_out=self.tc["dst_geoid"], ) shift_arr, _ = vt._vertical_transform(vt.epsg_in, vt.epsg_out) - GridWriter.write(self.tc['trans_fn'], shift_arr, proc_region) - - self.manual_vert_grid = self.tc['trans_fn'] + GridWriter.write(self.tc["trans_fn"], shift_arr, proc_region) + self.manual_vert_grid = self.tc["trans_fn"] def get_components(self): """Returns the components: @@ -176,7 +186,7 @@ def get_components(self): - Grid Path (str) or None """ - if self.tc['want_vertical'] and not self.manual_vert_grid: + if self.tc["want_vertical"] and not self.manual_vert_grid: self.set_vertical_transform() # Define Hub: NAD83(2011) 2D - maybe this should be wgs(transit) @@ -184,7 +194,7 @@ def get_components(self): # only build 2D transformer (for proj). We'll apply the vertical grid ourselves, # as pyproj seems finicky when it comes to this. - t_to_hub = Transformer.from_crs(self.tc['src_crs'], hub_crs, always_xy=True) - t_from_hub = Transformer.from_crs(hub_crs, self.tc['dst_crs'], always_xy=True) + t_to_hub = Transformer.from_crs(self.tc["src_crs"], hub_crs, always_xy=True) + t_from_hub = Transformer.from_crs(hub_crs, self.tc["dst_crs"], always_xy=True) return t_to_hub, t_from_hub, self.manual_vert_grid diff --git a/src/transformez/transform.py b/src/transformez/transform.py index 5970b7f..e443585 100644 --- a/src/transformez/transform.py +++ b/src/transformez/transform.py @@ -28,12 +28,25 @@ WGS84_EPSG = 4979 NAD83_EPSG = 6319 + class VerticalTransform: """Generate a vertical transformation grid using Transformez.""" - def __init__(self, region, nx, ny, epsg_in, epsg_out, - geoid_in=None, geoid_out=None, epoch_in=2010.0, epoch_out=2010.0, - decay_pixels=100, cache_dir=None, verbose=True): + def __init__( + self, + region, + nx, + ny, + epsg_in, + epsg_out, + geoid_in=None, + geoid_out=None, + epoch_in=2010.0, + epoch_out=2010.0, + decay_pixels=100, + cache_dir=None, + verbose=True, + ): self.region = region self.nx = nx @@ -77,15 +90,15 @@ def _get_native_ellipsoid(self, epsg, ref_type): if ref_type in ["surface", "global_tidal"]: # NOAA VDatum = NAD83, Global = WGS84 - region = Datums.SURFACES[epsg].get('region') - return NAD83_EPSG if region == 'usa' else WGS84_EPSG - elif ref_type == 'cdn': + region = Datums.SURFACES[epsg].get("region") + return NAD83_EPSG if region == "usa" else WGS84_EPSG + elif ref_type == "cdn": # Look up in definitions, default to NAD83 for US geoids - return Datums.CDN.get(epsg, {}).get('ellipsoid', NAD83_EPSG) - elif ref_type == 'htdp': + return Datums.CDN.get(epsg, {}).get("ellipsoid", NAD83_EPSG) + elif ref_type == "htdp": # If it's a Frame, it is its own native ellipsoid return epsg - return WGS84_EPSG # Default + return WGS84_EPSG # Default def fetch_grid(self, module_name, **kwargs): """Generic fetcher wrapper using the new fetchez API.""" @@ -95,7 +108,7 @@ def fetch_grid(self, module_name, **kwargs): region=self.region.to_list(), outdir=self.cache_dir, threads=2, - **kwargs + **kwargs, ) valid = [] @@ -104,25 +117,33 @@ def fetch_grid(self, module_name, **kwargs): if not os.path.exists(fn): continue - if fn.endswith('.zip'): - datatype = kwargs.get('datatype') + if fn.endswith(".zip"): + datatype = kwargs.get("datatype") fns_to_extract = [datatype] if datatype else None - extracted = fetchez.utils.p_f_unzip(fn, fns=fns_to_extract, outdir=self.cache_dir) - valid.extend([f for f in extracted if f.endswith(('.gtx', '.tif', '.grd', '.nc'))]) - - elif fn.endswith('.gz'): + extracted = fetchez.utils.p_f_unzip( + fn, fns=fns_to_extract, outdir=self.cache_dir + ) + valid.extend( + [ + f + for f in extracted + if f.endswith((".gtx", ".tif", ".grd", ".nc")) + ] + ) + + elif fn.endswith(".gz"): try: out_fn = os.path.splitext(fn)[0] if not os.path.exists(out_fn): logger.info(f"Decompressing {fn}...") - with gzip.open(fn, 'rb') as f_in: - with open(out_fn, 'wb') as f_out: + with gzip.open(fn, "rb") as f_in: + with open(out_fn, "wb") as f_out: shutil.copyfileobj(f_in, f_out) valid.append(out_fn) except Exception as e: logger.error(f"Failed to decompress {fn}: {e}") - elif fn.endswith(('.gtx', '.tif', '.grd', '.nc', '.mss')): + elif fn.endswith((".gtx", ".tif", ".grd", ".nc", ".mss")): valid.append(fn) return valid @@ -140,26 +161,34 @@ def fetch_grid_(self, module_name, **kwargs): fetchez.core.run_fetchez([fetcher], threads=2) valid = [] for r in fetcher.results: - fn = r['dst_fn'] + fn = r["dst_fn"] if not os.path.exists(fn): continue - if fn.endswith('.zip'): - extracted = fetchez.utils.p_f_unzip(fn, fns=[r['data_type']], outdir=self.cache_dir) - valid.extend([f for f in extracted if f.endswith(('.gtx', '.tif', '.grd', '.nc'))]) - elif fn.endswith('.gz'): + if fn.endswith(".zip"): + extracted = fetchez.utils.p_f_unzip( + fn, fns=[r["data_type"]], outdir=self.cache_dir + ) + valid.extend( + [ + f + for f in extracted + if f.endswith((".gtx", ".tif", ".grd", ".nc")) + ] + ) + elif fn.endswith(".gz"): try: out_fn = os.path.splitext(fn)[0] if not os.path.exists(out_fn): logger.info(f"Decompressing {fn}...") - with gzip.open(fn, 'rb') as f_in: - with open(out_fn, 'wb') as f_out: + with gzip.open(fn, "rb") as f_in: + with open(out_fn, "wb") as f_out: shutil.copyfileobj(f_in, f_out) valid.append(out_fn) except Exception as e: logger.error(f"Failed to decompress {fn}: {e}") - elif fn.endswith(('.gtx', '.tif', '.grd', '.nc', '.mss')): + elif fn.endswith((".gtx", ".tif", ".grd", ".nc", ".mss")): valid.append(fn) return valid @@ -168,21 +197,25 @@ def _get_grid(self, provider, name): return np.zeros((self.ny, self.nx)) if not provider: - provider = 'proj' + provider = "proj" - if 'geoid=' in name: - name = name.split('=')[1] + if "geoid=" in name: + name = name.split("=")[1] files = self.fetch_grid(provider, datatype=name, query=name) if not files: return np.zeros((self.ny, self.nx)) - if provider == 'seanoe' or provider == 'fes': + if provider == "seanoe" or provider == "fes": var_name = "lat_elevation" if "lat" in name.lower() else "msl_elevation" nc_path = f"netcdf:{files[0]}:{var_name}" - return GridEngine.load_and_interpolate([nc_path], self.region, self.nx, self.ny, decay_pixels=self.decay_pixels) + return GridEngine.load_and_interpolate( + [nc_path], self.region, self.nx, self.ny, decay_pixels=self.decay_pixels + ) - return GridEngine.load_and_interpolate(files, self.region, self.nx, self.ny, decay_pixels=self.decay_pixels) + return GridEngine.load_and_interpolate( + files, self.region, self.nx, self.ny, decay_pixels=self.decay_pixels + ) def _get_htdp_shift(self, epsg_from, epsg_to, epoch_from, epoch_to): """Calculate Frame Shift via HTDP.""" @@ -191,10 +224,13 @@ def _get_htdp_shift(self, epsg_from, epsg_to, epoch_from, epoch_to): return np.zeros((self.ny, self.nx)) from . import htdp + try: logger.info(f" [HTDP] Frame Shift: EPSG:{epsg_from} -> EPSG:{epsg_to}") tool = htdp.HTDP(verbose=False) - grid = tool.run_grid(self.region, self.nx, self.ny, epsg_from, epsg_to, epoch_from, epoch_to) + grid = tool.run_grid( + self.region, self.nx, self.ny, epsg_from, epsg_to, epoch_from, epoch_to + ) if np.any(grid): logger.info(f" [HTDP] Component Shift (Mean: {np.mean(grid):.3f}m)") return grid @@ -209,7 +245,7 @@ def _fetch_geoid_with_fallback(self, target_geoid): """ # Ordered list of preferred US geoids (Newest to Oldest) - us_geoids = ['g2018', 'g2012b', 'geoid09'] + us_geoids = ["g2018", "g2012b", "geoid09"] if target_geoid in us_geoids: start_idx = us_geoids.index(target_geoid) @@ -219,12 +255,14 @@ def _fetch_geoid_with_fallback(self, target_geoid): for g in geoids_to_try: geoid_def = Datums.GEOIDS.get(g, {}) - provider = geoid_def.get('provider', 'proj') + provider = geoid_def.get("provider", "proj") grid = self._get_grid(provider, g) if np.any(grid): if g != target_geoid and self.verbose: - logger.info(f" [Geoid Fallback] '{target_geoid}' lacks coverage here. Falling back to '{g}'.") + logger.info( + f" [Geoid Fallback] '{target_geoid}' lacks coverage here. Falling back to '{g}'." + ) return grid, g return np.zeros((self.ny, self.nx)), target_geoid @@ -238,15 +276,15 @@ def _get_vdatum_chain(self, datum_name, geoid_name): desc = [] # Tidal -> LMSL - if datum_name not in ['msl', '5714', 'lmsl']: - grid = self._get_grid('vdatum', datum_name) + if datum_name not in ["msl", "5714", "lmsl"]: + grid = self._get_grid("vdatum", datum_name) if np.isnan(grid).all() or (grid == 0).all(): return None, f"Missing Tidal Grid: {datum_name}" hydro_shift += grid desc.append(f"({datum_name}->LMSL)") # LMSL -> Ortho (TSS) - tss = self._get_grid('vdatum', 'tss') + tss = self._get_grid("vdatum", "tss") if np.isnan(tss).all() or (tss == 0).all(): return None, "Outside VDatum coverage (Missing TSS)" @@ -255,7 +293,7 @@ def _get_vdatum_chain(self, datum_name, geoid_name): # Ortho -> NAD83 (Geoid) # We fetch the geoid, but DO NOT add it to the shift yet! - actual_geoid = geoid_name if geoid_name else 'g2018' + actual_geoid = geoid_name if geoid_name else "g2018" geoid_grid, used_geoid = self._fetch_geoid_with_fallback(actual_geoid) desc.append(f"Geoid({used_geoid}->NAD83)") @@ -267,12 +305,18 @@ def _get_vdatum_chain(self, datum_name, geoid_name): if np.isnan(hydro_shift).any(): proxy_name = Datums.get_global_proxy(datum_name) if proxy_name: - logger.info(f"Partial VDatum coverage detected. Fetching {proxy_name.upper()} (FES) for offshore blending...") - global_shift, d_global = self._get_global_chain(proxy_name, model='fes2014') + logger.info( + f"Partial VDatum coverage detected. Fetching {proxy_name.upper()} (FES) for offshore blending..." + ) + global_shift, d_global = self._get_global_chain( + proxy_name, model="fes2014" + ) if global_shift is not None and np.any(global_shift): # We have valid FES data. We must align it to NAD83. - htdp_wgs_to_nad = self._get_htdp_shift(WGS84_EPSG, NAD83_EPSG, self.epoch_in, 2010.0) + htdp_wgs_to_nad = self._get_htdp_shift( + WGS84_EPSG, NAD83_EPSG, self.epoch_in, 2010.0 + ) fes_full = global_shift + htdp_wgs_to_nad hydro_shift = GridEngine.coastal_aware_composite( @@ -280,14 +324,18 @@ def _get_vdatum_chain(self, datum_name, geoid_name): global_grid=fes_full, decay_pixels=self.decay_pixels, buffer_pixels=10, - max_discontinuity=0.5 + max_discontinuity=0.5, ) desc.append(f"Blended w/ Global({proxy_name.upper()})") else: - hydro_shift = GridEngine.fill_nans(hydro_shift, decay_pixels=self.decay_pixels, buffer_pixels=10) + hydro_shift = GridEngine.fill_nans( + hydro_shift, decay_pixels=self.decay_pixels, buffer_pixels=10 + ) desc.append("Inland Hydro Decay") else: - hydro_shift = GridEngine.fill_nans(hydro_shift, decay_pixels=self.decay_pixels, buffer_pixels=10) + hydro_shift = GridEngine.fill_nans( + hydro_shift, decay_pixels=self.decay_pixels, buffer_pixels=10 + ) desc.append("Inland Hydro Decay") total_shift = hydro_shift + geoid_grid @@ -325,7 +373,7 @@ def _get_global_chain(self, datum_name, model="fes2014"): mss_name = model_def["grids"].get("mss") if mss_name: mss_grid = self._get_grid(provider, mss_name) - if provider == 'seanoe' or 'fes' in model.lower(): + if provider == "seanoe" or "fes" in model.lower(): mss_grid -= 0.70 desc.append("MSS->WGS84(TP_Corr)") else: @@ -351,34 +399,36 @@ def _step_to_hub(self, epsg, ref_type, geoid=None, epoch=None): chain_desc = "" if ref_type in ["surface", "global_tidal"]: - datum_name = Datums.SURFACES[epsg]['name'] - region_tag = Datums.SURFACES[epsg].get('region') + datum_name = Datums.SURFACES[epsg]["name"] + region_tag = Datums.SURFACES[epsg].get("region") - if region_tag == 'usa': + if region_tag == "usa": s, d = self._get_vdatum_chain(datum_name, geoid) if s is None: native_epsg = WGS84_EPSG proxy_name = Datums.get_global_proxy(datum_name) if proxy_name: - s, d = self._get_global_chain(proxy_name, model='fes2014') + s, d = self._get_global_chain(proxy_name, model="fes2014") # d = f"Global({proxy_name}) [Proxy] -> WGS84" chain_shift, chain_desc = s, d - elif region_tag == 'global': - chain_shift, chain_desc = self._get_global_chain(datum_name) + elif region_tag == "global": + chain_shift, chain_desc = self._get_global_chain(datum_name) - elif ref_type == 'cdn': - target_geoid = geoid if geoid else 'g2018' + elif ref_type == "cdn": + target_geoid = geoid if geoid else "g2018" chain_shift, used_geoid = self._fetch_geoid_with_fallback(target_geoid) chain_desc = f"Ortho(via {used_geoid}) -> Frame({native_epsg})" - elif ref_type == 'htdp': + elif ref_type == "htdp": chain_shift = np.zeros((self.ny, self.nx)) chain_desc = f"Frame({epsg})" if chain_shift is not None: if native_epsg != self.hub_epsg: - htdp_shift = self._get_htdp_shift(native_epsg, self.hub_epsg, epoch, self.epoch_out) + htdp_shift = self._get_htdp_shift( + native_epsg, self.hub_epsg, epoch, self.epoch_out + ) chain_shift += htdp_shift chain_desc += f" + Frame({native_epsg}->{self.hub_epsg})" return chain_shift, chain_desc @@ -397,16 +447,18 @@ def _step_from_hub(self, epsg, ref_type, geoid=None, epoch=None): htdp_shift = np.zeros((self.ny, self.nx)) if self.hub_epsg != native_epsg: - htdp_shift = self._get_htdp_shift(self.hub_epsg, native_epsg, self.epoch_in, epoch) - total_out += htdp_shift - desc_parts.append(f"Hub({self.hub_epsg}->{native_epsg})") + htdp_shift = self._get_htdp_shift( + self.hub_epsg, native_epsg, self.epoch_in, epoch + ) + total_out += htdp_shift + desc_parts.append(f"Hub({self.hub_epsg}->{native_epsg})") if ref_type in ["surface", "global_tidal"]: - datum_name = Datums.SURFACES[epsg]['name'] - region_tag = Datums.SURFACES[epsg].get('region') - chain_geoid = geoid if geoid else 'g2018' + datum_name = Datums.SURFACES[epsg]["name"] + region_tag = Datums.SURFACES[epsg].get("region") + chain_geoid = geoid if geoid else "g2018" - if region_tag == 'usa': + if region_tag == "usa": s, d = self._get_vdatum_chain(datum_name, chain_geoid) if s is None: proxy_name = Datums.get_global_proxy(datum_name) @@ -416,30 +468,34 @@ def _step_from_hub(self, epsg, ref_type, geoid=None, epoch=None): if desc_parts: desc_parts.pop() - s, d = self._get_global_chain(proxy_name, model='fes2014') + s, d = self._get_global_chain(proxy_name, model="fes2014") if s is not None: total_out -= s desc_parts.append(f"GlobalProxy({proxy_name})") else: - return np.zeros((self.ny, self.nx)), "FAILED Output Global Chain" + return np.zeros( + (self.ny, self.nx) + ), "FAILED Output Global Chain" else: return np.zeros((self.ny, self.nx)), "FAILED Output Chain" else: total_out -= s desc_parts.append(f"Native -> VDatum({datum_name})") - elif region_tag == 'global': + elif region_tag == "global": s, d = self._get_global_chain(datum_name) if s is not None: total_out -= s desc_parts.append(f"Native -> Global({datum_name})") - elif ref_type == 'cdn': - target_geoid = geoid if geoid else 'g2018' + elif ref_type == "cdn": + target_geoid = geoid if geoid else "g2018" geoid_grid, used_geoid = self._fetch_geoid_with_fallback(target_geoid) if not np.any(geoid_grid): - logger.warning(f"Geoid {target_geoid} (and fallbacks) not found/covered.") + logger.warning( + f"Geoid {target_geoid} (and fallbacks) not found/covered." + ) total_out -= geoid_grid desc_parts.append(f"Native -> Ortho(via {used_geoid})") @@ -459,7 +515,9 @@ def _vertical_transform(self, epsg_in, epsg_out): return total_shift, total_unc # Input -> Hub - grid_1, desc_1 = self._step_to_hub(self.epsg_in, self.ref_in, self.geoid_in, self.epoch_in) + grid_1, desc_1 = self._step_to_hub( + self.epsg_in, self.ref_in, self.geoid_in, self.epoch_in + ) if np.any(grid_1): logger.info(f" 1. {desc_1}") total_shift += grid_1 @@ -467,7 +525,9 @@ def _vertical_transform(self, epsg_in, epsg_out): logger.info(f" 1. {desc_1} (No Shift/Zero)") # Hub -> Output - grid_2, desc_2 = self._step_from_hub(self.epsg_out, self.ref_out, self.geoid_out, self.epoch_out) + grid_2, desc_2 = self._step_from_hub( + self.epsg_out, self.ref_out, self.geoid_out, self.epoch_out + ) if np.any(grid_2): logger.info(f" 2. {desc_2}") total_shift += grid_2 @@ -478,7 +538,9 @@ def _vertical_transform(self, epsg_in, epsg_out): mean_shift = np.nanmean(total_shift) min_shift = np.nanmin(total_shift) max_shift = np.nanmax(total_shift) - logger.info(f" => Total Shift Applied (Mean: {mean_shift:.3f}m | Min: {min_shift:.3f}m | Max: {max_shift:.3f}m)") + logger.info( + f" => Total Shift Applied (Mean: {mean_shift:.3f}m | Min: {min_shift:.3f}m | Max: {max_shift:.3f}m)" + ) else: logger.info(" => Total Shift Applied (Zero / Identity)") diff --git a/src/transformez/utils.py b/src/transformez/utils.py index 6d46a4b..d7bc197 100644 --- a/src/transformez/utils.py +++ b/src/transformez/utils.py @@ -18,8 +18,10 @@ import rasterio logger = logging.getLogger(__name__) -cmd_exists = lambda x: any(os.access(os.path.join(path, x), os.X_OK) - for path in os.environ['PATH'].split(os.pathsep)) +cmd_exists = lambda x: any( + os.access(os.path.join(path, x), os.X_OK) + for path in os.environ["PATH"].split(os.pathsep) +) def run_cmd(args): @@ -31,7 +33,7 @@ def run_cmd(args): args, shell=False if isinstance(args, list) else True, capture_output=True, - text=True + text=True, ) return result.stdout, result.returncode @@ -84,10 +86,7 @@ def query(self, x, y): cols = np.floor(cols_f).astype(int) rows = np.floor(rows_f).astype(int) - valid = ( - (rows >= 0) & (rows < self.height) & - (cols >= 0) & (cols < self.width) - ) + valid = (rows >= 0) & (rows < self.height) & (cols >= 0) & (cols < self.width) results = np.full_like(q_x, self.default_nodata, dtype=self.data.dtype) if np.any(valid): diff --git a/src/transformez/vdatum.py b/src/transformez/vdatum.py index cd90b6a..a1d75fe 100644 --- a/src/transformez/vdatum.py +++ b/src/transformez/vdatum.py @@ -23,10 +23,21 @@ class Vdatum: - def __init__(self, jar=None, ivert="navd88:m:height", overt="mhw:m:height", - ihorz="NAD83_2011", ohorz="NAD83_2011", region="4", fmt="txt", - xyzl="0,1,2", skip=0, delim="space", result_dir="result", - verbose=False): + def __init__( + self, + jar=None, + ivert="navd88:m:height", + overt="mhw:m:height", + ihorz="NAD83_2011", + ohorz="NAD83_2011", + region="4", + fmt="txt", + xyzl="0,1,2", + skip=0, + delim="space", + result_dir="result", + verbose=False, + ): self.jar = jar self.ivert = ivert self.overt = overt @@ -69,8 +80,8 @@ def vdatum_get_version(self): if self.jar is not None: out, _ = utils.run_cmd(f"java -jar {self.jar} -", verbose=self.verbose) for i in out.decode("utf-8").split("\n"): - if '- v' in i.strip(): - return i.strip().split('v')[-1] + if "- v" in i.strip(): + return i.strip().split("v")[-1] return None def vdatum_xyz(self, xyz): @@ -80,12 +91,13 @@ def vdatum_xyz(self, xyz): self.vdatum_locate_jar() if self.jar is not None: epoch_str = f"epoch:{self.epoch} " if self.epoch is not None else "" - vdc = (f"ihorz:{self.ihorz} ivert:{self.ivert} ohorz:{self.ohorz} overt:{self.overt} " - f"-nodata -pt:{xyz[0]},{xyz[1]},{xyz[2]} {epoch_str}region:{self.region}") + vdc = ( + f"ihorz:{self.ihorz} ivert:{self.ivert} ohorz:{self.ohorz} overt:{self.overt} " + f"-nodata -pt:{xyz[0]},{xyz[1]},{xyz[2]} {epoch_str}region:{self.region}" + ) out, _ = utils.run_cmd( - f"java -Djava.awt.headless=false -jar {self.jar} {vdc}", - verbose=False + f"java -Djava.awt.headless=false -jar {self.jar} {vdc}", verbose=False ) z = xyz[2] for i in out.split("\n"): @@ -115,9 +127,11 @@ def run_vdatum(self, src_fn): self.vdatum_locate_jar() if self.jar is not None: epoch_str = f"epoch:{self.epoch} " if self.epoch is not None else "" - vdc = (f"ihorz:{self.ihorz} ivert:{self.ivert} ohorz:{self.ohorz} overt:{self.overt} " - f"-nodata -file:txt:{self.delim},{self.xyzl},skip{self.skip}:{src_fn}:{self.result_dir} " - f"{epoch_str}region:{self.region}") + vdc = ( + f"ihorz:{self.ihorz} ivert:{self.ivert} ohorz:{self.ohorz} overt:{self.overt} " + f"-nodata -file:txt:{self.delim},{self.xyzl},skip{self.skip}:{src_fn}:{self.result_dir} " + f"{epoch_str}region:{self.region}" + ) return utils.run_cmd(f"java -jar {self.jar} {vdc}", verbose=self.verbose) else: return [], -1