|
11 | 11 | :license: MIT, see LICENSE for more details. |
12 | 12 | """ |
13 | 13 |
|
14 | | -import os |
15 | 14 | import sys |
16 | | -import argparse |
| 15 | +import click |
17 | 16 | import logging |
18 | | -import rasterio |
19 | | - |
20 | | -from . import __version__ |
21 | | -from .transform import VerticalTransform |
22 | | -from .definitions import Datums |
23 | | -from .grid_engine import plot_grid, GridWriter, GridEngine |
24 | | -from .htdp import HAS_HTDP, download_htdp |
25 | | - |
26 | | -from fetchez import spatial |
27 | | -from fetchez import utils |
28 | | -from fetchez.spatial import parse_region, fix_argparse_region |
29 | | - |
30 | | -logging.basicConfig(level=logging.INFO, format="[ %(levelname)s ] %(message)s", stream=sys.stderr) |
31 | | -logger = logging.getLogger("transformez") |
32 | | -logging.getLogger("fetchez").setLevel(logging.WARNING) |
33 | | - |
34 | | - |
35 | | -def parse_compound_datum(datum_arg): |
36 | | - """Parse string 'EPSG:GEOID' or 'NAME'.""" |
37 | | - |
38 | | - if not datum_arg: |
39 | | - return None, None |
40 | | - s = str(datum_arg) |
41 | | - if ':' in s: |
42 | | - parts = s.split(':') |
43 | | - # Handle "5703:geoid=g2012b" or "5703:g2012b" |
44 | | - geoid_part = parts[1] |
45 | | - if 'geoid=' in geoid_part: |
46 | | - geoid = geoid_part.split('=')[1] |
47 | | - else: |
48 | | - geoid = geoid_part |
49 | | - return Datums.get_vdatum_by_name(parts[0]), geoid |
50 | | - else: |
51 | | - return Datums.get_vdatum_by_name(s), None |
52 | | - |
53 | | - |
54 | | -def list_supported_datums(): |
55 | | - """Pretty print supported datums.""" |
56 | | - |
57 | | - print(f"\nTransformez v{__version__} - Supported Datums\n") |
58 | | - |
59 | | - print("--- Tidal Surfaces (Local & Global) ---") |
60 | | - for k, v in Datums.SURFACES.items(): |
61 | | - region = v.get("region", "global").upper() |
62 | | - print(f" {v['name']:<10} : {v['description']} [{region}]") |
| 17 | +from transformez import api |
63 | 18 |
|
64 | | - print("\n--- Ellipsoidal / Frame (HTDP) ---") |
65 | | - print(f" {'NAD83':<10} : North American Datum 1983 (EPSG:6319)") |
66 | | - print(f" {'WGS84':<10} : World Geodetic System 1984 (EPSG:4979)") |
67 | | - |
68 | | - print("\n--- Orthometric / Geoid-Based ---") |
69 | | - for k, v in Datums.CDN.items(): |
70 | | - print(f" {v['name']:<20} (Default Geoid: {v.get('default_geoid', 'None')})") |
71 | | - |
72 | | - print("\n--- Available Geoids ---") |
73 | | - print(f" {', '.join(Datums.GEOIDS.keys())}") |
74 | | - print("\n") |
| 19 | +logger = logging.getLogger(__name__) |
75 | 20 |
|
76 | 21 |
|
| 22 | +@click.group(name="transform") |
| 23 | +@click.version_option(package_name='transformez') |
77 | 24 | def transformez_cli(): |
78 | | - parser = argparse.ArgumentParser( |
79 | | - #description=f"{utils.CYAN}%(prog)s{utils.RESET} ({__version__}) :: Global Vertical Datum Transformer", |
80 | | - description=utils._cli_logo("transformez", "Global Vertical Datum Transformer", __version__), |
81 | | - epilog="Examples:\n" |
82 | | - " transformez -R -166/-164/63/64 -I mllw -O 4979 -E3s\n" |
83 | | - " transformez input_dem.tif -I mllw -O 5703:geoid=g2012b", |
84 | | - formatter_class=argparse.RawDescriptionHelpFormatter |
85 | | - ) |
86 | | - |
87 | | - parser.add_argument('input_file', nargs='?', help="Input DEM (GeoTIFF) to transform.") |
88 | | - |
89 | | - grp_loc = parser.add_argument_group('Location & Resolution (if no input file)') |
90 | | - grp_loc.add_argument('-R', '--region', nargs='?', help='Region (West/East/South/North).') |
91 | | - grp_loc.add_argument('-E', '--increment', nargs='?', help='Resolution (e.g. 1s, 0.0001).') |
92 | | - |
93 | | - grp_dat = parser.add_argument_group('Vertical Datums') |
94 | | - grp_dat.add_argument('-I', '--input-datum', '--vdatum-in', |
95 | | - help='Source Datum (e.g. "mllw", "5703", "4979").') |
96 | | - grp_dat.add_argument('-O', '--output-datum', '--vdatum-out', |
97 | | - help='Target Datum (e.g. "5703:geoid=g2012b", "4979").') |
98 | | - |
99 | | - grp_out = parser.add_argument_group('Output') |
100 | | - grp_out.add_argument('-o', '--output', help='Output filename (default: auto-named).') |
101 | | - grp_out.add_argument('--preview', action='store_true', help='Show plot of shift grid before saving.') |
102 | | - |
103 | | - grp_proc = parser.add_argument_group('Process') |
104 | | - grp_proc.add_argument('--decay-pixels', type=int, default=100, |
105 | | - help='Number of pixels to decay tidal shifts inland (default: 100). Set to 0 for strict VDatum boundaries.') |
106 | | - |
107 | | - grp_sys = parser.add_argument_group('System') |
108 | | - grp_sys.add_argument('--list-datums', action='store_true', help='List supported datums.') |
109 | | - grp_sys.add_argument('--download-htdp', action='store_true', help='Download the NGS HTDP software to the current directory.') |
110 | | - grp_sys.add_argument('--cache-dir', help='Override cache directory.') |
111 | | - grp_sys.add_argument('--verbose', action='store_true', help='Enable debug logging.') |
112 | | - grp_sys.add_argument( |
113 | | - "-v", "--version", action="version", version=f"%(prog)s {__version__}" |
114 | | - ) |
115 | | - |
116 | | - fixed_argv = fix_argparse_region(sys.argv[1:]) |
117 | | - args = parser.parse_args(fixed_argv) |
118 | | - |
119 | | - if args.download_htdp: |
120 | | - download_htdp() |
121 | | - sys.exit(0) |
122 | | - |
123 | | - if args.list_datums: |
124 | | - list_supported_datums() |
125 | | - sys.exit(0) |
126 | | - |
127 | | - # Validation & Setup |
128 | | - if args.verbose: |
129 | | - logger.setLevel(logging.DEBUG) |
130 | | - logging.getLogger('fetchez').setLevel(logging.INFO) |
131 | | - |
132 | | - if not HAS_HTDP: |
133 | | - logger.warning("=" * 60) |
134 | | - logger.warning("HTDP tool not found in PATH.") |
135 | | - logger.warning("Frame transformations (e.g., NAD83 <-> WGS84) will return zero shift.") |
136 | | - logger.warning("Run 'transformez --download-htdp' to get the software.") |
137 | | - logger.warning("=" * 60) |
138 | | - |
139 | | - cache_dir = args.cache_dir or os.path.join(os.getcwd(), "transformez_cache") |
140 | | - if not os.path.exists(cache_dir): |
141 | | - os.makedirs(cache_dir) |
142 | | - |
143 | | - # Parse Datums |
144 | | - epsg_in, geoid_in = parse_compound_datum(args.input_datum) |
145 | | - epsg_out, geoid_out = parse_compound_datum(args.output_datum) |
146 | | - |
147 | | - if not epsg_in or not epsg_out: |
148 | | - parser.print_help() |
149 | | - logger.error("Invalid datum specified. Use --list-datums to see options.") |
150 | | - sys.exit(1) |
151 | | - |
152 | | - # Mode Selection (DEM vs Grid) |
153 | | - target_dem = args.input_file |
154 | | - |
155 | | - if target_dem: |
156 | | - # Transform DEM |
157 | | - if not os.path.exists(target_dem): |
158 | | - parser.print_help() |
159 | | - logger.error(f"Input file not found: {target_dem}") |
160 | | - sys.exit(1) |
161 | | - |
162 | | - logger.info(f"Reading bounds from: {target_dem}") |
163 | | - with rasterio.open(target_dem) as src: |
164 | | - bounds = src.bounds |
165 | | - region_obj = spatial.Region(bounds.left, bounds.right, bounds.bottom, bounds.top) |
166 | | - nx, ny = src.width, src.height |
167 | | - |
168 | | - if not args.output: |
169 | | - base, ext = os.path.splitext(target_dem) |
170 | | - dst_fn = f"{base}_trans_{args.output_datum.replace(':','_')}{ext}" |
| 25 | + """Apply vertical datum transformations and generate shift grids.""" |
| 26 | + |
| 27 | + pass |
| 28 | + |
| 29 | +@transformez_cli.command("run") |
| 30 | +@click.argument("input_file", required=False) |
| 31 | +@click.option("-R", "--region", help="Bounding box or location string (if no input file).") |
| 32 | +@click.option("-E", "--increment", help="Resolution (e.g., 1s, 30m) (if no input file).") |
| 33 | +@click.option("-I", "--input-datum", required=True, help="Source Datum (e.g., 'mllw', '5703').") |
| 34 | +@click.option("-O", "--output-datum", required=True, help="Target Datum (e.g., '4979', '5703:g2012b').") |
| 35 | +@click.option("--out", "-o", help="Output filename (default: auto-named).") |
| 36 | +@click.option("--decay-pixels", type=int, default=100, help="Number of pixels to decay tidal shifts inland.") |
| 37 | +def transform_run(input_file, region, increment, input_datum, output_datum, out, decay_pixels): |
| 38 | + """Transform a raster's vertical datum or generate a standalone shift grid. |
| 39 | +
|
| 40 | + If an INPUT_FILE is provided, that specific raster is transformed in place. |
| 41 | + If no INPUT_FILE is provided, -R and -E must be used to generate a shift grid. |
| 42 | +
|
| 43 | + Examples: |
| 44 | + Transform a DEM : transformez run my_dem.tif -I mllw -O 5703 |
| 45 | + Generate a Grid : transformez run -R loc:"Miami" -E 1s -I mllw -O 4979 |
| 46 | + """ |
| 47 | + if input_file: |
| 48 | + click.secho(f"🚀 Transforming raster: {input_file}", fg="cyan", bold=True) |
| 49 | + click.echo(f" Shift: {input_datum} ➔ {output_datum}") |
| 50 | + |
| 51 | + result = api.transform_raster( |
| 52 | + input_raster=input_file, |
| 53 | + datum_in=input_datum, |
| 54 | + datum_out=output_datum, |
| 55 | + decay_pixels=decay_pixels, |
| 56 | + output_raster=out, |
| 57 | + verbose=True |
| 58 | + ) |
| 59 | + |
| 60 | + if result: |
| 61 | + click.secho(f"✅ Successfully transformed raster: {result}", fg="green", bold=True) |
171 | 62 | else: |
172 | | - dst_fn = args.output |
173 | | - |
174 | | - elif args.region: |
175 | | - # Generate Grid |
176 | | - if not args.increment: |
177 | | - parser.print_help() |
178 | | - logger.error("Increment (-E) is required when generating a grid from scratch.") |
| 63 | + click.secho("❌ Failed to transform raster.", fg="red") |
179 | 64 | sys.exit(1) |
180 | 65 |
|
181 | | - regions = parse_region(args.region) |
182 | | - region_obj = regions[0] |
183 | | - |
184 | | - # Parse Increment |
185 | | - try: |
186 | | - inc_val = utils.str2inc(args.increment) |
187 | | - nx = int(region_obj.width / inc_val) |
188 | | - ny = int(region_obj.height / inc_val) |
189 | | - except Exception as e: |
190 | | - parser.print_help() |
191 | | - logger.error(f"Invalid increment: {args.increment}: {e}") |
192 | | - sys.exit(1) |
193 | | - |
194 | | - if not args.output: |
195 | | - dst_fn = f"shift_{args.input_datum}_to_{args.output_datum.replace(':','_')}.tif" |
| 66 | + elif region and increment: |
| 67 | + click.secho(f"🚀 Generating vertical shift grid for region...", fg="cyan", bold=True) |
| 68 | + click.echo(f" Shift: {input_datum} ➔ {output_datum} @ {increment}") |
| 69 | + |
| 70 | + # Auto-generate an output name if one wasn't provided |
| 71 | + out_fn = out or f"shift_{input_datum}_to_{output_datum.replace(':', '_')}.tif" |
| 72 | + |
| 73 | + result = api.generate_grid( |
| 74 | + region=region, |
| 75 | + increment=increment, |
| 76 | + datum_in=input_datum, |
| 77 | + datum_out=output_datum, |
| 78 | + decay_pixels=decay_pixels, |
| 79 | + out_fn=out_fn, |
| 80 | + verbose=True |
| 81 | + ) |
| 82 | + |
| 83 | + if result is not None: |
| 84 | + click.secho(f"✅ Successfully generated shift grid: {out_fn}", fg="green", bold=True) |
196 | 85 | else: |
197 | | - dst_fn = args.output |
| 86 | + click.secho("❌ Failed to generate shift grid.", fg="red") |
| 87 | + sys.exit(1) |
198 | 88 |
|
199 | 89 | else: |
200 | | - parser.print_help() |
201 | | - logger.error("Either an input file OR a region (-R) is required.") |
| 90 | + click.secho("❌ Error: You must provide either an INPUT_FILE or both --region and --increment.", fg="red") |
202 | 91 | sys.exit(1) |
203 | 92 |
|
204 | | - vt = VerticalTransform( |
205 | | - region=region_obj, |
206 | | - nx=nx, ny=ny, |
207 | | - epsg_in=epsg_in, epsg_out=epsg_out, |
208 | | - geoid_in=geoid_in, geoid_out=geoid_out, |
209 | | - decay_pixels=args.decay_pixels, |
210 | | - cache_dir=cache_dir, |
211 | | - verbose=args.verbose |
212 | | - ) |
213 | | - |
214 | | - logger.info(f"Computing Shift: {args.input_datum} -> {args.output_datum}") |
215 | | - shift_array, _ = vt._vertical_transform(vt.epsg_in, vt.epsg_out) |
216 | | - |
217 | | - if shift_array is None: |
218 | | - logger.error("Transformation failed (No coverage found).") |
219 | | - sys.exit(1) |
220 | 93 |
|
221 | | - if args.preview: |
222 | | - plot_grid(shift_array, region_obj, title=f"{args.input_datum} -> {args.output_datum}") |
| 94 | +@transformez_cli.command("list") |
| 95 | +def transform_list(): |
| 96 | + """List all supported vertical datums, EPSG codes, and geoids.""" |
| 97 | + try: |
| 98 | + from transformez.definitions import Datums |
223 | 99 |
|
224 | | - if target_dem: |
225 | | - logger.info(f"Applying shift to {target_dem}...") |
226 | | - GridEngine.apply_vertical_shift(target_dem, shift_array, dst_fn) |
227 | | - else: |
228 | | - logger.info(f"Writing shift {dst_fn}...") |
229 | | - GridWriter.write(dst_fn, shift_array, region_obj) |
| 100 | + click.secho("\n🌊 Supported Tidal Surfaces:", fg="cyan", bold=True) |
| 101 | + # For tidal datums, the user types the dictionary key (e.g., 'mllw') |
| 102 | + for key, v in Datums.SURFACES.items(): |
| 103 | + region_str = v.get('region', 'global').upper() |
| 104 | + click.echo(f" {key:<12} : {v.get('name', key):<30} [{region_str}]") |
| 105 | + |
| 106 | + click.secho("\n🌐 Ellipsoidal / Frame Datums (EPSG):", fg="cyan", bold=True) |
| 107 | + # For ellipsoidal, explicitly list the EPSG codes |
| 108 | + click.echo(f" {'4979':<12} : WGS84 - World Geodetic System 1984") |
| 109 | + click.echo(f" {'6319':<12} : NAD83 - North American Datum 1983") |
| 110 | + |
| 111 | + click.secho("\n🏔️ Orthometric / Geoid-Based (EPSG):", fg="cyan", bold=True) |
| 112 | + # For orthometric, the key in Datums.CDN is typically the EPSG code (e.g., '5703') |
| 113 | + for epsg_key, v in Datums.CDN.items(): |
| 114 | + # Fallback to the key if 'epsg' isn't explicitly defined in the dict |
| 115 | + epsg_code = v.get('epsg', epsg_key) |
| 116 | + geoid_str = v.get('default_geoid', 'None') |
| 117 | + click.echo(f" {str(epsg_code):<12} : {v.get('name', 'Unknown'):<30} (Default Geoid: {geoid_str})") |
| 118 | + |
| 119 | + click.secho("\n🌍 Available Geoids:", fg="cyan", bold=True) |
| 120 | + click.echo(f" {', '.join(Datums.GEOIDS.keys())}") |
| 121 | + |
| 122 | + click.secho("\n💡 Pro-Tip:", fg="yellow", bold=True, nl=False) |
| 123 | + click.echo(" Combine an EPSG and a specific Geoid using a colon (e.g., -O 5703:g2012b)\n") |
230 | 124 |
|
231 | | - logger.info(f"Success: {dst_fn}") |
| 125 | + except ImportError: |
| 126 | + click.secho("❌ Error: Could not load Transformez datum definitions.", fg="red") |
232 | 127 |
|
233 | 128 | if __name__ == '__main__': |
234 | 129 | transformez_cli() |
0 commit comments