Skip to content

Enable Option for High Res Mesh in Ocean #894

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
+465 −45
Open

Enable Option for High Res Mesh in Ocean #894

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
0f9fdbb
Enable coarse mesh resolution over the open ocean by using geojson file.
andrewdnolan Jun 20, 2024
af99278
Enable high resolution in ocean
alexolinhager Jan 16, 2025
9a5d0cd
Create option for high res in ocean for gis mesh
alexolinhager Feb 18, 2025
b3fc57e
Move geojson bnds setup to greenland mesh.py
alexolinhager Feb 19, 2025
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
83 changes: 46 additions & 37 deletions compass/landice/mesh.py
View file
Original file line number Diff line number Diff line change
@@ -265,12 +265,10 @@ def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None,
low_bed = section.getfloat('low_bed')
high_bed = section.getfloat('high_bed')

# convert km to m
cull_distance = section.getfloat('cull_distance') * 1.e3

# Cell spacing function based on union of masks
if section.get('use_bed') == 'True':
logger.info('Using bed elevation for spacing.')

if flood_fill_iStart is not None and flood_fill_jStart is not None:
logger.info('calling gridded_flood_fill to find \
bedTopography <= low_bed connected to the ocean.')
@@ -290,18 +288,30 @@ def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None,
k = 0.05 # This works well, but could try other values
spacing_bed = min_spac + (max_spac - min_spac) / (1.0 + np.exp(
-k * (bed - np.mean([high_bed, low_bed]))))
# We only want bed topography to influence spacing within high_dist_bed
# from the ice margin. In the region between high_dist_bed and
# low_dist_bed, use a linear ramp to damp influence of bed topo.
spacing_bed[dist_to_grounding_line >= low_dist_bed] = (
(1.0 - (dist_to_grounding_line[
dist_to_grounding_line >= low_dist_bed] -
low_dist_bed) / (high_dist_bed - low_dist_bed)) *
spacing_bed[dist_to_grounding_line >= low_dist_bed] +
(dist_to_grounding_line[dist_to_grounding_line >=
low_dist_bed] - low_dist_bed) /
(high_dist_bed - low_dist_bed) * max_spac)
spacing_bed[dist_to_grounding_line >= high_dist_bed] = max_spac

# If max_res_in_ocn is true, use the minimum cell spacing for
# all ocean cells. Important for ocean-coupled runs where
# resolving fjords and ocean bathymetry is necessary for
# thermal forcing parameterizations. Otherwize, telescope out
# to coarse resolution with distance from the ice front.
if section.get('max_res_in_ocn') == 'True':
spacing_bed[np.logical_and(bed < 0, thk == 0)] = min_spac
print("Maximizing resolution in ocean. {} ocean cells \
".format(np.sum(np.logical_and(bed < 0, thk == 0))))
else:
# We only want bed topography to influence spacing within
# high_dist_bed from the ice margin. In the region
# between high_dist_bed and low_dist_bed, use a linear
# ramp to damp influence of bed topo.
spacing_bed[dist_to_grounding_line >= low_dist_bed] = (
(1.0 - (dist_to_grounding_line[
dist_to_grounding_line >= low_dist_bed] -
low_dist_bed) / (high_dist_bed - low_dist_bed)) *
spacing_bed[dist_to_grounding_line >= low_dist_bed] +
(dist_to_grounding_line[dist_to_grounding_line >=
low_dist_bed] - low_dist_bed) /
(high_dist_bed - low_dist_bed) * max_spac)
spacing_bed[dist_to_grounding_line >= high_dist_bed] = max_spac
if flood_fill_iStart is not None and flood_fill_jStart is not None:
spacing_bed[low_bed_mask == 0] = max_spac
# Do one more flood fill to eliminate isolated pockets
@@ -366,17 +376,6 @@ def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None,
for width in [spacing_bed, spacing_speed, spacing_edge, spacing_gl]:
cell_width = np.minimum(cell_width, width)

# Set large cell_width in areas we are going to cull anyway (speeds up
# whole process). Use 3x the cull_distance to avoid this affecting
# cell size in the final mesh. There may be a more rigorous way to set
# that distance.
if dist_to_edge is not None:
mask = np.logical_and(
thk == 0.0, dist_to_edge > (3. * cull_distance))
logger.info('Setting cell_width in outer regions to max_spac '
f'for {mask.sum()} cells')
cell_width[mask] = max_spac

return cell_width


@@ -518,7 +517,7 @@ def get_dist_to_edge_and_gl(self, thk, topg, x, y,


def build_cell_width(self, section_name, gridded_dataset,
flood_fill_start=[None, None]):
flood_fill_start=[None, None], calc_geom_bnds=True):
"""
Determine MPAS mesh cell size based on user-defined density function.
@@ -543,6 +542,11 @@ def build_cell_width(self, section_name, gridded_dataset,
fill. Most cases will use ``[None, None]``, which will just start the
flood fill in the center of the gridded dataset.
calc_geom_bnds : logical
Option to calculate geom_points and geom_edges needed for jigsaw within
build_cell_width. Default is to perform calculation, but the user may
opt out if these are determined elsewhere (e.g., using a geojson file)
Returns
-------
cell_width : numpy.ndarray
@@ -581,15 +585,16 @@ def build_cell_width(self, section_name, gridded_dataset,

f.close()

# Get bounds defined by user, or use bound of gridded dataset
bnds = [np.min(x1), np.max(x1), np.min(y1), np.max(y1)]
bnds_options = ['x_min', 'x_max', 'y_min', 'y_max']
for index, option in enumerate(bnds_options):
bnd = section.get(option)
if bnd != 'None':
bnds[index] = float(bnd)
# If necessary, get bounds defined by user or use bound of gridded dataset
if calc_geom_bnds:
bnds = [np.min(x1), np.max(x1), np.min(y1), np.max(y1)]
bnds_options = ['x_min', 'x_max', 'y_min', 'y_max']
for index, option in enumerate(bnds_options):
bnd = section.get(option)
if bnd != 'None':
bnds[index] = float(bnd)

geom_points, geom_edges = set_rectangular_geom_points_and_edges(*bnds)
geom_points, geom_edges = set_rectangular_geom_points_and_edges(*bnds)

# Remove ice not connected to the ice sheet.
flood_mask = gridded_flood_fill(thk)
@@ -611,8 +616,12 @@ def build_cell_width(self, section_name, gridded_dataset,
flood_fill_iStart=flood_fill_start[0],
flood_fill_jStart=flood_fill_start[1])

return (cell_width.astype('float64'), x1.astype('float64'),
y1.astype('float64'), geom_points, geom_edges, flood_mask)
if not calc_geom_bnds:
return (cell_width.astype('float64'), x1.astype('float64'),
y1.astype('float64'), flood_mask)
else:
return (cell_width.astype('float64'), x1.astype('float64'),
y1.astype('float64'), geom_points, geom_edges, flood_mask)


def build_mali_mesh(self, cell_width, x1, y1, geom_points,
342 changes: 342 additions & 0 deletions compass/landice/tests/greenland/gis_contShelfExtent_PS.geojson
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,342 @@
{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"properties": {
"name": "westCentralGreenland1",
"tags": "westCentralGreenland",
"object": "region",
"component": "landice",
"author": "http://icesat4.gsfc.nasa.gov/cryo_data/ant_grn_drainage_systems.php"
},
"geometry": {
"type": "Polygon",
"coordinates": [
[
[
-588895.3500545282,
-2114024.5158124366
],
[
-600330.202835335,
-2052276.3107960785
],
[
-558601.8188636905,
-1995223.4835278299
],
[
-538010.2551799342,
-1931066.8713195215
],
[
-511710.09378407826,
-1884183.9749182125
],
[
-542584.1962922566,
-1860170.784078517
],
[
-555162.5343511449,
-1821292.2846237745
],
[
-535723.2846237727,
-1786987.726281353
],
[
-579175.7251908389,
-1736674.374045801
],
[
-589467.0926935658,
-1673782.6837513635
],
[
-582606.1810250811,
-1623469.331515813
],
[
-556306.0196292251,
-1605173.567066521
],
[
-587180.1221374033,
-1573155.9792802609
],
[
-597471.4896401304,
-1539994.9062159217
],
[
-648928.327153762,
-1498829.436205017
],
[
-651215.2977099229,
-1431363.8047982552
],
[
-651215.2977100034,
-1059731.0894221598
],
[
-635177.8481473451,
-1040280.8406720696
],
[
-522001.46128683415,
-960247.8702290629
],
[
-350478.6695747058,
-839038.4307524698
],
[
-199538.6128680486,
-736124.7557251963
],
[
-41737.64449291128,
-642358.962922575
],
[
45167.236641222284,
-640071.9923664134
],
[
76041.33914940123,
-634354.5659760092
],
[
224694.42529989238,
-632067.595419848
],
[
252138.07197382947,
-640071.9923664138
],
[
299020.9683751388,
-670946.0948745939
],
[
347047.35005452833,
-654937.3009814639
],
[
375634.48200654646,
-661798.212649949
],
[
408795.5550708885,
-700676.7121046956
],
[
455678.45147219795,
-697246.2562704541
],
[
494556.95092694415,
-731550.8146128787
],
[
640923.0665212895,
-788725.0785169315
],
[
718680.0654308012,
-861908.1363141248
],
[
734688.8593239472,
-914508.4591058589
],
[
778141.2998910325,
-930517.252999008
],
[
830741.6226826613,
-980830.6052344592
],
[
876481.0338058897,
-1051726.6924754628
],
[
878768.0043620515,
-1250693.130861505
],
[
862759.2104689217,
-1509120.8037077426
],
[
878768.0043620509,
-1628043.2726281346
],
[
878768.0043620516,
-2032265.3184296673
],
[
879339.7470010923,
-2104876.633587793
],
[
858757.0119956397,
-2120885.427480922
],
[
863330.9531079635,
-2166624.8386041536
],
[
831313.3653217041,
-2210077.279171221
],
[
808443.6597600909,
-2290121.248636872
],
[
831313.3653217062,
-2415904.629225754
],
[
771852.1308615117,
-2537114.0687023154
],
[
694095.1319520242,
-2612584.097055645
],
[
604903.2802617289,
-2621731.979280291
],
[
490554.75245365844,
-2724645.6543075596
],
[
481406.87022901414,
-2788680.8298800867
],
[
433380.48854962486,
-2845855.0937841283
],
[
389928.04798255826,
-2896168.4460196844
],
[
309884.0785169081,
-2909890.2693566554
],
[
264144.6673936795,
-2944194.82769908
],
[
234414.05016358144,
-3005943.0327154496
],
[
245848.90294438964,
-3058543.355507171
],
[
197822.5212650006,
-3184326.736096079
],
[
158944.02181025772,
-3353562.5572520755
],
[
120065.5223555122,
-3385580.1450383463
],
[
-53744.239912762794,
-3381006.203926024
],
[
-94909.70992366974,
-3330692.851690453
],
[
-170379.73827699898,
-3296388.293348022
],
[
-243562.79607416524,
-3230066.1472193203
],
[
-252710.67829881026,
-3186613.706652242
],
[
-309884.9422028462,
-3127152.4721920304
],
[
-309884.94220284483,
-3063117.2966194954
],
[
-346476.47110142827,
-3040247.591057881
],
[
-410511.6466739453,
-2857289.9465649365
],
[
-447103.17557252746,
-2809263.5648855423
],
[
-490555.6161395937,
-2754376.2715376625
],
[
-465398.94002181716,
-2708636.8604144296
],
[
-474546.82224646193,
-2646888.655398066
],
[
-517999.26281352824,
-2587427.4209378655
],
[
-506564.4100327197,
-2447922.217012012
],
[
-598356.8312565115,
-2370130.1091510067
],
[
-604904.1439476594,
-2290121.2486368697
],
[
-529434.115594332,
-2173485.750272634
]
]
]
}
}
]
}
8 changes: 8 additions & 0 deletions ...ce/tests/greenland/greenland_only_outline_45km_buffer_latlon_singlepart.geojson
View file

Large diffs are not rendered by default.

77 changes: 69 additions & 8 deletions compass/landice/tests/greenland/mesh.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
import json
import os

import jigsawpy
import numpy as np
import xarray as xr
from mpas_tools.scrip.from_mpas import scrip_from_mpas
@@ -31,7 +33,7 @@ def __init__(self, test_case):
Parameters
----------
test_case : compass.TestCase
The test case this step belongs to
The test case this step belongs to/
"""
super().__init__(test_case=test_case, name='mesh', cpus_per_task=128,
@@ -53,6 +55,10 @@ def __init__(self, test_case):
self.add_input_file(filename='greenland_2km_2024_01_29.epsg3413.nc',
target='greenland_2km_2024_01_29.epsg3413.nc',
database='')
self.add_input_file(filename='greenland_only_outline_45km_buffer_latlon_singlepart.geojson', # noqa: E501
package='compass.landice.tests.greenland',
target='greenland_only_outline_45km_buffer_latlon_singlepart.geojson', # noqa: E501
database=None)

# no setup() method is needed

@@ -80,18 +86,29 @@ def run(self):
source_gridded_dataset_2km = 'greenland_2km_2024_01_29.epsg3413.nc'

logger.info('calling build_cell_width')
cell_width, x1, y1, geom_points, geom_edges, floodMask = \
build_cell_width(
self, section_name=section_name,
gridded_dataset=source_gridded_dataset_2km,
flood_fill_start=[100, 700])
if section_gis.get("define_bnds_by_geojson") == 'True':
geom_points, geom_edges = set_geojson_geom_points_and_edges(self)

cell_width, x1, y1, floodMask = \
build_cell_width(self,
section_name=section_name,
gridded_dataset=source_gridded_dataset_2km,
flood_fill_start=[100, 700],
calc_geom_bnds=False)
else:
cell_width, x1, y1, geom_points, geom_edges, floodMask = \
build_cell_width(
self, section_name=section_name,
gridded_dataset=source_gridded_dataset_2km,
flood_fill_start=[100, 700])

# Now build the base mesh and perform the standard interpolation
build_mali_mesh(
self, cell_width, x1, y1, geom_points, geom_edges,
mesh_name=self.mesh_filename, section_name=section_name,
gridded_dataset=source_gridded_dataset_1km,
projection=src_proj, geojson_file=None)
gridded_dataset=source_gridded_dataset_1km, projection=src_proj,
geojson_file="greenland_only_outline_45km_buffer_latlon_singlepart.geojson", # noqa: E501
)

