diff --git a/conda_package/docs/landice/adjust_bed_to_haf.rst b/conda_package/docs/landice/adjust_bed_to_haf.rst new file mode 100644 index 000000000..18a51de31 --- /dev/null +++ b/conda_package/docs/landice/adjust_bed_to_haf.rst @@ -0,0 +1,151 @@ +.. _landice_adjust_bed_to_haf: + +************************************ +Adjusting Bed to Height Above Flotation +************************************ + +The ``adjust_bed_to_haf`` command-line tool adjusts bed topography within +grounding line regions to achieve a specified height above flotation (HAF). + +Height Above Flotation +====================== + +Height above flotation (HAF) is a measure of how far an ice column is from +hydrostatic equilibrium with ocean water. It is defined as: + +.. math:: + + \text{HAF} = z_{\text{surface}} - z_{\text{flotation}} + +where :math:`z_{\text{surface}}` is the ice surface elevation and +:math:`z_{\text{flotation}}` is the surface elevation at which the ice would +be in flotation equilibrium. + +For ice in flotation equilibrium: + +.. math:: + + \rho_{\text{ice}} \cdot H = \rho_{\text{ocean}} \cdot d + +where :math:`H` is ice thickness, :math:`d` is the draft (depth below sea +level), and :math:`\rho_{\text{ice}}` and :math:`\rho_{\text{ocean}}` are ice +and ocean water densities, respectively. + +Usage +===== + +Basic usage: + +.. code-block:: bash + + adjust_bed_to_haf \ + --mesh input_mesh.nc \ + --geojson grounding_line.geojson \ + --output output_mesh.nc \ + --target-haf 10.0 + +This will adjust the bed topography in cells within the grounding line polygons +to achieve a height above flotation of 10 meters. + +Command-Line Options +==================== + +Required Arguments +------------------ + +``-m, --mesh MESH_FILE`` + MALI mesh file in NetCDF format containing ice thickness and bed topography + +``-g, --geojson GEOJSON_FILE`` + GeoJSON file containing grounding line delineations as polygon geometries + +Optional Arguments +------------------ + +``-o, --output OUTPUT_FILE`` + Output mesh file. If not specified, the input mesh file is modified in place. + +``--target-haf HAF`` + Target height above flotation in meters (default: 10.0) + +``--sea-level LEVEL`` + Sea level in meters (default: 0.0) + +``--rho-ice DENSITY`` + Ice density in kg/m³ (default: 910.0) + +``--rho-ocean DENSITY`` + Ocean water density in kg/m³ (default: 1028.0) + +``--thickness-var VARNAME`` + Name of thickness variable in mesh file (default: 'thickness') + +``--bed-var VARNAME`` + Name of bed topography variable in mesh file (default: 'bedTopography') + +Example +======= + +The following example adjusts bed topography within Thwaites Glacier grounding +line regions to achieve a 15-meter height above flotation: + +.. code-block:: bash + + adjust_bed_to_haf \ + --mesh AIS_4to20km_mesh.nc \ + --geojson Thwaites_GL_2014_pinning_points.geojson \ + --output AIS_4to20km_mesh_adjusted.nc \ + --target-haf 15.0 \ + --rho-ice 918.0 \ + --rho-ocean 1028.0 + +The tool will: + +1. Load the grounding line polygons from the GeoJSON file +2. Identify all mesh cells that fall within these polygons +3. Calculate the current height above flotation in these regions +4. Adjust the bed topography to achieve the target HAF while preserving ice thickness +5. Update the output file with metadata about the changes + +Technical Details +================= + +The tool calculates the required bed elevation :math:`b` to achieve a target +HAF given ice thickness :math:`H`: + +.. math:: + + b = z_{\text{sea}} + \text{HAF}_{\text{target}} - H \cdot \frac{\rho_{\text{ice}}}{\rho_{\text{ocean}}} + +This ensures that with the existing ice thickness, the ice surface will be at +the specified height above the flotation level. + +Notes +===== + +- The GeoJSON file should use WGS 84 (EPSG:4326) coordinates (longitude, latitude) +- The tool converts mesh coordinates from radians to degrees for comparison with GeoJSON +- If no cells are found within the grounding line polygons, check that the coordinate systems are compatible +- The tool updates the file's ``history`` and ``comment`` global attributes to document the modifications +- Ice thickness is not modified; only bed topography is adjusted + +Dependencies +============ + +Required Python packages: + +- netCDF4 +- numpy +- shapely + +Optional (for better performance): + +- geopandas (falls back to json + shapely if not available) + +References +========== + +Wild, C. T., Alley, K. E., Muto, A., Pettit, E. C., Scambos, T. A., & Truffer, M. (2022). +Thwaites Glacier 2014 and 2019/20 grounding line positions. +U.S. Antarctic Program (USAP) Data Center. +doi: 10.15784/601499 diff --git a/conda_package/mpas_tools/landice/README_adjust_bed_to_haf.md b/conda_package/mpas_tools/landice/README_adjust_bed_to_haf.md new file mode 100644 index 000000000..798cf8d36 --- /dev/null +++ b/conda_package/mpas_tools/landice/README_adjust_bed_to_haf.md @@ -0,0 +1,145 @@ +# Adjust Bed to Height Above Flotation + +This tool adjusts bed topography within grounding line regions to achieve a specified height above flotation (HAF). + +## Quick Start + +```bash +adjust_bed_to_haf \ + --mesh input_mesh.nc \ + --geojson grounding_line.geojson \ + --output output_mesh.nc \ + --target-haf 10.0 +``` + +## What is Height Above Flotation? + +Height above flotation (HAF) measures how far ice is from hydrostatic equilibrium with ocean water: + +- **HAF = 0**: Ice is exactly at flotation (neither grounded nor floating freely) +- **HAF > 0**: Ice surface is above the flotation level (typical for grounded ice) +- **HAF < 0**: Ice surface is below flotation level (should not occur in steady state) + +## How It Works + +The tool: + +1. Loads grounding line polygons from a GeoJSON file +2. Identifies mesh cells within these polygons +3. Calculates required bed elevation to achieve target HAF +4. Updates the `bedTopography` field in the mesh file +5. Preserves ice thickness (only bed is modified) + +## Physics + +For a given ice thickness H and target HAF, the required bed elevation is: + +``` +bed = sea_level + HAF_target - H × (ρ_ice / ρ_ocean) +``` + +Where: +- ρ_ice = 910 kg/m³ (default) +- ρ_ocean = 1028 kg/m³ (default) + +## Input File Requirements + +### Mesh File (NetCDF) +Must contain: +- `lonCell`, `latCell`: Cell center coordinates (in radians) +- `thickness`: Ice thickness (m) +- `bedTopography`: Bed elevation (m, positive up) + +### GeoJSON File +Must contain: +- Polygon or MultiLineString geometries +- Coordinates in WGS 84 (longitude, latitude in degrees) + +## Example: Thwaites Glacier + +```bash +adjust_bed_to_haf \ + --mesh /path/to/AIS_4to20km_mesh.nc \ + --geojson /path/to/Thwaites_GL_2014_pinning_points.geojson \ + --output AIS_4to20km_mesh_haf15m.nc \ + --target-haf 15.0 \ + --rho-ice 918.0 +``` + +## Common Options + +- `--target-haf 10.0`: Set target height above flotation (meters) +- `--sea-level 0.0`: Set sea level (meters) +- `--rho-ice 910.0`: Set ice density (kg/m³) +- `--rho-ocean 1028.0`: Set ocean water density (kg/m³) +- `--output file.nc`: Save to new file (otherwise modifies in place) + +## Troubleshooting + +### "No cells found within grounding line polygons" + +This usually means: +1. Coordinate system mismatch (mesh in radians vs. GeoJSON in degrees) +2. Grounding line doesn't overlap with mesh domain +3. GeoJSON uses wrong projection + +**Solution**: Verify that your GeoJSON is in WGS 84 (EPSG:4326) coordinates. + +### Variable not found + +The tool looks for `thickness` and `bedTopography` by default. If your mesh uses different names: + +```bash +adjust_bed_to_haf \ + --thickness-var myThickness \ + --bed-var myBedTopo \ + ... +``` + +## Use Cases + +### 1. Initialize grounding line regions +Set grounding line areas to a specific HAF for model initialization: + +```bash +adjust_bed_to_haf --mesh init.nc --geojson gl.geojson --target-haf 10.0 +``` + +### 2. Sensitivity experiments +Test model response to different HAF values: + +```bash +for haf in 5 10 15 20; do + adjust_bed_to_haf \ + --mesh base.nc \ + --geojson gl.geojson \ + --output mesh_haf${haf}m.nc \ + --target-haf $haf +done +``` + +### 3. Adjust pinning points +Modify bed topography at pinning points to control ice sheet stability: + +```bash +adjust_bed_to_haf \ + --mesh mesh.nc \ + --geojson pinning_points.geojson \ + --target-haf 5.0 \ + --output mesh_weak_pinning.nc +``` + +## Dependencies + +**Required:** +- netCDF4 +- numpy +- shapely + +**Optional:** +- geopandas (recommended for better performance) + +## See Also + +- [Full documentation](../../docs/landice/adjust_bed_to_haf.rst) +- [MPAS-Tools landice module](https://mpas-dev.github.io/MPAS-Tools/stable/landice.html) diff --git a/conda_package/mpas_tools/landice/adjust_bed_to_haf.py b/conda_package/mpas_tools/landice/adjust_bed_to_haf.py new file mode 100755 index 000000000..4f80e6c86 --- /dev/null +++ b/conda_package/mpas_tools/landice/adjust_bed_to_haf.py @@ -0,0 +1,426 @@ +#!/usr/bin/env python3 +""" +Script to adjust bed topography within grounding line regions to achieve +a specified height above flotation (HAF). + +Height above flotation is defined as: + HAF = ice_surface - flotation_surface +where flotation_surface is the surface elevation at which ice would be +in hydrostatic equilibrium with ocean water. + +For a given ice thickness and target HAF, the required bed elevation is: + bed = ice_surface - thickness +where ice_surface is calculated to achieve the target HAF. +""" + +import argparse +import sys +from datetime import datetime + +import netCDF4 +import numpy as np +from shapely.geometry import Point, shape + +try: + import geopandas as gpd + HAS_GEOPANDAS = True +except ImportError: + HAS_GEOPANDAS = False + import json + + +def load_grounding_line_geojson(geojson_file): + """ + Load grounding line polygons from a GeoJSON file. + + Parameters + ---------- + geojson_file : str + Path to the GeoJSON file containing grounding line delineations + + Returns + ------- + geometries : list + List of shapely geometry objects + """ + if HAS_GEOPANDAS: + # Use geopandas if available + gdf = gpd.read_file(geojson_file) + geometries = gdf.geometry.tolist() + else: + # Fall back to json + shapely + with open(geojson_file, 'r') as f: + data = json.load(f) + + geometries = [] + for feature in data['features']: + geom = shape(feature['geometry']) + geometries.append(geom) + + return geometries + + +def find_cells_in_polygon(lon_cell, lat_cell, geometries): + """ + Find mesh cells that fall within any of the provided polygons. + + Parameters + ---------- + lon_cell : ndarray + Longitude of cell centers (in degrees) + lat_cell : ndarray + Latitude of cell centers (in degrees) + geometries : list + List of shapely geometry objects + + Returns + ------- + mask : ndarray (bool) + Boolean mask indicating which cells are inside the polygons + """ + n_cells = len(lon_cell) + mask = np.zeros(n_cells, dtype=bool) + + print(f'Checking {n_cells} cells against {len(geometries)} geometries...') + + for i, (lon, lat) in enumerate(zip(lon_cell, lat_cell)): + if i % 10000 == 0: + print(f' Processed {i}/{n_cells} cells...') + + point = Point(lon, lat) + for geom in geometries: + if geom.contains(point) or geom.intersects(point): + mask[i] = True + break + + print(f'Found {np.sum(mask)} cells within polygons') + return mask + + +def calculate_flotation_thickness(bed_elevation, sea_level=0.0, + rho_ice=910.0, rho_ocean=1028.0): + """ + Calculate the ice thickness at which ice would be in flotation. + + Parameters + ---------- + bed_elevation : ndarray + Bed topography (positive up, m) + sea_level : float + Sea level (m), default 0.0 + rho_ice : float + Ice density (kg/m^3), default 910.0 + rho_ocean : float + Ocean water density (kg/m^3), default 1028.0 + + Returns + ------- + flotation_thickness : ndarray + Ice thickness at flotation (m) + """ + # For ice to float: rho_ice * thickness = rho_ocean * draft + # where draft = sea_level - bed_elevation + # Therefore: thickness_flotation = (rho_ocean / rho_ice) * draft + + draft = sea_level - bed_elevation + flotation_thickness = (rho_ocean / rho_ice) * draft + + # Thickness must be positive + flotation_thickness = np.maximum(flotation_thickness, 0.0) + + return flotation_thickness + + +def calculate_required_bed(thickness, target_haf, sea_level=0.0, + rho_ice=910.0, rho_ocean=1028.0): + """ + Calculate the bed elevation required to achieve a target height above + flotation for a given ice thickness. + + Parameters + ---------- + thickness : ndarray + Ice thickness (m) + target_haf : float + Target height above flotation (m) + sea_level : float + Sea level (m), default 0.0 + rho_ice : float + Ice density (kg/m^3), default 910.0 + rho_ocean : float + Ocean water density (kg/m^3), default 1028.0 + + Returns + ------- + bed_elevation : ndarray + Required bed elevation (m) + """ + # At flotation: + # ice_surface_flotation = sea_level + thickness * (1 - rho_ice/rho_ocean) + # + # With target HAF: + # ice_surface = ice_surface_flotation + target_haf + # + # Since ice_surface = bed + thickness: + # bed = ice_surface - thickness + # = ice_surface_flotation + target_haf - thickness + # = sea_level + thickness * (1 - rho_ice/rho_ocean) + target_haf - thickness + # = sea_level + target_haf - thickness * rho_ice/rho_ocean + + bed_elevation = sea_level + target_haf - thickness * (rho_ice / rho_ocean) + + return bed_elevation + + +def adjust_bed_to_haf(): + """ + Main function to adjust bed topography to achieve target height above + flotation within grounding line regions. + """ + + parser = argparse.ArgumentParser( + description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter + ) + parser.add_argument( + '-m', '--mesh', + dest='mesh_file', + required=True, + help='MALI mesh file (NetCDF format)' + ) + parser.add_argument( + '-g', '--geojson', + dest='geojson_file', + required=True, + help='GeoJSON file containing grounding line delineations' + ) + parser.add_argument( + '-o', '--output', + dest='output_file', + help='Output mesh file. If not specified, modifies input file in place.' + ) + parser.add_argument( + '--target-haf', + dest='target_haf', + type=float, + default=10.0, + help='Target height above flotation in meters (default: 10.0)' + ) + parser.add_argument( + '--sea-level', + dest='sea_level', + type=float, + default=0.0, + help='Sea level in meters (default: 0.0)' + ) + parser.add_argument( + '--rho-ice', + dest='rho_ice', + type=float, + default=910.0, + help='Ice density in kg/m^3 (default: 910.0)' + ) + parser.add_argument( + '--rho-ocean', + dest='rho_ocean', + type=float, + default=1028.0, + help='Ocean water density in kg/m^3 (default: 1028.0)' + ) + parser.add_argument( + '--thickness-var', + dest='thickness_var', + default='thickness', + help='Name of thickness variable in mesh file (default: thickness)' + ) + parser.add_argument( + '--bed-var', + dest='bed_var', + default='bedTopography', + help='Name of bed topography variable in mesh file (default: bedTopography)' + ) + + args = parser.parse_args() + + # Validate inputs + if not HAS_GEOPANDAS: + print('Warning: geopandas not available, using json + shapely instead') + + print('\n** Adjusting bed topography to achieve target height above flotation **') + print(f'Input mesh file: {args.mesh_file}') + print(f'Grounding line GeoJSON: {args.geojson_file}') + print(f'Target HAF: {args.target_haf} m') + print(f'Sea level: {args.sea_level} m') + print(f'Ice density: {args.rho_ice} kg/m^3') + print(f'Ocean density: {args.rho_ocean} kg/m^3') + print() + + # Load grounding line polygons + print('Loading grounding line geometries...') + geometries = load_grounding_line_geojson(args.geojson_file) + print(f'Loaded {len(geometries)} geometry features') + print() + + # Open mesh file + print('Opening mesh file...') + if args.output_file: + # Copy to output file + import shutil + shutil.copy2(args.mesh_file, args.output_file) + mesh_file = args.output_file + print(f'Copied input to output file: {args.output_file}') + else: + mesh_file = args.mesh_file + print('Modifying mesh file in place') + + with netCDF4.Dataset(mesh_file, 'r+') as mesh: + + # Read mesh coordinates + print('Reading mesh coordinates...') + # Coordinates are typically in radians, convert to degrees + lon_cell = np.degrees(mesh.variables['lonCell'][:]) + lat_cell = np.degrees(mesh.variables['latCell'][:]) + n_cells = len(lon_cell) + print(f'Mesh has {n_cells} cells') + print() + + # Find cells within grounding line polygons + print('Identifying cells within grounding line regions...') + mask = find_cells_in_polygon(lon_cell, lat_cell, geometries) + print() + + if np.sum(mask) == 0: + print('ERROR: No cells found within grounding line polygons!') + print('Check that:') + print(' 1. GeoJSON and mesh use compatible coordinate systems') + print(' 2. GeoJSON geometries overlap with mesh extent') + sys.exit(1) + + # Read thickness and bed topography + print('Reading ice thickness and bed topography...') + if args.thickness_var not in mesh.variables: + print(f'ERROR: Variable "{args.thickness_var}" not found in mesh file') + print(f'Available variables: {list(mesh.variables.keys())}') + sys.exit(1) + + if args.bed_var not in mesh.variables: + print(f'ERROR: Variable "{args.bed_var}" not found in mesh file') + print(f'Available variables: {list(mesh.variables.keys())}') + sys.exit(1) + + thickness = mesh.variables[args.thickness_var][:] + bed_topo = mesh.variables[args.bed_var][:] + + # Handle potential time dimension + if len(thickness.shape) > 1: + # Assume time is first dimension + thickness = thickness[0, :] + bed_topo = bed_topo[0, :] + has_time_dim = True + else: + has_time_dim = False + + print(f'Thickness range: [{np.min(thickness):.2f}, {np.max(thickness):.2f}] m') + print(f'Bed topography range: [{np.min(bed_topo):.2f}, {np.max(bed_topo):.2f}] m') + print() + + # Calculate current HAF in the region + print('Calculating current height above flotation...') + flotation_thickness = calculate_flotation_thickness( + bed_topo, + sea_level=args.sea_level, + rho_ice=args.rho_ice, + rho_ocean=args.rho_ocean + ) + ice_surface = bed_topo + thickness + flotation_surface = args.sea_level + flotation_thickness * ( + 1.0 - args.rho_ice / args.rho_ocean + ) + current_haf = ice_surface - flotation_surface + + print(f'Current HAF in region (mean): {np.mean(current_haf[mask]):.2f} m') + print(f'Current HAF in region (min): {np.min(current_haf[mask]):.2f} m') + print(f'Current HAF in region (max): {np.max(current_haf[mask]):.2f} m') + print() + + # Calculate new bed elevation + print('Calculating new bed topography...') + new_bed = calculate_required_bed( + thickness[mask], + args.target_haf, + sea_level=args.sea_level, + rho_ice=args.rho_ice, + rho_ocean=args.rho_ocean + ) + + print(f'New bed range in region: [{np.min(new_bed):.2f}, {np.max(new_bed):.2f}] m') + bed_change = new_bed - bed_topo[mask] + print(f'Bed change (mean): {np.mean(bed_change):.2f} m') + print(f'Bed change (min): {np.min(bed_change):.2f} m') + print(f'Bed change (max): {np.max(bed_change):.2f} m') + print() + + # Update bed topography + print('Updating bed topography in mesh file...') + bed_topo[mask] = new_bed + + if has_time_dim: + # Write back with time dimension + mesh.variables[args.bed_var][0, :] = bed_topo + else: + mesh.variables[args.bed_var][:] = bed_topo + + # Update global attributes + timestamp = datetime.now().strftime('%Y-%m-%d %H:%M:%S') + if 'history' in mesh.ncattrs(): + history = mesh.getncattr('history') + history = f'{timestamp}: adjust_bed_to_haf.py\n{history}' + else: + history = f'{timestamp}: adjust_bed_to_haf.py' + mesh.setncattr('history', history) + + comment = ( + f'Bed topography adjusted within grounding line regions ' + f'to achieve target HAF = {args.target_haf} m. ' + f'Modified {np.sum(mask)} cells using grounding line data from ' + f'{args.geojson_file}.' + ) + if 'comment' in mesh.ncattrs(): + existing_comment = mesh.getncattr('comment') + comment = f'{comment} {existing_comment}' + mesh.setncattr('comment', comment) + + print('Successfully updated mesh file!') + print() + + # Verify the result + print('Verifying updated HAF...') + flotation_thickness = calculate_flotation_thickness( + bed_topo, + sea_level=args.sea_level, + rho_ice=args.rho_ice, + rho_ocean=args.rho_ocean + ) + ice_surface = bed_topo + thickness + flotation_surface = args.sea_level + flotation_thickness * ( + 1.0 - args.rho_ice / args.rho_ocean + ) + new_haf = ice_surface - flotation_surface + + print(f'New HAF in region (mean): {np.mean(new_haf[mask]):.2f} m') + print(f'New HAF in region (min): {np.min(new_haf[mask]):.2f} m') + print(f'New HAF in region (max): {np.max(new_haf[mask]):.2f} m') + print() + + if np.allclose(new_haf[mask], args.target_haf, rtol=1e-6): + print(f'✓ Successfully achieved target HAF = {args.target_haf} m') + else: + print(f'Note: HAF values may differ slightly from target due to ' + f'numerical precision') + + print() + print('Done!') + + +if __name__ == '__main__': + adjust_bed_to_haf() diff --git a/conda_package/pyproject.toml b/conda_package/pyproject.toml index 6cf26b4c8..9b53ac30b 100644 --- a/conda_package/pyproject.toml +++ b/conda_package/pyproject.toml @@ -164,6 +164,7 @@ create_landice_grid_from_generic_mpas_grid = "mpas_tools.landice.create:create_f define_landice_cull_mask = "mpas_tools.landice.cull:define_cull_mask" interpolate_to_mpasli_grid = "mpas_tools.landice.interpolate:interpolate_to_mpasli_grid" mark_domain_boundaries_dirichlet = "mpas_tools.landice.boundary:mark_domain_boundaries_dirichlet" +adjust_bed_to_haf = "mpas_tools.landice.adjust_bed_to_haf:adjust_bed_to_haf" add_critical_land_blockages_to_mask = "mpas_tools.ocean.coastline_alteration:main_add_critical_land_blockages" add_land_locked_cells_to_mask = "mpas_tools.ocean.coastline_alteration:main_add_land_locked_cells_to_mask" widen_transect_edge_masks = "mpas_tools.ocean.coastline_alteration:main_widen_transect_edge_masks" diff --git a/examples/adjust_bed_to_haf_example.py b/examples/adjust_bed_to_haf_example.py new file mode 100755 index 000000000..dda7f2b6a --- /dev/null +++ b/examples/adjust_bed_to_haf_example.py @@ -0,0 +1,195 @@ +#!/usr/bin/env python3 +""" +Example script demonstrating how to use the adjust_bed_to_haf tool +programmatically (as a Python module) rather than from the command line. + +This is useful if you want to integrate the functionality into a larger +workflow or perform additional processing. +""" + +import sys +import numpy as np +import netCDF4 +from shapely.geometry import Point, shape + +# Import the functions from the landice module +# (assuming mpas_tools is installed) +try: + from mpas_tools.landice.adjust_bed_to_haf import ( + load_grounding_line_geojson, + find_cells_in_polygon, + calculate_flotation_thickness, + calculate_required_bed + ) +except ImportError: + print("Error: mpas_tools not installed or not in PYTHONPATH") + print("Install with: cd conda_package && pip install -e .") + sys.exit(1) + + +def example_workflow(): + """ + Example workflow for adjusting bed topography to achieve target HAF. + """ + + # Define input files + mesh_file = "/Users/trhille/Documents/ISMIP7/Antarctica/mesh/4to20km/AIS_4to20km_r03_20260910_ASE_extracted.nc" + geojson_file = "/Users/trhille/Documents/ISMIP7/Antarctica/wild_grounding_lines/Thwaites_GLs_2014_201920/Thwaites_GL_2014_pinning_points.geojson" + output_file = "output_mesh_programmatic.nc" + + # Define parameters + target_haf = 12.0 # meters + rho_ice = 910.0 # kg/m^3 + rho_ocean = 1028.0 # kg/m^3 + sea_level = 0.0 # meters + + print("="*60) + print("Adjust Bed to HAF - Programmatic Example") + print("="*60) + print() + + # Step 1: Load grounding line geometries + print("Step 1: Loading grounding line geometries...") + geometries = load_grounding_line_geojson(geojson_file) + print(f" Loaded {len(geometries)} geometry features") + print() + + # Step 2: Open mesh file and get coordinates + print("Step 2: Opening mesh file...") + with netCDF4.Dataset(mesh_file, 'r') as mesh: + lon_cell = np.degrees(mesh.variables['lonCell'][:]) + lat_cell = np.degrees(mesh.variables['latCell'][:]) + thickness = mesh.variables['thickness'][:] + bed_topo = mesh.variables['bedTopography'][:] + + # Handle potential time dimension + if len(thickness.shape) > 1: + thickness = thickness[0, :] + bed_topo = bed_topo[0, :] + + print(f" Mesh has {len(lon_cell)} cells") + print(f" Thickness range: [{np.min(thickness):.1f}, {np.max(thickness):.1f}] m") + print(f" Bed range: [{np.min(bed_topo):.1f}, {np.max(bed_topo):.1f}] m") + print() + + # Step 3: Find cells within grounding line polygons + print("Step 3: Finding cells within grounding line regions...") + mask = find_cells_in_polygon(lon_cell, lat_cell, geometries) + n_cells_in_region = np.sum(mask) + print(f" Found {n_cells_in_region} cells within grounding line") + print(f" ({100.0 * n_cells_in_region / len(mask):.2f}% of total mesh)") + print() + + if n_cells_in_region == 0: + print("ERROR: No cells found in region!") + return + + # Step 4: Calculate current HAF + print("Step 4: Calculating current height above flotation...") + flotation_thickness = calculate_flotation_thickness( + bed_topo, + sea_level=sea_level, + rho_ice=rho_ice, + rho_ocean=rho_ocean + ) + ice_surface = bed_topo + thickness + flotation_surface = sea_level + flotation_thickness * (1.0 - rho_ice / rho_ocean) + current_haf = ice_surface - flotation_surface + + print(f" Current HAF in region:") + print(f" Mean: {np.mean(current_haf[mask]):.2f} m") + print(f" Min: {np.min(current_haf[mask]):.2f} m") + print(f" Max: {np.max(current_haf[mask]):.2f} m") + print(f" Std: {np.std(current_haf[mask]):.2f} m") + print() + + # Step 5: Calculate new bed elevation + print(f"Step 5: Calculating new bed for target HAF = {target_haf} m...") + new_bed = calculate_required_bed( + thickness[mask], + target_haf, + sea_level=sea_level, + rho_ice=rho_ice, + rho_ocean=rho_ocean + ) + + bed_change = new_bed - bed_topo[mask] + print(f" Required bed change:") + print(f" Mean: {np.mean(bed_change):.2f} m") + print(f" Min: {np.min(bed_change):.2f} m") + print(f" Max: {np.max(bed_change):.2f} m") + print(f" Std: {np.std(bed_change):.2f} m") + print() + + # Step 6: Update mesh and save + print("Step 6: Updating mesh file...") + import shutil + shutil.copy2(mesh_file, output_file) + + with netCDF4.Dataset(output_file, 'r+') as mesh: + bed_var = mesh.variables['bedTopography'] + + # Read current bed + if len(bed_var.shape) > 1: + current_bed = bed_var[0, :] + else: + current_bed = bed_var[:] + + # Update bed in region + current_bed[mask] = new_bed + + # Write back + if len(bed_var.shape) > 1: + bed_var[0, :] = current_bed + else: + bed_var[:] = current_bed + + # Update metadata + from datetime import datetime + timestamp = datetime.now().strftime('%Y-%m-%d %H:%M:%S') + history = f'{timestamp}: Programmatic adjust_bed_to_haf example\n' + if 'history' in mesh.ncattrs(): + history += mesh.getncattr('history') + mesh.setncattr('history', history) + + print(f" Saved to: {output_file}") + print() + + # Step 7: Verify result + print("Step 7: Verifying result...") + with netCDF4.Dataset(output_file, 'r') as mesh: + new_bed_topo = mesh.variables['bedTopography'][:] + if len(new_bed_topo.shape) > 1: + new_bed_topo = new_bed_topo[0, :] + + new_flotation_thickness = calculate_flotation_thickness( + new_bed_topo, + sea_level=sea_level, + rho_ice=rho_ice, + rho_ocean=rho_ocean + ) + new_ice_surface = new_bed_topo + thickness + new_flotation_surface = sea_level + new_flotation_thickness * (1.0 - rho_ice / rho_ocean) + verified_haf = new_ice_surface - new_flotation_surface + + print(f" Verified HAF in region:") + print(f" Mean: {np.mean(verified_haf[mask]):.2f} m") + print(f" Min: {np.min(verified_haf[mask]):.2f} m") + print(f" Max: {np.max(verified_haf[mask]):.2f} m") + print() + + if np.allclose(verified_haf[mask], target_haf, rtol=1e-6): + print("✓ SUCCESS: Target HAF achieved!") + else: + diff = np.abs(verified_haf[mask] - target_haf) + print(f" Max difference from target: {np.max(diff):.6f} m") + print(" (small differences expected due to numerical precision)") + + print() + print("="*60) + print("Example completed successfully!") + print("="*60) + + +if __name__ == '__main__': + example_workflow() diff --git a/landice/mesh_tools_li/README_adjust_bed_to_haf.md b/landice/mesh_tools_li/README_adjust_bed_to_haf.md new file mode 100644 index 000000000..50f958ee4 --- /dev/null +++ b/landice/mesh_tools_li/README_adjust_bed_to_haf.md @@ -0,0 +1,124 @@ +# adjust_bed_to_haf.py + +Adjust bed topography within grounding line regions to achieve a specified height above flotation (HAF). + +## Quick Start + +```bash +python adjust_bed_to_haf.py \ + --mesh input_mesh.nc \ + --geojson grounding_line.geojson \ + --projection ais-bedmap2 \ + --output output_mesh.nc \ + --target-haf 10.0 +``` + +## Usage + +``` +python adjust_bed_to_haf.py [-h] -m FILENAME -g FILENAME -p PROJECTION [-o FILENAME] + [--target-haf TARGET_HAF] [--sea-level SEA_LEVEL] + [--rho-ice RHO_ICE] [--rho-ocean RHO_OCEAN] + [--thickness-var THICKNESS_VAR] [--bed-var BED_VAR] +``` + +### Required Arguments + +- `-m`, `--mesh FILENAME`: MALI mesh file (NetCDF format) +- `-g`, `--geojson FILENAME`: GeoJSON file containing grounding line delineations +- `-p`, `--projection PROJECTION`: Projection of the MALI mesh (see Available Projections below) + +### Optional Arguments + +- `-o`, `--output FILENAME`: Output mesh file (default: modifies input in place) +- `--target-haf TARGET_HAF`: Target height above flotation in meters (default: 10.0) +- `--sea-level SEA_LEVEL`: Sea level in meters (default: 0.0) +- `--rho-ice RHO_ICE`: Ice density in kg/m³ (default: 910.0) +- `--rho-ocean RHO_OCEAN`: Ocean water density in kg/m³ (default: 1028.0) +- `--thickness-var THICKNESS_VAR`: Name of thickness variable (default: 'thickness') +- `--bed-var BED_VAR`: Name of bed topography variable (default: 'bedTopography') +- `--buffer-distance DISTANCE`: Buffer distance for LineString geometries in degrees (default: 0.01, ~1km at poles). Use 0 for Polygon geometries. +- `--thickness-threshold THRESHOLD`: Minimum ice thickness (m) for cells to be modified (default: 1.0). Cells with thickness below this are skipped, preventing unrealistic bed adjustments in ice-free areas. + +## Available Projections + +The script supports the following MALI mesh projections: + +- **ais-bedmap2**: Antarctic BEDMAP2 projection (WGS84 ellipsoid) + - Standard for Antarctica ice sheet models + - Polar stereographic with standard parallel at -71°S +- **ais-bedmap2-sphere**: BEDMAP2 projection on sphere + - Use for coupled MALI-SeaLevelModel simulations +- **gis-bamber**: Greenland Bamber projection +- **gis-gimp**: Greenland GIMP projection +- **latlon**: Standard latitude/longitude (WGS84) + +The script automatically detects the GeoJSON file's CRS and transforms the mesh coordinates to match for accurate spatial queries. + +## Example + +```bash +python adjust_bed_to_haf.py \ + --mesh AIS_4to20km_mesh.nc \ + --geojson Thwaites_GL_2014_pinning_points.geojson \ + --projection ais-bedmap2 \ + --output AIS_4to20km_mesh_adjusted.nc \ + --target-haf 15.0 +``` + +## What is Height Above Flotation? + +Height above flotation (HAF) measures how far ice is from hydrostatic equilibrium: +- **HAF = 0**: Ice at flotation (neutral buoyancy) +- **HAF > 0**: Ice grounded above flotation level +- **HAF < 0**: Ice below flotation (uncommon) + +## Physics + +The tool calculates required bed elevation using: + +``` +bed = sea_level + HAF_target - thickness × (ρ_ice / ρ_ocean) +``` + +This ensures the ice surface achieves the target HAF given the existing thickness. + +## Input Requirements + +### Mesh File +Must contain: +- `xCell`, `yCell`: Cell coordinates in the specified projection (meters) +- `thickness`: Ice thickness (m) +- `bedTopography`: Bed elevation (m) + +### GeoJSON File +Must contain: +- Polygon, LineString, or MultiLineString geometries +- CRS/projection information (typically WGS 84, EPSG:4326) +- The script will auto-detect the GeoJSON CRS and transform coordinates as needed + +**Note:** If your GeoJSON contains LineString or MultiLineString geometries (e.g., grounding lines), the script automatically buffers them to create polygons. Adjust `--buffer-distance` if needed (default: 0.01 degrees, ~1 km). + +## Dependencies + +**Required:** +- xarray +- numpy +- shapely +- pyproj + +**Optional:** +- geopandas (recommended for better performance and CRS detection) + +## Notes + +- The script automatically transforms mesh coordinates to match the GeoJSON CRS for accurate spatial matching +- Ice thickness is preserved; only bed topography is adjusted +- The output file includes updated `history` and `comment` attributes +- If no output file is specified, the input file is modified in place +- The mesh projection must be correctly specified using the `--projection` argument + +## References + +Wild, C. T., et al. (2022). Thwaites Glacier 2014 and 2019/20 grounding line positions. +U.S. Antarctic Program (USAP) Data Center. doi: 10.15784/601499 diff --git a/landice/mesh_tools_li/README_map_regional_to_global.md b/landice/mesh_tools_li/README_map_regional_to_global.md new file mode 100644 index 000000000..17aa73a0f --- /dev/null +++ b/landice/mesh_tools_li/README_map_regional_to_global.md @@ -0,0 +1,202 @@ +# map_regional_to_global_mesh.py + +Map cells from a regional mesh back to a global mesh based on exact coordinate matching. + +## Purpose + +This script is designed for workflows where: +1. A regional mesh was extracted from a larger global mesh (e.g., using compass subdomain extractor) +2. Variables in the regional mesh were modified +3. The modified values need to be copied back to the corresponding cells in the global mesh +4. All other global mesh cells must remain unchanged + +The script uses **exact coordinate matching** (not interpolation) since the regional mesh cells are a subset of the global mesh cells with identical coordinates. + +## Quick Start + +```bash +python map_regional_to_global_mesh.py \ + --regional regional_mesh.nc \ + --global global_mesh.nc \ + --output updated_global_mesh.nc \ + --vars bedTopography thickness +``` + +## Usage + +``` +python map_regional_to_global_mesh.py [-h] -r FILENAME -g FILENAME -o FILENAME + -v VARIABLES [VARIABLES ...] + [--coord-type {spherical,planar}] + [--tolerance TOLERANCE] + [--verify-only] +``` + +### Required Arguments + +- `-r`, `--regional FILENAME`: Regional mesh file (NetCDF format) +- `-g`, `--global FILENAME`: Global mesh file (NetCDF format) +- `-o`, `--output FILENAME`: Output file (must differ from global file) +- `-v`, `--vars VARIABLES`: Variable(s) to copy from regional to global mesh + +### Optional Arguments + +- `--coord-type {spherical,planar}`: Coordinate type (default: spherical) + - `spherical`: Use lonCell/latCell (for spherical meshes) + - `planar`: Use xCell/yCell (for planar meshes) +- `--tolerance TOLERANCE`: Coordinate matching tolerance (default: 1e-10) +- `--verify-only`: Only verify mapping without copying data + +## Example Workflow + +### 1. Extract regional mesh from global mesh +```bash +# Using compass subdomain extractor or similar tool +compass extract-subdomain \ + --mesh global_mesh.nc \ + --output regional_mesh.nc \ + --region amundsen_sea +``` + +### 2. Modify regional mesh +```bash +# Make modifications to the regional mesh +python adjust_bed_to_haf.py \ + --mesh regional_mesh.nc \ + --geojson grounding_line.geojson \ + --target-haf 10.0 +``` + +### 3. Map back to global mesh +```bash +# Copy modified variables back to global mesh +python map_regional_to_global_mesh.py \ + --regional regional_mesh.nc \ + --global global_mesh.nc \ + --output global_mesh_updated.nc \ + --vars bedTopography +``` + +## Use Cases + +### Copy modified bed topography +```bash +python map_regional_to_global_mesh.py \ + -r ASE_mesh_modified.nc \ + -g AIS_mesh.nc \ + -o AIS_mesh_updated.nc \ + -v bedTopography +``` + +### Copy multiple variables +```bash +python map_regional_to_global_mesh.py \ + -r ASE_mesh_modified.nc \ + -g AIS_mesh.nc \ + -o AIS_mesh_updated.nc \ + -v bedTopography thickness surfaceSpeed temperature +``` + +### Verify mapping without copying +```bash +python map_regional_to_global_mesh.py \ + -r ASE_mesh.nc \ + -g AIS_mesh.nc \ + -o dummy_output.nc \ + -v bedTopography \ + --verify-only +``` + +### Use planar coordinates +```bash +python map_regional_to_global_mesh.py \ + -r regional_planar.nc \ + -g global_planar.nc \ + -o global_updated.nc \ + -v thickness \ + --coord-type planar +``` + +## How It Works + +1. **Load both meshes**: Opens regional and global mesh files +2. **Extract coordinates**: Gets lonCell/latCell (spherical) or xCell/yCell (planar) +3. **Build spatial index**: Creates a KDTree for fast nearest-neighbor search +4. **Build mapping**: Uses KDTree to efficiently find matching global cells for all regional cells +5. **Verify mapping**: Checks that all matches are within tolerance +6. **Copy variables**: Updates specified variables in global mesh for matched cells only +7. **Write output**: Saves updated global mesh to new file + +**Performance**: Uses scipy's KDTree for O(n log m) spatial search instead of O(n*m) brute force. Typical mapping time: seconds instead of hours for large meshes. + +## Coordinate Matching + +The script uses exact coordinate matching with a small tolerance (default 1e-10) to account for floating-point precision: + +- **Spherical**: Matches based on (lonCell, latCell) in radians +- **Planar**: Matches based on (xCell, yCell) in meters + +For each regional cell, the script finds the global cell with minimum distance. If the distance exceeds the tolerance, the regional cell is marked as unmatched. + +## Variable Handling + +The script automatically handles variables with different dimensions: + +- **1D** `(nCells)`: Direct cell-to-cell copy +- **2D** `(nCells, nVertLevels)`: Copies all vertical levels +- **2D** `(Time, nCells)`: Copies all time levels +- **3D** `(Time, nCells, nVertLevels)`: Copies all time and vertical levels + +## Error Checking + +The script includes several safety checks: + +- Verifies output file differs from global file (prevents accidental overwrite) +- Checks that all requested variables exist in both meshes +- Reports unmatched regional cells +- Verifies coordinate matching is within tolerance +- Updates file metadata with history and comments + +## Output + +The script reports: +- Number of cells in each mesh +- Number of matched cells +- Number of unmatched regional cells (if any) +- Coordinate matching errors (mean and max) +- Variables successfully copied + +## Notes + +- The output file must differ from the global input file (safety feature) +- Unmatched regional cells are skipped (with warning) +- The global mesh is loaded into memory, so ensure sufficient RAM for large meshes +- Original global file is never modified (read-only) +- All cells in global mesh not matching regional cells remain unchanged + +## Dependencies + +Required Python packages: +- xarray +- numpy +- scipy (for KDTree spatial indexing) + +## Common Issues + +### "No matching cells found" +- Check that regional mesh was actually extracted from this global mesh +- Verify coordinate type (spherical vs planar) +- Try increasing tolerance + +### "WARNING: Some matches exceed tolerance" +- Regional and global meshes may have been modified with different precision +- Try increasing tolerance with `--tolerance 1e-6` + +### Memory issues +- For very large global meshes, ensure sufficient RAM +- Consider processing variables one at a time + +## See Also + +- `adjust_bed_to_haf.py`: Tool for modifying bed topography +- `compass extract-subdomain`: Tool for extracting regional meshes diff --git a/landice/mesh_tools_li/adjust_bed_to_haf.py b/landice/mesh_tools_li/adjust_bed_to_haf.py new file mode 100755 index 000000000..334b12b5b --- /dev/null +++ b/landice/mesh_tools_li/adjust_bed_to_haf.py @@ -0,0 +1,573 @@ +#!/usr/bin/env python +''' +Adjust bed topography within grounding line regions to achieve a specified +height above flotation (HAF). + +Height above flotation is defined as: + HAF = ice_surface - flotation_surface +where flotation_surface is the surface elevation at which ice would be in +hydrostatic equilibrium with ocean water. + +For a given ice thickness and target HAF, the required bed elevation is +calculated from the flotation condition. The tool loads grounding line +delineations from a GeoJSON file, identifies mesh cells within those regions, +and adjusts the bed topography to achieve the target HAF while preserving +ice thickness. + +The modified mesh is written to the required output file (--output), or if +not specified, the input file is modified in place. + +Trevor Hillebrand, 2026 +''' + +import sys +import numpy as np +import xarray as xr +from argparse import ArgumentParser, RawDescriptionHelpFormatter +from datetime import datetime +from shapely.geometry import Point, shape +from pyproj import Transformer, CRS + +# Try to import mpas_tools for proper NetCDF writing +try: + from mpas_tools.io import write_netcdf + HAS_MPAS_TOOLS = True +except ImportError: + HAS_MPAS_TOOLS = False + import warnings + warnings.warn('mpas_tools not available. Output may not be in proper MPAS format.') + +try: + import geopandas as gpd + HAS_GEOPANDAS = True +except ImportError: + HAS_GEOPANDAS = False + import json + +# Define available projections (from mpas_tools.landice.projections) +PROJECTIONS = { + 'gis-bamber': ( + '+proj=stere +lat_ts=71.0 +lat_0=90 +lon_0=321.0 +k_0=1.0 ' + '+x_0=800000.0 +y_0=3400000.0 +ellps=WGS84' + ), + 'gis-gimp': ( + '+proj=stere +lat_ts=70.0 +lat_0=90 +lon_0=315.0 +k_0=1.0 +x_0=0.0 ' + '+y_0=0.0 +ellps=WGS84' + ), + 'ais-bedmap2': ( + '+proj=stere +lat_ts=-71.0 +lat_0=-90 +lon_0=0.0 +k_0=1.0 +x_0=0.0 ' + '+y_0=0.0 +ellps=WGS84' + ), + 'ais-bedmap2-sphere': ( + '+proj=stere +lat_ts=-71.0 +lat_0=-90 +lon_0=0.0 +k_0=1.0 +x_0=0.0 ' + '+y_0=0.0 +ellps=sphere' + ), + 'latlon': '+proj=longlat +ellps=WGS84', +} + + +def parse_args(): + parser = ArgumentParser(description=__doc__, + formatter_class=RawDescriptionHelpFormatter) + parser.add_argument('-m', '--mesh', dest='mesh_file', required=True, + metavar='FILENAME', + help='MALI mesh file (NetCDF format)') + parser.add_argument('-g', '--geojson', dest='geojson_file', required=True, + metavar='FILENAME', + help='GeoJSON file containing grounding line delineations') + parser.add_argument('-o', '--output', dest='output_file', + metavar='FILENAME', + help='Output mesh file. If not specified, modifies input ' + 'file in place.') + parser.add_argument('--target-haf', dest='target_haf', type=float, + default=10.0, + help='Target height above flotation in meters ' + '(default: 10.0)') + parser.add_argument('--sea-level', dest='sea_level', type=float, + default=0.0, + help='Sea level in meters (default: 0.0)') + parser.add_argument('--rho-ice', dest='rho_ice', type=float, + default=910.0, + help='Ice density in kg/m^3 (default: 910.0)') + parser.add_argument('--rho-ocean', dest='rho_ocean', type=float, + default=1028.0, + help='Ocean water density in kg/m^3 (default: 1028.0)') + parser.add_argument('--thickness-var', dest='thickness_var', + default='thickness', + help='Name of thickness variable in mesh file ' + '(default: thickness)') + parser.add_argument('--bed-var', dest='bed_var', default='bedTopography', + help='Name of bed topography variable in mesh file ' + '(default: bedTopography)') + parser.add_argument('-p', '--projection', dest='projection', + choices=list(PROJECTIONS.keys()), + required=True, + help='Projection of the MALI mesh. Available: ' + + ', '.join(PROJECTIONS.keys())) + parser.add_argument('--buffer-distance', dest='buffer_distance', + type=float, default=0.01, + help='Buffer distance for LineString geometries in ' + 'degrees (default: 0.01, ~1km). Use 0 for Polygons.') + parser.add_argument('--thickness-threshold', dest='thickness_threshold', + type=float, default=1.0, + help='Minimum ice thickness (m) for cells to be modified. ' + 'Cells with thickness below this are skipped ' + '(default: 1.0)') + + return parser.parse_args() + + +def load_grounding_line_geojson(geojson_file, buffer_distance=0.0): + ''' + Load grounding line geometries from a GeoJSON file and extract CRS. + + Parameters + ---------- + geojson_file : str + Path to the GeoJSON file containing grounding line delineations + buffer_distance : float + Buffer distance to apply to LineString geometries (in CRS units) + + Returns + ------- + geometries : list + List of shapely geometry objects (buffered if LineStrings) + crs : pyproj.CRS or None + CRS of the GeoJSON file + ''' + if HAS_GEOPANDAS: + # Use geopandas if available + gdf = gpd.read_file(geojson_file) + geometries = gdf.geometry.tolist() + crs = gdf.crs if gdf.crs is not None else CRS.from_epsg(4326) + else: + # Fall back to json + shapely + with open(geojson_file, 'r') as f: + data = json.load(f) + + geometries = [] + for feature in data['features']: + geom = shape(feature['geometry']) + geometries.append(geom) + + # Try to extract CRS from GeoJSON + if 'crs' in data and 'properties' in data['crs']: + crs_name = data['crs']['properties'].get('name', '') + if 'EPSG' in crs_name or 'epsg' in crs_name: + # Extract EPSG code + epsg_code = int(crs_name.split(':')[-1]) + crs = CRS.from_epsg(epsg_code) + else: + # Default to WGS84 + crs = CRS.from_epsg(4326) + else: + # No CRS specified, assume WGS84 (standard for GeoJSON) + crs = CRS.from_epsg(4326) + + # Buffer LineString/MultiLineString geometries to create polygons + if buffer_distance > 0: + buffered_geometries = [] + for geom in geometries: + geom_type = geom.geom_type + if geom_type in ['LineString', 'MultiLineString']: + # Buffer the line to create a polygon + buffered_geom = geom.buffer(buffer_distance) + buffered_geometries.append(buffered_geom) + else: + # Keep polygons as-is + buffered_geometries.append(geom) + geometries = buffered_geometries + + return geometries, crs + + +def transform_mesh_coords(x_cell, y_cell, mesh_proj_str, target_crs): + ''' + Transform mesh coordinates from mesh projection to target CRS. + + Parameters + ---------- + x_cell : ndarray + X coordinates of cell centers in mesh projection + y_cell : ndarray + Y coordinates of cell centers in mesh projection + mesh_proj_str : str + Proj4 string defining the mesh projection + target_crs : pyproj.CRS + Target coordinate reference system + + Returns + ------- + x_transformed : ndarray + Transformed x coordinates + y_transformed : ndarray + Transformed y coordinates + ''' + mesh_crs = CRS.from_proj4(mesh_proj_str) + transformer = Transformer.from_crs(mesh_crs, target_crs, always_xy=True) + x_transformed, y_transformed = transformer.transform(x_cell, y_cell) + return x_transformed, y_transformed + + +def find_cells_in_polygon(x_cell, y_cell, geometries): + ''' + Find mesh cells that fall within any of the provided polygons. + + Parameters + ---------- + x_cell : ndarray + X coordinates of cell centers (in same CRS as geometries) + y_cell : ndarray + Y coordinates of cell centers (in same CRS as geometries) + geometries : list + List of shapely geometry objects + + Returns + ------- + mask : ndarray (bool) + Boolean mask indicating which cells are inside the polygons + ''' + n_cells = len(x_cell) + mask = np.zeros(n_cells, dtype=bool) + + print(f'Checking {n_cells} cells against {len(geometries)} geometries...') + + for i, (x, y) in enumerate(zip(x_cell, y_cell)): + if i % 10000 == 0: + print(f' Processed {i}/{n_cells} cells...') + + point = Point(x, y) + for geom in geometries: + if geom.contains(point) or geom.intersects(point): + mask[i] = True + break + + print(f'Found {np.sum(mask)} cells within polygons') + return mask + + +def calculate_flotation_thickness(bed_elevation, sea_level=0.0, + rho_ice=910.0, rho_ocean=1028.0): + ''' + Calculate the ice thickness at which ice would be in flotation. + + Parameters + ---------- + bed_elevation : ndarray + Bed topography (positive up, m) + sea_level : float + Sea level (m), default 0.0 + rho_ice : float + Ice density (kg/m^3), default 910.0 + rho_ocean : float + Ocean water density (kg/m^3), default 1028.0 + + Returns + ------- + flotation_thickness : ndarray + Ice thickness at flotation (m) + ''' + # For ice to float: rho_ice * thickness = rho_ocean * draft + # where draft = sea_level - bed_elevation + # Therefore: thickness_flotation = (rho_ocean / rho_ice) * draft + + draft = sea_level - bed_elevation + flotation_thickness = (rho_ocean / rho_ice) * draft + + # Thickness must be positive + flotation_thickness = np.maximum(flotation_thickness, 0.0) + + return flotation_thickness + + +def calculate_required_bed(thickness, target_haf, sea_level=0.0, + rho_ice=910.0, rho_ocean=1028.0): + ''' + Calculate the bed elevation required to achieve a target height above + flotation for a given ice thickness. + + Parameters + ---------- + thickness : ndarray + Ice thickness (m) + target_haf : float + Target height above flotation (m) + sea_level : float + Sea level (m), default 0.0 + rho_ice : float + Ice density (kg/m^3), default 910.0 + rho_ocean : float + Ocean water density (kg/m^3), default 1028.0 + + Returns + ------- + bed_elevation : ndarray + Required bed elevation (m) + ''' + # At flotation: + # ice_surface_flotation = sea_level + thickness * (1 - rho_ice/rho_ocean) + # + # With target HAF: + # ice_surface = ice_surface_flotation + target_haf + # + # Since ice_surface = bed + thickness: + # bed = ice_surface - thickness + # = ice_surface_flotation + target_haf - thickness + # = sea_level + thickness * (1 - rho_ice/rho_ocean) + target_haf + # - thickness + # = sea_level + target_haf - thickness * rho_ice/rho_ocean + + bed_elevation = sea_level + target_haf - thickness * (rho_ice / rho_ocean) + + return bed_elevation + + +def main(): + ''' + Main function to adjust bed topography to achieve target height above + flotation within grounding line regions. + ''' + + args = parse_args() + + # Validate inputs + if not HAS_GEOPANDAS: + print('Warning: geopandas not available, using json + shapely instead') + + print('\n** Adjusting bed topography to achieve target height above ' + 'flotation **') + print(f'Input mesh file: {args.mesh_file}') + print(f'Grounding line GeoJSON: {args.geojson_file}') + print(f'Mesh projection: {args.projection}') + print(f'Target HAF: {args.target_haf} m') + print(f'Sea level: {args.sea_level} m') + print(f'Ice density: {args.rho_ice} kg/m^3') + print(f'Ocean density: {args.rho_ocean} kg/m^3') + print() + + # Load grounding line polygons and get CRS + print('Loading grounding line geometries...') + geometries, geojson_crs = load_grounding_line_geojson( + args.geojson_file, args.buffer_distance + ) + print(f'Loaded {len(geometries)} geometry features') + # Report geometry types + geom_types = set(geom.geom_type for geom in geometries) + print(f'Geometry types: {", ".join(geom_types)}') + print(f'GeoJSON CRS: {geojson_crs}') + if args.buffer_distance > 0: + print(f'Applied buffer distance: {args.buffer_distance} degrees ' + f'(~{args.buffer_distance * 111:.1f} km)') + print() + + # Open mesh file with xarray + print('Opening mesh file...') + ds = xr.open_dataset(args.mesh_file) + + # Read mesh coordinates + print('Reading mesh coordinates...') + x_cell = ds['xCell'].values + y_cell = ds['yCell'].values + n_cells = len(x_cell) + print(f'Mesh has {n_cells} cells') + print(f'Mesh coordinate range:') + print(f' X: [{np.min(x_cell):.1f}, {np.max(x_cell):.1f}] m') + print(f' Y: [{np.min(y_cell):.1f}, {np.max(y_cell):.1f}] m') + print() + + # Transform mesh coordinates to GeoJSON CRS + print('Transforming mesh coordinates to GeoJSON CRS...') + mesh_proj_str = PROJECTIONS[args.projection] + x_transformed, y_transformed = transform_mesh_coords( + x_cell, y_cell, mesh_proj_str, geojson_crs + ) + print(f'Transformed coordinate range:') + print(f' X: [{np.min(x_transformed):.6f}, {np.max(x_transformed):.6f}]') + print(f' Y: [{np.min(y_transformed):.6f}, {np.max(y_transformed):.6f}]') + print() + + # Find cells within grounding line polygons + print('Identifying cells within grounding line regions...') + mask = find_cells_in_polygon(x_transformed, y_transformed, geometries) + print() + + if np.sum(mask) == 0: + print('ERROR: No cells found within grounding line polygons!') + print('Check that:') + print(' 1. Mesh projection is correct (specified: {})'.format( + args.projection)) + print(' 2. GeoJSON geometries overlap with mesh extent') + print(' 3. GeoJSON CRS was correctly detected') + sys.exit(1) + + # Read thickness and bed topography + print('Reading ice thickness and bed topography...') + if args.thickness_var not in ds.variables: + print(f'ERROR: Variable "{args.thickness_var}" not found in mesh file') + print(f'Available variables: {list(ds.variables.keys())}') + sys.exit(1) + + if args.bed_var not in ds.variables: + print(f'ERROR: Variable "{args.bed_var}" not found in mesh file') + print(f'Available variables: {list(ds.variables.keys())}') + sys.exit(1) + + thickness = ds[args.thickness_var].values + bed_topo = ds[args.bed_var].values + + # Handle potential time dimension + if len(thickness.shape) > 1: + # Assume time is first dimension + thickness = thickness[0, :] + bed_topo = bed_topo[0, :] + has_time_dim = True + else: + has_time_dim = False + + print(f'Thickness range: [{np.min(thickness):.2f}, ' + f'{np.max(thickness):.2f}] m') + print(f'Bed topography range: [{np.min(bed_topo):.2f}, ' + f'{np.max(bed_topo):.2f}] m') + print() + + # Apply thickness threshold to mask + print(f'Applying thickness threshold ({args.thickness_threshold} m)...') + ice_mask = thickness >= args.thickness_threshold + print(f' Cells in grounding line region: {np.sum(mask)}') + print(f' Cells with ice >= {args.thickness_threshold} m: {np.sum(ice_mask)}') + + # Combine spatial mask with thickness mask + mask = mask & ice_mask + print(f' Cells to modify (region AND ice present): {np.sum(mask)}') + print() + + if np.sum(mask) == 0: + print('ERROR: No cells with ice found within grounding line polygons!') + print('Check that:') + print(f' 1. Thickness threshold ({args.thickness_threshold} m) is appropriate') + print(' 2. Ice is present in the grounding line region') + sys.exit(1) + + # Calculate current HAF in the region + print('Calculating current height above flotation...') + flotation_thickness = calculate_flotation_thickness( + bed_topo, + sea_level=args.sea_level, + rho_ice=args.rho_ice, + rho_ocean=args.rho_ocean + ) + ice_surface = bed_topo + thickness + flotation_surface = args.sea_level + flotation_thickness * ( + 1.0 - args.rho_ice / args.rho_ocean + ) + current_haf = ice_surface - flotation_surface + + print(f'Current HAF in region (mean): {np.mean(current_haf[mask]):.2f} m') + print(f'Current HAF in region (min): {np.min(current_haf[mask]):.2f} m') + print(f'Current HAF in region (max): {np.max(current_haf[mask]):.2f} m') + print() + + # Calculate new bed elevation + print('Calculating new bed topography...') + new_bed = calculate_required_bed( + thickness[mask], + args.target_haf, + sea_level=args.sea_level, + rho_ice=args.rho_ice, + rho_ocean=args.rho_ocean + ) + + print(f'New bed range in region: [{np.min(new_bed):.2f}, ' + f'{np.max(new_bed):.2f}] m') + bed_change = new_bed - bed_topo[mask] + print(f'Bed change (mean): {np.mean(bed_change):.2f} m') + print(f'Bed change (min): {np.min(bed_change):.2f} m') + print(f'Bed change (max): {np.max(bed_change):.2f} m') + print() + + # Update bed topography in dataset + print('Updating bed topography...') + if has_time_dim: + ds[args.bed_var].values[0, mask] = new_bed + else: + ds[args.bed_var].values[mask] = new_bed + + # Update comment attribute + comment = ( + f'Bed topography adjusted within grounding line regions ' + f'to achieve target HAF = {args.target_haf} m. ' + f'Modified {np.sum(mask)} cells with ice thickness >= ' + f'{args.thickness_threshold} m using grounding line data from ' + f'{args.geojson_file}.' + ) + if 'comment' in ds.attrs: + existing_comment = ds.attrs['comment'] + comment = f'{comment} {existing_comment}' + ds.attrs['comment'] = comment + + # Write output + output_file = args.output_file if args.output_file else args.mesh_file + print(f'Writing to {output_file}...') + + if HAS_MPAS_TOOLS: + # Use mpas_tools for proper MPAS format + # write_netcdf will automatically add history with command line + write_netcdf(ds, output_file, format='NETCDF3_64BIT_DATA') + else: + # Fallback to xarray if mpas_tools not available + # Manually update history + timestamp = datetime.now().strftime('%Y-%m-%d %H:%M:%S') + command_line = ' '.join(sys.argv) + if 'history' in ds.attrs: + history = ds.attrs['history'] + history = f'{timestamp}: {command_line}\n{history}' + else: + history = f'{timestamp}: {command_line}' + ds.attrs['history'] = history + + encoding = {var: {'_FillValue': None} for var in ds.data_vars} + ds.to_netcdf(output_file, encoding=encoding) + ds.close() + + print('Successfully updated mesh file!') + print() + + # Verify the result + print('Verifying updated HAF...') + ds_verify = xr.open_dataset(output_file) + + verified_bed = ds_verify[args.bed_var].values + if has_time_dim: + verified_bed = verified_bed[0, :] + + flotation_thickness = calculate_flotation_thickness( + verified_bed, + sea_level=args.sea_level, + rho_ice=args.rho_ice, + rho_ocean=args.rho_ocean + ) + ice_surface = verified_bed + thickness + flotation_surface = args.sea_level + flotation_thickness * ( + 1.0 - args.rho_ice / args.rho_ocean + ) + new_haf = ice_surface - flotation_surface + + print(f'New HAF in region (mean): {np.mean(new_haf[mask]):.2f} m') + print(f'New HAF in region (min): {np.min(new_haf[mask]):.2f} m') + print(f'New HAF in region (max): {np.max(new_haf[mask]):.2f} m') + print() + + if np.allclose(new_haf[mask], args.target_haf, rtol=1e-6): + print(f'✓ Successfully achieved target HAF = {args.target_haf} m') + else: + print(f'Note: HAF values may differ slightly from target due to ' + f'numerical precision') + + ds_verify.close() + + print() + print('Done!') + + +if __name__ == '__main__': + main() diff --git a/landice/mesh_tools_li/convert_to_budd_mu.py b/landice/mesh_tools_li/convert_to_budd_mu.py new file mode 100644 index 000000000..2126a10c1 --- /dev/null +++ b/landice/mesh_tools_li/convert_to_budd_mu.py @@ -0,0 +1,189 @@ +from argparse import ArgumentParser, RawDescriptionHelpFormatter +from mpas_tools.io import write_netcdf +import xarray as xr +import numpy as np +from scipy.spatial import cKDTree + + +parser = ArgumentParser(description=__doc__, + formatter_class=RawDescriptionHelpFormatter) +parser.add_argument('-i', '--input_file', dest='input_file', required=True, + metavar='FILENAME', + help='MALI file (NetCDF format) containing muFriction' + ' for Weertman friction law') +parser.add_argument('-o', '--output_file', dest='output_file', required=True, + metavar='FILENAME', + help='Destination file for converted muFriction' + ' for Budd friction law') +args = parser.parse_args() + +input_filename = args.input_file +output_filename = args.output_file +assert input_filename != output_filename, \ + "Input file and output file must have different names!" + +mu_friction_max = 100 # Replace with your desired threshold +fill_method = "idw" # Either "idw" or "nearest" +idw_neighbors = 8 +idw_power = 2.0 + +def fill_large_values( + field, + x_cell, + y_cell, + maximum, + method="idw", + neighbors=8, + power=2.0, +): + """ + Replace values greater than `maximum` using spatial interpolation. + + Parameters + ---------- + field : xr.DataArray + Field containing an ``nCells`` dimension. + x_cell, y_cell : xr.DataArray + MPAS cell-center coordinates. + maximum : float + Values greater than this threshold will be replaced. + method : {"idw", "nearest"} + Spatial filling method. + neighbors : int + Number of neighbors used for IDW. + power : float + Distance exponent used for IDW. + """ + if "nCells" not in field.dims: + raise ValueError("field must contain an 'nCells' dimension") + + if method not in {"idw", "nearest"}: + raise ValueError("method must be either 'idw' or 'nearest'") + + # Put nCells last, allowing the function to handle optional Time or + # other leading dimensions. + original_dims = field.dims + field_work = field.transpose( + *[dim for dim in field.dims if dim != "nCells"], + "nCells", + ) + + values = np.asarray(field_work.values).copy() + original_shape = values.shape + values_2d = values.reshape(-1, original_shape[-1]) + + coordinates = np.column_stack( + [np.asarray(x_cell.values), np.asarray(y_cell.values)] + ) + + for row in values_2d: + target_mask = np.isfinite(row) & (row > maximum) + donor_mask = np.isfinite(row) & (row <= maximum) + + if not np.any(target_mask): + continue + + if not np.any(donor_mask): + raise ValueError( + "No valid muFriction cells are available for interpolation" + ) + + donor_values = row[donor_mask] + donor_coordinates = coordinates[donor_mask] + target_coordinates = coordinates[target_mask] + + tree = cKDTree(donor_coordinates) + + if method == "nearest": + _, indices = tree.query(target_coordinates, k=1) + row[target_mask] = donor_values[indices] + + else: + neighbor_count = min(neighbors, donor_values.size) + distances, indices = tree.query( + target_coordinates, + k=neighbor_count, + ) + + # Ensure two-dimensional arrays when k=1. + distances = np.atleast_2d(distances) + indices = np.atleast_2d(indices) + + if distances.shape[0] != target_coordinates.shape[0]: + distances = distances.T + indices = indices.T + + # A tiny distance floor prevents division by zero. + weights = 1.0 / np.maximum(distances, 1.0e-12) ** power + neighbor_values = donor_values[indices] + + row[target_mask] = np.sum( + weights * neighbor_values, axis=1 + ) / np.sum(weights, axis=1) + + filled = xr.DataArray( + values_2d.reshape(original_shape), + dims=field_work.dims, + coords=field_work.coords, + attrs=field.attrs, + name=field.name, + ) + + return filled.transpose(*original_dims) + + +# ---------------------------------------------------------------------- +# Load and process the data +# ---------------------------------------------------------------------- + +ds = xr.open_dataset(input_filename) + +thickness = ds["thickness"] +bed_topography = ds["bedTopography"] +mu_friction = ds["muFriction"] + +# Compute max(0, -(rho_ocean / rho_ice) * bedTopography). +rho_ocean = 1028.0 +rho_ice = 910.0 +gravity = 9.81 + +inner_term = -(rho_ocean / rho_ice) * bed_topography +inner_max = xr.where(inner_term > 0.0, inner_term, 0.0) + +denominator = 1.e-3 * ( # account for kPa units in Albany + rho_ice * gravity * thickness + - rho_ice * gravity * inner_max +) + +denominator = xr.where( + denominator > 1.0e-12, + denominator, + 1.0e-12, +) + +mu_friction /= denominator +# First fill anomalously large muFriction values. +mu_friction_filled = fill_large_values( + field=mu_friction, + x_cell=ds["xCell"], + y_cell=ds["yCell"], + maximum=mu_friction_max, + method=fill_method, + neighbors=idw_neighbors, + power=idw_power, +) + +# Normalize the filled field and update the dataset. +ds["muFriction"] = mu_friction_filled +ds["muFriction"].attrs.update(mu_friction.attrs) + +#ds.to_netcdf( +# output_filename, +# engine="netcdf4", +# format="NETCDF3_64BIT_OFFSET" +#) + +write_netcdf(ds, output_filename, format='NETCDF3_64BIT_DATA') +ds.close() + +print(f"Processing complete. Output saved to {output_filename!r}.") diff --git a/landice/mesh_tools_li/map_regional_to_global_mesh.py b/landice/mesh_tools_li/map_regional_to_global_mesh.py new file mode 100755 index 000000000..cccccc62a --- /dev/null +++ b/landice/mesh_tools_li/map_regional_to_global_mesh.py @@ -0,0 +1,410 @@ +#!/usr/bin/env python +''' +Map cells from a regional mesh back to a global mesh based on exact coordinate +matching. + +This script is designed for workflows where a regional mesh has been extracted +from a larger global mesh (e.g., using compass subdomain extractor), modified, +and then needs to be mapped back to the global mesh. Since the regional mesh +cells are a subset of the global mesh cells with identical coordinates, this +script creates an exact coordinate-based mapping rather than using interpolation. + +The script: +1. Loads both regional and global meshes +2. Finds cells in the global mesh that match regional mesh coordinates exactly +3. Copies specified variables from regional to global mesh for matched cells +4. Leaves all other global mesh cells unchanged + +Trevor Hillebrand, 2026 +''' + +import sys +import numpy as np +import xarray as xr +from argparse import ArgumentParser, RawDescriptionHelpFormatter +from datetime import datetime +from scipy.spatial import KDTree + +# Try to import mpas_tools for proper NetCDF writing +try: + from mpas_tools.io import write_netcdf + HAS_MPAS_TOOLS = True +except ImportError: + HAS_MPAS_TOOLS = False + import warnings + warnings.warn('mpas_tools not available. Output may not be in proper MPAS format.') + + +def parse_args(): + parser = ArgumentParser(description=__doc__, + formatter_class=RawDescriptionHelpFormatter) + parser.add_argument('-r', '--regional', dest='regional_file', required=True, + metavar='FILENAME', + help='Regional mesh file (NetCDF format)') + parser.add_argument('-g', '--global', dest='global_file', required=True, + metavar='FILENAME', + help='Global mesh file (NetCDF format)') + parser.add_argument('-o', '--output', dest='output_file', required=True, + metavar='FILENAME', + help='Output mesh file (must differ from global file)') + parser.add_argument('-v', '--vars', dest='variables', required=True, + nargs='+', + help='Variable(s) to copy from regional to global mesh') + parser.add_argument('--coord-type', dest='coord_type', + choices=['spherical', 'planar'], default='spherical', + help='Coordinate type: spherical (lon/lat) or planar ' + '(x/y) (default: spherical)') + parser.add_argument('--tolerance', dest='tolerance', type=float, + default=1e-10, + help='Coordinate matching tolerance (default: 1e-10)') + parser.add_argument('--verify-only', dest='verify_only', action='store_true', + help='Only verify mapping without copying data') + + return parser.parse_args() + + +def build_coordinate_mapping(regional_coords, global_coords, tolerance=1e-10): + ''' + Build a mapping from regional mesh cell indices to global mesh cell indices + based on exact coordinate matching using a KDTree for fast spatial search. + + Parameters + ---------- + regional_coords : ndarray, shape (n_regional, 2) + Regional mesh coordinates (x, y) or (lon, lat) + global_coords : ndarray, shape (n_global, 2) + Global mesh coordinates (x, y) or (lon, lat) + tolerance : float + Maximum distance for considering coordinates as matching + + Returns + ------- + mapping : dict + Dictionary mapping regional cell index to global cell index + unmatched_regional : list + List of regional cell indices that have no match in global mesh + ''' + n_regional = regional_coords.shape[0] + n_global = global_coords.shape[0] + + print(f'Building coordinate mapping...') + print(f' Regional mesh: {n_regional} cells') + print(f' Global mesh: {n_global} cells') + print(f' Tolerance: {tolerance}') + + # Build KDTree for fast spatial search + print(f' Building spatial index (KDTree)...') + tree = KDTree(global_coords) + print(f' Spatial index built!') + + # Query all regional coordinates at once (vectorized) + print(f' Querying nearest neighbors...') + distances, indices = tree.query(regional_coords, k=1) + print(f' Query complete!') + + # Build mapping dictionary + mapping = {} + unmatched_regional = [] + + for i in range(n_regional): + if distances[i] <= tolerance: + mapping[i] = indices[i] + else: + unmatched_regional.append(i) + + print(f' Matched cells: {len(mapping)} / {n_regional}') + if unmatched_regional: + print(f' WARNING: {len(unmatched_regional)} regional cells have no ' + f'match in global mesh') + + return mapping, unmatched_regional + + +def verify_mapping(mapping, regional_coords, global_coords, tolerance): + ''' + Verify that the mapping is correct by checking coordinate differences. + + Parameters + ---------- + mapping : dict + Regional to global cell index mapping + regional_coords : ndarray + Regional mesh coordinates + global_coords : ndarray + Global mesh coordinates + tolerance : float + Matching tolerance + + Returns + ------- + max_error : float + Maximum coordinate error + ''' + print('Verifying mapping...') + + errors = [] + for regional_idx, global_idx in mapping.items(): + regional_coord = regional_coords[regional_idx, :] + global_coord = global_coords[global_idx, :] + error = np.sqrt(np.sum((regional_coord - global_coord)**2)) + errors.append(error) + + errors = np.array(errors) + max_error = np.max(errors) + mean_error = np.mean(errors) + + print(f' Mean coordinate error: {mean_error:.2e}') + print(f' Max coordinate error: {max_error:.2e}') + print(f' Tolerance: {tolerance:.2e}') + + if max_error > tolerance: + print(' WARNING: Some matches exceed tolerance!') + return False + else: + print(' ✓ All matches within tolerance') + return True + + +def copy_variables(regional_ds, global_ds, mapping, variables): + ''' + Copy specified variables from regional to global mesh using the mapping. + + Parameters + ---------- + regional_ds : xarray.Dataset + Regional mesh dataset + global_ds : xarray.Dataset + Global mesh dataset + mapping : dict + Regional to global cell index mapping + variables : list + List of variable names to copy + + Returns + ------- + global_ds : xarray.Dataset + Modified global mesh dataset + ''' + print(f'Copying variables: {", ".join(variables)}') + + for var in variables: + print(f' Processing {var}...') + + if var not in regional_ds.variables: + print(f' ERROR: Variable "{var}" not found in regional mesh') + continue + + if var not in global_ds.variables: + print(f' ERROR: Variable "{var}" not found in global mesh') + continue + + regional_var = regional_ds[var].values + + # Get variable shape + var_shape = regional_var.shape + print(f' Regional shape: {var_shape}') + print(f' Global shape: {global_ds[var].shape}') + + # Make a copy of the global variable data to modify + global_var = global_ds[var].values.copy() + + # Handle different dimensionalities + if len(var_shape) == 1: + # 1D variable (nCells) + for regional_idx, global_idx in mapping.items(): + global_var[global_idx] = regional_var[regional_idx] + + elif len(var_shape) == 2: + # 2D variable, could be (Time, nCells) or (nCells, nVertLevels) + dim_names = regional_ds[var].dims + + if 'nCells' in dim_names: + cell_dim = dim_names.index('nCells') + + if cell_dim == 0: + # (nCells, nVertLevels) + for regional_idx, global_idx in mapping.items(): + global_var[global_idx, :] = regional_var[regional_idx, :] + else: + # (Time, nCells) or (nVertLevels, nCells) + for regional_idx, global_idx in mapping.items(): + global_var[:, global_idx] = regional_var[:, regional_idx] + else: + print(f' WARNING: Cannot determine cell dimension for {var}') + continue + + elif len(var_shape) == 3: + # 3D variable (Time, nCells, nVertLevels) + dim_names = regional_ds[var].dims + + if 'nCells' in dim_names: + cell_dim = dim_names.index('nCells') + + if cell_dim == 1: + # (Time, nCells, nVertLevels) + for regional_idx, global_idx in mapping.items(): + global_var[:, global_idx, :] = \ + regional_var[:, regional_idx, :] + else: + print(f' WARNING: Unexpected cell dimension position ' + f'for {var}') + continue + else: + print(f' WARNING: Cannot determine cell dimension for {var}') + continue + + else: + print(f' WARNING: Variable {var} has {len(var_shape)} dimensions,' + f' not currently supported') + continue + + # Explicitly update the global dataset with modified data + # Using indexing assignment to ensure the data is actually written + global_ds[var].data = global_var + print(f' ✓ Copied {var}') + + return global_ds + + +def main(): + ''' + Main function to map regional mesh data to global mesh. + ''' + + args = parse_args() + + # Validate that output file differs from global file + if args.output_file == args.global_file: + print('ERROR: Output file must differ from global file') + print(' (to avoid overwriting the original global mesh)') + sys.exit(1) + + print('\n** Mapping regional mesh to global mesh **') + print(f'Regional mesh: {args.regional_file}') + print(f'Global mesh: {args.global_file}') + print(f'Output file: {args.output_file}') + print(f'Variables to copy: {", ".join(args.variables)}') + print(f'Coordinate type: {args.coord_type}') + print() + + # Load meshes + print('Loading regional mesh...') + regional_ds = xr.open_dataset(args.regional_file) + n_regional = regional_ds.dims['nCells'] + print(f' Regional mesh has {n_regional} cells') + + print('Loading global mesh...') + global_ds = xr.open_dataset(args.global_file) + n_global = global_ds.dims['nCells'] + print(f' Global mesh has {n_global} cells') + print() + + # Get coordinates based on type + if args.coord_type == 'spherical': + print('Using spherical coordinates (lonCell, latCell)...') + regional_coords = np.column_stack([ + regional_ds['lonCell'].values, + regional_ds['latCell'].values + ]) + global_coords = np.column_stack([ + global_ds['lonCell'].values, + global_ds['latCell'].values + ]) + else: + print('Using planar coordinates (xCell, yCell)...') + regional_coords = np.column_stack([ + regional_ds['xCell'].values, + regional_ds['yCell'].values + ]) + global_coords = np.column_stack([ + global_ds['xCell'].values, + global_ds['yCell'].values + ]) + + print() + + # Build coordinate mapping + mapping, unmatched = build_coordinate_mapping( + regional_coords, global_coords, args.tolerance + ) + + if not mapping: + print('ERROR: No matching cells found!') + sys.exit(1) + + print() + + # Verify mapping + mapping_ok = verify_mapping(mapping, regional_coords, global_coords, + args.tolerance) + if not mapping_ok: + print('WARNING: Mapping verification failed, but continuing...') + + print() + + # Report unmatched cells if any + if unmatched: + print(f'WARNING: {len(unmatched)} regional cells have no match in ' + f'global mesh') + print('These cells will be skipped.') + print() + + # If verify-only mode, stop here + if args.verify_only: + print('Verify-only mode: stopping before copying data') + return + + # Copy variables + print('Copying variables from regional to global mesh...') + global_ds = copy_variables(regional_ds, global_ds, mapping, args.variables) + print() + + # Update comment attribute + comment = ( + f'Variables {", ".join(args.variables)} updated from regional mesh ' + f'{args.regional_file}. {len(mapping)} cells modified.' + ) + if 'comment' in global_ds.attrs: + existing_comment = global_ds.attrs['comment'] + comment = f'{comment} {existing_comment}' + global_ds.attrs['comment'] = comment + + # Write output + print(f'Writing output to {args.output_file}...') + + if HAS_MPAS_TOOLS: + # Use mpas_tools for proper MPAS format + # write_netcdf will automatically add history with command line + write_netcdf(global_ds, args.output_file, format="NETCDF3_64BIT_DATA") + else: + # Fallback to xarray if mpas_tools not available + # Manually update history + timestamp = datetime.now().strftime('%Y-%m-%d %H:%M:%S') + command_line = ' '.join(sys.argv) + if 'history' in global_ds.attrs: + history = global_ds.attrs['history'] + history = f'{timestamp}: {command_line}\n{history}' + else: + history = f'{timestamp}: {command_line}' + global_ds.attrs['history'] = history + + encoding = {var: {'_FillValue': None} for var in global_ds.data_vars} + global_ds.to_netcdf(args.output_file, encoding=encoding) + + print('Done!') + print() + print(f'Summary:') + print(f' Regional cells: {n_regional}') + print(f' Global cells: {n_global}') + print(f' Matched cells: {len(mapping)}') + print(f' Unmatched regional cells: {len(unmatched)}') + print(f' Variables copied: {len(args.variables)}') + print(f' Output: {args.output_file}') + + # Clean up + regional_ds.close() + global_ds.close() + + +if __name__ == '__main__': + main() diff --git a/test_adjust_bed.sh b/test_adjust_bed.sh new file mode 100755 index 000000000..a255d9863 --- /dev/null +++ b/test_adjust_bed.sh @@ -0,0 +1,14 @@ +#!/bin/bash +# Test script for adjust_bed_to_haf.py + +# Test with the provided files using the standalone script +python3 landice/mesh_tools_li/adjust_bed_to_haf.py \ + --mesh /Users/trhille/Documents/ISMIP7/Antarctica/mesh/4to20km/AIS_4to20km_r03_20260910_ASE_extracted.nc \ + --geojson /Users/trhille/Documents/ISMIP7/Antarctica/wild_grounding_lines/Thwaites_GLs_2014_201920/Thwaites_GL_2014_pinning_points.geojson \ + --projection ais-bedmap2 \ + --output test_output.nc \ + --target-haf 10.0 + +echo "" +echo "Test completed. Output saved to test_output.nc" +echo "You can compare the original and modified bed topography using ncview or similar tools."