From b57823fac38697f04725f3bf5fc0b94d980a153c Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Wed, 8 Oct 2025 14:47:38 -0700 Subject: [PATCH 01/17] Version 0 of ECCO to MALI interpolation --- .../mesh_tools_li/interpolate_ecco_to_mali.py | 401 ++++++++++++++++++ 1 file changed, 401 insertions(+) create mode 100755 landice/mesh_tools_li/interpolate_ecco_to_mali.py diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py new file mode 100755 index 000000000..3576beefc --- /dev/null +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -0,0 +1,401 @@ +#!/usr/bin/env python + +import os +import shutil +import xarray as xr +import subprocess +import numpy as np +import pandas as pd +from scipy import stats as st +from mpas_tools.logging import check_call +from mpas_tools.scrip.from_mpas import scrip_from_mpas +from mpas_tools.ocean.depth import compute_zmid +from pyproj import Proj, Transformer +from argparse import ArgumentParser +import time +import cftime +import glob + +class eccoToMaliInterp: + + def __init__(self): + print("Gathering Information ... ") + parser = ArgumentParser( + prog='interpolate_ecco_to_mali.py', + description='Interpolates ECCO reanalysis data onto a MALI mesh') + parser.add_argument("--eccoDir", dest="eccoDir", required=True, help="Directory where ECCO files are stored. Should be only netcdf files in that directory", metavar="FILENAME") + parser.add_argument("--maliMesh", dest="maliMesh", required=True, help="MALI mesh to be interpolated onto", metavar="FILENAME") + parser.add_argument("--tileNums", type=int, nargs='+', dest="tileNums", required=True, help="ECCO grid regional tile numbers to pull data from") + parser.add_argument("--eccoVars", type=str, nargs='+', dest="eccoVars", required=True, help="ECCO variables to interpolate, current options are any combination of 'THETA', 'SALT', 'SIarea'") + parser.add_argument("--maxDepth", type=int, dest="maxDepth", default=1500, help="Maximum depth in meters to include in the interpolation") + parser.add_argument("--ntasks", dest="ntasks", type=str, help="Number of processors to use with ESMF_RegridWeightGen", default='128') + parser.add_argument("--method", dest="method", type=str, help="Remapping method, either 'bilinear' or 'neareststod'", default='conserve') + parser.add_argument("--outFile", dest="outputFile", help="Desired name of output file", metavar="FILENAME", default="ecco_to_mali_remapped.nc") + parser.add_argument("--meshVars", dest="includeMeshVars", help="whether to include mesh variables in resulting MALI forcing file", action='store_true') + args = parser.parse_args() + self.options = args + + self.mapping_file_name = f'ecco_to_mali_mapping_{self.options.method}.nc' + + def merge_MITgcm_files_to_unstruct_grid(self): + ncfiles = glob.glob(os.path.join(self.options.eccoDir, "*.nc")) + tIter = 0 + saltBool = 0 + thetaBool = 0 + siaBool = 0 + + #loop through MITgcm tile files and create unstructured arrays for each variable in eccoVars + for tile in self.options.tileNums: + print(f"Processing Tile {tile}") + o3dmBool = 0 + for file in ncfiles: + print(f"Processing File {file}") + paddedTile = str(tile).zfill(4) + print(f"Padded Tile: {paddedTile}") + if paddedTile in file: + print("Padded tile is in file") + if 'SALT' in file and 'SALT' in self.options.eccoVars: + validFile = self.options.eccoDir + '/SALT.' + paddedTile + '.nc' + ds_ecco = xr.open_dataset(validFile) + + # only open depth variable once, regardless of which variables we are interested in + if saltBool == 0 and thetaBool == 0: + z = ds_ecco['dep'].values + + if tIter == 0: + sal = ds_ecco['SALT'][:].values + nt,nz,nx,ny = np.shape(sal) + sal = sal.reshape(nz, nt, nx * ny) + else: + sl = ds_ecco['SALT'][:].values + nt,nz,nx,ny = np.shape(sl) + sl = sl.reshape(nz, nt, nx * ny) + sal = np.concatenate((sal,sl), axis=2) + + # ECCO 'land' variable is equivalent to MALI's orig3dOceanMask, despite the counterintuitive name. + # Only need to call this once per tile, but needs to be using a theta/salinity file because it includes depth + # No need to call this if only interpolating SIarea + if tIter == 0 and o3dmBool == 0: + orig3dOceanMask = ds_ecco['land'].values + nz,nx,ny = np.shape(orig3dOceanMask) + orig3dOceanMask = orig3dOceanMask.reshape(nz, nx * ny) + od3mBool = 1 + elif tIter > 0 and o3dmBool == 0: + o3dm = ds_ecco['land'].values + nz,nx,ny = np.shape(o3dm) + o3dm = o3dm.reshape(nz, nx * ny) + orig3dOceanMask = np.concatenate((orig3dOceanMask, o3dm), axis=1) + o3dmBool = 1 + + # set MALI invalid value + boolZ,boolXY = np.where(orig3dOceanMask == 0) + for iZ,iXY in zip(boolZ,boolXY): + sal[iZ,:,iXY] = 1e36 + + saltBool = 1 + + # Repeat for potential temperature if needed + elif 'THETA' in file and 'THETA' in self.options.eccoVars: + validFile = self.options.eccoDir + '/THETA.' + paddedTile + '.nc' + ds_ecco = xr.open_dataset(validFile) + + if saltBool == 0 and thetaBool == 0: + z = ds_ecco['dep'].values + + if tIter == 0: + theta = ds_ecco['THETA'][:].values + nt,nz,nx,ny = np.shape(theta) + theta = theta.reshape(nz, nt, nx * ny) + else: + th = ds_ecco['THETA'][:].values + nt,nz,nx,ny = np.shape(th) + th = th.reshape(nz, nt, nx * ny) + theta = np.concatenate((theta,th), axis=2) + + # ECCO 'land' variable is equivalent to MALI's orig3dOceanMask, despite the counterintuitive name. + # Only need to call this once per tile, but needs to be using a theta/salinity file because it includes depth + # No need to call this if only interpolating SIarea + if tIter == 0 and o3dmBool == 0: + orig3dOceanMask = ds_ecco['land'].values + nz,nx,ny = np.shape(orig3dOceanMask) + orig3dOceanMask = orig3dOceanMask.reshape(nz, nx * ny) + od3mBool = 1 + + elif tIter > 0 and o3dmBool == 0: + o3dm = ds_ecco['land'].values + nz,nx,ny = np.shape(o3dm) + o3dm = o3dm.reshape(nz, nx * ny) + orig3dOceanMask = np.concatenate((orig3dOceanMask, o3dm), axis=1) + o3dmBool = 1 + + # set MALI invalid value + boolZ,boolXY = np.where(orig3dOceanMask == 0) + for iZ,iXY in zip(boolZ,boolXY): + theta[iZ,:,iXY] = 1e36 + + thetaBool = 1 + + # Repeat for sea ice fraction if needed + elif 'SIarea' in file and 'SIarea' in self.options.eccoVars: + siaBool = 1 + validFile = self.options.eccoDir + '/SIarea.' + paddedTile + '.nc' + ds_ecco = xr.open_dataset(validFile) + if tIter == 0: + sia = ds_ecco['SIarea'][:].values + nt,nx,ny = np.shape(sia) + sia = sia.reshape(nt, nx * ny) + else: + si = ds_ecco['SIarea'][:].values + nt,nx,ny = np.shape(si) + si = si.reshape(nt, nx * ny) + sia = np.concatenate((sia,si), axis=1) + else : + raise ValueError(f"eccoDir: {self.options.eccoDir}, paddedTile: {paddedTile}") + + if ('SALT' in self.options.eccoVars and saltBool == 0): + raise ValueError(f"SALT netcdf file for tile {tile} not found in {self.options.eccoDir}") + + if ('THETA' in self.options.eccoVars and thetaBool == 0): + raise ValueError(f"THETA netcdf file for tile {tile} not found in {self.options.eccoDir}") + + if ('SIarea' in self.options.eccoVars and siaBool == 0): + raise ValueError(f"SIarea netcdf file for tile {tile} not found in {self.options.eccoDir}") + + + # Doesn't matter which file these come from, just that there is only one set per tile, so just use most + # recent valid file in each tile + + if tIter == 0: + if thetaBool == 1 : + validFile = self.options.eccoDir + '/THETA.' + paddedTile + '.nc' + elif saltBool == 1 : + validFile = self.options.eccoDir + '/SALT.' + paddedTile + '.nc' + elif siaBool == 1: + validFile = self.options.eccoDir + '/SIarea.' + paddedTile + '.nc' + else : + raise ValueError("No variable bool found") + + ds_ecco = xr.open_dataset(validFile) + lon = ds_ecco['lon'][:].values.flatten() + lat = ds_ecco['lat'][:].values.flatten() + else: + ds_ecco = xr.open_dataset(validFile) + ln = ds_ecco['lon'][:].values.flatten() + lon = np.concatenate((lon, ln)) + lt = ds_ecco['lat'][:].values.flatten() + lat = np.concatenate((lat,lt)) + + tIter = tIter + 1 + + # combine into new netcdf. Need to load depth and time from original ecco files. Assuming time/depth is consistent + # across all files, just load info from most recent + + # Establish dataset, and fill in standard dimensions. other dimensions are added on an as-needed basis + ds_unstruct = xr.Dataset() + ds_unstruct.expand_dims(["nCells", "Time"]) + + # Translate ECCO time to MALI xtime format as save to dataset + time = ds_ecco['tim'].values + xtime = [pd.to_datetime(dt).strftime("%Y-%m-%d_%H:%M:%S").ljust(64) for dt in time] + + DA_xtime = xr.DataArray(np.array(xtime).astype('S64'), dims=("Time")) + DA_xtime.encoding.update({"char_dim_name":"StrLen"}) + DA_xtime.attrs['long_name'] = "model time, with format 'YYYY-MM-DD_HH:MM:SS'" + DA_xtime.attrs['standard_name'] = 'time' + ds_unstruct['xtime'] = DA_xtime + + # save unstructured variables with MALI/MPAS-O names + if 'z' in locals(): + ds_unstruct.expand_dims("nISMIP6OceanLayers") + ds_unstruct.expand_dims("TWO") + z = z[z <= self.options.maxDepth] + + # create bnds array + zBnds = np.zeros((len(z), 2)) + zBnds[0, 0] = 0.0 + for i in range(1, len(z)): + zBnds[i, 0] = 0.5 * (z[i - 1] + z[i]) + for i in range(0, len(z) - 1): + zBnds[i, 1] = 0.5 * (z[i] + z[i + 1]) + zBnds[-1, 1] = z[-1] + (z[-1] - zBnds[-1, 0]) + + DA_z = xr.DataArray(-z.astype('float64'), dims=("nISMIP6OceanLayers")) #MALI wants negative depth coordinates + DA_z.attrs['long_name'] = 'depth' + DA_z.attrs['units'] = 'm' + ds_unstruct['ismip6shelfMelt_zOcean'] = DA_z + + DA_zBnds = xr.DataArray(-zBnds.astype('float64'), dims=("nISMIP6OceanLayers", "TWO")) #MALI wants negative depth coordinates + DA_zBnds.attrs['long_name'] = 'Bounds for ismip6 ocean layers' + DA_zBnds.attrs['units'] = 'm' + ds_unstruct['ismip6shelfMelt_zBndsOcean'] = DA_zBnds + + if 'theta' in locals(): + print(f"z shape: {z.shape}") + print(f"theta shape: {theta.shape}") + theta = theta[0:len(z),:,:] + print(f"new theta shape: {theta.shape}") + DA_theta = xr.DataArray(theta.astype('float64'), dims=("nISMIP6OceanLayers", "Time", "nCells")) + DA_theta.attrs['long_name'] = 'Potential_Temperature' + DA_theta.attrs['units'] = 'degC' + ds_unstruct['oceanTemperature'] = DA_theta + + if 'sal' in locals(): + sal = sal[0:len(z),:,:] + DA_salt = xr.DataArray(sal.astype('float64'), dims=("nISMIP6OceanLayers", "Time", "nCells")) + DA_salt.attrs['long_name'] = 'Salinity' + DA_salt.attrs['units'] = 'psu' + ds_unstruct['oceanSalinity'] = DA_salt + + if 'sia' in locals(): + DA_sia = xr.DataArray(sia.astype('float64'), dims=("Time", "nCells")) + DA_sia.attrs['long_name'] = 'Sea ice fractional ice-covered area [0 to 1]' + DA_sia.attrs['units'] = 'm2/m2' + ds_unstruct['iceAreaCell'] = DA_sia + + if 'orig3dOceanMask' in locals(): + orig3dOceanMask = orig3dOceanMask[0:len(z),:] + DA_o3dm = xr.DataArray(orig3dOceanMask.astype('int32'), dims=("nISMIP6OceanLayers", "nCells")) + DA_o3dm.attrs['long_name'] = ("3D mask of original valid ocean data. Because it is 3d," + " it can include the ocean domain inside ice-shelf cavities." + " orig3dOceanMask is altered online to include ice-dammed inland seas") + ds_unstruct['orig3dOceanMask'] = DA_o3dm + + DA_lon = xr.DataArray(lon.astype('float64'), dims=("nCells")) + DA_lon.attrs['long_name'] = 'longitude' + DA_lon.attrs['standard_name'] = 'longitude' + DA_lon.attrs['units'] = 'degree_east' + ds_unstruct['lon'] = DA_lon # Keep original lon/lat names of coordinates for projection transformer + + DA_lat = xr.DataArray(lat.astype('float64'), dims=("nCells")) + DA_lat.attrs['long_name'] = 'latitude' + DA_lat.attrs['standard_name'] = 'latitude' + DA_lat.attrs['units'] = 'degree_north' + ds_unstruct['lat'] = DA_lat + + self.ecco_unstruct = "ecco_combined_unstructured.nc" + ds_unstruct.to_netcdf(self.ecco_unstruct) + ds_unstruct.close() + + def remap_ecco_to_mali(self): + + print("Creating Scrip Files") + ecco_scripfile = "tmp_ecco_scrip.nc" + mali_scripfile = "tmp_mali_scrip.nc" + ecco_unstruct = self.ecco_unstruct + file, ext = os.path.splitext(ecco_unstruct) + ecco_unstruct_xy = file + '.tmp' + ext + + #convert ecco mesh to polar stereographic + EM = xr.open_dataset(ecco_unstruct) + lon = EM['lon'].values + lat = EM['lat'].values + #use gis-gimp projection, defined in mpas_tools.landice.projections + proj_string = ( + "+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" + ) + transformer = Transformer.from_crs("EPSG:4326", proj_string, always_xy=True) + x,y = transformer.transform(lon,lat) + + EM['x'] = xr.DataArray(x,dims=("nCells")) + EM['y'] = xr.DataArray(y,dims=("nCells")) + EM.to_netcdf(ecco_unstruct_xy) + EM.close() + + args = ['create_scrip_file_from_planar_rectangular_grid', + '-i', ecco_unstruct_xy, + '-s', ecco_scripfile, + '-p', 'gis-gimp', + '-r', '1'] + check_call(args) + + scrip_from_mpas(self.options.maliMesh, mali_scripfile) + + #create mapping file + print("Creating Mapping File") + args = ['srun', + '-n', self.options.ntasks, 'ESMF_RegridWeightGen', + '--source', ecco_scripfile, + '--destination', mali_scripfile, + '--weight', self.mapping_file_name, + '--method', self.options.method, + '--netcdf4', + "--dst_regional", "--src_regional", '--ignore_unmapped'] + check_call(args) + + #remap + print("Start remapping ... ") + + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", ecco_unstruct]) + subprocess.run(["ncpdq", "-O", "-a", "Time,nISMIP6OceanLayers,nCells", ecco_unstruct, ecco_unstruct]) + + print("ncremap ...") + st = time.time() + # remap the input data + args_remap = ["ncremap", + "-i", ecco_unstruct_xy, + "-o", self.options.outputFile, + "-m", self.mapping_file_name] + check_call(args_remap) + + subprocess.run(["ncpdq", "-O", "-a", "Time,nCells,nISMIP6OceanLayers", self.options.outputFile, self.options.outputFile]) + nd = time.time() + tm = nd - st + print(f"ncremap completed: {tm} seconds") + print("Remapping completed") + + # Transfer MALI mesh variables if needed + if self.options.includeMeshVars: + ds_mali = xr.open_dataset(self.options.maliMesh, decode_times=False, decode_cf=False) + ds_out = xr.open_dataset(self.options.outputFile, decode_times=False, decode_cf=False) + ds_out['angleEdge'] = ds_mali['angleEdge'] + ds_out['areaCell'] = ds_mali['areaCell'] + ds_out['areaTriangle'] = ds_mali['areaTriangle'] + ds_out['cellsOnCell'] = ds_mali['cellsOnCell'] + ds_out['cellsOnEdge'] = ds_mali['cellsOnEdge'] + ds_out['cellsOnVertex'] = ds_mali['cellsOnVertex'] + ds_out['dcEdge'] = ds_mali['dcEdge'] + ds_out['dvEdge'] = ds_mali['dvEdge'] + ds_out['edgesOnCell'] = ds_mali['edgesOnCell'] + ds_out['edgesOnEdge'] = ds_mali['edgesOnEdge'] + ds_out['edgesOnVertex'] = ds_mali['edgesOnVertex'] + ds_out['indexToCellID'] = ds_mali['indexToCellID'] + ds_out['indexToEdgeID'] = ds_mali['indexToEdgeID'] + ds_out['indexToVertexID'] = ds_mali['indexToVertexID'] + ds_out['kiteAreasOnVertex'] = ds_mali['kiteAreasOnVertex'] + ds_out['latCell'] = ds_mali['latCell'] + ds_out['latEdge'] = ds_mali['latEdge'] + ds_out['latVertex'] = ds_mali['latVertex'] + ds_out['lonCell'] = ds_mali['lonCell'] + ds_out['lonEdge'] = ds_mali['lonEdge'] + ds_out['lonVertex'] = ds_mali['lonVertex'] + ds_out['meshDensity'] = ds_mali['meshDensity'] + ds_out['nEdgesOnCell'] = ds_mali['nEdgesOnCell'] + ds_out['nEdgesOnEdge'] = ds_mali['nEdgesOnEdge'] + ds_out['verticesOnCell'] = ds_mali['verticesOnCell'] + ds_out['verticesOnEdge'] = ds_mali['verticesOnEdge'] + ds_out['weightsOnEdge'] = ds_mali['weightsOnEdge'] + ds_out['xCell'] = ds_mali['xCell'] + ds_out['xEdge'] = ds_mali['xEdge'] + ds_out['xVertex'] = ds_mali['xVertex'] + ds_out['yCell'] = ds_mali['yCell'] + ds_out['yEdge'] = ds_mali['yEdge'] + ds_out['yVertex'] = ds_mali['yVertex'] + ds_out['zCell'] = ds_mali['zCell'] + ds_out['zEdge'] = ds_mali['zEdge'] + ds_out['zVertex'] = ds_mali['zVertex'] + ds_out.to_netcdf(self.options.outputFile, mode='w') + ds_out.close() + ds_mali.close() + +def main(): + run = eccoToMaliInterp() + + run.merge_MITgcm_files_to_unstruct_grid() + + run.remap_ecco_to_mali() + +if __name__ == "__main__": + main() + + From 6c9d90d3315bc73343040d668b8f6b82e2a9b16b Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Thu, 9 Oct 2025 14:38:17 -0700 Subject: [PATCH 02/17] Reorder to remap each tile Rearranges script so that remapping happens for each individual tile. Allows ECCO arrays to remain gridded and improves efficiency of ESMF_RegridWeightGen. --- .../mesh_tools_li/interpolate_ecco_to_mali.py | 378 ++++++++---------- 1 file changed, 157 insertions(+), 221 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index 3576beefc..3ec4a6a31 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -39,7 +39,6 @@ def __init__(self): def merge_MITgcm_files_to_unstruct_grid(self): ncfiles = glob.glob(os.path.join(self.options.eccoDir, "*.nc")) - tIter = 0 saltBool = 0 thetaBool = 0 siaBool = 0 @@ -62,38 +61,22 @@ def merge_MITgcm_files_to_unstruct_grid(self): if saltBool == 0 and thetaBool == 0: z = ds_ecco['dep'].values - if tIter == 0: - sal = ds_ecco['SALT'][:].values - nt,nz,nx,ny = np.shape(sal) - sal = sal.reshape(nz, nt, nx * ny) - else: - sl = ds_ecco['SALT'][:].values - nt,nz,nx,ny = np.shape(sl) - sl = sl.reshape(nz, nt, nx * ny) - sal = np.concatenate((sal,sl), axis=2) + sal = ds_ecco['SALT'][:].values # ECCO 'land' variable is equivalent to MALI's orig3dOceanMask, despite the counterintuitive name. # Only need to call this once per tile, but needs to be using a theta/salinity file because it includes depth # No need to call this if only interpolating SIarea - if tIter == 0 and o3dmBool == 0: + if o3dmBool == 0: orig3dOceanMask = ds_ecco['land'].values - nz,nx,ny = np.shape(orig3dOceanMask) - orig3dOceanMask = orig3dOceanMask.reshape(nz, nx * ny) od3mBool = 1 - elif tIter > 0 and o3dmBool == 0: - o3dm = ds_ecco['land'].values - nz,nx,ny = np.shape(o3dm) - o3dm = o3dm.reshape(nz, nx * ny) - orig3dOceanMask = np.concatenate((orig3dOceanMask, o3dm), axis=1) - o3dmBool = 1 # set MALI invalid value - boolZ,boolXY = np.where(orig3dOceanMask == 0) - for iZ,iXY in zip(boolZ,boolXY): - sal[iZ,:,iXY] = 1e36 + boolZ,boolX,boolY = np.where(orig3dOceanMask == 0) + for iZ,iX,iY in zip(boolZ,boolX,boolY): + sal[:,iZ,iX,iY] = 1e36 saltBool = 1 - + # Repeat for potential temperature if needed elif 'THETA' in file and 'THETA' in self.options.eccoVars: validFile = self.options.eccoDir + '/THETA.' + paddedTile + '.nc' @@ -102,53 +85,30 @@ def merge_MITgcm_files_to_unstruct_grid(self): if saltBool == 0 and thetaBool == 0: z = ds_ecco['dep'].values - if tIter == 0: - theta = ds_ecco['THETA'][:].values - nt,nz,nx,ny = np.shape(theta) - theta = theta.reshape(nz, nt, nx * ny) - else: - th = ds_ecco['THETA'][:].values - nt,nz,nx,ny = np.shape(th) - th = th.reshape(nz, nt, nx * ny) - theta = np.concatenate((theta,th), axis=2) + theta = ds_ecco['THETA'][:].values # ECCO 'land' variable is equivalent to MALI's orig3dOceanMask, despite the counterintuitive name. # Only need to call this once per tile, but needs to be using a theta/salinity file because it includes depth # No need to call this if only interpolating SIarea - if tIter == 0 and o3dmBool == 0: + if o3dmBool == 0: orig3dOceanMask = ds_ecco['land'].values - nz,nx,ny = np.shape(orig3dOceanMask) - orig3dOceanMask = orig3dOceanMask.reshape(nz, nx * ny) od3mBool = 1 - elif tIter > 0 and o3dmBool == 0: - o3dm = ds_ecco['land'].values - nz,nx,ny = np.shape(o3dm) - o3dm = o3dm.reshape(nz, nx * ny) - orig3dOceanMask = np.concatenate((orig3dOceanMask, o3dm), axis=1) - o3dmBool = 1 - # set MALI invalid value - boolZ,boolXY = np.where(orig3dOceanMask == 0) - for iZ,iXY in zip(boolZ,boolXY): - theta[iZ,:,iXY] = 1e36 + boolZ,boolX,boolY = np.where(orig3dOceanMask == 0) + for iZ,iX,iY in zip(boolZ,boolX,boolY): + theta[:,iZ,iX,iY] = 1e36 thetaBool = 1 # Repeat for sea ice fraction if needed elif 'SIarea' in file and 'SIarea' in self.options.eccoVars: - siaBool = 1 validFile = self.options.eccoDir + '/SIarea.' + paddedTile + '.nc' ds_ecco = xr.open_dataset(validFile) - if tIter == 0: - sia = ds_ecco['SIarea'][:].values - nt,nx,ny = np.shape(sia) - sia = sia.reshape(nt, nx * ny) - else: - si = ds_ecco['SIarea'][:].values - nt,nx,ny = np.shape(si) - si = si.reshape(nt, nx * ny) - sia = np.concatenate((sia,si), axis=1) + + sia = ds_ecco['SIarea'][:].values + + siaBool = 1 else : raise ValueError(f"eccoDir: {self.options.eccoDir}, paddedTile: {paddedTile}") @@ -165,184 +125,162 @@ def merge_MITgcm_files_to_unstruct_grid(self): # Doesn't matter which file these come from, just that there is only one set per tile, so just use most # recent valid file in each tile - if tIter == 0: - if thetaBool == 1 : - validFile = self.options.eccoDir + '/THETA.' + paddedTile + '.nc' - elif saltBool == 1 : - validFile = self.options.eccoDir + '/SALT.' + paddedTile + '.nc' - elif siaBool == 1: - validFile = self.options.eccoDir + '/SIarea.' + paddedTile + '.nc' - else : - raise ValueError("No variable bool found") - - ds_ecco = xr.open_dataset(validFile) - lon = ds_ecco['lon'][:].values.flatten() - lat = ds_ecco['lat'][:].values.flatten() - else: - ds_ecco = xr.open_dataset(validFile) - ln = ds_ecco['lon'][:].values.flatten() - lon = np.concatenate((lon, ln)) - lt = ds_ecco['lat'][:].values.flatten() - lat = np.concatenate((lat,lt)) - - tIter = tIter + 1 + lon = ds_ecco['lon'][:].values.flatten() + lat = ds_ecco['lat'][:].values.flatten() - # combine into new netcdf. Need to load depth and time from original ecco files. Assuming time/depth is consistent - # across all files, just load info from most recent - - # Establish dataset, and fill in standard dimensions. other dimensions are added on an as-needed basis - ds_unstruct = xr.Dataset() - ds_unstruct.expand_dims(["nCells", "Time"]) - - # Translate ECCO time to MALI xtime format as save to dataset - time = ds_ecco['tim'].values - xtime = [pd.to_datetime(dt).strftime("%Y-%m-%d_%H:%M:%S").ljust(64) for dt in time] - - DA_xtime = xr.DataArray(np.array(xtime).astype('S64'), dims=("Time")) - DA_xtime.encoding.update({"char_dim_name":"StrLen"}) - DA_xtime.attrs['long_name'] = "model time, with format 'YYYY-MM-DD_HH:MM:SS'" - DA_xtime.attrs['standard_name'] = 'time' - ds_unstruct['xtime'] = DA_xtime - - # save unstructured variables with MALI/MPAS-O names - if 'z' in locals(): - ds_unstruct.expand_dims("nISMIP6OceanLayers") - ds_unstruct.expand_dims("TWO") - z = z[z <= self.options.maxDepth] - - # create bnds array - zBnds = np.zeros((len(z), 2)) - zBnds[0, 0] = 0.0 - for i in range(1, len(z)): - zBnds[i, 0] = 0.5 * (z[i - 1] + z[i]) - for i in range(0, len(z) - 1): - zBnds[i, 1] = 0.5 * (z[i] + z[i + 1]) - zBnds[-1, 1] = z[-1] + (z[-1] - zBnds[-1, 0]) + # combine into new netcdf. Need to load depth and time from original ecco files. Assuming time/depth is consistent + # across all files, just load info from most recent - DA_z = xr.DataArray(-z.astype('float64'), dims=("nISMIP6OceanLayers")) #MALI wants negative depth coordinates - DA_z.attrs['long_name'] = 'depth' - DA_z.attrs['units'] = 'm' - ds_unstruct['ismip6shelfMelt_zOcean'] = DA_z - - DA_zBnds = xr.DataArray(-zBnds.astype('float64'), dims=("nISMIP6OceanLayers", "TWO")) #MALI wants negative depth coordinates - DA_zBnds.attrs['long_name'] = 'Bounds for ismip6 ocean layers' - DA_zBnds.attrs['units'] = 'm' - ds_unstruct['ismip6shelfMelt_zBndsOcean'] = DA_zBnds - - if 'theta' in locals(): - print(f"z shape: {z.shape}") - print(f"theta shape: {theta.shape}") - theta = theta[0:len(z),:,:] - print(f"new theta shape: {theta.shape}") - DA_theta = xr.DataArray(theta.astype('float64'), dims=("nISMIP6OceanLayers", "Time", "nCells")) - DA_theta.attrs['long_name'] = 'Potential_Temperature' - DA_theta.attrs['units'] = 'degC' - ds_unstruct['oceanTemperature'] = DA_theta - - if 'sal' in locals(): - sal = sal[0:len(z),:,:] - DA_salt = xr.DataArray(sal.astype('float64'), dims=("nISMIP6OceanLayers", "Time", "nCells")) - DA_salt.attrs['long_name'] = 'Salinity' - DA_salt.attrs['units'] = 'psu' - ds_unstruct['oceanSalinity'] = DA_salt + # Establish dataset, and fill in standard dimensions. other dimensions are added on an as-needed basis + ds_tile = xr.Dataset() + ds_tile.expand_dims(["nCells", "Time"]) + + # Translate ECCO time to MALI xtime format as save to dataset + time = ds_ecco['tim'].values + xtime = [pd.to_datetime(dt).strftime("%Y-%m-%d_%H:%M:%S").ljust(64) for dt in time] + + DA_xtime = xr.DataArray(np.array(xtime).astype('S64'), dims=("Time")) + DA_xtime.encoding.update({"char_dim_name":"StrLen"}) + DA_xtime.attrs['long_name'] = "model time, with format 'YYYY-MM-DD_HH:MM:SS'" + DA_xtime.attrs['standard_name'] = 'time' + ds_tile['xtime'] = DA_xtime + + # save unstructured variables with MALI/MPAS-O names + if 'z' in locals(): + ds_tile.expand_dims("nISMIP6OceanLayers") + ds_tile.expand_dims("TWO") + z = z[z <= self.options.maxDepth] + + # create bnds array + zBnds = np.zeros((len(z), 2)) + zBnds[0, 0] = 0.0 + for i in range(1, len(z)): + zBnds[i, 0] = 0.5 * (z[i - 1] + z[i]) + for i in range(0, len(z) - 1): + zBnds[i, 1] = 0.5 * (z[i] + z[i + 1]) + zBnds[-1, 1] = z[-1] + (z[-1] - zBnds[-1, 0]) - if 'sia' in locals(): - DA_sia = xr.DataArray(sia.astype('float64'), dims=("Time", "nCells")) - DA_sia.attrs['long_name'] = 'Sea ice fractional ice-covered area [0 to 1]' - DA_sia.attrs['units'] = 'm2/m2' - ds_unstruct['iceAreaCell'] = DA_sia + DA_z = xr.DataArray(-z.astype('float64'), dims=("nISMIP6OceanLayers")) #MALI wants negative depth coordinates + DA_z.attrs['long_name'] = 'depth' + DA_z.attrs['units'] = 'm' + ds_tile['ismip6shelfMelt_zOcean'] = DA_z - if 'orig3dOceanMask' in locals(): - orig3dOceanMask = orig3dOceanMask[0:len(z),:] - DA_o3dm = xr.DataArray(orig3dOceanMask.astype('int32'), dims=("nISMIP6OceanLayers", "nCells")) - DA_o3dm.attrs['long_name'] = ("3D mask of original valid ocean data. Because it is 3d," - " it can include the ocean domain inside ice-shelf cavities." - " orig3dOceanMask is altered online to include ice-dammed inland seas") - ds_unstruct['orig3dOceanMask'] = DA_o3dm - - DA_lon = xr.DataArray(lon.astype('float64'), dims=("nCells")) - DA_lon.attrs['long_name'] = 'longitude' - DA_lon.attrs['standard_name'] = 'longitude' - DA_lon.attrs['units'] = 'degree_east' - ds_unstruct['lon'] = DA_lon # Keep original lon/lat names of coordinates for projection transformer + DA_zBnds = xr.DataArray(-zBnds.astype('float64'), dims=("nISMIP6OceanLayers", "TWO")) #MALI wants negative depth coordinates + DA_zBnds.attrs['long_name'] = 'Bounds for ismip6 ocean layers' + DA_zBnds.attrs['units'] = 'm' + ds_tile['ismip6shelfMelt_zBndsOcean'] = DA_zBnds + + if 'theta' in locals(): + print(f"z shape: {z.shape}") + print(f"theta shape: {theta.shape}") + theta = theta[:,0:len(z),:,:] + print(f"new theta shape: {theta.shape}") + DA_theta = xr.DataArray(theta.astype('float64'), dims=("Time","nISMIP6OceanLayers", "X", "Y")) + DA_theta.attrs['long_name'] = 'Potential_Temperature' + DA_theta.attrs['units'] = 'degC' + ds_tile['oceanTemperature'] = DA_theta + + if 'sal' in locals(): + sal = sal[:,0:len(z),:,:] + DA_salt = xr.DataArray(sal.astype('float64'), dims=("Time","nISMIP6OceanLayers", "X", "Y")) + DA_salt.attrs['long_name'] = 'Salinity' + DA_salt.attrs['units'] = 'psu' + ds_tile['oceanSalinity'] = DA_salt + + if 'sia' in locals(): + DA_sia = xr.DataArray(sia.astype('float64'), dims=("Time", "X", "Y")) + DA_sia.attrs['long_name'] = 'Sea ice fractional ice-covered area [0 to 1]' + DA_sia.attrs['units'] = 'm2/m2' + ds_tile['iceAreaCell'] = DA_sia - DA_lat = xr.DataArray(lat.astype('float64'), dims=("nCells")) - DA_lat.attrs['long_name'] = 'latitude' - DA_lat.attrs['standard_name'] = 'latitude' - DA_lat.attrs['units'] = 'degree_north' - ds_unstruct['lat'] = DA_lat - - self.ecco_unstruct = "ecco_combined_unstructured.nc" - ds_unstruct.to_netcdf(self.ecco_unstruct) - ds_unstruct.close() + if 'orig3dOceanMask' in locals(): + orig3dOceanMask = orig3dOceanMask[0:len(z),:] + DA_o3dm = xr.DataArray(orig3dOceanMask.astype('int32'), dims=("nISMIP6OceanLayers", "X", "Y")) + DA_o3dm.attrs['long_name'] = ("3D mask of original valid ocean data. Because it is 3d," + " it can include the ocean domain inside ice-shelf cavities." + " orig3dOceanMask is altered online to include ice-dammed inland seas") + ds_tile['orig3dOceanMask'] = DA_o3dm - def remap_ecco_to_mali(self): + DA_lon = xr.DataArray(lon.astype('float64'), dims=("X")) + DA_lon.attrs['long_name'] = 'longitude' + DA_lon.attrs['standard_name'] = 'longitude' + DA_lon.attrs['units'] = 'degree_east' + ds_tile['lon'] = DA_lon # Keep original lon/lat names of coordinates for projection transformer - print("Creating Scrip Files") - ecco_scripfile = "tmp_ecco_scrip.nc" - mali_scripfile = "tmp_mali_scrip.nc" - ecco_unstruct = self.ecco_unstruct - file, ext = os.path.splitext(ecco_unstruct) - ecco_unstruct_xy = file + '.tmp' + ext + DA_lat = xr.DataArray(lat.astype('float64'), dims=("Y")) + DA_lat.attrs['long_name'] = 'latitude' + DA_lat.attrs['standard_name'] = 'latitude' + DA_lat.attrs['units'] = 'degree_north' + ds_tile['lat'] = DA_lat + + #convert ecco mesh to polar stereographic + #use gis-gimp projection, defined in mpas_tools.landice.projections + proj_string = ( + "+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" + ) + transformer = Transformer.from_crs("EPSG:4326", proj_string, always_xy=True) + x,y = transformer.transform(lon,lat) + + DA_x['x'] = xr.DataArray(x.astype('float64'), dims=("X")) + ds_tile['x'] = DA_x + + DA_y['y'] = xr.DataArray(y.astype('float64'), dims=("Y")) + ds_tile['y'] = DA_y - #convert ecco mesh to polar stereographic - EM = xr.open_dataset(ecco_unstruct) - lon = EM['lon'].values - lat = EM['lat'].values - #use gis-gimp projection, defined in mpas_tools.landice.projections - proj_string = ( - "+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" - ) - transformer = Transformer.from_crs("EPSG:4326", proj_string, always_xy=True) - x,y = transformer.transform(lon,lat) - - EM['x'] = xr.DataArray(x,dims=("nCells")) - EM['y'] = xr.DataArray(y,dims=("nCells")) - EM.to_netcdf(ecco_unstruct_xy) - EM.close() + ecco_tile = f"ecco_combined_tile{tile}.nc" + ds_tile.to_netcdf(ecco_tile) + ds_tile.close() - args = ['create_scrip_file_from_planar_rectangular_grid', - '-i', ecco_unstruct_xy, - '-s', ecco_scripfile, - '-p', 'gis-gimp', - '-r', '1'] - check_call(args) - - scrip_from_mpas(self.options.maliMesh, mali_scripfile) + print("Creating Scrip Files") + filename = os.path.basename(ecco_tile) + base, ext = os.path.splitext(filename) + ecco_scripfile = base + '.scrip' + ext + + filename = os.path.basename(self.options.maliMesh) + base, ext = os.path.splitext(filename) + mali_scripfile = base + '.scrip' + ext + + args = ['create_scrip_file_from_planar_rectangular_grid', + '-i', ecco_tile, + '-s', ecco_scripfile, + '-p', 'gis-gimp', + '-r', '2'] + check_call(args) + + scrip_from_mpas(self.options.maliMesh, mali_scripfile) - #create mapping file - print("Creating Mapping File") - args = ['srun', - '-n', self.options.ntasks, 'ESMF_RegridWeightGen', - '--source', ecco_scripfile, - '--destination', mali_scripfile, - '--weight', self.mapping_file_name, - '--method', self.options.method, - '--netcdf4', - "--dst_regional", "--src_regional", '--ignore_unmapped'] - check_call(args) + #create mapping file + print("Creating Mapping File") + args = ['srun', + '-n', self.options.ntasks, 'ESMF_RegridWeightGen', + '--source', ecco_scripfile, + '--destination', mali_scripfile, + '--weight', self.mapping_file_name, + '--method', self.options.method, + '--netcdf4', + "--dst_regional", "--src_regional", '--ignore_unmapped'] + check_call(args) - #remap - print("Start remapping ... ") + #remap + print("Start remapping ... ") - subprocess.run(["ncatted", "-a", "_FillValue,,d,,", ecco_unstruct]) - subprocess.run(["ncpdq", "-O", "-a", "Time,nISMIP6OceanLayers,nCells", ecco_unstruct, ecco_unstruct]) + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", ecco_tile]) - print("ncremap ...") - st = time.time() - # remap the input data - args_remap = ["ncremap", - "-i", ecco_unstruct_xy, - "-o", self.options.outputFile, - "-m", self.mapping_file_name] - check_call(args_remap) + print("ncremap ...") + st = time.time() + # remap the input data + args = ["ncremap", + "-i", ecco_tile, + "-o", self.options.outputFile, + "-m", self.mapping_file_name] + check_call(args) - subprocess.run(["ncpdq", "-O", "-a", "Time,nCells,nISMIP6OceanLayers", self.options.outputFile, self.options.outputFile]) - nd = time.time() - tm = nd - st - print(f"ncremap completed: {tm} seconds") - print("Remapping completed") + subprocess.run(["ncpdq", "-O", "-a", "Time,nCells,nISMIP6OceanLayers", self.options.outputFile, self.options.outputFile]) + nd = time.time() + tm = nd - st + print(f"ncremap completed: {tm} seconds") + print("Remapping completed") # Transfer MALI mesh variables if needed if self.options.includeMeshVars: @@ -393,8 +331,6 @@ def main(): run.merge_MITgcm_files_to_unstruct_grid() - run.remap_ecco_to_mali() - if __name__ == "__main__": main() From becc717d2301212901dc95c4d399daf3c06402bc Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Thu, 16 Oct 2025 13:49:33 -0700 Subject: [PATCH 03/17] Revert "Reorder to remap each tile" This reverts commit bf4de92db9d8f40500493c15cf33ed145166cb95. --- .../mesh_tools_li/interpolate_ecco_to_mali.py | 378 ++++++++++-------- 1 file changed, 221 insertions(+), 157 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index 3ec4a6a31..3576beefc 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -39,6 +39,7 @@ def __init__(self): def merge_MITgcm_files_to_unstruct_grid(self): ncfiles = glob.glob(os.path.join(self.options.eccoDir, "*.nc")) + tIter = 0 saltBool = 0 thetaBool = 0 siaBool = 0 @@ -61,22 +62,38 @@ def merge_MITgcm_files_to_unstruct_grid(self): if saltBool == 0 and thetaBool == 0: z = ds_ecco['dep'].values - sal = ds_ecco['SALT'][:].values + if tIter == 0: + sal = ds_ecco['SALT'][:].values + nt,nz,nx,ny = np.shape(sal) + sal = sal.reshape(nz, nt, nx * ny) + else: + sl = ds_ecco['SALT'][:].values + nt,nz,nx,ny = np.shape(sl) + sl = sl.reshape(nz, nt, nx * ny) + sal = np.concatenate((sal,sl), axis=2) # ECCO 'land' variable is equivalent to MALI's orig3dOceanMask, despite the counterintuitive name. # Only need to call this once per tile, but needs to be using a theta/salinity file because it includes depth # No need to call this if only interpolating SIarea - if o3dmBool == 0: + if tIter == 0 and o3dmBool == 0: orig3dOceanMask = ds_ecco['land'].values + nz,nx,ny = np.shape(orig3dOceanMask) + orig3dOceanMask = orig3dOceanMask.reshape(nz, nx * ny) od3mBool = 1 + elif tIter > 0 and o3dmBool == 0: + o3dm = ds_ecco['land'].values + nz,nx,ny = np.shape(o3dm) + o3dm = o3dm.reshape(nz, nx * ny) + orig3dOceanMask = np.concatenate((orig3dOceanMask, o3dm), axis=1) + o3dmBool = 1 # set MALI invalid value - boolZ,boolX,boolY = np.where(orig3dOceanMask == 0) - for iZ,iX,iY in zip(boolZ,boolX,boolY): - sal[:,iZ,iX,iY] = 1e36 + boolZ,boolXY = np.where(orig3dOceanMask == 0) + for iZ,iXY in zip(boolZ,boolXY): + sal[iZ,:,iXY] = 1e36 saltBool = 1 - + # Repeat for potential temperature if needed elif 'THETA' in file and 'THETA' in self.options.eccoVars: validFile = self.options.eccoDir + '/THETA.' + paddedTile + '.nc' @@ -85,30 +102,53 @@ def merge_MITgcm_files_to_unstruct_grid(self): if saltBool == 0 and thetaBool == 0: z = ds_ecco['dep'].values - theta = ds_ecco['THETA'][:].values + if tIter == 0: + theta = ds_ecco['THETA'][:].values + nt,nz,nx,ny = np.shape(theta) + theta = theta.reshape(nz, nt, nx * ny) + else: + th = ds_ecco['THETA'][:].values + nt,nz,nx,ny = np.shape(th) + th = th.reshape(nz, nt, nx * ny) + theta = np.concatenate((theta,th), axis=2) # ECCO 'land' variable is equivalent to MALI's orig3dOceanMask, despite the counterintuitive name. # Only need to call this once per tile, but needs to be using a theta/salinity file because it includes depth # No need to call this if only interpolating SIarea - if o3dmBool == 0: + if tIter == 0 and o3dmBool == 0: orig3dOceanMask = ds_ecco['land'].values + nz,nx,ny = np.shape(orig3dOceanMask) + orig3dOceanMask = orig3dOceanMask.reshape(nz, nx * ny) od3mBool = 1 + elif tIter > 0 and o3dmBool == 0: + o3dm = ds_ecco['land'].values + nz,nx,ny = np.shape(o3dm) + o3dm = o3dm.reshape(nz, nx * ny) + orig3dOceanMask = np.concatenate((orig3dOceanMask, o3dm), axis=1) + o3dmBool = 1 + # set MALI invalid value - boolZ,boolX,boolY = np.where(orig3dOceanMask == 0) - for iZ,iX,iY in zip(boolZ,boolX,boolY): - theta[:,iZ,iX,iY] = 1e36 + boolZ,boolXY = np.where(orig3dOceanMask == 0) + for iZ,iXY in zip(boolZ,boolXY): + theta[iZ,:,iXY] = 1e36 thetaBool = 1 # Repeat for sea ice fraction if needed elif 'SIarea' in file and 'SIarea' in self.options.eccoVars: + siaBool = 1 validFile = self.options.eccoDir + '/SIarea.' + paddedTile + '.nc' ds_ecco = xr.open_dataset(validFile) - - sia = ds_ecco['SIarea'][:].values - - siaBool = 1 + if tIter == 0: + sia = ds_ecco['SIarea'][:].values + nt,nx,ny = np.shape(sia) + sia = sia.reshape(nt, nx * ny) + else: + si = ds_ecco['SIarea'][:].values + nt,nx,ny = np.shape(si) + si = si.reshape(nt, nx * ny) + sia = np.concatenate((sia,si), axis=1) else : raise ValueError(f"eccoDir: {self.options.eccoDir}, paddedTile: {paddedTile}") @@ -125,162 +165,184 @@ def merge_MITgcm_files_to_unstruct_grid(self): # Doesn't matter which file these come from, just that there is only one set per tile, so just use most # recent valid file in each tile - lon = ds_ecco['lon'][:].values.flatten() - lat = ds_ecco['lat'][:].values.flatten() + if tIter == 0: + if thetaBool == 1 : + validFile = self.options.eccoDir + '/THETA.' + paddedTile + '.nc' + elif saltBool == 1 : + validFile = self.options.eccoDir + '/SALT.' + paddedTile + '.nc' + elif siaBool == 1: + validFile = self.options.eccoDir + '/SIarea.' + paddedTile + '.nc' + else : + raise ValueError("No variable bool found") - # combine into new netcdf. Need to load depth and time from original ecco files. Assuming time/depth is consistent - # across all files, just load info from most recent - - # Establish dataset, and fill in standard dimensions. other dimensions are added on an as-needed basis - ds_tile = xr.Dataset() - ds_tile.expand_dims(["nCells", "Time"]) - - # Translate ECCO time to MALI xtime format as save to dataset - time = ds_ecco['tim'].values - xtime = [pd.to_datetime(dt).strftime("%Y-%m-%d_%H:%M:%S").ljust(64) for dt in time] - - DA_xtime = xr.DataArray(np.array(xtime).astype('S64'), dims=("Time")) - DA_xtime.encoding.update({"char_dim_name":"StrLen"}) - DA_xtime.attrs['long_name'] = "model time, with format 'YYYY-MM-DD_HH:MM:SS'" - DA_xtime.attrs['standard_name'] = 'time' - ds_tile['xtime'] = DA_xtime - - # save unstructured variables with MALI/MPAS-O names - if 'z' in locals(): - ds_tile.expand_dims("nISMIP6OceanLayers") - ds_tile.expand_dims("TWO") - z = z[z <= self.options.maxDepth] + ds_ecco = xr.open_dataset(validFile) + lon = ds_ecco['lon'][:].values.flatten() + lat = ds_ecco['lat'][:].values.flatten() + else: + ds_ecco = xr.open_dataset(validFile) + ln = ds_ecco['lon'][:].values.flatten() + lon = np.concatenate((lon, ln)) + lt = ds_ecco['lat'][:].values.flatten() + lat = np.concatenate((lat,lt)) + + tIter = tIter + 1 + + # combine into new netcdf. Need to load depth and time from original ecco files. Assuming time/depth is consistent + # across all files, just load info from most recent + + # Establish dataset, and fill in standard dimensions. other dimensions are added on an as-needed basis + ds_unstruct = xr.Dataset() + ds_unstruct.expand_dims(["nCells", "Time"]) + + # Translate ECCO time to MALI xtime format as save to dataset + time = ds_ecco['tim'].values + xtime = [pd.to_datetime(dt).strftime("%Y-%m-%d_%H:%M:%S").ljust(64) for dt in time] + + DA_xtime = xr.DataArray(np.array(xtime).astype('S64'), dims=("Time")) + DA_xtime.encoding.update({"char_dim_name":"StrLen"}) + DA_xtime.attrs['long_name'] = "model time, with format 'YYYY-MM-DD_HH:MM:SS'" + DA_xtime.attrs['standard_name'] = 'time' + ds_unstruct['xtime'] = DA_xtime + + # save unstructured variables with MALI/MPAS-O names + if 'z' in locals(): + ds_unstruct.expand_dims("nISMIP6OceanLayers") + ds_unstruct.expand_dims("TWO") + z = z[z <= self.options.maxDepth] + + # create bnds array + zBnds = np.zeros((len(z), 2)) + zBnds[0, 0] = 0.0 + for i in range(1, len(z)): + zBnds[i, 0] = 0.5 * (z[i - 1] + z[i]) + for i in range(0, len(z) - 1): + zBnds[i, 1] = 0.5 * (z[i] + z[i + 1]) + zBnds[-1, 1] = z[-1] + (z[-1] - zBnds[-1, 0]) - # create bnds array - zBnds = np.zeros((len(z), 2)) - zBnds[0, 0] = 0.0 - for i in range(1, len(z)): - zBnds[i, 0] = 0.5 * (z[i - 1] + z[i]) - for i in range(0, len(z) - 1): - zBnds[i, 1] = 0.5 * (z[i] + z[i + 1]) - zBnds[-1, 1] = z[-1] + (z[-1] - zBnds[-1, 0]) + DA_z = xr.DataArray(-z.astype('float64'), dims=("nISMIP6OceanLayers")) #MALI wants negative depth coordinates + DA_z.attrs['long_name'] = 'depth' + DA_z.attrs['units'] = 'm' + ds_unstruct['ismip6shelfMelt_zOcean'] = DA_z + + DA_zBnds = xr.DataArray(-zBnds.astype('float64'), dims=("nISMIP6OceanLayers", "TWO")) #MALI wants negative depth coordinates + DA_zBnds.attrs['long_name'] = 'Bounds for ismip6 ocean layers' + DA_zBnds.attrs['units'] = 'm' + ds_unstruct['ismip6shelfMelt_zBndsOcean'] = DA_zBnds + + if 'theta' in locals(): + print(f"z shape: {z.shape}") + print(f"theta shape: {theta.shape}") + theta = theta[0:len(z),:,:] + print(f"new theta shape: {theta.shape}") + DA_theta = xr.DataArray(theta.astype('float64'), dims=("nISMIP6OceanLayers", "Time", "nCells")) + DA_theta.attrs['long_name'] = 'Potential_Temperature' + DA_theta.attrs['units'] = 'degC' + ds_unstruct['oceanTemperature'] = DA_theta + + if 'sal' in locals(): + sal = sal[0:len(z),:,:] + DA_salt = xr.DataArray(sal.astype('float64'), dims=("nISMIP6OceanLayers", "Time", "nCells")) + DA_salt.attrs['long_name'] = 'Salinity' + DA_salt.attrs['units'] = 'psu' + ds_unstruct['oceanSalinity'] = DA_salt - DA_z = xr.DataArray(-z.astype('float64'), dims=("nISMIP6OceanLayers")) #MALI wants negative depth coordinates - DA_z.attrs['long_name'] = 'depth' - DA_z.attrs['units'] = 'm' - ds_tile['ismip6shelfMelt_zOcean'] = DA_z + if 'sia' in locals(): + DA_sia = xr.DataArray(sia.astype('float64'), dims=("Time", "nCells")) + DA_sia.attrs['long_name'] = 'Sea ice fractional ice-covered area [0 to 1]' + DA_sia.attrs['units'] = 'm2/m2' + ds_unstruct['iceAreaCell'] = DA_sia - DA_zBnds = xr.DataArray(-zBnds.astype('float64'), dims=("nISMIP6OceanLayers", "TWO")) #MALI wants negative depth coordinates - DA_zBnds.attrs['long_name'] = 'Bounds for ismip6 ocean layers' - DA_zBnds.attrs['units'] = 'm' - ds_tile['ismip6shelfMelt_zBndsOcean'] = DA_zBnds - - if 'theta' in locals(): - print(f"z shape: {z.shape}") - print(f"theta shape: {theta.shape}") - theta = theta[:,0:len(z),:,:] - print(f"new theta shape: {theta.shape}") - DA_theta = xr.DataArray(theta.astype('float64'), dims=("Time","nISMIP6OceanLayers", "X", "Y")) - DA_theta.attrs['long_name'] = 'Potential_Temperature' - DA_theta.attrs['units'] = 'degC' - ds_tile['oceanTemperature'] = DA_theta - - if 'sal' in locals(): - sal = sal[:,0:len(z),:,:] - DA_salt = xr.DataArray(sal.astype('float64'), dims=("Time","nISMIP6OceanLayers", "X", "Y")) - DA_salt.attrs['long_name'] = 'Salinity' - DA_salt.attrs['units'] = 'psu' - ds_tile['oceanSalinity'] = DA_salt - - if 'sia' in locals(): - DA_sia = xr.DataArray(sia.astype('float64'), dims=("Time", "X", "Y")) - DA_sia.attrs['long_name'] = 'Sea ice fractional ice-covered area [0 to 1]' - DA_sia.attrs['units'] = 'm2/m2' - ds_tile['iceAreaCell'] = DA_sia + if 'orig3dOceanMask' in locals(): + orig3dOceanMask = orig3dOceanMask[0:len(z),:] + DA_o3dm = xr.DataArray(orig3dOceanMask.astype('int32'), dims=("nISMIP6OceanLayers", "nCells")) + DA_o3dm.attrs['long_name'] = ("3D mask of original valid ocean data. Because it is 3d," + " it can include the ocean domain inside ice-shelf cavities." + " orig3dOceanMask is altered online to include ice-dammed inland seas") + ds_unstruct['orig3dOceanMask'] = DA_o3dm - if 'orig3dOceanMask' in locals(): - orig3dOceanMask = orig3dOceanMask[0:len(z),:] - DA_o3dm = xr.DataArray(orig3dOceanMask.astype('int32'), dims=("nISMIP6OceanLayers", "X", "Y")) - DA_o3dm.attrs['long_name'] = ("3D mask of original valid ocean data. Because it is 3d," - " it can include the ocean domain inside ice-shelf cavities." - " orig3dOceanMask is altered online to include ice-dammed inland seas") - ds_tile['orig3dOceanMask'] = DA_o3dm + DA_lon = xr.DataArray(lon.astype('float64'), dims=("nCells")) + DA_lon.attrs['long_name'] = 'longitude' + DA_lon.attrs['standard_name'] = 'longitude' + DA_lon.attrs['units'] = 'degree_east' + ds_unstruct['lon'] = DA_lon # Keep original lon/lat names of coordinates for projection transformer - DA_lon = xr.DataArray(lon.astype('float64'), dims=("X")) - DA_lon.attrs['long_name'] = 'longitude' - DA_lon.attrs['standard_name'] = 'longitude' - DA_lon.attrs['units'] = 'degree_east' - ds_tile['lon'] = DA_lon # Keep original lon/lat names of coordinates for projection transformer + DA_lat = xr.DataArray(lat.astype('float64'), dims=("nCells")) + DA_lat.attrs['long_name'] = 'latitude' + DA_lat.attrs['standard_name'] = 'latitude' + DA_lat.attrs['units'] = 'degree_north' + ds_unstruct['lat'] = DA_lat + + self.ecco_unstruct = "ecco_combined_unstructured.nc" + ds_unstruct.to_netcdf(self.ecco_unstruct) + ds_unstruct.close() - DA_lat = xr.DataArray(lat.astype('float64'), dims=("Y")) - DA_lat.attrs['long_name'] = 'latitude' - DA_lat.attrs['standard_name'] = 'latitude' - DA_lat.attrs['units'] = 'degree_north' - ds_tile['lat'] = DA_lat - - #convert ecco mesh to polar stereographic - #use gis-gimp projection, defined in mpas_tools.landice.projections - proj_string = ( - "+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" - ) - transformer = Transformer.from_crs("EPSG:4326", proj_string, always_xy=True) - x,y = transformer.transform(lon,lat) - - DA_x['x'] = xr.DataArray(x.astype('float64'), dims=("X")) - ds_tile['x'] = DA_x - - DA_y['y'] = xr.DataArray(y.astype('float64'), dims=("Y")) - ds_tile['y'] = DA_y + def remap_ecco_to_mali(self): - ecco_tile = f"ecco_combined_tile{tile}.nc" - ds_tile.to_netcdf(ecco_tile) - ds_tile.close() + print("Creating Scrip Files") + ecco_scripfile = "tmp_ecco_scrip.nc" + mali_scripfile = "tmp_mali_scrip.nc" + ecco_unstruct = self.ecco_unstruct + file, ext = os.path.splitext(ecco_unstruct) + ecco_unstruct_xy = file + '.tmp' + ext - print("Creating Scrip Files") - filename = os.path.basename(ecco_tile) - base, ext = os.path.splitext(filename) - ecco_scripfile = base + '.scrip' + ext - - filename = os.path.basename(self.options.maliMesh) - base, ext = os.path.splitext(filename) - mali_scripfile = base + '.scrip' + ext - - args = ['create_scrip_file_from_planar_rectangular_grid', - '-i', ecco_tile, - '-s', ecco_scripfile, - '-p', 'gis-gimp', - '-r', '2'] - check_call(args) - - scrip_from_mpas(self.options.maliMesh, mali_scripfile) + #convert ecco mesh to polar stereographic + EM = xr.open_dataset(ecco_unstruct) + lon = EM['lon'].values + lat = EM['lat'].values + #use gis-gimp projection, defined in mpas_tools.landice.projections + proj_string = ( + "+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" + ) + transformer = Transformer.from_crs("EPSG:4326", proj_string, always_xy=True) + x,y = transformer.transform(lon,lat) + + EM['x'] = xr.DataArray(x,dims=("nCells")) + EM['y'] = xr.DataArray(y,dims=("nCells")) + EM.to_netcdf(ecco_unstruct_xy) + EM.close() + + args = ['create_scrip_file_from_planar_rectangular_grid', + '-i', ecco_unstruct_xy, + '-s', ecco_scripfile, + '-p', 'gis-gimp', + '-r', '1'] + check_call(args) + + scrip_from_mpas(self.options.maliMesh, mali_scripfile) - #create mapping file - print("Creating Mapping File") - args = ['srun', - '-n', self.options.ntasks, 'ESMF_RegridWeightGen', - '--source', ecco_scripfile, - '--destination', mali_scripfile, - '--weight', self.mapping_file_name, - '--method', self.options.method, - '--netcdf4', - "--dst_regional", "--src_regional", '--ignore_unmapped'] - check_call(args) + #create mapping file + print("Creating Mapping File") + args = ['srun', + '-n', self.options.ntasks, 'ESMF_RegridWeightGen', + '--source', ecco_scripfile, + '--destination', mali_scripfile, + '--weight', self.mapping_file_name, + '--method', self.options.method, + '--netcdf4', + "--dst_regional", "--src_regional", '--ignore_unmapped'] + check_call(args) - #remap - print("Start remapping ... ") + #remap + print("Start remapping ... ") - subprocess.run(["ncatted", "-a", "_FillValue,,d,,", ecco_tile]) + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", ecco_unstruct]) + subprocess.run(["ncpdq", "-O", "-a", "Time,nISMIP6OceanLayers,nCells", ecco_unstruct, ecco_unstruct]) - print("ncremap ...") - st = time.time() - # remap the input data - args = ["ncremap", - "-i", ecco_tile, - "-o", self.options.outputFile, - "-m", self.mapping_file_name] - check_call(args) + print("ncremap ...") + st = time.time() + # remap the input data + args_remap = ["ncremap", + "-i", ecco_unstruct_xy, + "-o", self.options.outputFile, + "-m", self.mapping_file_name] + check_call(args_remap) - subprocess.run(["ncpdq", "-O", "-a", "Time,nCells,nISMIP6OceanLayers", self.options.outputFile, self.options.outputFile]) - nd = time.time() - tm = nd - st - print(f"ncremap completed: {tm} seconds") - print("Remapping completed") + subprocess.run(["ncpdq", "-O", "-a", "Time,nCells,nISMIP6OceanLayers", self.options.outputFile, self.options.outputFile]) + nd = time.time() + tm = nd - st + print(f"ncremap completed: {tm} seconds") + print("Remapping completed") # Transfer MALI mesh variables if needed if self.options.includeMeshVars: @@ -331,6 +393,8 @@ def main(): run.merge_MITgcm_files_to_unstruct_grid() + run.remap_ecco_to_mali() + if __name__ == "__main__": main() From bb750964d9be21c50dead6be7494c9b408a70121 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Thu, 30 Oct 2025 14:37:55 -0700 Subject: [PATCH 04/17] Create scrip file Incorporates ECCO grid file to create custom scrip file for ECCO format --- .../mesh_tools_li/interpolate_ecco_to_mali.py | 248 +++++++++++++----- 1 file changed, 181 insertions(+), 67 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index 3576beefc..1cd440dc2 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -48,12 +48,122 @@ def merge_MITgcm_files_to_unstruct_grid(self): for tile in self.options.tileNums: print(f"Processing Tile {tile}") o3dmBool = 0 + gridBool = 0 for file in ncfiles: print(f"Processing File {file}") paddedTile = str(tile).zfill(4) print(f"Padded Tile: {paddedTile}") if paddedTile in file: print("Padded tile is in file") + + if gridBool == 0: + # start processing with grid files + gridFile = self.options.eccoDir + '/GRID.' + paddedTile + '.nc' + grd = xr.open_dataset(gridFile) + + # identify cell centers and corners + # lat/lon centers + XC = grd['XC'][:].values + YC = grd['YC'][:].values + + # lat/lon corners + XG = grd['XG'][:].values + YG = grd['YG'][:].values + DXG = grd['DXG'][:].values + + Xse_corner = np.zeros((XC.shape)) + Xsw_corner = np.zeros((XC.shape)) + Xnw_corner = np.zeros((XC.shape)) + Xne_corner = np.zeros((XC.shape)) + + Yse_corner = np.zeros((YC.shape)) + Ysw_corner = np.zeros((YC.shape)) + Ynw_corner = np.zeros((YC.shape)) + Yne_corner = np.zeros((YC.shape)) + + nx,ny = XC.shape + for i in np.arange(0,nx-1): + for ii in np.arange(0,ny-1): + Xse_corner[i,ii] = XG[i,ii] + Xsw_corner[i,ii]= XG[i+1,ii] + Xnw_corner[i,ii]= XG[i+1,ii+1] + Xne_corner[i,ii] = XG[i,ii+1] + + Yse_corner[i,ii]= YG[i,ii] + Ysw_corner[i,ii]= YG[i+1,ii] + Ynw_corner[i,ii]= YG[i+1,ii+1] + Yne_corner[i,ii]= YG[i,ii+1] + + # Approximate corners on leftmost x edge + Xse_corner[nx,ii] = XG[nx,ii] + Xsw_corner[nx,ii] = DXG[nx,ii] + XG[nx,ii] + Xnw_corner[nx,ii] = DXG[nx,ii+1] + XG[nx,ii+1] + Xne_corner[nx,ii] = XG[nx,ii+1] + + Yse_corner[nx,ii] = YG[nx,ii] + Ysw_corner[nx,ii] = YG[nx,ii] + Ynw_corner[nx,ii] = YG[nx,ii+1] + Yne_corner[nx,ii] = YG[nx,ii+1] + + # Approximate corners on uppermost y edge + Xse_corner[i,ny] = XG[i,ny] + Xsw_corner[i,ny] = XG[i+1,ny] + Xnw_corner[i,ny] = XG[i+1,ny] + Xne_corner[i,ny] = XG[i,ny] + + Yse_corner[i,ny] = YG[i,ny] + Ysw_corner[i,ny] = YG[i+1,ny] + Ynw_corner[i,ny] = YG[i+1,ny] + DYG[i,ny] + Yne_corner[i,ny] = YG[i,ny] + DYG[i-1,ny] + + #Approximate corners of uppermost/leftmost cell + Xse_corner[nx,ny] = XG[nx,ny] + Xsw_corner[nx,ny] = XG[nx,ny] + DXG[nx,ny] + Xnw_corner[nx,ny] = XG[nx,ny] + DXG[nx,ny] + Xne_corner[nx,ny] = XG[nx,ny] + + Yse_corner[nx,ny] = YG[nx,ny] + Ysw_corner[nx,ny] = YG[nx,ny] + Ynw_corner[nx,ny] = YG[nx,ny] + DYG[nx,ny] + Yne_corner[nx,ny] = YG[nx,ny] + DYG[nx-1,ny] + + # concatenate and combine tiles + if tIter == 0: + lon = XC.flatten() + lat = YC.flatten() + lon_corner1 = Xse_corner.flatten() + lon_corner2 = Xsw_corner.flatten() + lon_corner3 = Xnw_corner.flatten() + lon_corner4 = Xne_corner.flatten() + lat_corner1 = Yse_corner.flatten() + lat_corner2 = Ysw_corner.flatten() + lat_corner3 = Ynw_corner.flatten() + lat_corner4 = Yne_corner.flatten() + else: + ln = XC.flatten() + lon = np.concatenate((lon, ln)) + lt = YC.flatten() + lat = np.concatenate((lat,lt)) + lnc1 = Xse_corner.flatten() + lon_corner1 = np.concatenate((lon_corner1,lnc1)) + lnc2 = Xsw_corner.flatten() + lon_corner2 = np.concatenate((lon_corner2,lnc2)) + lnc3 = Xnw_corner.flatten() + lon_corner3 = np.concatenate((lon_corner3,lnc3)) + lnc4 = Xne_corner.flatten() + lon_corner4 = np.concatenate((lon_corner4,lnc4)) + + ltc1 = Yse_corner.flatten() + lat_corner1 = np.concatenate((lat_corner1,ltc1)) + ltc2 = Ysw_corner.flatten() + lat_corner2 = np.concatenate((lat_corner2,ltc2)) + ltc3 = Ynw_corner.flatten() + lat_corner3 = np.concatenate((lat_corner3,ltc3)) + ltc4 = Yne_corner.flatten() + lat_corner4 = np.concatenate((lat_corner4,ltc4)) + + gridBool = 1 + if 'SALT' in file and 'SALT' in self.options.eccoVars: validFile = self.options.eccoDir + '/SALT.' + paddedTile + '.nc' ds_ecco = xr.open_dataset(validFile) @@ -64,10 +174,12 @@ def merge_MITgcm_files_to_unstruct_grid(self): if tIter == 0: sal = ds_ecco['SALT'][:].values + # omit last row/column to be consistent with lat/lon nt,nz,nx,ny = np.shape(sal) sal = sal.reshape(nz, nt, nx * ny) else: sl = ds_ecco['SALT'][:].values + # omit last row/column to be consistent with lat/lon nt,nz,nx,ny = np.shape(sl) sl = sl.reshape(nz, nt, nx * ny) sal = np.concatenate((sal,sl), axis=2) @@ -77,11 +189,13 @@ def merge_MITgcm_files_to_unstruct_grid(self): # No need to call this if only interpolating SIarea if tIter == 0 and o3dmBool == 0: orig3dOceanMask = ds_ecco['land'].values + # omit last row/column to be consistent with lat/lon nz,nx,ny = np.shape(orig3dOceanMask) orig3dOceanMask = orig3dOceanMask.reshape(nz, nx * ny) od3mBool = 1 elif tIter > 0 and o3dmBool == 0: o3dm = ds_ecco['land'].values + # omit last row/column to be consistent with lat/lon nz,nx,ny = np.shape(o3dm) o3dm = o3dm.reshape(nz, nx * ny) orig3dOceanMask = np.concatenate((orig3dOceanMask, o3dm), axis=1) @@ -117,17 +231,18 @@ def merge_MITgcm_files_to_unstruct_grid(self): # No need to call this if only interpolating SIarea if tIter == 0 and o3dmBool == 0: orig3dOceanMask = ds_ecco['land'].values + # omit last row/column to be consistent with lat/lon nz,nx,ny = np.shape(orig3dOceanMask) orig3dOceanMask = orig3dOceanMask.reshape(nz, nx * ny) od3mBool = 1 - elif tIter > 0 and o3dmBool == 0: o3dm = ds_ecco['land'].values + # omit last row/column to be consistent with lat/lon nz,nx,ny = np.shape(o3dm) o3dm = o3dm.reshape(nz, nx * ny) orig3dOceanMask = np.concatenate((orig3dOceanMask, o3dm), axis=1) o3dmBool = 1 - + # set MALI invalid value boolZ,boolXY = np.where(orig3dOceanMask == 0) for iZ,iXY in zip(boolZ,boolXY): @@ -149,8 +264,6 @@ def merge_MITgcm_files_to_unstruct_grid(self): nt,nx,ny = np.shape(si) si = si.reshape(nt, nx * ny) sia = np.concatenate((sia,si), axis=1) - else : - raise ValueError(f"eccoDir: {self.options.eccoDir}, paddedTile: {paddedTile}") if ('SALT' in self.options.eccoVars and saltBool == 0): raise ValueError(f"SALT netcdf file for tile {tile} not found in {self.options.eccoDir}") @@ -160,30 +273,6 @@ def merge_MITgcm_files_to_unstruct_grid(self): if ('SIarea' in self.options.eccoVars and siaBool == 0): raise ValueError(f"SIarea netcdf file for tile {tile} not found in {self.options.eccoDir}") - - - # Doesn't matter which file these come from, just that there is only one set per tile, so just use most - # recent valid file in each tile - - if tIter == 0: - if thetaBool == 1 : - validFile = self.options.eccoDir + '/THETA.' + paddedTile + '.nc' - elif saltBool == 1 : - validFile = self.options.eccoDir + '/SALT.' + paddedTile + '.nc' - elif siaBool == 1: - validFile = self.options.eccoDir + '/SIarea.' + paddedTile + '.nc' - else : - raise ValueError("No variable bool found") - - ds_ecco = xr.open_dataset(validFile) - lon = ds_ecco['lon'][:].values.flatten() - lat = ds_ecco['lat'][:].values.flatten() - else: - ds_ecco = xr.open_dataset(validFile) - ln = ds_ecco['lon'][:].values.flatten() - lon = np.concatenate((lon, ln)) - lt = ds_ecco['lat'][:].values.flatten() - lat = np.concatenate((lat,lt)) tIter = tIter + 1 @@ -195,9 +284,11 @@ def merge_MITgcm_files_to_unstruct_grid(self): ds_unstruct.expand_dims(["nCells", "Time"]) # Translate ECCO time to MALI xtime format as save to dataset + dates = [] time = ds_ecco['tim'].values xtime = [pd.to_datetime(dt).strftime("%Y-%m-%d_%H:%M:%S").ljust(64) for dt in time] - + dates.append(xtime) + print(f"xtime: {np.array(xtime).astype('S64').shape}") DA_xtime = xr.DataArray(np.array(xtime).astype('S64'), dims=("Time")) DA_xtime.encoding.update({"char_dim_name":"StrLen"}) DA_xtime.attrs['long_name'] = "model time, with format 'YYYY-MM-DD_HH:MM:SS'" @@ -272,50 +363,71 @@ def merge_MITgcm_files_to_unstruct_grid(self): DA_lat.attrs['units'] = 'degree_north' ds_unstruct['lat'] = DA_lat + DA_lat = xr.DataArray(lat.astype('float64'), dims=("nCells")) + DA_lat.attrs['long_name'] = 'latitude' + DA_lat.attrs['standard_name'] = 'latitude' + DA_lat.attrs['units'] = 'degree_north' + ds_unstruct['lat'] = DA_lat + self.ecco_unstruct = "ecco_combined_unstructured.nc" ds_unstruct.to_netcdf(self.ecco_unstruct) ds_unstruct.close() + print("ECCO restructuring complete") + + print("Creating ECCO scrip file") + # Create scrip file for combined unstructured grid + grid_corner_lon = np.zeros((len(lon_corner1),4)) + grid_corner_lat = np.zeros((len(lat_corner1),4)) + grid_center_lon = np.zeros((len(lon),)) + grid_center_lat = np.zeros((len(lat),)) + grid_corner_lon[:,0] = lon_corner1 + grid_corner_lon[:,1] = lon_corner2 + grid_corner_lon[:,2] = lon_corner3 + grid_corner_lon[:,3] = lon_corner4 + grid_corner_lat[:,0] = lat_corner1 + grid_corner_lat[:,1] = lat_corner2 + grid_corner_lat[:,2] = lat_corner3 + grid_corner_lat[:,3] = lat_corner4 + grid_center_lon = lon + grid_center_lat = lat + grid_imask = np.ones((len(lon),)) # assume using all points + grid_dims = np.array(lon.shape) - def remap_ecco_to_mali(self): + scrip = xr.Dataset() + scrip['grid_corner_lon'] = xr.DataArray(grid_corner_lon.astype('float64'), + dims=('grid_size','grid_corners'), + attrs={'units':'degrees'}) + scrip['grid_corner_lat'] = xr.DataArray(grid_corner_lat.astype('float64'), + dims=('grid_size','grid_corners'), + attrs={'units':'degrees'}) + scrip['grid_center_lon'] = xr.DataArray(grid_center_lon.astype('float64'), + dims=('grid_size'), + attrs={'units':'degrees'}) + scrip['grid_center_lat'] = xr.DataArray(grid_center_lat.astype('float64'), + dims=('grid_size'), + attrs={'units':'degrees'}) + scrip['grid_imask'] = xr.DataArray(grid_imask.astype('int32'), + dims=('grid_size'), + attrs={'units':'unitless'}) + scrip['grid_dims'] = xr.DataArray(grid_dims, + dims=('grid_rank'), + attrs={'units':'unitless'}) + self.ecco_scripfile = 'ecco_combined_unstructured.scrip.nc' + scrip.to_netcdf(self.ecco_scripfile) + scrip.close() + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", self.ecco_scripfile]) + print("ECCO scrip file complete") - print("Creating Scrip Files") - ecco_scripfile = "tmp_ecco_scrip.nc" - mali_scripfile = "tmp_mali_scrip.nc" - ecco_unstruct = self.ecco_unstruct - file, ext = os.path.splitext(ecco_unstruct) - ecco_unstruct_xy = file + '.tmp' + ext - - #convert ecco mesh to polar stereographic - EM = xr.open_dataset(ecco_unstruct) - lon = EM['lon'].values - lat = EM['lat'].values - #use gis-gimp projection, defined in mpas_tools.landice.projections - proj_string = ( - "+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" - ) - transformer = Transformer.from_crs("EPSG:4326", proj_string, always_xy=True) - x,y = transformer.transform(lon,lat) - - EM['x'] = xr.DataArray(x,dims=("nCells")) - EM['y'] = xr.DataArray(y,dims=("nCells")) - EM.to_netcdf(ecco_unstruct_xy) - EM.close() - - args = ['create_scrip_file_from_planar_rectangular_grid', - '-i', ecco_unstruct_xy, - '-s', ecco_scripfile, - '-p', 'gis-gimp', - '-r', '1'] - check_call(args) - + def remap_ecco_to_mali(self): + print("Creating MALI scrip file") + mali_scripfile = "mali.scrip.nc" scrip_from_mpas(self.options.maliMesh, mali_scripfile) #create mapping file print("Creating Mapping File") args = ['srun', '-n', self.options.ntasks, 'ESMF_RegridWeightGen', - '--source', ecco_scripfile, + '--source', self.ecco_scripfile, '--destination', mali_scripfile, '--weight', self.mapping_file_name, '--method', self.options.method, @@ -326,14 +438,14 @@ def remap_ecco_to_mali(self): #remap print("Start remapping ... ") - subprocess.run(["ncatted", "-a", "_FillValue,,d,,", ecco_unstruct]) - subprocess.run(["ncpdq", "-O", "-a", "Time,nISMIP6OceanLayers,nCells", ecco_unstruct, ecco_unstruct]) + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", self.ecco_unstruct]) + subprocess.run(["ncpdq", "-O", "-a", "Time,nISMIP6OceanLayers,nCells", self.ecco_unstruct, self.ecco_unstruct]) print("ncremap ...") st = time.time() # remap the input data args_remap = ["ncremap", - "-i", ecco_unstruct_xy, + "-i", self.ecco_unstruct, "-o", self.options.outputFile, "-m", self.mapping_file_name] check_call(args_remap) @@ -347,7 +459,7 @@ def remap_ecco_to_mali(self): # Transfer MALI mesh variables if needed if self.options.includeMeshVars: ds_mali = xr.open_dataset(self.options.maliMesh, decode_times=False, decode_cf=False) - ds_out = xr.open_dataset(self.options.outputFile, decode_times=False, decode_cf=False) + ds_out = xr.open_dataset(self.options.outputFile, decode_times=False, decode_cf=False, mode='r+') ds_out['angleEdge'] = ds_mali['angleEdge'] ds_out['areaCell'] = ds_mali['areaCell'] ds_out['areaTriangle'] = ds_mali['areaTriangle'] @@ -384,9 +496,11 @@ def remap_ecco_to_mali(self): ds_out['zCell'] = ds_mali['zCell'] ds_out['zEdge'] = ds_mali['zEdge'] ds_out['zVertex'] = ds_mali['zVertex'] - ds_out.to_netcdf(self.options.outputFile, mode='w') + outFileWithMeshVars = self.options.outputFile[0:-3] + '.withMeshVars.nc' + ds_out.to_netcdf(outFileWithMeshVars) ds_out.close() ds_mali.close() + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", outFileWithMeshVars]) def main(): run = eccoToMaliInterp() From 72b067a182e06425159220279479dcadcf961411 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Tue, 4 Nov 2025 10:34:33 -0800 Subject: [PATCH 05/17] Hardcode edge cell corners 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).. --- .../mesh_tools_li/interpolate_ecco_to_mali.py | 412 +++++++++--------- 1 file changed, 208 insertions(+), 204 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index 1cd440dc2..878e39443 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -10,7 +10,6 @@ from mpas_tools.logging import check_call from mpas_tools.scrip.from_mpas import scrip_from_mpas from mpas_tools.ocean.depth import compute_zmid -from pyproj import Proj, Transformer from argparse import ArgumentParser import time import cftime @@ -25,7 +24,6 @@ def __init__(self): description='Interpolates ECCO reanalysis data onto a MALI mesh') parser.add_argument("--eccoDir", dest="eccoDir", required=True, help="Directory where ECCO files are stored. Should be only netcdf files in that directory", metavar="FILENAME") parser.add_argument("--maliMesh", dest="maliMesh", required=True, help="MALI mesh to be interpolated onto", metavar="FILENAME") - parser.add_argument("--tileNums", type=int, nargs='+', dest="tileNums", required=True, help="ECCO grid regional tile numbers to pull data from") parser.add_argument("--eccoVars", type=str, nargs='+', dest="eccoVars", required=True, help="ECCO variables to interpolate, current options are any combination of 'THETA', 'SALT', 'SIarea'") parser.add_argument("--maxDepth", type=int, dest="maxDepth", default=1500, help="Maximum depth in meters to include in the interpolation") parser.add_argument("--ntasks", dest="ntasks", type=str, help="Number of processors to use with ESMF_RegridWeightGen", default='128') @@ -36,16 +34,16 @@ def __init__(self): self.options = args self.mapping_file_name = f'ecco_to_mali_mapping_{self.options.method}.nc' + self.tileNums = [14, 27] # Hardcoding ECCO tiles for now. Indexing is dependent on tile order, do not change def merge_MITgcm_files_to_unstruct_grid(self): ncfiles = glob.glob(os.path.join(self.options.eccoDir, "*.nc")) - tIter = 0 saltBool = 0 thetaBool = 0 siaBool = 0 #loop through MITgcm tile files and create unstructured arrays for each variable in eccoVars - for tile in self.options.tileNums: + for tile in self.tileNums: print(f"Processing Tile {tile}") o3dmBool = 0 gridBool = 0 @@ -56,114 +54,6 @@ def merge_MITgcm_files_to_unstruct_grid(self): if paddedTile in file: print("Padded tile is in file") - if gridBool == 0: - # start processing with grid files - gridFile = self.options.eccoDir + '/GRID.' + paddedTile + '.nc' - grd = xr.open_dataset(gridFile) - - # identify cell centers and corners - # lat/lon centers - XC = grd['XC'][:].values - YC = grd['YC'][:].values - - # lat/lon corners - XG = grd['XG'][:].values - YG = grd['YG'][:].values - DXG = grd['DXG'][:].values - - Xse_corner = np.zeros((XC.shape)) - Xsw_corner = np.zeros((XC.shape)) - Xnw_corner = np.zeros((XC.shape)) - Xne_corner = np.zeros((XC.shape)) - - Yse_corner = np.zeros((YC.shape)) - Ysw_corner = np.zeros((YC.shape)) - Ynw_corner = np.zeros((YC.shape)) - Yne_corner = np.zeros((YC.shape)) - - nx,ny = XC.shape - for i in np.arange(0,nx-1): - for ii in np.arange(0,ny-1): - Xse_corner[i,ii] = XG[i,ii] - Xsw_corner[i,ii]= XG[i+1,ii] - Xnw_corner[i,ii]= XG[i+1,ii+1] - Xne_corner[i,ii] = XG[i,ii+1] - - Yse_corner[i,ii]= YG[i,ii] - Ysw_corner[i,ii]= YG[i+1,ii] - Ynw_corner[i,ii]= YG[i+1,ii+1] - Yne_corner[i,ii]= YG[i,ii+1] - - # Approximate corners on leftmost x edge - Xse_corner[nx,ii] = XG[nx,ii] - Xsw_corner[nx,ii] = DXG[nx,ii] + XG[nx,ii] - Xnw_corner[nx,ii] = DXG[nx,ii+1] + XG[nx,ii+1] - Xne_corner[nx,ii] = XG[nx,ii+1] - - Yse_corner[nx,ii] = YG[nx,ii] - Ysw_corner[nx,ii] = YG[nx,ii] - Ynw_corner[nx,ii] = YG[nx,ii+1] - Yne_corner[nx,ii] = YG[nx,ii+1] - - # Approximate corners on uppermost y edge - Xse_corner[i,ny] = XG[i,ny] - Xsw_corner[i,ny] = XG[i+1,ny] - Xnw_corner[i,ny] = XG[i+1,ny] - Xne_corner[i,ny] = XG[i,ny] - - Yse_corner[i,ny] = YG[i,ny] - Ysw_corner[i,ny] = YG[i+1,ny] - Ynw_corner[i,ny] = YG[i+1,ny] + DYG[i,ny] - Yne_corner[i,ny] = YG[i,ny] + DYG[i-1,ny] - - #Approximate corners of uppermost/leftmost cell - Xse_corner[nx,ny] = XG[nx,ny] - Xsw_corner[nx,ny] = XG[nx,ny] + DXG[nx,ny] - Xnw_corner[nx,ny] = XG[nx,ny] + DXG[nx,ny] - Xne_corner[nx,ny] = XG[nx,ny] - - Yse_corner[nx,ny] = YG[nx,ny] - Ysw_corner[nx,ny] = YG[nx,ny] - Ynw_corner[nx,ny] = YG[nx,ny] + DYG[nx,ny] - Yne_corner[nx,ny] = YG[nx,ny] + DYG[nx-1,ny] - - # concatenate and combine tiles - if tIter == 0: - lon = XC.flatten() - lat = YC.flatten() - lon_corner1 = Xse_corner.flatten() - lon_corner2 = Xsw_corner.flatten() - lon_corner3 = Xnw_corner.flatten() - lon_corner4 = Xne_corner.flatten() - lat_corner1 = Yse_corner.flatten() - lat_corner2 = Ysw_corner.flatten() - lat_corner3 = Ynw_corner.flatten() - lat_corner4 = Yne_corner.flatten() - else: - ln = XC.flatten() - lon = np.concatenate((lon, ln)) - lt = YC.flatten() - lat = np.concatenate((lat,lt)) - lnc1 = Xse_corner.flatten() - lon_corner1 = np.concatenate((lon_corner1,lnc1)) - lnc2 = Xsw_corner.flatten() - lon_corner2 = np.concatenate((lon_corner2,lnc2)) - lnc3 = Xnw_corner.flatten() - lon_corner3 = np.concatenate((lon_corner3,lnc3)) - lnc4 = Xne_corner.flatten() - lon_corner4 = np.concatenate((lon_corner4,lnc4)) - - ltc1 = Yse_corner.flatten() - lat_corner1 = np.concatenate((lat_corner1,ltc1)) - ltc2 = Ysw_corner.flatten() - lat_corner2 = np.concatenate((lat_corner2,ltc2)) - ltc3 = Ynw_corner.flatten() - lat_corner3 = np.concatenate((lat_corner3,ltc3)) - ltc4 = Yne_corner.flatten() - lat_corner4 = np.concatenate((lat_corner4,ltc4)) - - gridBool = 1 - if 'SALT' in file and 'SALT' in self.options.eccoVars: validFile = self.options.eccoDir + '/SALT.' + paddedTile + '.nc' ds_ecco = xr.open_dataset(validFile) @@ -172,30 +62,34 @@ def merge_MITgcm_files_to_unstruct_grid(self): if saltBool == 0 and thetaBool == 0: z = ds_ecco['dep'].values - if tIter == 0: + if tile == 14: sal = ds_ecco['SALT'][:].values - # omit last row/column to be consistent with lat/lon + # omit rows/columns where no corner information in scrip file + sal = sal[:,:,:,1:-1] nt,nz,nx,ny = np.shape(sal) - sal = sal.reshape(nz, nt, nx * ny) - else: + sal = sal.reshape(nt, nz, nx * ny) + elif tile == 27: sl = ds_ecco['SALT'][:].values - # omit last row/column to be consistent with lat/lon + # omit rows/columns where no corner information in scrip file + sl = sl[:,:,0:-1,0:-1] nt,nz,nx,ny = np.shape(sl) - sl = sl.reshape(nz, nt, nx * ny) + sl = sl.reshape(nt, nz, nx * ny) sal = np.concatenate((sal,sl), axis=2) # ECCO 'land' variable is equivalent to MALI's orig3dOceanMask, despite the counterintuitive name. # Only need to call this once per tile, but needs to be using a theta/salinity file because it includes depth # No need to call this if only interpolating SIarea - if tIter == 0 and o3dmBool == 0: + if tile == 14 and o3dmBool == 0: orig3dOceanMask = ds_ecco['land'].values - # omit last row/column to be consistent with lat/lon + # omit rows/columns where no corner information in scrip file + orig3dOceanMask = orig3dOceanMask[:,:,1:-1] nz,nx,ny = np.shape(orig3dOceanMask) orig3dOceanMask = orig3dOceanMask.reshape(nz, nx * ny) od3mBool = 1 - elif tIter > 0 and o3dmBool == 0: + elif tile == 27 and o3dmBool == 0: o3dm = ds_ecco['land'].values - # omit last row/column to be consistent with lat/lon + # omit rows/columns where no corner information in scrip file + o3dm = o3dm[:,0:-1,0:-1] nz,nx,ny = np.shape(o3dm) o3dm = o3dm.reshape(nz, nx * ny) orig3dOceanMask = np.concatenate((orig3dOceanMask, o3dm), axis=1) @@ -204,7 +98,7 @@ def merge_MITgcm_files_to_unstruct_grid(self): # set MALI invalid value boolZ,boolXY = np.where(orig3dOceanMask == 0) for iZ,iXY in zip(boolZ,boolXY): - sal[iZ,:,iXY] = 1e36 + sal[:,iZ,iXY] = 1e36 saltBool = 1 @@ -216,28 +110,34 @@ def merge_MITgcm_files_to_unstruct_grid(self): if saltBool == 0 and thetaBool == 0: z = ds_ecco['dep'].values - if tIter == 0: + if tile == 14: theta = ds_ecco['THETA'][:].values + # omit rows/columns where no corner information in scrip file + theta = theta[:,:,:,1:-1] nt,nz,nx,ny = np.shape(theta) - theta = theta.reshape(nz, nt, nx * ny) - else: + theta = theta.reshape(nt, nz, nx * ny) + elif tile == 27: th = ds_ecco['THETA'][:].values + # omit rows/columns where no corner information in scrip file + th = th[:,:,0:-1,0:-1] nt,nz,nx,ny = np.shape(th) - th = th.reshape(nz, nt, nx * ny) + th = th.reshape(nt, nz, nx * ny) theta = np.concatenate((theta,th), axis=2) # ECCO 'land' variable is equivalent to MALI's orig3dOceanMask, despite the counterintuitive name. # Only need to call this once per tile, but needs to be using a theta/salinity file because it includes depth # No need to call this if only interpolating SIarea - if tIter == 0 and o3dmBool == 0: + if tile == 14 and o3dmBool == 0: orig3dOceanMask = ds_ecco['land'].values - # omit last row/column to be consistent with lat/lon + # omit rows/columns where no corner information in scrip file + orig3dOceanMask = orig3dOceanMask[:,:,1:-1] nz,nx,ny = np.shape(orig3dOceanMask) orig3dOceanMask = orig3dOceanMask.reshape(nz, nx * ny) od3mBool = 1 - elif tIter > 0 and o3dmBool == 0: + elif tile == 27 and o3dmBool == 0: o3dm = ds_ecco['land'].values - # omit last row/column to be consistent with lat/lon + # omit rows/columns where no corner information in scrip file + o3dm = o3dm[:,0:-1,0:-1] nz,nx,ny = np.shape(o3dm) o3dm = o3dm.reshape(nz, nx * ny) orig3dOceanMask = np.concatenate((orig3dOceanMask, o3dm), axis=1) @@ -246,7 +146,7 @@ def merge_MITgcm_files_to_unstruct_grid(self): # set MALI invalid value boolZ,boolXY = np.where(orig3dOceanMask == 0) for iZ,iXY in zip(boolZ,boolXY): - theta[iZ,:,iXY] = 1e36 + theta[:,iZ,iXY] = 1e36 thetaBool = 1 @@ -255,12 +155,16 @@ def merge_MITgcm_files_to_unstruct_grid(self): siaBool = 1 validFile = self.options.eccoDir + '/SIarea.' + paddedTile + '.nc' ds_ecco = xr.open_dataset(validFile) - if tIter == 0: + if tile == 14: sia = ds_ecco['SIarea'][:].values + # omit rows/columns where no corner information in scrip file + sia = sia[:,:,1:-1] nt,nx,ny = np.shape(sia) sia = sia.reshape(nt, nx * ny) - else: + elif tile == 27: si = ds_ecco['SIarea'][:].values + # omit rows/columns where no corner information in scrip file + si = si[:,0:-1,0:-1] nt,nx,ny = np.shape(si) si = si.reshape(nt, nx * ny) sia = np.concatenate((sia,si), axis=1) @@ -274,8 +178,6 @@ def merge_MITgcm_files_to_unstruct_grid(self): if ('SIarea' in self.options.eccoVars and siaBool == 0): raise ValueError(f"SIarea netcdf file for tile {tile} not found in {self.options.eccoDir}") - tIter = tIter + 1 - # combine into new netcdf. Need to load depth and time from original ecco files. Assuming time/depth is consistent # across all files, just load info from most recent @@ -288,8 +190,7 @@ def merge_MITgcm_files_to_unstruct_grid(self): time = ds_ecco['tim'].values xtime = [pd.to_datetime(dt).strftime("%Y-%m-%d_%H:%M:%S").ljust(64) for dt in time] dates.append(xtime) - print(f"xtime: {np.array(xtime).astype('S64').shape}") - DA_xtime = xr.DataArray(np.array(xtime).astype('S64'), dims=("Time")) + DA_xtime = xr.DataArray(np.array(xtime,dtype=object), dims=("Time")) DA_xtime.encoding.update({"char_dim_name":"StrLen"}) DA_xtime.attrs['long_name'] = "model time, with format 'YYYY-MM-DD_HH:MM:SS'" DA_xtime.attrs['standard_name'] = 'time' @@ -323,16 +224,16 @@ def merge_MITgcm_files_to_unstruct_grid(self): if 'theta' in locals(): print(f"z shape: {z.shape}") print(f"theta shape: {theta.shape}") - theta = theta[0:len(z),:,:] + theta = theta[:,0:len(z),:] print(f"new theta shape: {theta.shape}") - DA_theta = xr.DataArray(theta.astype('float64'), dims=("nISMIP6OceanLayers", "Time", "nCells")) + DA_theta = xr.DataArray(theta.astype('float64'), dims=("Time", "nISMIP6OceanLayers", "nCells")) DA_theta.attrs['long_name'] = 'Potential_Temperature' DA_theta.attrs['units'] = 'degC' ds_unstruct['oceanTemperature'] = DA_theta if 'sal' in locals(): - sal = sal[0:len(z),:,:] - DA_salt = xr.DataArray(sal.astype('float64'), dims=("nISMIP6OceanLayers", "Time", "nCells")) + sal = sal[:,0:len(z),:] + DA_salt = xr.DataArray(sal.astype('float64'), dims=("Time", "nISMIP6OceanLayers", "nCells")) DA_salt.attrs['long_name'] = 'Salinity' DA_salt.attrs['units'] = 'psu' ds_unstruct['oceanSalinity'] = DA_salt @@ -351,73 +252,11 @@ def merge_MITgcm_files_to_unstruct_grid(self): " orig3dOceanMask is altered online to include ice-dammed inland seas") ds_unstruct['orig3dOceanMask'] = DA_o3dm - DA_lon = xr.DataArray(lon.astype('float64'), dims=("nCells")) - DA_lon.attrs['long_name'] = 'longitude' - DA_lon.attrs['standard_name'] = 'longitude' - DA_lon.attrs['units'] = 'degree_east' - ds_unstruct['lon'] = DA_lon # Keep original lon/lat names of coordinates for projection transformer - - DA_lat = xr.DataArray(lat.astype('float64'), dims=("nCells")) - DA_lat.attrs['long_name'] = 'latitude' - DA_lat.attrs['standard_name'] = 'latitude' - DA_lat.attrs['units'] = 'degree_north' - ds_unstruct['lat'] = DA_lat - - DA_lat = xr.DataArray(lat.astype('float64'), dims=("nCells")) - DA_lat.attrs['long_name'] = 'latitude' - DA_lat.attrs['standard_name'] = 'latitude' - DA_lat.attrs['units'] = 'degree_north' - ds_unstruct['lat'] = DA_lat - self.ecco_unstruct = "ecco_combined_unstructured.nc" ds_unstruct.to_netcdf(self.ecco_unstruct) ds_unstruct.close() print("ECCO restructuring complete") - print("Creating ECCO scrip file") - # Create scrip file for combined unstructured grid - grid_corner_lon = np.zeros((len(lon_corner1),4)) - grid_corner_lat = np.zeros((len(lat_corner1),4)) - grid_center_lon = np.zeros((len(lon),)) - grid_center_lat = np.zeros((len(lat),)) - grid_corner_lon[:,0] = lon_corner1 - grid_corner_lon[:,1] = lon_corner2 - grid_corner_lon[:,2] = lon_corner3 - grid_corner_lon[:,3] = lon_corner4 - grid_corner_lat[:,0] = lat_corner1 - grid_corner_lat[:,1] = lat_corner2 - grid_corner_lat[:,2] = lat_corner3 - grid_corner_lat[:,3] = lat_corner4 - grid_center_lon = lon - grid_center_lat = lat - grid_imask = np.ones((len(lon),)) # assume using all points - grid_dims = np.array(lon.shape) - - scrip = xr.Dataset() - scrip['grid_corner_lon'] = xr.DataArray(grid_corner_lon.astype('float64'), - dims=('grid_size','grid_corners'), - attrs={'units':'degrees'}) - scrip['grid_corner_lat'] = xr.DataArray(grid_corner_lat.astype('float64'), - dims=('grid_size','grid_corners'), - attrs={'units':'degrees'}) - scrip['grid_center_lon'] = xr.DataArray(grid_center_lon.astype('float64'), - dims=('grid_size'), - attrs={'units':'degrees'}) - scrip['grid_center_lat'] = xr.DataArray(grid_center_lat.astype('float64'), - dims=('grid_size'), - attrs={'units':'degrees'}) - scrip['grid_imask'] = xr.DataArray(grid_imask.astype('int32'), - dims=('grid_size'), - attrs={'units':'unitless'}) - scrip['grid_dims'] = xr.DataArray(grid_dims, - dims=('grid_rank'), - attrs={'units':'unitless'}) - self.ecco_scripfile = 'ecco_combined_unstructured.scrip.nc' - scrip.to_netcdf(self.ecco_scripfile) - scrip.close() - subprocess.run(["ncatted", "-a", "_FillValue,,d,,", self.ecco_scripfile]) - print("ECCO scrip file complete") - def remap_ecco_to_mali(self): print("Creating MALI scrip file") mali_scripfile = "mali.scrip.nc" @@ -501,10 +340,175 @@ def remap_ecco_to_mali(self): ds_out.close() ds_mali.close() subprocess.run(["ncatted", "-a", "_FillValue,,d,,", outFileWithMeshVars]) + + def create_ECCO_scrip_file(self): + + print("Creating ECCO scrip file") + + for tile in self.tileNums : + # start processing with grid files + paddedTile = str(tile).zfill(4) + gridFile = self.options.eccoDir + '/GRID.' + paddedTile + '.nc' + grd = xr.open_dataset(gridFile) + + # identify cell centers and corners + # lat/lon centers + XC = grd['XC'][:].values + YC = grd['YC'][:].values + + # lat/lon corners + XG = grd['XG'][:].values + YG = grd['YG'][:].values + + Xse_corner = np.zeros((XC.shape)) + Xsw_corner = np.zeros((XC.shape)) + Xnw_corner = np.zeros((XC.shape)) + Xne_corner = np.zeros((XC.shape)) + + Yse_corner = np.zeros((YC.shape)) + Ysw_corner = np.zeros((YC.shape)) + Ynw_corner = np.zeros((YC.shape)) + Yne_corner = np.zeros((YC.shape)) + + nx,ny = XC.shape + for i in np.arange(0,nx-1): + for ii in np.arange(0,ny-1): + Xse_corner[i,ii] = XG[i,ii] + Xsw_corner[i,ii] = XG[i+1,ii] + Xnw_corner[i,ii] = XG[i+1,ii+1] + Xne_corner[i,ii] = XG[i,ii+1] + + Yse_corner[i,ii] = YG[i,ii] + Ysw_corner[i,ii] = YG[i+1,ii] + Ynw_corner[i,ii] = YG[i+1,ii+1] + Yne_corner[i,ii] = YG[i,ii+1] + + if tile == 14: + # Stitch together corner info from neighboring tile + gridFile = self.options.eccoDir + '/GRID.0027.nc' + grd27 = xr.open_dataset(gridFile) + xg27 = grd27['XG'][:].values + yg27 = grd27['YG'][:].values + xg = np.flip(xg27[1:,0]) + yg = np.flip(yg27[1:,0]) + + for i in np.arange(0,ny-1): + Xse_corner[nx-1,i] = XG[nx-1,i] + Xne_corner[nx-1,i] = XG[nx-1,i+1] + Xnw_corner[nx-1,i] = xg[i] + Xsw_corner[nx-1,i] = xg[i-1] + + Yse_corner[nx-1,i] = YG[nx-1,i] + Yne_corner[nx-1,i] = YG[nx-1,i+1] + Ynw_corner[nx-1,i] = yg[i] + Ysw_corner[nx-1,i] = yg[i-1] + + # Remove edge cells without corner information + XC = XC[:,1:-1] + YC = YC[:,1:-1] + Xse_corner = Xse_corner[:,1:-1] + Xne_corner = Xne_corner[:,1:-1] + Xnw_corner = Xnw_corner[:,1:-1] + Xsw_corner = Xsw_corner[:,1:-1] + Yse_corner = Yse_corner[:,1:-1] + Yne_corner = Yne_corner[:,1:-1] + Ynw_corner = Ynw_corner[:,1:-1] + Ysw_corner = Ysw_corner[:,1:-1] + lon = XC.flatten() + lat = YC.flatten() + lon_corner1 = Xse_corner.flatten() + lon_corner2 = Xsw_corner.flatten() + lon_corner3 = Xnw_corner.flatten() + lon_corner4 = Xne_corner.flatten() + lat_corner1 = Yse_corner.flatten() + lat_corner2 = Ysw_corner.flatten() + lat_corner3 = Ynw_corner.flatten() + lat_corner4 = Yne_corner.flatten() + + elif tile == 27 : + # Remove edge cells without corner information + XC = XC[0:-1,0:-1] + YC = YC[0:-1,0:-1] + Xse_corner = Xse_corner[0:-1,0:-1] + Xne_corner = Xne_corner[0:-1,0:-1] + Xnw_corner = Xnw_corner[0:-1,0:-1] + Xsw_corner = Xsw_corner[0:-1,0:-1] + Yse_corner = Yse_corner[0:-1,0:-1] + Yne_corner = Yne_corner[0:-1,0:-1] + Ynw_corner = Ynw_corner[0:-1,0:-1] + Ysw_corner = Ysw_corner[0:-1,0:-1] + + ln = XC.flatten() + lon = np.concatenate((lon, ln)) + lt = YC.flatten() + lat = np.concatenate((lat,lt)) + lnc1 = Xse_corner.flatten() + lon_corner1 = np.concatenate((lon_corner1,lnc1)) + lnc2 = Xsw_corner.flatten() + lon_corner2 = np.concatenate((lon_corner2,lnc2)) + lnc3 = Xnw_corner.flatten() + lon_corner3 = np.concatenate((lon_corner3,lnc3)) + lnc4 = Xne_corner.flatten() + lon_corner4 = np.concatenate((lon_corner4,lnc4)) + + ltc1 = Yse_corner.flatten() + lat_corner1 = np.concatenate((lat_corner1,ltc1)) + ltc2 = Ysw_corner.flatten() + lat_corner2 = np.concatenate((lat_corner2,ltc2)) + ltc3 = Ynw_corner.flatten() + lat_corner3 = np.concatenate((lat_corner3,ltc3)) + ltc4 = Yne_corner.flatten() + lat_corner4 = np.concatenate((lat_corner4,ltc4)) + + # Create scrip file for combined unstructured grid + grid_corner_lon = np.zeros((len(lon_corner1),4)) + grid_corner_lat = np.zeros((len(lat_corner1),4)) + grid_center_lon = np.zeros((len(lon),)) + grid_center_lat = np.zeros((len(lat),)) + grid_corner_lon[:,0] = lon_corner1 + grid_corner_lon[:,1] = lon_corner2 + grid_corner_lon[:,2] = lon_corner3 + grid_corner_lon[:,3] = lon_corner4 + grid_corner_lat[:,0] = lat_corner1 + grid_corner_lat[:,1] = lat_corner2 + grid_corner_lat[:,2] = lat_corner3 + grid_corner_lat[:,3] = lat_corner4 + grid_center_lon = lon + grid_center_lat = lat + grid_imask = np.ones((len(lon),)) # assume using all points + grid_dims = np.array(lon.shape) + + scrip = xr.Dataset() + scrip['grid_corner_lon'] = xr.DataArray(grid_corner_lon.astype('float64'), + dims=('grid_size','grid_corners'), + attrs={'units':'degrees'}) + scrip['grid_corner_lat'] = xr.DataArray(grid_corner_lat.astype('float64'), + dims=('grid_size','grid_corners'), + attrs={'units':'degrees'}) + scrip['grid_center_lon'] = xr.DataArray(grid_center_lon.astype('float64'), + dims=('grid_size'), + attrs={'units':'degrees'}) + scrip['grid_center_lat'] = xr.DataArray(grid_center_lat.astype('float64'), + dims=('grid_size'), + attrs={'units':'degrees'}) + scrip['grid_imask'] = xr.DataArray(grid_imask.astype('int32'), + dims=('grid_size'), + attrs={'units':'unitless'}) + scrip['grid_dims'] = xr.DataArray(grid_dims.astype('int32'), + dims=('grid_rank'), + attrs={'units':'unitless'}) + self.ecco_scripfile = 'ecco_combined_unstructured.scrip.nc' + scrip.to_netcdf(self.ecco_scripfile) + scrip.close() + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", self.ecco_scripfile]) + print("ECCO scrip file complete") + def main(): run = eccoToMaliInterp() + run.create_ECCO_scrip_file() + run.merge_MITgcm_files_to_unstruct_grid() run.remap_ecco_to_mali() From fc89af2a6f541b4cc2158f3dea173674213475fd Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Thu, 4 Dec 2025 11:10:31 -0800 Subject: [PATCH 06/17] Fix xtime in output file Addresses formating issue with xtime variable in initial output file. Still need to fix formatting in 'meshVars' output file. --- landice/mesh_tools_li/interpolate_ecco_to_mali.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index 878e39443..bc18da1a7 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -186,16 +186,13 @@ def merge_MITgcm_files_to_unstruct_grid(self): ds_unstruct.expand_dims(["nCells", "Time"]) # Translate ECCO time to MALI xtime format as save to dataset - dates = [] time = ds_ecco['tim'].values - xtime = [pd.to_datetime(dt).strftime("%Y-%m-%d_%H:%M:%S").ljust(64) for dt in time] - dates.append(xtime) - DA_xtime = xr.DataArray(np.array(xtime,dtype=object), dims=("Time")) - DA_xtime.encoding.update({"char_dim_name":"StrLen"}) - DA_xtime.attrs['long_name'] = "model time, with format 'YYYY-MM-DD_HH:MM:SS'" - DA_xtime.attrs['standard_name'] = 'time' + xtime_str = [pd.to_datetime(dt).strftime("%Y-%m-%d_%H:%M:%S").ljust(64) for dt in time] + xtime = np.array(xtime_str, dtype = np.dtype(('S',64))) + DA_xtime = xr.DataArray(xtime, dims=["Time"]) + DA_xtime.encoding.update({"char_dim_name": "StrLen"}) ds_unstruct['xtime'] = DA_xtime - + # save unstructured variables with MALI/MPAS-O names if 'z' in locals(): ds_unstruct.expand_dims("nISMIP6OceanLayers") @@ -335,6 +332,8 @@ def remap_ecco_to_mali(self): ds_out['zCell'] = ds_mali['zCell'] ds_out['zEdge'] = ds_mali['zEdge'] ds_out['zVertex'] = ds_mali['zVertex'] + # <>: Creating a new netcdf file like this changes the format of xtime. Need to add + # encoding update to maintain original formatting outFileWithMeshVars = self.options.outputFile[0:-3] + '.withMeshVars.nc' ds_out.to_netcdf(outFileWithMeshVars) ds_out.close() From 0fd822482aab407ccc8a3d25219a641c6fa4b1c6 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Thu, 4 Dec 2025 13:04:49 -0800 Subject: [PATCH 07/17] Add time dimension to orig3dOceanMask --- landice/mesh_tools_li/interpolate_ecco_to_mali.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index bc18da1a7..08cabe1c5 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -243,7 +243,8 @@ def merge_MITgcm_files_to_unstruct_grid(self): if 'orig3dOceanMask' in locals(): orig3dOceanMask = orig3dOceanMask[0:len(z),:] - DA_o3dm = xr.DataArray(orig3dOceanMask.astype('int32'), dims=("nISMIP6OceanLayers", "nCells")) + orig3dOceanMask = np.tile(orig3dOceanMask[np.newaxis, :, :], (nt, 1, 1)) + DA_o3dm = xr.DataArray(orig3dOceanMask.astype('int32'), dims=("Time", "nISMIP6OceanLayers", "nCells")) DA_o3dm.attrs['long_name'] = ("3D mask of original valid ocean data. Because it is 3d," " it can include the ocean domain inside ice-shelf cavities." " orig3dOceanMask is altered online to include ice-dammed inland seas") From 373979e1bf739c7c21535f8f51ec88fd7c02943e Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Wed, 10 Dec 2025 09:10:40 -0800 Subject: [PATCH 08/17] Ensure Time in unlimited dimension --- landice/mesh_tools_li/interpolate_ecco_to_mali.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index 08cabe1c5..c42e942eb 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -251,7 +251,7 @@ def merge_MITgcm_files_to_unstruct_grid(self): ds_unstruct['orig3dOceanMask'] = DA_o3dm self.ecco_unstruct = "ecco_combined_unstructured.nc" - ds_unstruct.to_netcdf(self.ecco_unstruct) + ds_unstruct.to_netcdf(self.ecco_unstruct, unlimited_dims=['Time']) ds_unstruct.close() print("ECCO restructuring complete") @@ -336,7 +336,7 @@ def remap_ecco_to_mali(self): # <>: Creating a new netcdf file like this changes the format of xtime. Need to add # encoding update to maintain original formatting outFileWithMeshVars = self.options.outputFile[0:-3] + '.withMeshVars.nc' - ds_out.to_netcdf(outFileWithMeshVars) + ds_out.to_netcdf(outFileWithMeshVars, unlimited_dims=['Time']) ds_out.close() ds_mali.close() subprocess.run(["ncatted", "-a", "_FillValue,,d,,", outFileWithMeshVars]) From 1708452ad237db3d8481454e5ac219f4fad8d6c6 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Thu, 11 Dec 2025 08:46:35 -0800 Subject: [PATCH 09/17] Invalidate bad ocean cells after remapping Conservative remapping averages defined and undefined ocean cells. Use orig3dOceanMask to identify these cells and treat as undefined --- .../mesh_tools_li/interpolate_ecco_to_mali.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index c42e942eb..bb08f0177 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -244,7 +244,7 @@ def merge_MITgcm_files_to_unstruct_grid(self): if 'orig3dOceanMask' in locals(): orig3dOceanMask = orig3dOceanMask[0:len(z),:] orig3dOceanMask = np.tile(orig3dOceanMask[np.newaxis, :, :], (nt, 1, 1)) - DA_o3dm = xr.DataArray(orig3dOceanMask.astype('int32'), dims=("Time", "nISMIP6OceanLayers", "nCells")) + DA_o3dm = xr.DataArray(orig3dOceanMask.astype('float64'), dims=("Time", "nISMIP6OceanLayers", "nCells")) DA_o3dm.attrs['long_name'] = ("3D mask of original valid ocean data. Because it is 3d," " it can include the ocean domain inside ice-shelf cavities." " orig3dOceanMask is altered online to include ice-dammed inland seas") @@ -293,6 +293,22 @@ def remap_ecco_to_mali(self): print(f"ncremap completed: {tm} seconds") print("Remapping completed") + # During conservative remapping, some valid ocean cells are averaged with invalid ocean cells. Use + # float version of orig3dOceanMask to identify these cells and remove. Resave orig3dOceanMask as int32 + # so MALI recognizes it + ds_out = xr.open_dataset(self.options.outputFile, decode_times=False, decode_cf=False, mode='r+') + o3dm = ds_out['orig3dOceanMask'] + temp = ds_out['oceanTemperature'] + sal = ds_out['oceanSalinity'] + + mask = o3dm < 1.0 + ds_out['orig3dOceanMask'] = o3dm.where(~mask, 0).astype('int32') + ds_out['oceanTemperature'] = temp.where(~mask, 1e36).astype('float64') + ds_out['oceanSalinity'] = sal.where(~mask, 1e36).astype('float64') + + ds_out.to_netcdf('test_remapped_file.nc', 'w') + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", 'test_remapped_file.nc']) + # Transfer MALI mesh variables if needed if self.options.includeMeshVars: ds_mali = xr.open_dataset(self.options.maliMesh, decode_times=False, decode_cf=False) From bc38f7cc804a3211413acba89649f61207929bf8 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Thu, 11 Dec 2025 09:19:04 -0800 Subject: [PATCH 10/17] Synthesize netcdf file names Creates a consistent naming scheme so that only one output file is created and the remaining temporary files are removed --- .../mesh_tools_li/interpolate_ecco_to_mali.py | 33 +++++++++++++------ 1 file changed, 23 insertions(+), 10 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index bb08f0177..aa48af368 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -250,7 +250,7 @@ def merge_MITgcm_files_to_unstruct_grid(self): " orig3dOceanMask is altered online to include ice-dammed inland seas") ds_unstruct['orig3dOceanMask'] = DA_o3dm - self.ecco_unstruct = "ecco_combined_unstructured.nc" + self.ecco_unstruct = "ecco_combined_unstructured.tmp.nc" ds_unstruct.to_netcdf(self.ecco_unstruct, unlimited_dims=['Time']) ds_unstruct.close() print("ECCO restructuring complete") @@ -281,13 +281,14 @@ def remap_ecco_to_mali(self): print("ncremap ...") st = time.time() # remap the input data + remappedOutFile = self.options.outputFile[0:-3] + '.tmp.remapped.nc' args_remap = ["ncremap", "-i", self.ecco_unstruct, - "-o", self.options.outputFile, + "-o", remappedOutFile, "-m", self.mapping_file_name] check_call(args_remap) - subprocess.run(["ncpdq", "-O", "-a", "Time,nCells,nISMIP6OceanLayers", self.options.outputFile, self.options.outputFile]) + subprocess.run(["ncpdq", "-O", "-a", "Time,nCells,nISMIP6OceanLayers", remappedOutFile, remappedOutFile]) nd = time.time() tm = nd - st print(f"ncremap completed: {tm} seconds") @@ -296,7 +297,7 @@ def remap_ecco_to_mali(self): # During conservative remapping, some valid ocean cells are averaged with invalid ocean cells. Use # float version of orig3dOceanMask to identify these cells and remove. Resave orig3dOceanMask as int32 # so MALI recognizes it - ds_out = xr.open_dataset(self.options.outputFile, decode_times=False, decode_cf=False, mode='r+') + ds_out = xr.open_dataset(remappedOutFile, decode_times=False, decode_cf=False, mode='r+') o3dm = ds_out['orig3dOceanMask'] temp = ds_out['oceanTemperature'] sal = ds_out['oceanSalinity'] @@ -306,8 +307,16 @@ def remap_ecco_to_mali(self): ds_out['oceanTemperature'] = temp.where(~mask, 1e36).astype('float64') ds_out['oceanSalinity'] = sal.where(~mask, 1e36).astype('float64') - ds_out.to_netcdf('test_remapped_file.nc', 'w') - subprocess.run(["ncatted", "-a", "_FillValue,,d,,", 'test_remapped_file.nc']) + # Save final output file if not adding mesh variables + if self.options.includeMeshVars: + altRemappedOutFile = self.options.outputFile[0:-3] + '.tmp.altRemapped.nc' + else : + altRemappedOutFile = self.options.outputFile + + # <>: 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(altRemappedOutFile, 'w') + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", altRemappedOutFile]) # Transfer MALI mesh variables if needed if self.options.includeMeshVars: @@ -351,12 +360,16 @@ def remap_ecco_to_mali(self): ds_out['zVertex'] = ds_mali['zVertex'] # <>: Creating a new netcdf file like this changes the format of xtime. Need to add # encoding update to maintain original formatting - outFileWithMeshVars = self.options.outputFile[0:-3] + '.withMeshVars.nc' - ds_out.to_netcdf(outFileWithMeshVars, unlimited_dims=['Time']) + ds_out.to_netcdf(self.options.outputFile, unlimited_dims=['Time']) ds_out.close() ds_mali.close() - subprocess.run(["ncatted", "-a", "_FillValue,,d,,", outFileWithMeshVars]) - + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", self.options.outputFile]) + + # Remove temporary output files + for filename in os.listdir("."): + if "tmp" in filename and os.path.isfile(filename): + os.remove(filename) + def create_ECCO_scrip_file(self): print("Creating ECCO scrip file") From 1cb7841250fc06da3a6b382b28f24a81eff4d13a Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Thu, 11 Dec 2025 11:20:03 -0800 Subject: [PATCH 11/17] Account for rounding error in valid cell masking --- landice/mesh_tools_li/interpolate_ecco_to_mali.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index aa48af368..d23b7e4d6 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -302,10 +302,10 @@ def remap_ecco_to_mali(self): temp = ds_out['oceanTemperature'] sal = ds_out['oceanSalinity'] - mask = o3dm < 1.0 - ds_out['orig3dOceanMask'] = o3dm.where(~mask, 0).astype('int32') - ds_out['oceanTemperature'] = temp.where(~mask, 1e36).astype('float64') - ds_out['oceanSalinity'] = sal.where(~mask, 1e36).astype('float64') + mask = o3dm > 0.9999 #Account for rounding error. Good cells are not exactly 1 + ds_out['orig3dOceanMask'] = o3dm.where(mask, 0).astype('int32') + ds_out['oceanTemperature'] = temp.where(mask, 1e36).astype('float64') + ds_out['oceanSalinity'] = sal.where(mask, 1e36).astype('float64') # Save final output file if not adding mesh variables if self.options.includeMeshVars: From e4605a34124cb02979d2daf82596d819fee34655 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Mon, 15 Dec 2025 08:37:00 -0800 Subject: [PATCH 12/17] Define icebergFjordMask Extrapolates iceCellArea to fill ocean cells, and defines icebergFjordMask wherever iceCellArea is above a threshold value --- .../mesh_tools_li/interpolate_ecco_to_mali.py | 112 ++++++++++++++---- 1 file changed, 89 insertions(+), 23 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index d23b7e4d6..45d273d2f 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -15,6 +15,11 @@ import cftime import glob +""" +<< TO DO >>: +1) Fix xtime after in final output files +2) Mask out Jakobshaven/Ilulissat for icebergFjordMask +""" class eccoToMaliInterp: def __init__(self): @@ -30,6 +35,9 @@ def __init__(self): parser.add_argument("--method", dest="method", type=str, help="Remapping method, either 'bilinear' or 'neareststod'", default='conserve') parser.add_argument("--outFile", dest="outputFile", help="Desired name of output file", metavar="FILENAME", default="ecco_to_mali_remapped.nc") parser.add_argument("--meshVars", dest="includeMeshVars", help="whether to include mesh variables in resulting MALI forcing file", action='store_true') + parser.add_argument("--iceAreaThresh", dest="iceAreaThresh", type=float, help="Threshold sea ice fraction used to define icebergFjordMask.", default=0.5) + parser.add_argument("--extrapIters", dest="extrapIters", type=str, help="Number of sea ice extrapolation iterations when regridding", default='50') + parser.add_argument("--keepTmpFiles", dest="keepTmpFiles", help="Whether to keep temporary netcdf files created throughout script. Useful for debugging", action="store_true") args = parser.parse_args() self.options = args @@ -239,15 +247,17 @@ def merge_MITgcm_files_to_unstruct_grid(self): DA_sia = xr.DataArray(sia.astype('float64'), dims=("Time", "nCells")) DA_sia.attrs['long_name'] = 'Sea ice fractional ice-covered area [0 to 1]' DA_sia.attrs['units'] = 'm2/m2' - ds_unstruct['iceAreaCell'] = DA_sia + # create separate netcdf file for SIA so it can be extrapolated during remapping + ds_sia = xr.Dataset() + ds_sia['iceAreaCell'] = DA_sia + self.siaFile = 'sia.tmp.nc' + ds_sia.to_netcdf(self.siaFile) + ds_sia.close() if 'orig3dOceanMask' in locals(): orig3dOceanMask = orig3dOceanMask[0:len(z),:] orig3dOceanMask = np.tile(orig3dOceanMask[np.newaxis, :, :], (nt, 1, 1)) DA_o3dm = xr.DataArray(orig3dOceanMask.astype('float64'), dims=("Time", "nISMIP6OceanLayers", "nCells")) - DA_o3dm.attrs['long_name'] = ("3D mask of original valid ocean data. Because it is 3d," - " it can include the ocean domain inside ice-shelf cavities." - " orig3dOceanMask is altered online to include ice-dammed inland seas") ds_unstruct['orig3dOceanMask'] = DA_o3dm self.ecco_unstruct = "ecco_combined_unstructured.tmp.nc" @@ -261,7 +271,7 @@ def remap_ecco_to_mali(self): scrip_from_mpas(self.options.maliMesh, mali_scripfile) #create mapping file - print("Creating Mapping File") + print("Creating Primary Mapping File") args = ['srun', '-n', self.options.ntasks, 'ESMF_RegridWeightGen', '--source', self.ecco_scripfile, @@ -292,7 +302,64 @@ def remap_ecco_to_mali(self): nd = time.time() tm = nd - st print(f"ncremap completed: {tm} seconds") - print("Remapping completed") + print("Primary Remapping completed") + + # repeat remapping for sia, including extrapolation of sia to fill domain + 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) + + #remap + print("Start SIA remapping ... ") + + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", self.siaFile]) + + print("ncremap of SIA...") + st = time.time() + # remap the input data + siaRemappedOutFile = self.options.outputFile[0:-3] + '.tmp.sia.remapped.nc' + args_remap = ["ncremap", + "-i", self.siaFile, + "-o", siaRemappedOutFile, + "-m", sia_mapping_file] + check_call(args_remap) + + nd = time.time() + tm = nd - st + print(f"ncremap completed: {tm} seconds") + print("SIA Remapping completed") + + print("Applying final edits to remapped file ...") + # define icebergFjordMask from SIA + ds_sia = xr.open_dataset(siaRemappedOutFile) + iac = ds_sia['iceAreaCell'].values + icebergFjordMask = np.zeros(iac.shape) + # Find cells where sea ice fraction exceeds threshold + mask = iac > self.options.iceAreaThresh + icebergFjordMask[mask] = 1 + ifm = xr.DataArray(icebergFjordMask.astype('int32'), dims=("Time", "nCells")) + sia = xr.DataArray(iac.astype('float32'), dims=("Time", "nCells")) # During conservative remapping, some valid ocean cells are averaged with invalid ocean cells. Use # float version of orig3dOceanMask to identify these cells and remove. Resave orig3dOceanMask as int32 @@ -303,25 +370,23 @@ def remap_ecco_to_mali(self): sal = ds_out['oceanSalinity'] mask = o3dm > 0.9999 #Account for rounding error. Good cells are not exactly 1 - ds_out['orig3dOceanMask'] = o3dm.where(mask, 0).astype('int32') + ds_out['orig3dOceanMask'] = xr.DataArray(mask.astype('int32'), dims=("Time", "nCells", "nISMIP6OceanLayers")) ds_out['oceanTemperature'] = temp.where(mask, 1e36).astype('float64') ds_out['oceanSalinity'] = sal.where(mask, 1e36).astype('float64') + # Add sia and icebergBergMask into main file + ds_out['iceAreaCell'] = sia + ds_out['icebergFjordMask'] = ifm # Save final output file if not adding mesh variables - if self.options.includeMeshVars: - altRemappedOutFile = self.options.outputFile[0:-3] + '.tmp.altRemapped.nc' - else : - altRemappedOutFile = self.options.outputFile - - # <>: 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(altRemappedOutFile, 'w') - subprocess.run(["ncatted", "-a", "_FillValue,,d,,", altRemappedOutFile]) - - # Transfer MALI mesh variables if needed - if self.options.includeMeshVars: + if not self.options.includeMeshVars: + # <>: 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]) + + else: + # Transfer MALI mesh variables if needed ds_mali = xr.open_dataset(self.options.maliMesh, decode_times=False, decode_cf=False) - ds_out = xr.open_dataset(self.options.outputFile, decode_times=False, decode_cf=False, mode='r+') ds_out['angleEdge'] = ds_mali['angleEdge'] ds_out['areaCell'] = ds_mali['areaCell'] ds_out['areaTriangle'] = ds_mali['areaTriangle'] @@ -366,9 +431,10 @@ def remap_ecco_to_mali(self): subprocess.run(["ncatted", "-a", "_FillValue,,d,,", self.options.outputFile]) # Remove temporary output files - for filename in os.listdir("."): - if "tmp" in filename and os.path.isfile(filename): - os.remove(filename) + if not self.options.keepTmpFiles : + for filename in os.listdir("."): + if "tmp" in filename and os.path.isfile(filename): + os.remove(filename) def create_ECCO_scrip_file(self): From daac57eea134b150c2ca01d079479c855578d22d Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Mon, 15 Dec 2025 12:11:11 -0800 Subject: [PATCH 13/17] Add geojson to mask icebergFjordMask Adds an option to use a geojson file to define a region where icebergFjordMask is permanently 1 --- .../mesh_tools_li/interpolate_ecco_to_mali.py | 26 ++++++++++++++++--- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index 45d273d2f..0fdd667af 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -14,11 +14,14 @@ import time import cftime import glob +import json +from shapely.geometry import shape +from shapely.prepared import prep +from shapely.geometry import Point """ << TO DO >>: 1) Fix xtime after in final output files -2) Mask out Jakobshaven/Ilulissat for icebergFjordMask """ class eccoToMaliInterp: @@ -30,6 +33,7 @@ def __init__(self): parser.add_argument("--eccoDir", dest="eccoDir", required=True, help="Directory where ECCO files are stored. Should be only netcdf files in that directory", metavar="FILENAME") parser.add_argument("--maliMesh", dest="maliMesh", required=True, help="MALI mesh to be interpolated onto", metavar="FILENAME") parser.add_argument("--eccoVars", type=str, nargs='+', dest="eccoVars", required=True, help="ECCO variables to interpolate, current options are any combination of 'THETA', 'SALT', 'SIarea'") + parser.add_argument("--geojson", dest="geojson", help="Optional geojson used to permanently define icebergFjordMask within a region") parser.add_argument("--maxDepth", type=int, dest="maxDepth", default=1500, help="Maximum depth in meters to include in the interpolation") parser.add_argument("--ntasks", dest="ntasks", type=str, help="Number of processors to use with ESMF_RegridWeightGen", default='128') parser.add_argument("--method", dest="method", type=str, help="Remapping method, either 'bilinear' or 'neareststod'", default='conserve') @@ -352,19 +356,33 @@ def remap_ecco_to_mali(self): print("Applying final edits to remapped file ...") # define icebergFjordMask from SIA + ds_out = xr.open_dataset(remappedOutFile, decode_times=False, decode_cf=False, mode='r+') ds_sia = xr.open_dataset(siaRemappedOutFile) iac = ds_sia['iceAreaCell'].values icebergFjordMask = np.zeros(iac.shape) # Find cells where sea ice fraction exceeds threshold - mask = iac > self.options.iceAreaThresh - icebergFjordMask[mask] = 1 + seaice_mask = iac > self.options.iceAreaThresh + icebergFjordMask[seaice_mask] = 1 + + # Use geojson to add any additional cells to icebergFjordMask + if hasattr(self.options, "geojson"): + with open(self.options.geojson) as f: + geo = json.load(f) + polygon = prep(shape(geo["features"][0]["geometry"])) + + lonCell = ds_out['lon'][:].values.ravel() - 360 # Convert to geojson longitude convention + latCell = ds_out['lat'][:].values.ravel() + + points = [Point(xy) for xy in zip(lonCell, latCell)] + geo_mask = np.array([polygon.covers(p) for p in points]).astype('bool') + icebergFjordMask[:,geo_mask] = 1 + ifm = xr.DataArray(icebergFjordMask.astype('int32'), dims=("Time", "nCells")) sia = xr.DataArray(iac.astype('float32'), dims=("Time", "nCells")) # During conservative remapping, some valid ocean cells are averaged with invalid ocean cells. Use # float version of orig3dOceanMask to identify these cells and remove. Resave orig3dOceanMask as int32 # so MALI recognizes it - ds_out = xr.open_dataset(remappedOutFile, decode_times=False, decode_cf=False, mode='r+') o3dm = ds_out['orig3dOceanMask'] temp = ds_out['oceanTemperature'] sal = ds_out['oceanSalinity'] From 88454d3934b509d74ad443907547dd2aa904669d Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Mon, 15 Dec 2025 16:38:49 -0800 Subject: [PATCH 14/17] Fix xtime in final output file Moves saving of xtime time down to where final output file is created in order to preserve its format --- .../mesh_tools_li/interpolate_ecco_to_mali.py | 20 +++++++++---------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index 0fdd667af..557c3f374 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -19,10 +19,6 @@ from shapely.prepared import prep from shapely.geometry import Point -""" -<< TO DO >>: -1) Fix xtime after in final output files -""" class eccoToMaliInterp: def __init__(self): @@ -197,14 +193,10 @@ def merge_MITgcm_files_to_unstruct_grid(self): ds_unstruct = xr.Dataset() ds_unstruct.expand_dims(["nCells", "Time"]) - # Translate ECCO time to MALI xtime format as save to dataset + # Save ECCO time variable as datetime string. Will be converted into xtime format when saving final output file time = ds_ecco['tim'].values - xtime_str = [pd.to_datetime(dt).strftime("%Y-%m-%d_%H:%M:%S").ljust(64) for dt in time] - xtime = np.array(xtime_str, dtype = np.dtype(('S',64))) - DA_xtime = xr.DataArray(xtime, dims=["Time"]) - DA_xtime.encoding.update({"char_dim_name": "StrLen"}) - ds_unstruct['xtime'] = DA_xtime - + self.xtime_str = [pd.to_datetime(dt).strftime("%Y-%m-%d_%H:%M:%S").ljust(64) for dt in time] + # save unstructured variables with MALI/MPAS-O names if 'z' in locals(): ds_unstruct.expand_dims("nISMIP6OceanLayers") @@ -394,6 +386,12 @@ def remap_ecco_to_mali(self): # Add sia and icebergBergMask into main file ds_out['iceAreaCell'] = sia ds_out['icebergFjordMask'] = ifm + + # save xtime correctly + xtime = np.array(self.xtime_str, dtype = np.dtype(('S',64))) + DA_xtime = xr.DataArray(xtime, dims=["Time"]) + DA_xtime.encoding.update({"char_dim_name": "StrLen"}) + ds_out['xtime'] = DA_xtime # Save final output file if not adding mesh variables if not self.options.includeMeshVars: From 049c309ef3399a42f97ce55dd2887460c32df927 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Tue, 23 Dec 2025 15:10:45 -0800 Subject: [PATCH 15/17] netcdf3 64-bit offset Creates a copy of the output file with netcdf3 64-bit offset formatting to be consistent with MALI requirements --- landice/mesh_tools_li/interpolate_ecco_to_mali.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index 557c3f374..3dc10ad3d 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -452,6 +452,11 @@ def remap_ecco_to_mali(self): if "tmp" in filename and os.path.isfile(filename): os.remove(filename) + # 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' + subprocess.run(["ncks", "-6", self.options.outputFile, nc64offset]) + def create_ECCO_scrip_file(self): print("Creating ECCO scrip file") From fb4f4024d93637d8f9a14f4ea8257e006159f1ce Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Tue, 23 Dec 2025 15:25:32 -0800 Subject: [PATCH 16/17] Subtract 1 month from ECCO time 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. --- landice/mesh_tools_li/interpolate_ecco_to_mali.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index 3dc10ad3d..b2235f8d0 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -194,9 +194,11 @@ def merge_MITgcm_files_to_unstruct_grid(self): ds_unstruct.expand_dims(["nCells", "Time"]) # Save ECCO time variable as datetime string. Will be converted into xtime format when saving final output file + # 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).strftime("%Y-%m-%d_%H:%M:%S").ljust(64) for dt in time] - + self.xtime_str = [(pd.to_datetime(dt) - pd.to_datetime(dt)).strftime("%Y-%m-%d_%H:%M:%S").ljust(64) for dt in time] + # save unstructured variables with MALI/MPAS-O names if 'z' in locals(): ds_unstruct.expand_dims("nISMIP6OceanLayers") From 970414149675e2ca0e77f923a471042815c07829 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Tue, 23 Dec 2025 15:39:10 -0800 Subject: [PATCH 17/17] No orig3dOceanMask where TS is invalid 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 --- landice/mesh_tools_li/interpolate_ecco_to_mali.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/landice/mesh_tools_li/interpolate_ecco_to_mali.py b/landice/mesh_tools_li/interpolate_ecco_to_mali.py index b2235f8d0..c1c62d5ed 100755 --- a/landice/mesh_tools_li/interpolate_ecco_to_mali.py +++ b/landice/mesh_tools_li/interpolate_ecco_to_mali.py @@ -381,7 +381,9 @@ def remap_ecco_to_mali(self): temp = ds_out['oceanTemperature'] sal = ds_out['oceanSalinity'] - mask = o3dm > 0.9999 #Account for rounding error. Good cells are not exactly 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) ds_out['orig3dOceanMask'] = xr.DataArray(mask.astype('int32'), dims=("Time", "nCells", "nISMIP6OceanLayers")) ds_out['oceanTemperature'] = temp.where(mask, 1e36).astype('float64') ds_out['oceanSalinity'] = sal.where(mask, 1e36).astype('float64')