Skip to content

Interpolate ECCO ASTE reanalysis to MALI grid#691

Open
alexolinhager wants to merge 19 commits intoMPAS-Dev:masterfrom
alexolinhager:interpolate_ecco_to_mali
Open

Interpolate ECCO ASTE reanalysis to MALI grid#691
alexolinhager wants to merge 19 commits intoMPAS-Dev:masterfrom
alexolinhager:interpolate_ecco_to_mali

Conversation

@alexolinhager
Copy link
Contributor

@alexolinhager alexolinhager commented Nov 19, 2025

This PR introduces landice/mesh_tools_li/interpolate_ecco_to_mali.py, a script that interpolates ECCO ASTE R1 reanalysis data (https://arcticdata.io/data/10.18739/A2CV4BS5K/) onto a MALI mesh. Relevant ECCO files are currently stored at: /global/cfs/cdirs/fanssie/users/ahager/gisMeshSetup/forcingFiles/ASTE on Perlmutter, although these should eventually be moved to a permanent location. Currently supported fields are ocean temperature, ocean salinity, and sea ice fraction. ECCO grid files are unconventional for MITgcm simulations and do not contain adequate grid cell corner information for tile edge cells (necessary for creating scrip file and remapping). Therefore, this information needed to be hardcoded into interpolate_ecco_to_mali.py, using grid corner information from adjacent MITgcm tiles. Currently supported tiles are 14 and 27, which cover the western and northern coasts of Greenland.

Script can be tested with (on permutter):
python /global/cfs/cdirs/fanssie/users/ahager/MPAS-Tools/landice/mesh_tools_li/interpolate_ecco_to_mali.py --eccoDir /global/cfs/cdirs/fanssie/users/ahager/gisMeshSetup/forcingFiles/ASTE --maliMesh /global/cfs/cdirs/fanssie/MALI_input_files/UummannaqDisko_1to10km_r01/optimization/FO/UummannaqDisko_1to10km_r01_20251030_FO.nc --eccoVars SALT THETA SIarea --meshVars --geojson /global/cfs/cdirs/fanssie/users/ahager/gisMeshSetup/forcingFiles/UummannaqDisko/geojson/Ilulissat_Basin.geojson

Rearranges script so that remapping happens for each individual tile.
Allows ECCO arrays to remain gridded and improves efficiency of
ESMF_RegridWeightGen.
This reverts commit bf4de92db9d8f40500493c15cf33ed145166cb95.
Incorporates ECCO grid file to create custom scrip file for ECCO format
Hardcodes ECCO tile numbers and explicitly defines edge cells along
boundary for the two tiles. Currently only tiles 14 and 27 are
supported. Not enough information exists in ecco grid files to define
boundary cells edges without hardcoding (unconventional MITgcm grid file
format)..
Addresses formating issue with xtime variable in initial output file.
Still need to fix formatting in 'meshVars' output file.
Conservative remapping averages defined and undefined ocean cells. Use
orig3dOceanMask to identify these cells and treat as undefined
Creates a consistent naming scheme so that only one output file is
created and the remaining temporary files are removed
Extrapolates iceCellArea to fill ocean cells, and defines
icebergFjordMask wherever iceCellArea is above a threshold value
Adds an option to use a geojson file to define a region where
icebergFjordMask is permanently 1
@alexolinhager alexolinhager force-pushed the interpolate_ecco_to_mali branch from f13ce7c to daac57e Compare December 15, 2025 20:59
Moves saving of xtime time down to where final output file is created in
order to preserve its format
Creates a copy of the output file with netcdf3 64-bit offset
formatting to be consistent with MALI requirements
Monthly averaged ECCO data is output at the beginning of following
month. Need to subtract a month from the time stamp in order to force
MALI with the correct monthly ocean conditions.
Previously, a few invalid TS cells were designated as valid in the
orig3dOceanMask after remmapping, which caused issues when extrapolating
ocean properties. This commit add a check to remove this cells if they
occur
Changes default mapping method to bilinear. Better for remapping coarse
resolution ECCO data onto finer MALI mesh
Copy link
Contributor

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

This PR adds a new script to interpolate ECCO ASTE R1 reanalysis data (ocean temperature, salinity, sea ice fraction) onto a MALI mesh. It handles the unconventional ECCO grid format by hardcoding corner information for tiles 14 and 27 (covering western and northern Greenland coasts).

Changes:

  • New interpolate_ecco_to_mali.py script that merges MITgcm tile files into an unstructured grid, creates SCRIP files, and remaps ECCO data onto MALI meshes using ESMF regridding.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

# <<NOTE>>: Creating a new netcdf file like this changes the format of xtime. Need to add
# encoding update to maintain original formatting
ds_out.to_netcdf(self.options.outputFile, unlimited_dims=['Time'])
subprocess.run(["ncatted", "-a", "_FillValue,,d,,", altRemappedOutFile])
orig3dOceanMask = orig3dOceanMask[:,:,1:-1]
nz,nx,ny = np.shape(orig3dOceanMask)
orig3dOceanMask = orig3dOceanMask.reshape(nz, nx * ny)
od3mBool = 1

# Account for rounding error. Good cells are not exactly 1 and temp/sal invalid values may be averaged out.
# It is otherwise possible to get some invalid temp/sal values align with good orig3dOceanMask values after remapping
mask = (o3dm > 0.9999) & (temp < 100) & (salt < 100)
# Subtract one month from each time step. Ecco to output every month as an average of the previous month. We
# want to force MALI with ECCO data from the same month (avoid 1 month lag)
time = ds_ecco['tim'].values
self.xtime_str = [(pd.to_datetime(dt) - pd.to_datetime(dt)).strftime("%Y-%m-%d_%H:%M:%S").ljust(64) for dt in time]
subprocess.run(["ncpdq", "-O", "-a", "Time,nISMIP6OceanLayers,nCells", self.ecco_unstruct, self.ecco_unstruct])

print("ncremap ...")
st = time.time()
icebergFjordMask[seaice_mask] = 1

# Use geojson to add any additional cells to icebergFjordMask
if hasattr(self.options, "geojson"):
Comment on lines +306 to +329
ds_siaScrip = xr.open_dataset(self.ecco_scripfile)
ds_unstruct = xr.open_dataset(self.ecco_unstruct)
o3dm = ds_unstruct['orig3dOceanMask'][1,1,:].values
valid = (o3dm > 0.9999).astype('int32') # Account for rounding error
ds_siaScrip['grid_imask'] = xr.DataArray(valid.astype('int32'),
dims=('grid_size'),
attrs={'units':'unitless'})
sia_scripfile = 'ecco_sia.scrip.nc'
ds_siaScrip.to_netcdf(sia_scripfile)
ds_siaScrip.close()

print("Creating SIA mapping File")
sia_mapping_file = 'ecco_to_mali_mapping_bilinear_sia.nc'
args = ['srun',
'-n', self.options.ntasks, 'ESMF_RegridWeightGen',
'--source', sia_scripfile,
'--destination', mali_scripfile,
'--weight', sia_mapping_file,
'--method', 'bilinear',
'--extrap_method', 'creep',
'--extrap_num_levels', self.options.extrapIters,
'--netcdf4',
"--dst_regional", "--src_regional", '--ignore_unmapped']
check_call(args)

# Convert to netcdf3 64-bit offset. Much faster to do this with ncks than create the netcdf3 originally
# with xarray.to_netcdf.
nc64offset = self.options.outFile[0:-3] + '.64bitOffset.nc'
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants