Runonff mapping fixes

[alperaltuntas]

This is the last of a series of PRs to implement a robust runoff grid selection and runoff to ocean mapping generation features. The other PRs were:

ESMCI/visualCaseGen#26
https://github.com/NCAR/mom6_bathy/pull/46

Changes include:

• misc minor refactorings.
• ability to specify runoff grid
• a robust implementation of runoff to ocean mapping.

[alperaltuntas]

I've attempted to resolve the conflicts, which I think is due to the merge of #159, but I am unable to run the CrocoGallery examples to test my conflict resolution. @manishvenu Has the CrocoGallery been updated for the latest changes? I am interested in running three_boundary_from_t232.ipynb in particular.

[manishvenu]

I've attempted to resolve the conflicts, which I think is due to the merge of #159, but I am unable to run the CrocoGallery examples to test my conflict resolution. @manishvenu Has the CrocoGallery been updated for the latest changes? I am interested in running three_boundary_from_t232.ipynb in particular.

The CrocoGallery is not updated. However that specific notebook shouldn't have issues. What is your error?

[alperaltuntas]

When I run:

case.configure_forcings(
    date_range = ["2020-01-01 00:00:00", "2020-02-01 00:00:00"],
    boundaries = ["south", "north","east"], 
    rmax = 250.0,
    fold = 250.0
)

I get:

INFO:GLORYS:This data access method retuns a script at path /glade/derecho/scratch/altuntas/tmp/tmpqhs46fs5/get_glorys_data.sh to run to get access data 
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Input In [14], in <cell line: 0>()
----> 1 case.configure_forcings(
      2     date_range = ["2020-01-01 00:00:00", "2020-02-01 00:00:00"],
      3     boundaries = ["south", "north","east"], 
      4     rmax = 250.0,
      5     fold = 250.0
      6 )

File /glade/work/altuntas/croc/CrocoDash/CrocoDash/case.py:436, in Case.configure_forcings(self, date_range, boundaries, tidal_constituents, tpxo_elevation_filepath, tpxo_velocity_filepath, product_name, function_name, product_info, too_much_data, rmax, fold, chl_processed_filepath, global_river_nutrients_filepath, marbl_ic_filepath)
    431 self.forcing_product_name = product_name.lower()
    432 if (
    433     ProductRegistry.product_exists(product_name)
    434     and ProductRegistry.product_is_of_type(product_name,ForcingProduct)
    435 ):
--> 436     self.configure_initial_and_boundary_conditions(
    437         date_range=date_range,
    438         boundaries=boundaries,
    439         product_name=product_name,
    440         function_name=function_name,
    441         too_much_data=too_much_data,
    442     )
    443 else:
    444     raise ValueError("Product / Data Path is not supported quite yet")

File /glade/work/altuntas/croc/CrocoDash/CrocoDash/case.py:610, in Case.configure_initial_and_boundary_conditions(self, date_range, boundaries, product_name, function_name, too_much_data)
    607 spec.loader.exec_module(self.driver)
    609 if not self._too_much_data:
--> 610     self.driver.main(
    611         regrid_dataset_piecewise=False, merge_piecewise_dataset=False
    612     )
    613 else:
    614     print(
    615         f"Extract Forcings workflow was called, please go to the extract forcings path: {self.extract_forcings_path} and run the driver script there."
    616     )

File /glade/derecho/scratch/altuntas/crocodile_2025/croc_input/c.e30.C_JRA.nwa232.004/glorys/extract_forcings/driver.py:41, in main(get_dataset_piecewise, regrid_dataset_piecewise, merge_piecewise_dataset)
     39 # Call get_dataset_piecewise
     40 if get_dataset_piecewise:
---> 41     gdp.get_dataset_piecewise(
     42         product_name=config["forcing"]["product_name"],
     43         function_name=config["forcing"]["function_name"],
     44         product_information=config["forcing"]["information"],
     45         date_format=config["dates"]["format"],
     46         start_date=config["dates"]["start"],
     47         end_date=config["dates"]["end"],
     48         hgrid_path=config["paths"]["hgrid_path"],
     49         step_days=int(config["general"]["step"]),
     50         output_dir=config["paths"]["raw_dataset_path"],
     51         boundary_number_conversion=config["general"]["boundary_number_conversion"],
     52         run_initial_condition=config["general"]["run_initial_condition"],
     53         run_boundary_conditions=config["general"]["run_boundary_conditions"],
     54         preview=config["general"]["preview"],
     55     )
     57 # Call regrid_dataset_piecewise
     58 if regrid_dataset_piecewise:

File /glade/derecho/scratch/altuntas/crocodile_2025/croc_input/c.e30.C_JRA.nwa232.004/glorys/extract_forcings/code/get_dataset_piecewise.py:79, in get_dataset_piecewise(product_name, function_name, product_information, date_format, start_date, end_date, hgrid_path, step_days, output_dir, boundary_number_conversion, run_initial_condition, run_boundary_conditions, preview)
     77 # Get lat,lon information for each boundary
     78 hgrid = xr.open_dataset(hgrid_path)
---> 79 boundary_info = Grid.get_bounding_boxes_of_rectangular_grid(hgrid)
     81 # Set up date range, pd.date_range is exclusive of the end_date
     82 dates = (
     83     pd.date_range(start=start_date, end=end_date, freq=f"{step_days}D")
     84     .to_pydatetime()
     85     .tolist()
     86 )

File /glade/work/altuntas/croc/CrocoDash/CrocoDash/visualCaseGen/external/mom6_bathy/mom6_bathy/grid.py:385, in Grid.get_bounding_boxes_of_rectangular_grid(cls, hgrid)
    383 else:
    384     grid_check = Grid.from_supergrid_ds(hgrid)
--> 385     assert grid_check.is_rectangular()
    386     assert not Grid.is_cyclic_x(hgrid)
    388 init_result = {
    389     "lon_min": float(hgrid.x.min()),
    390     "lon_max": float(hgrid.x.max()),
    391     "lat_min": float(hgrid.y.min()),
    392     "lat_max": float(hgrid.y.max()),
    393 }

AssertionError: 

[manishvenu]

For some reason this grid is failing the is_rectangular check, which is super weird, I'll try running this.

[manishvenu]

Oh actually, the is_rectangular check is not required for this, we just need a check for an actual rectangle, and change the is_rectangular to is_rectilinear.

Do you mind removing the lines in your mom6_bathy? I can open a PR tday to make these changes.

[alperaltuntas]

Oh actually, the is_rectangular check is not required for this, we just need a check for an actual rectangle, and change the is_rectangular to is_rectilinear.

Do you mind removing the lines in your mom6_bathy? I can open a PR tday to make these changes.

Removing the line fixed that issue. But I encountered another error in the process_forcings() step:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [17], in <cell line: 0>()
----> 1 case.process_forcings()

File /glade/work/altuntas/croc/CrocoDash/CrocoDash/case.py:734, in Case.process_forcings(self, process_initial_condition, process_tides, process_velocity_tracers, process_bgc, process_chl, process_river_nutrients)
    729 if not self._configure_forcings_called:
    730     raise RuntimeError(
    731         "configure_forcings() must be called before process_forcings()."
    732     )
--> 734 self.process_initial_and_boundary_conditions(
    735     process_initial_condition, process_velocity_tracers
    736 )
    737 if self.configured_bgc and process_bgc:
    738     self.process_bgc_iron_forcing()

File /glade/work/altuntas/croc/CrocoDash/CrocoDash/case.py:1037, in Case.process_initial_and_boundary_conditions(self, process_initial_condition, process_velocity_tracers)
   1032     print(
   1033         f"Please make sure to execute large_data_workflow as described in {self.extract_forcings_path}"
   1034     )
   1036 if process_initial_condition or process_velocity_tracers:
-> 1037     self.driver.main(
   1038         get_dataset_piecewise=False,
   1039         regrid_dataset_piecewise=True,
   1040         merge_piecewise_dataset=True,
   1041     )

File /glade/derecho/scratch/altuntas/crocodile_2025/croc_input/c.e30.C_JRA.nwa232.004/glorys/extract_forcings/driver.py:59, in main(get_dataset_piecewise, regrid_dataset_piecewise, merge_piecewise_dataset)
     57 # Call regrid_dataset_piecewise
     58 if regrid_dataset_piecewise:
---> 59     rdp.regrid_dataset_piecewise(
     60         config["paths"]["raw_dataset_path"],
     61         config["file_regex"]["raw_dataset_pattern"],
     62         config["dates"]["format"],
     63         config["dates"]["start"],
     64         config["dates"]["end"],
     65         config["paths"]["hgrid_path"],
     66         config["paths"]["bathymetry_path"],
     67         config["forcing"]["information"],
     68         config["paths"]["regridded_dataset_path"],
     69         config["general"]["boundary_number_conversion"],
     70         config["general"]["run_initial_condition"],
     71         config["general"]["run_boundary_conditions"],
     72         config["paths"]["vgrid_path"],
     73         config["general"]["preview"],
     74     )
     76 # Call merge_dataset_piecewise
     77 if merge_piecewise_dataset:

