Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions src/transformez/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -61,6 +64,7 @@ def _find_proj_lib():

return None


target_proj_lib = _find_proj_lib()

if "PROJ_LIB" in os.environ:
Expand Down
57 changes: 32 additions & 25 deletions src/transformez/api.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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:
Expand All @@ -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()

Expand Down Expand Up @@ -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)
Expand All @@ -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.

Expand All @@ -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)
Expand All @@ -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)
Expand Down
72 changes: 52 additions & 20 deletions src/transformez/cli.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,13 +15,12 @@
import click
import logging
from transformez import api
from transformez import grid_engine

logger = logging.getLogger(__name__)


@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."""

Expand All @@ -30,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.
Expand All @@ -58,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
Expand All @@ -81,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)


Expand All @@ -107,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)
Expand All @@ -119,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()
Loading
Loading