# Create scrip file for the newly generated mesh
logger.info('creating scrip file for destination mesh')
@@ -158,3 +175,47 @@ def run(self):
ds["observedThicknessTendencyUncertainty"] = dHdtErr
# Write the data to disk
ds.to_netcdf(self.mesh_filename, 'a')


def set_geojson_geom_points_and_edges(self):
self.add_input_file(filename='gis_contShelfExtent.geojson',
package='compass.landice.tests.greenland',
target='gis_contShelfExtent.geojson',
database=None)

with open('gis_contShelfExtent.geojson', 'r') as meshMarginFile:
geojson_data = json.load(meshMarginFile)

x = []
y = []
start_edge = []
end_edge = []
ct = 0

for feature in geojson_data['features']:
geometry = feature['geometry']
for coord in geometry['coordinates']:
for sub_coord in coord:
x.append(sub_coord[0])
y.append(sub_coord[1])

start_edge.append(ct)
end_edge.append(ct + 1)
ct = ct + 1
x = np.array(x)
y = np.array(y)
start_edge = np.array(start_edge)
end_edge = np.array(end_edge)

start_edge[-1] = ct - 1
end_edge[-1] = 0

geom_points = np.array([((x[i], y[i]), 0)
for i in range(len(y))],
dtype=jigsawpy.jigsaw_msh_t.VERT2_t)

geom_edges = np.array([((start_edge[i], end_edge[i]), 0)
for i in range(len(start_edge))],
dtype=jigsawpy.jigsaw_msh_t.EDGE2_t)

return (geom_points, geom_edges)