File /glade/derecho/scratch/altuntas/crocodile_2025/croc_input/c.e30.C_JRA.nwa232.004/glorys/extract_forcings/code/regrid_dataset_piecewise.py:222, in regrid_dataset_piecewise(folder, input_dataset_regex, date_format, start_date, end_date, hgrid_path, bathymetry, dataset_varnames, output_folder, boundary_number_conversion, run_initial_condition, run_boundary_conditions, vgrid_path, preview)
    218     logger.info(
    219         f"Initial condition files already exist. They will be skipped."
    220     )
    221 else:
--> 222     expt.setup_initial_condition(
    223         file_path, dataset_varnames, arakawa_grid=None
    224     )
    225 if (expt.mom_input_dir / "init_eta_filled.nc").exists():
    226     logger.info(
    227         f"Initial condition filled files already exist. They will be skipped."
    228     )

File /glade/work/altuntas/croc/CrocoDash/CrocoDash/rm6/regional_mom6/regional_mom6.py:1091, in experiment.setup_initial_condition(self, raw_ic_path, varnames, arakawa_grid, vcoord_type, rotational_method, regridding_method)
   1085     ic_raw[varnames["tracers"]["temp"]].attrs["units"] = "degrees Celsius"
   1086 # NaNs might be here from the land mask of the model that the IC has come from.
   1087 # If they're not removed then the coastlines from this other grid will be retained!
   1088 # The land mask comes from the bathymetry file, so we don't need NaNs
   1089 # to tell MOM6 where the land is.
   1090 ic_raw_tracers = (
-> 1091     ic_raw_tracers.interpolate_na(
   1092         reprocessed_var_map["tracer_x_coord"], method="linear"
   1093     )
   1094     .ffill(reprocessed_var_map["tracer_x_coord"])
   1095     .bfill(reprocessed_var_map["tracer_x_coord"])
   1096     .ffill(reprocessed_var_map["tracer_y_coord"])
   1097     .bfill(reprocessed_var_map["tracer_y_coord"])
   1098     .ffill(reprocessed_var_map["depth_coord"])
   1099 )
   1101 ic_raw_u = (
   1102     ic_raw_u.interpolate_na(reprocessed_var_map["u_x_coord"], method="linear")
   1103     .ffill(reprocessed_var_map["u_x_coord"])
   (...)
   1107     .ffill(reprocessed_var_map["depth_coord"])
   1108 )
   1110 ic_raw_v = (
   1111     ic_raw_v.interpolate_na(reprocessed_var_map["v_x_coord"], method="linear")
   1112     .ffill(reprocessed_var_map["v_x_coord"])
   (...)
   1116     .ffill(reprocessed_var_map["depth_coord"])
   1117 )

File /glade/work/altuntas/miniconda3/envs/CrocoDash/lib/python3.11/site-packages/xarray/core/dataset.py:6588, in Dataset.interpolate_na(self, dim, method, limit, use_coordinate, max_gap, **kwargs)
   6471 """Fill in NaNs by interpolating according to different methods.
   6472 
   6473 Parameters
   (...)
   6584     D        (x) float64 5.0 3.0 1.0 -1.0 4.0
   6585 """
   6586 from xarray.core.missing import _apply_over_vars_with_dim, interp_na
-> 6588 new = _apply_over_vars_with_dim(
   6589     interp_na,
   6590     self,
   6591     dim=dim,
   6592     method=method,
   6593     limit=limit,
   6594     use_coordinate=use_coordinate,
   6595     max_gap=max_gap,
   6596     **kwargs,
   6597 )
   6598 return new

File /glade/work/altuntas/miniconda3/envs/CrocoDash/lib/python3.11/site-packages/xarray/core/missing.py:214, in _apply_over_vars_with_dim(func, self, dim, **kwargs)
    212 for name, var in self.data_vars.items():
    213     if dim in var.dims:
--> 214         ds[name] = func(var, dim=dim, **kwargs)
    215     else:
    216         ds[name] = var

File /glade/work/altuntas/miniconda3/envs/CrocoDash/lib/python3.11/site-packages/xarray/core/missing.py:357, in interp_na(self, dim, use_coordinate, method, limit, max_gap, keep_attrs, **kwargs)
    355     warnings.filterwarnings("ignore", "overflow", RuntimeWarning)
    356     warnings.filterwarnings("ignore", "invalid value", RuntimeWarning)
--> 357     arr = apply_ufunc(
    358         interpolator,
    359         self,
    360         index,
    361         input_core_dims=[[dim], [dim]],
    362         output_core_dims=[[dim]],
    363         output_dtypes=[self.dtype],
    364         dask="parallelized",
    365         vectorize=True,
    366         keep_attrs=keep_attrs,
    367     ).transpose(*self.dims)
    369 if limit is not None:
    370     arr = arr.where(valids)

File /glade/work/altuntas/miniconda3/envs/CrocoDash/lib/python3.11/site-packages/xarray/core/computation.py:1267, in apply_ufunc(func, input_core_dims, output_core_dims, exclude_dims, vectorize, join, dataset_join, dataset_fill_value, keep_attrs, kwargs, dask, output_dtypes, output_sizes, meta, dask_gufunc_kwargs, on_missing_core_dim, *args)
   1265 # feed DataArray apply_variable_ufunc through apply_dataarray_vfunc
   1266 elif any(isinstance(a, DataArray) for a in args):
-> 1267     return apply_dataarray_vfunc(
   1268         variables_vfunc,
   1269         *args,
   1270         signature=signature,
   1271         join=join,
   1272         exclude_dims=exclude_dims,
   1273         keep_attrs=keep_attrs,
   1274     )
   1275 # feed Variables directly through apply_variable_ufunc
   1276 elif any(isinstance(a, Variable) for a in args):

File /glade/work/altuntas/miniconda3/envs/CrocoDash/lib/python3.11/site-packages/xarray/core/computation.py:315, in apply_dataarray_vfunc(func, signature, join, exclude_dims, keep_attrs, *args)
    310 result_coords, result_indexes = build_output_coords_and_indexes(
    311     args, signature, exclude_dims, combine_attrs=keep_attrs
    312 )
    314 data_vars = [getattr(a, "variable", a) for a in args]
--> 315 result_var = func(*data_vars)
    317 out: tuple[DataArray, ...] | DataArray
    318 if signature.num_outputs > 1:

File /glade/work/altuntas/miniconda3/envs/CrocoDash/lib/python3.11/site-packages/xarray/core/computation.py:768, in apply_variable_ufunc(func, signature, exclude_dims, dask, output_dtypes, vectorize, keep_attrs, dask_gufunc_kwargs, *args)
    766             for axis, dim in enumerate(core_dims, start=-len(core_dims)):
    767                 if len(data.chunks[axis]) != 1:
--> 768                     raise ValueError(
    769                         f"dimension {dim} on {n}th function argument to "
    770                         "apply_ufunc with dask='parallelized' consists of "
    771                         "multiple chunks, but is also a core dimension. To "
    772                         "fix, either rechunk into a single array chunk along "
    773                         f"this dimension, i.e., ``.chunk(dict({dim}=-1))``, or "
    774                         "pass ``allow_rechunk=True`` in ``dask_gufunc_kwargs`` "
    775                         "but beware that this may significantly increase memory usage."
    776                     )
    777     dask_gufunc_kwargs["allow_rechunk"] = True
    779 output_sizes = dask_gufunc_kwargs.pop("output_sizes", {})

ValueError: dimension longitude on 0th function argument to apply_ufunc with dask='parallelized' consists of multiple chunks, but is also a core dimension. To fix, either rechunk into a single array chunk along this dimension, i.e., ``.chunk(dict(longitude=-1))``, or pass ``allow_rechunk=True`` in ``dask_gufunc_kwargs`` but beware that this may significantly increase memory usage.

[manishvenu]

Ah shoot, this is a known issue. This pr should help: CROCODILE-CESM/regional-mom6#64

I think we now use this in our specific checkout in regional mom6 and the workshop branch.

[manishvenu]

Looks great!

[manishvenu]

I feel like there was a reason we write the cice grid in init, but I can't remember why, so this probably works!

[manishvenu]

These xml changes should be in _update_forcing_variables I think (though it doesn't matter much because I'll be moving this to ForcingConfigRegistry in the new PR)

[github-actions]

📄 Preview your docs here:
👉 https://CROCODILE-CESM.github.io/CrocoDash/pr-166/index.html