From 2933a05931f094cdb87be892e6228de8d5ad2c67 Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Wed, 24 Apr 2024 15:48:34 -0700 Subject: [PATCH 01/12] Built out `Ismip6GrISForcing` testgroup. Still a work in progress. --- compass/landice/__init__.py | 2 + .../tests/ismip6_GrIS_forcing/__init__.py | 29 ++++ .../create_mapping_files.py | 140 ++++++++++++++++++ .../tests/ismip6_GrIS_forcing/experiments.yml | 35 +++++ .../tests/ismip6_GrIS_forcing/file_finders.py | 137 +++++++++++++++++ .../forcing_gen/__init__.py | 94 ++++++++++++ .../forcing_gen/forcing_gen.cfg | 11 ++ .../ismip6_GrIS_forcing/process_forcing.py | 99 +++++++++++++ 8 files changed, 547 insertions(+) create mode 100644 compass/landice/tests/ismip6_GrIS_forcing/__init__.py create mode 100644 compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py create mode 100644 compass/landice/tests/ismip6_GrIS_forcing/experiments.yml create mode 100644 compass/landice/tests/ismip6_GrIS_forcing/file_finders.py create mode 100644 compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py create mode 100644 compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/forcing_gen.cfg create mode 100644 compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py diff --git a/compass/landice/__init__.py b/compass/landice/__init__.py index 97b33a010c..2b531dfdb5 100644 --- a/compass/landice/__init__.py +++ b/compass/landice/__init__.py @@ -10,6 +10,7 @@ from compass.landice.tests.humboldt import Humboldt from compass.landice.tests.hydro_radial import HydroRadial from compass.landice.tests.ismip6_forcing import Ismip6Forcing +from compass.landice.tests.ismip6_GrIS_forcing import Ismip6GrISForcing from compass.landice.tests.ismip6_run import Ismip6Run from compass.landice.tests.isunnguata_sermia import IsunnguataSermia from compass.landice.tests.kangerlussuaq import Kangerlussuaq @@ -43,6 +44,7 @@ def __init__(self): self.add_test_group(Humboldt(mpas_core=self)) self.add_test_group(HydroRadial(mpas_core=self)) self.add_test_group(Ismip6Forcing(mpas_core=self)) + self.add_test_group(Ismip6GrISForcing(mpas_core=self)) self.add_test_group(Ismip6Run(mpas_core=self)) self.add_test_group(IsunnguataSermia(mpas_core=self)) self.add_test_group(Kangerlussuaq(mpas_core=self)) diff --git a/compass/landice/tests/ismip6_GrIS_forcing/__init__.py b/compass/landice/tests/ismip6_GrIS_forcing/__init__.py new file mode 100644 index 0000000000..52b9386d36 --- /dev/null +++ b/compass/landice/tests/ismip6_GrIS_forcing/__init__.py @@ -0,0 +1,29 @@ +import yaml + +from compass.landice.tests.ismip6_GrIS_forcing.forcing_gen import ForcingGen +from compass.testgroup import TestGroup + + +class Ismip6GrISForcing(TestGroup): + """ + A test group for processing the forcing for ISMIP6 GrIS projections + """ + def __init__(self, mpas_core): + """ + mpas_core : compass.landice.Landice + the MPAS core that this test group belongs to + """ + super().__init__(mpas_core=mpas_core, name='ismip6_GrIS_forcing') + + self.add_test_case(ForcingGen(test_group=self)) + + with open("experiments.yml", "r") as f: + experiments = yaml.safe_load(f) + + # should I filter through the experiments dictionary to make sure + # everything is valid.... + + self.experiments = experiments + + def __read_and_validate_experiment_file(self): + pass diff --git a/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py b/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py new file mode 100644 index 0000000000..e84b3cd9ad --- /dev/null +++ b/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py @@ -0,0 +1,140 @@ +import glob +import os +import shutil + +import xarray as xr +from mpas_tools.logging import check_call +from mpas_tools.scrip.from_mpas import scrip_from_mpas +from pyproj import Proj +from pyremap.descriptor import ProjectionGridDescriptor + +from compass.step import Step + + +class CreateMappingFiles(Step): + """ + """ + + def __init__(self, test_case): + """ + """ + name = "create_mapping_files" + subdir = "mapping_files" + + super().__init__(test_case=test_case, name=name, subdir=subdir, + cpus_per_task=128, min_cpus_per_task=1) + + # initalize the attributes with empty values for now, attributes will + # be propely in `setup` method so user specified config options can + # be parsed + self.mali_mesh = None + + def setup(self): + """ + """ + + # parse user specified parameters from the config + config = self.config + section = config["ISMIP6_GrIS_Forcing"] + + # read the mesh path/name from the config file + mali_mesh = section.get("MALI_mesh_fp") + + # make sure the mesh path/name is valid and accessible + if not os.path.exists(mali_mesh): + raise FileNotFoundError(f"{mali_mesh}") + + # add mesh name as an attribute and an input file to the step + self.mali_mesh = mali_mesh + self.add_input_file(filename=self.mali_mesh) + + # loop over mapping file `test_case` attributes and add the current + # steps `work_dir` to the absolute path. This will ensure other steps + # will be able to access the mapping files created here + for attr in ["mali_mesh_scrip", + "ismip6_GrIS_scrip", + "remapping_weights"]: + # get the attribute values + file_name = getattr(self.test_case, attr) + # combine filename w/ the testcases work dir + full_path = os.path.join(self.work_dir, file_name) + # update the attribute value to include the full path + setattr(self.test_case, attr, full_path) + + # Add the testcase attributes as output files; since they'll be + # generated as part of this step + self.add_output_file(filename=self.test_case.mali_mesh_scrip) + self.add_output_file(filename=self.test_case.ismip6_GrIS_scrip) + self.add_output_file(filename=self.test_case.remapping_weights) + + def run(self): + """ + """ + + # unpack the filepaths stored as attributes of the test_case + mali_mesh_scrip = self.test_case.mali_mesh_scrip + ismip6_GrIS_scrip = self.test_case.ismip6_GrIS_scrip + remapping_weights = self.test_case.remapping_weights + + # create scrip file describing MALI mesh that forcing variables + # will be interpolated onto + scrip_from_mpas(self.mali_mesh, mali_mesh_scrip) + + # get a list of the projection experiments requesed. + # (i.e. experiments names that use ISMIP6 forcing) + proj_exprs = list(filter(lambda x: "Exp" in x, + self.test_case.experiments.keys())) + + if proj_exprs: + # b/c all GrIS forcing files are on the same grid, it doesn't + # matter what expr we use; so just use the first suitable candiate + expr = proj_exprs[0] + + expr_params = self.test_case.experiments[expr] + + GCM = expr_params["GCM"] + scenario = expr_params["Scenario"] + variable = "thermal_forcing" + # get the filepath for any given forcing file (since all files are on + # the same grid), in order to generate a scrip file for the grid + forcing_fp = self.test_case.findForcingFiles(GCM, scenario, variable) + + # WGS 84 / NSIDC Sea Ice Polar Stereographic North Projection + epsg_3413 = Proj("EPSG:3413") + + # open up forcing grid; get x and y coordinates as numpy arrays + with xr.open_dataset(forcing_fp) as ds: + x = ds.x.values + y = ds.y.values + + # create pyremap `MeshDescriptor` obj for the ISMIP6 GrIS forcing grid + forcingDescriptor = ProjectionGridDescriptor.create( + projection=epsg_3413, x=x, y=y, meshName=forcing_fp) + + # write scrip file describing ISMIP6 GrIS forcing grid + forcingDescriptor.to_scrip(ismip6_GrIS_scrip) + + # generate mapping weights between ismip6 forcing data and the MALI + # mesh, using `ESMF_RegridWeightGen` + args = ['srun', '-n', '128', 'ESMF_RegridWeightGen', + '-s', ismip6_GrIS_scrip, + '-d', mali_mesh_scrip, + '-w', remapping_weights, + '--method', 'conserve', + "--netcdf4", + # "--no_log", + "--dst_regional", "--src_regional", '--ignore_unmapped'] + + check_call(args, logger=self.logger) + + # `ESMF_RegridWeightGen` will generate a log file for each processor, + # which really clutters up the directory. So, if log files are created + # move them to their own `logs` directory + log_files = glob.glob("PET*.RegridWeightGen.Log") + + # check if glob returns an empty list + if log_files: + # make the log directory + os.mkdir("logs/") + # copy the log files to the log directory; using list comprehension + _ = [shutil.move(f, "logs/") for f in log_files] diff --git a/compass/landice/tests/ismip6_GrIS_forcing/experiments.yml b/compass/landice/tests/ismip6_GrIS_forcing/experiments.yml new file mode 100644 index 0000000000..e5b959fbac --- /dev/null +++ b/compass/landice/tests/ismip6_GrIS_forcing/experiments.yml @@ -0,0 +1,35 @@ +ctrl: + Scenario: ctrl + GCM: RACMO + start: 1960 + end: 1989 + +hist: + Scenario: ctrl + GCM: RACMO + start: 1989 + end: 2015 + +Exp05: + Scenario: rcp8.5 + GCM: MIROC5 + start: 2015 + end: 2100 + +Exp06: + Scenario: rcp8.5 + GCM: NorESM1-M + start: 2015 + end: 2100 + +Exp07: + Scenario: rcp2.6 + GCM: MIROC5 + start: 2015 + end: 2100 + +Exp08: + Scenario: rcp8.5 + GCM: HadGEM2-ES + start: 2015 + end: 2100 diff --git a/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py b/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py new file mode 100644 index 0000000000..7802874711 --- /dev/null +++ b/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py @@ -0,0 +1,137 @@ +import glob +import os + +import xarray as xr + +# create mapping dictionary of ISMIP6 variables to MALI variable names +{"thermal_forcing": "ismip6_2dThermalForcing", + "basin_runoff": "basin_runoff"} + + +def check_year(fn, start=2015, end=2100): + + # strip the file extension from filename and + # convert the year at end of str to an int + fn_year = int(os.path.splitext(fn)[0].split("-")[-1]) + + # check if year is within range + return start <= fn_year <= end + + +class oceanFileFinder: + + def __init__(self, archive_fp, version="v4"): + + self.version = version + + # file strucutre within Ocean_Forcing directory for navigating + file_struct = f"Ocean_Forcing/Melt_Implementation/{version}" + + # file path to directory containing forcings files from each GCM + dir2GCMs = os.path.join(archive_fp, file_struct) + + # should do some checking here that all the fps work... + self.dir2GCMs = dir2GCMs + + def get_filename(self, GCM, scenario, variable): + + # convert from variable name within NetCDF file, to variable name as it + # appears in the filename + if variable == "thermal_forcing": + fn_var = "oceanThermalForcing" + elif variable == "basin_runoff": + fn_var = "basinRunoff" + else: + raise ValueError(f"invalid varibale name: {variable}") + + # within filename scenario have for removed + scenario_no_dot = "".join(scenario.split(".")) + + fp = (f"{GCM.lower()}_{scenario}/" + f"MAR3.9_{GCM}_{scenario_no_dot}_{fn_var}_{self.version}.nc") + + # check that file exists!!! + return os.path.join(self.dir2GCMs, fp) + + +class atmosphereFileFinder: + + def __init__(self, archive_fp, version="v1", workdir="./"): + + self.version = version + + if os.path.exists(workdir): + self.workdir = workdir + else: + print("gotta do error checking") + + # file strucutre within Atmosphere_Forcing directory for navigating + file_struct = f"Atmosphere_Forcing/aSMB_observed/{version}" + + # file path to directory containing forcings files from each GCM + dir2GCMs = os.path.join(archive_fp, file_struct) + + # should do some checking here that all the fps work... + self.dir2GCMs = dir2GCMs + + def get_filename(self, GCM, scenario, variable, start=2015, end=2100): + """ + """ + if GCM == "NorESM1-M": + GCM = "NorESM1" + # get a list of all the yearly files within the period of intrest + yearly_files = self.__find_yearly_files(GCM, scenario, variable, + start, end) + # still need to make the output filename to write combined files to + out_fn = f"MAR3.9_{GCM}_{scenario}_{variable}_{start}--{end}.nc" + # relative to the workdir, which we've already checked if if existed + # make the filepath for the nested GCM/scenario/var direcotry struct. + top_fp = os.path.join(self.workdir, f"{GCM}-{scenario}/{variable}") + + # make the nested directory strucure; if needed + if not os.path.exists(top_fp): + os.makedirs(top_fp, exist_ok=True) + + out_fp = os.path.join(top_fp, out_fn) + + # only combine the files if they don't exist yet + if not os.path.exists(out_fp): + # with output filename, combine the files to a single + self.__combine_files(yearly_files, out_fp) + + return out_fp, yearly_files + + def __find_yearly_files(self, GCM, scenario, variable, start, end): + """ + """ + # within filename scenario have for removed + scenario_no_dot = "".join(scenario.split(".")) + + # create filename template to be globbed + fp = (f"{GCM}-{scenario_no_dot}/{variable}/" + f"{variable}_MARv3.9-yearly-{GCM}-{scenario_no_dot}-*.nc") + + # glob the files + all_files = glob.glob(os.path.join(self.dir2GCMs, fp)) + # all directories should have same number of files; so check here to + # make sure things worked properly + + # get the files within the desired period + files = sorted(f for f in all_files if check_year(f, start, end)) + + # make sure the length of fns returned matches the length (in years) + # of the period + period = (end - start) + 1 + + assert len(files) == period + return files + + def __combine_files(self, files, out_fn): + """ + """ + ds = xr.open_mfdataset(files, concat_dim="time", combine="nested", + data_vars='minimal', coords='minimal', + compat="broadcast_equals", + combine_attrs="override") + + ds.to_netcdf(out_fn) diff --git a/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py b/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py new file mode 100644 index 0000000000..49f0c95ff2 --- /dev/null +++ b/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py @@ -0,0 +1,94 @@ +from compass.landice.tests.ismip6_GrIS_forcing.create_mapping_files import ( + CreateMappingFiles, +) +from compass.landice.tests.ismip6_GrIS_forcing.file_finders import ( + atmosphereFileFinder, + oceanFileFinder, +) +from compass.landice.tests.ismip6_GrIS_forcing.process_forcing import ( + ProcessForcing, +) +from compass.testcase import TestCase + + +class ForcingGen(TestCase): + """ + A TestCase for remapping (i.e. interpolating) ISMIP6 GrIS forcing files + onto a MALI mesh + + Attributes + ---------- + + mali_mesh_scrip : str + filepath to the scrip file describing the MALI mesh + + ismip6_GrIS_scrip : str + filepath to the scrip file describing the grid the ISMIP6 GrIS + forcing data is on. (Note: All forcing files are on the same grid, + so only need one scrip file for all forcing files) + + remapping_weights : str + filepath to the `ESMF_RegirdWeightGen` generated mapping file + + experiments : dict + + """ + + def __init__(self, test_group): + """ + Parameters + ---------- + test_group : compass.landice.tests... + The test group that this test case belongs to + """ + name = "forcing_gen" + super().__init__(test_group=test_group, name=name) + + # filenames for remapping files, stored at the testcase level so they + # will be accesible to all steps in the testcase + self.mali_mesh_scrip = "mali_mesh.scrip.nc" + self.ismip6_GrIS_scrip = "ismip6_GrIS.scrip.nc" + self.remapping_weights = "ismip6_2_mali.weights.nc" + # place holder for file finders that will be initialized in `configure` + self.__atmFF = None + self.__ocnFF = None + + # precusssory step that builds scrip file for mali mesh, + # and generates a common weights file to be used in remapping + self.add_step(CreateMappingFiles(test_case=self)) + + # step that deals with racmo, do all I need to do is remap and average? + + # add steps that re-maps and processes downscaled GCM data for each + # experiment. + self.add_step(ProcessForcing(test_case=self)) + + def configure(self): + + config = self.config + # add ouputdir path to the remapping files + + # get the list of requested experiments + expr2run = config.getlist("ISMIP6_GrIS_Forcing", "experiments") + # get the dictionary of experiments, as defined in the yaml file + all_exprs = self.test_group.experiments + # get subset of dictionaries, based on requested expriments + self.experiments = {e: all_exprs[e] for e in expr2run} + + archive_fp = config.get("ISMIP6_GrIS_Forcing", "archive_fp") + + # initalize the oceanFileFinder + self.__ocnFF = oceanFileFinder(archive_fp) + self.__atmFF = atmosphereFileFinder(archive_fp) # , workdir=workdir) + + def findForcingFiles(self, GCM, scenario, variable): + """ + """ + + if variable in ["basin_runoff", "thermal_forcing"]: + forcing_fp = self.__ocnFF.get_filename(GCM, scenario, variable) + + if variable in ["aSMB", "aST", "dSMBdz", "dSTdz"]: + forcing_fp = self.__atmFF.get_filename(GCM, scenario, variable) + + return forcing_fp diff --git a/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/forcing_gen.cfg b/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/forcing_gen.cfg new file mode 100644 index 0000000000..3759c8a433 --- /dev/null +++ b/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/forcing_gen.cfg @@ -0,0 +1,11 @@ +[ISMIP6_GrIS_Forcing] + +# full path to MALI mesh forcing data will be regridded onto +MALI_mesh_fp = /global/cfs/cdirs/fanssie/MALI_input_files/GIS_1to10km_r02/GIS_1to10km_r02_20230202_relaxed.nc + +# Path to GrIS folder of the ISMIP6-Forcing-Ghub archive +archive_fp = /global/cfs/cdirs/fanssie/standard_datasets/ISMIP6-Forcing-Ghub/GrIS + +# list of experiments to generate forcing for. +# See `experiments.yml` for a list of supported experiments +experiments = ctrl, Exp05 diff --git a/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py b/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py new file mode 100644 index 0000000000..79898a9690 --- /dev/null +++ b/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py @@ -0,0 +1,99 @@ +import os + +import xarray as xr +from mpas_tools.logging import check_call + +from compass.step import Step + + +def datetime_2_xtime(ds, var="time"): + return xr.apply_ufunc(lambda t: t.strftime("%Y-%m-%d_%H:%M:%S").ljust(64), + ds[var], vectorize=True, + output_dtypes=["S"]) + + +class ProcessForcing(Step): + """ + + """ + + def __init__(self, test_case): + """ + """ + + name = "process_forcing" + + super().__init__(test_case=test_case, name=name, subdir=None, + cpus_per_task=1, min_cpus_per_task=1) + # read and store the experiment dictionary + + # initalize the FileFinders are store them as dict? + + def run(self): + + # list of variables to process + atmosphere_vars = ["aSMB", "aST", "dSMBdz", "dSTdz"] + ocean_vars = ["basin_runoff", "thermal_forcing"] + + # loop over experiments passed via config file + for expr_name, expr_dict in self.test_case.experiments.items(): + # print to logger which expriment we are on (i/n) + + GCM = expr_dict["GCM"] + scenario = expr_dict["Scenario"] + start = expr_dict["start"] + end = expr_dict["end"] + + # need a special way to treat the RACMO datasets + if expr_name in ["ctrl", "hist"]: + continue + + # loop over ocean variables seperately + for var in ocean_vars: + var_fp = self.test_case.findForcingFiles(GCM, scenario, var) + + remapped_fp = f"gis_{var}_{GCM}_{scenario}_{start}-{end}.nc" + + self.remap_variable(var_fp, remapped_fp) + print(var_fp, os.path.exists(var_fp)) + + # loop over atmosphere variables + for variable in atmosphere_vars: + pass + + def process_variable(self, variables): + # variable object be initialized in the setup? + + # remap using the landice framework function + + # open the remapped file for post processing + # remapped_ds = xr.open_dataset(self.out_fn) + + # if var == "SMB": + # SMB is not a var in ISMIP6 forcing files, need to use `SMB_ref` + # and the processed `aSMB` to produce a `SMB` forcing + + # rename varibales + # self.rename_variables(remapped_ds) + + # convert time to xtime + # self.create_xtime(remapped_ds) + pass + + def remap_variable(self, input_file, output_file): + """ + """ + + # remap the forcing file onto the MALI mesh + args = ["ncremap", + "-i", input_file, + "-o", output_file, + "-m", self.test_case.remapping_weights] + + check_call(args, logger=self.logger) + + def rename_variables(self, ds, var): + pass + + def create_xtime(self, ds): + pass From e33b7cd8047d5b22fcda45173a890d22bbbcfeaa Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Thu, 23 May 2024 09:47:31 -0700 Subject: [PATCH 02/12] Add first pass at racmo climatology calculation. --- .../tests/ismip6_GrIS_forcing/__init__.py | 15 +- .../create_mapping_files.py | 147 +++++++++++++----- .../forcing_gen/__init__.py | 19 ++- .../forcing_gen/forcing_gen.cfg | 11 ++ .../ismip6_GrIS_forcing/process_forcing.py | 19 ++- .../ref_smb_climatology.py | 104 +++++++++++++ 6 files changed, 263 insertions(+), 52 deletions(-) create mode 100644 compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py diff --git a/compass/landice/tests/ismip6_GrIS_forcing/__init__.py b/compass/landice/tests/ismip6_GrIS_forcing/__init__.py index 52b9386d36..232d8428fc 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/__init__.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/__init__.py @@ -1,3 +1,5 @@ +import os + import yaml from compass.landice.tests.ismip6_GrIS_forcing.forcing_gen import ForcingGen @@ -17,13 +19,18 @@ def __init__(self, mpas_core): self.add_test_case(ForcingGen(test_group=self)) - with open("experiments.yml", "r") as f: + # open, read, and validate the experiment file + self._read_and_validate_experiment_file() + + def _read_and_validate_experiment_file(self): + + # get the filepath to current module, needed of opening experiment file + module_path = os.path.dirname(os.path.realpath(__file__)) + + with open(f"{module_path}/experiments.yml", "r") as f: experiments = yaml.safe_load(f) # should I filter through the experiments dictionary to make sure # everything is valid.... self.experiments = experiments - - def __read_and_validate_experiment_file(self): - pass diff --git a/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py b/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py index e84b3cd9ad..3ea737e08d 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py @@ -28,6 +28,7 @@ def __init__(self, test_case): # be propely in `setup` method so user specified config options can # be parsed self.mali_mesh = None + self.racmo_grid = None def setup(self): """ @@ -35,10 +36,11 @@ def setup(self): # parse user specified parameters from the config config = self.config - section = config["ISMIP6_GrIS_Forcing"] + forcing_section = config["ISMIP6_GrIS_Forcing"] + smb_ref_section = config["smb_ref_climatology"] # read the mesh path/name from the config file - mali_mesh = section.get("MALI_mesh_fp") + mali_mesh = forcing_section.get("MALI_mesh_fp") # make sure the mesh path/name is valid and accessible if not os.path.exists(mali_mesh): @@ -48,12 +50,34 @@ def setup(self): self.mali_mesh = mali_mesh self.add_input_file(filename=self.mali_mesh) + # + racmo_directory = smb_ref_section.get("racmo_directory") + # this filename should probably just be hardcoded..... + racmo_grid_fn = smb_ref_section.get("racmo_grid_fn") + # combine the directory and filename + racmo_grid_fp = os.path.join(racmo_directory, racmo_grid_fn) + + # make sure the combined filename exists + if not os.path.exists(racmo_grid_fp): + # check if the parent directory exists + if not os.path.exists(racmo_directory): + raise FileNotFoundError(f"{racmo_directory} does not exist") + # the parent directory exists but the forcing file does not + else: + raise FileNotFoundError(f"{racmo_grid_fp} does not exist") + + # add the racmo grid as an attribute and as an input file + self.racmo_grid = racmo_grid_fp + self.add_input_file(filename=racmo_grid_fp) + # loop over mapping file `test_case` attributes and add the current # steps `work_dir` to the absolute path. This will ensure other steps # will be able to access the mapping files created here for attr in ["mali_mesh_scrip", - "ismip6_GrIS_scrip", - "remapping_weights"]: + "racmo_gis_scrip", + "ismip6_gis_scrip", + "racmo_2_mali_weights", + "ismip6_2_mali_weights"]: # get the attribute values file_name = getattr(self.test_case, attr) # combine filename w/ the testcases work dir @@ -64,8 +88,10 @@ def setup(self): # Add the testcase attributes as output files; since they'll be # generated as part of this step self.add_output_file(filename=self.test_case.mali_mesh_scrip) - self.add_output_file(filename=self.test_case.ismip6_GrIS_scrip) - self.add_output_file(filename=self.test_case.remapping_weights) + self.add_output_file(filename=self.test_case.racmo_gis_scrip) + self.add_output_file(filename=self.test_case.ismip6_gis_scrip) + self.add_output_file(filename=self.test_case.racmo_2_mali_weights) + self.add_output_file(filename=self.test_case.ismip6_2_mali_weights) def run(self): """ @@ -73,53 +99,60 @@ def run(self): # unpack the filepaths stored as attributes of the test_case mali_mesh_scrip = self.test_case.mali_mesh_scrip - ismip6_GrIS_scrip = self.test_case.ismip6_GrIS_scrip - remapping_weights = self.test_case.remapping_weights + racmo_gis_scrip = self.test_case.racmo_gis_scrip + ismip6_gis_scrip = self.test_case.ismip6_gis_scrip + racmo_2_mali_weights = self.test_case.racmo_2_mali_weights + ismip6_2_mali_weights = self.test_case.ismip6_2_mali_weights # create scrip file describing MALI mesh that forcing variables # will be interpolated onto scrip_from_mpas(self.mali_mesh, mali_mesh_scrip) - # get a list of the projection experiments requesed. - # (i.e. experiments names that use ISMIP6 forcing) - proj_exprs = list(filter(lambda x: "Exp" in x, - self.test_case.experiments.keys())) - - if proj_exprs: - # b/c all GrIS forcing files are on the same grid, it doesn't - # matter what expr we use; so just use the first suitable candiate - expr = proj_exprs[0] - - expr_params = self.test_case.experiments[expr] - - GCM = expr_params["GCM"] - scenario = expr_params["Scenario"] - variable = "thermal_forcing" - # get the filepath for any given forcing file (since all files are on - # the same grid), in order to generate a scrip file for the grid - forcing_fp = self.test_case.findForcingFiles(GCM, scenario, variable) - # WGS 84 / NSIDC Sea Ice Polar Stereographic North Projection epsg_3413 = Proj("EPSG:3413") + # make scrip file for racmo grid and mapping weights for racmo->mpas + self.make_forcing_descriptor_and_weights(self.racmo_grid, + racmo_gis_scrip, + mali_mesh_scrip, + epsg_3413, + racmo_2_mali_weights) + + # finding the right ismip6 forcing file is complicated, + # so let's call a helper function to do that. + ismip6_forcing_fp = self._find_ismip6_forcing_files() + + # make scrip file for ismip6 grid and mapping weights for ismip6->mpas + self.make_forcing_descriptor_and_weights(ismip6_forcing_fp, + ismip6_gis_scrip, + mali_mesh_scrip, + epsg_3413, + ismip6_2_mali_weights) + + def make_forcing_descriptor_and_weights(self, + forcing_file_fp, + forcing_scrip_fp, + mali_mesh_scrip_fp, + forcing_proj, + weights_file_fp): + # open up forcing grid; get x and y coordinates as numpy arrays - with xr.open_dataset(forcing_fp) as ds: + with xr.open_dataset(forcing_file_fp) as ds: x = ds.x.values y = ds.y.values - # create pyremap `MeshDescriptor` obj for the ISMIP6 GrIS forcing grid + # create pyremap `MeshDescriptor` obj for the forcing grid forcingDescriptor = ProjectionGridDescriptor.create( - projection=epsg_3413, x=x, y=y, meshName=forcing_fp) + projection=forcing_proj, x=x, y=y, meshName=forcing_file_fp) - # write scrip file describing ISMIP6 GrIS forcing grid - forcingDescriptor.to_scrip(ismip6_GrIS_scrip) + # write scrip file describing forcing grid + forcingDescriptor.to_scrip(forcing_scrip_fp) - # generate mapping weights between ismip6 forcing data and the MALI - # mesh, using `ESMF_RegridWeightGen` - args = ['srun', '-n', '128', 'ESMF_RegridWeightGen', - '-s', ismip6_GrIS_scrip, - '-d', mali_mesh_scrip, - '-w', remapping_weights, + # generate mapping weights between forcing data and the MALI mesh + args = ['srun', '-n', '128', 'ESMF_RegridWeigh:tGen', + '-s', forcing_scrip_fp, + '-d', mali_mesh_scrip_fp, + '-w', weights_file_fp, '--method', 'conserve', "--netcdf4", # "--no_log", @@ -134,7 +167,41 @@ def run(self): # check if glob returns an empty list if log_files: + # split the output weights filename so that only the basename w/ + # no extension is left, in order to clean up the log files + base_name = os.path.splitext(os.path.basename(weights_file_fp))[0] + dir_name = f"{base_name}.Logs/" + + # check if there is already a Log files directory, if so wipe it + if os.path.exists(dir_name): + shutil.rmtree(dir_name) + # make the log directory - os.mkdir("logs/") + os.mkdir(dir_name) + # copy the log files to the log directory; using list comprehension - _ = [shutil.move(f, "logs/") for f in log_files] + _ = [shutil.move(f, dir_name) for f in log_files] + + def _find_ismip6_forcing_files(self): + """ + """ + # get a list of all projection experiments requesed. + # (i.e. experiments names that use ISMIP6 forcing) + proj_exprs = list(filter(lambda x: "Exp" in x, + self.test_case.experiments.keys())) + + if proj_exprs: + # b/c all GrIS forcing files are on the same grid, it doesn't + # matter what expr we use; so just use the first suitable candiate + expr = proj_exprs[0] + + expr_params = self.test_case.experiments[expr] + GCM = expr_params["GCM"] + scenario = expr_params["Scenario"] + variable = "thermal_forcing" + + # get the filepath for any given forcing file (since all files are on + # the same grid), in order to generate a scrip file for the grid + forcing_fp = self.test_case.findForcingFiles(GCM, scenario, variable) + + return forcing_fp diff --git a/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py b/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py index 49f0c95ff2..3a69eb8b6c 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py @@ -8,6 +8,9 @@ from compass.landice.tests.ismip6_GrIS_forcing.process_forcing import ( ProcessForcing, ) +from compass.landice.tests.ismip6_GrIS_forcing.ref_smb_climatology import ( + SMBRefClimatology, +) from compass.testcase import TestCase @@ -44,11 +47,17 @@ def __init__(self, test_group): name = "forcing_gen" super().__init__(test_group=test_group, name=name) - # filenames for remapping files, stored at the testcase level so they - # will be accesible to all steps in the testcase + # NOTE: scrip and weight filenames are stored at the testcase level + # so they will be accesible to all steps in the testcase + + # filenames for scrip files (i.e. grid descriptors) self.mali_mesh_scrip = "mali_mesh.scrip.nc" - self.ismip6_GrIS_scrip = "ismip6_GrIS.scrip.nc" - self.remapping_weights = "ismip6_2_mali.weights.nc" + self.racmo_gis_scrip = "racmo_gis.scrip.nc" + self.ismip6_gis_scrip = "ismip6_GrIS.scrip.nc" + # filenames of remapping files + self.racmo_2_mali_weights = "racmo_2_mali.weights.nc" + self.ismip6_2_mali_weights = "ismip6_2_mali.weights.nc" + # place holder for file finders that will be initialized in `configure` self.__atmFF = None self.__ocnFF = None @@ -59,6 +68,8 @@ def __init__(self, test_group): # step that deals with racmo, do all I need to do is remap and average? + self.add_step(SMBRefClimatology(test_case=self)) + # add steps that re-maps and processes downscaled GCM data for each # experiment. self.add_step(ProcessForcing(test_case=self)) diff --git a/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/forcing_gen.cfg b/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/forcing_gen.cfg index 3759c8a433..becc05ce41 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/forcing_gen.cfg +++ b/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/forcing_gen.cfg @@ -1,3 +1,14 @@ +[smb_ref_climatology] +# path to racmo data... or could just have it on lcrc +racmo_directory = /global/cfs/cdirs/fanssie/standard_datasets/RACMO2p3 +# should this be included here, or just hardcoded w/in the step? +racmo_grid_fn = Icemask_Topo_Iceclasses_lon_lat_average_1km_GrIS.nc +# see above about hardcoding +racmo_smb_fn = smb_rec.1958-2019.BN_RACMO2.3p2_FGRN055_GrIS.MM.nc +# +climatology_start = 1960 +climatology_end = 1989 + [ISMIP6_GrIS_Forcing] # full path to MALI mesh forcing data will be regridded onto diff --git a/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py b/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py index 79898a9690..fc73f9f18c 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py @@ -12,6 +12,10 @@ def datetime_2_xtime(ds, var="time"): output_dtypes=["S"]) +{"thermal_forcing": "ismip6_2dThermalForcing", + "basin_runoff": "ismip6Runoff"} + + class ProcessForcing(Step): """ @@ -54,7 +58,10 @@ def run(self): remapped_fp = f"gis_{var}_{GCM}_{scenario}_{start}-{end}.nc" - self.remap_variable(var_fp, remapped_fp) + self.remap_variable(var_fp, + remapped_fp, + self.test_case.ismip6_2_mali_weights) + print(var_fp, os.path.exists(var_fp)) # loop over atmosphere variables @@ -80,7 +87,7 @@ def process_variable(self, variables): # self.create_xtime(remapped_ds) pass - def remap_variable(self, input_file, output_file): + def remap_variable(self, input_file, output_file, weights_file): """ """ @@ -88,11 +95,15 @@ def remap_variable(self, input_file, output_file): args = ["ncremap", "-i", input_file, "-o", output_file, - "-m", self.test_case.remapping_weights] + "-m", weights_file] check_call(args, logger=self.logger) - def rename_variables(self, ds, var): + def rename_variable_and_trim_dataset(self, ds, var): + + # drop unnecessary variables + ds = ds.drop_vars(["lat_vertices", "area", + "lon_vertices", "lat", "lon"]) pass def create_xtime(self, ds): diff --git a/compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py b/compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py new file mode 100644 index 0000000000..dee13a7a8c --- /dev/null +++ b/compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py @@ -0,0 +1,104 @@ +import os + +import xarray as xr +from mpas_tools.logging import check_call + +from compass.step import Step + + +class SMBRefClimatology(Step): + """ + """ + + def __init__(self, test_case): + """ + """ + name = "smb_ref_climatology" + subdir = None + + super().__init__(test_case=test_case, name=name, subdir=subdir) + + # initalize the attributes with empty values for now, attributes will + # be propely in `setup` method so user specified config options can + # be parsed + self.racmo_smb_fp = None + self.smb_ref_climatology = None + + def setup(self): + """ + """ + + # parse user specified parameters from the config + config = self.config + smb_ref_section = config["smb_ref_climatology"] + + # + racmo_directory = smb_ref_section.get("racmo_directory") + # this filename should probably just be hardcoded..... + racmo_smb_fn = smb_ref_section.get("racmo_smb_fn") + # combine the directory and filename + racmo_smb_fp = os.path.join(racmo_directory, racmo_smb_fn) + + # make sure the combined filename exists + if not os.path.exists(racmo_smb_fp): + # check if the parent directory exists + if not os.path.exists(racmo_directory): + raise FileNotFoundError(f"{racmo_directory} does not exist") + # the parent directory exists but the forcing file does not + else: + raise FileNotFoundError(f"{racmo_smb_fp} does not exist") + + # add the racmo smb as an attribute and as an input file + self.racmo_smb = racmo_smb_fp + self.add_input_file(filename=racmo_smb_fp) + + # get the start and end dates for the climatological mean + clima_start = smb_ref_section.getint("climatology_start") + clima_end = smb_ref_section.getint("climatology_end") + + # make a descriptive filename based on climatology period + clima_fn = f"racmo_climatology_{clima_start}--{clima_end}.nc" + + self.smb_ref_climatology = clima_fn + self.add_output_file(filename=clima_fn) + + def run(self): + """ + """ + + # parse user specified parameters from the config + config = self.config + smb_ref_section = config["smb_ref_climatology"] + # get the start and end dates for the climatological mean + clima_start = smb_ref_section.getint("climatology_start") + clima_end = smb_ref_section.getint("climatology_end") + + # remap the gridded racmo data onto the mpas grid + self.remap_variable(self.racmo_smb, + self.smb_ref_climatology, + self.test_case.racmo_2_mali_weights) + + ds = xr.open_dataset(self.smb_ref_climatology, decode_times=False) + + s_idx = ((clima_start - 1958) * 12) - 1 + e_idx = ((clima_end - 1958) * 12) - 1 + + climatology = ds.SMB_rec.isel(time=slice(s_idx, e_idx)).mean("time") + + climatology.to_netcdf(self.smb_ref_climatology, "w") + + # trim unused variables + + # take temporal average + + def remap_variable(self, input_file, output_file, weights_file): + """ + """ + + # remap the forcing file onto the MALI mesh + args = ["ncremap", + "-i", input_file, + "-o", output_file, + "-m", weights_file] + + check_call(args, logger=self.logger) From 9ac3954c103fb21be7849090976be1e725b1b8e9 Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Tue, 11 Jun 2024 15:23:10 -0700 Subject: [PATCH 03/12] Add `smb_ref_climatology` as testcase attribute. Ensure's the climatology can be found by multiple steps, in order process the SMB anomoly files provided by ISMIP6. --- .../forcing_gen/__init__.py | 4 +- .../ref_smb_climatology.py | 41 +++++++++++++------ 2 files changed, 31 insertions(+), 14 deletions(-) diff --git a/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py b/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py index 3a69eb8b6c..78a721df31 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py @@ -57,7 +57,8 @@ def __init__(self, test_group): # filenames of remapping files self.racmo_2_mali_weights = "racmo_2_mali.weights.nc" self.ismip6_2_mali_weights = "ismip6_2_mali.weights.nc" - + # filename of the reference climatology + self.smb_ref_climatology = None # place holder for file finders that will be initialized in `configure` self.__atmFF = None self.__ocnFF = None @@ -67,7 +68,6 @@ def __init__(self, test_group): self.add_step(CreateMappingFiles(test_case=self)) # step that deals with racmo, do all I need to do is remap and average? - self.add_step(SMBRefClimatology(test_case=self)) # add steps that re-maps and processes downscaled GCM data for each diff --git a/compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py b/compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py index dee13a7a8c..e9b6b718c1 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py @@ -1,6 +1,7 @@ import os import xarray as xr +from mpas_tools.io import write_netcdf from mpas_tools.logging import check_call from compass.step import Step @@ -58,9 +59,13 @@ def setup(self): # make a descriptive filename based on climatology period clima_fn = f"racmo_climatology_{clima_start}--{clima_end}.nc" + # add the steps workdir to the filename to make it a full path + clima_fp = os.path.join(self.work_dir, clima_fn) - self.smb_ref_climatology = clima_fn - self.add_output_file(filename=clima_fn) + # set the testcase attribute and add it as an output file for this + # step so the climatology will by useable by other steps + self.test_case.smb_ref_climatology = clima_fp + self.add_output_file(filename=self.test_case.smb_ref_climatology) def run(self): """ @@ -75,21 +80,32 @@ def run(self): # remap the gridded racmo data onto the mpas grid self.remap_variable(self.racmo_smb, - self.smb_ref_climatology, + self.test_case.smb_ref_climatology, self.test_case.racmo_2_mali_weights) - ds = xr.open_dataset(self.smb_ref_climatology, decode_times=False) + ds = xr.open_dataset(self.test_case.smb_ref_climatology, + decode_times=False) + # find indices of climatology start/end (TO DO: make more robust) s_idx = ((clima_start - 1958) * 12) - 1 e_idx = ((clima_end - 1958) * 12) - 1 - climatology = ds.SMB_rec.isel(time=slice(s_idx, e_idx)).mean("time") - - climatology.to_netcdf(self.smb_ref_climatology, "w") - - # trim unused variables - - # take temporal average + # calculate climatology + ds["SMB_rec"] = ds.SMB_rec.isel(time=slice(s_idx, e_idx)).mean("time") + # rename variables to match MALI/MPAS conventiosn + ds = ds.rename(SMB_rec="sfcMassBal", ncol="nCells") + # drop unused dimensions + ds = ds.drop_dims(["time", "nv"]) + # drop un-needed varibales + ds = ds.drop_vars(["area", "lat", "lon"]) + # convert `sfcMassBal` to MPAS units + ds["sfcMassBal"] /= (60 * 60 * 24 * 365) / 12. + # add a units attribute to `sfcMassBal` + ds["sfcMassBal"].attrs["units"] = "kg m-2 s-1" + # expand sfcMassBal dimension to match what MALI expects + ds["sfcMassBal"] = ds.sfcMassBal.expand_dims("Time") + # write the file + write_netcdf(ds, self.test_case.smb_ref_climatology) def remap_variable(self, input_file, output_file, weights_file): """ @@ -99,6 +115,7 @@ def remap_variable(self, input_file, output_file, weights_file): args = ["ncremap", "-i", input_file, "-o", output_file, - "-m", weights_file] + "-m", weights_file, + "-v", "SMB_rec"] check_call(args, logger=self.logger) From be69711b840b100364e09c83190a566e2eea9b58 Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Tue, 11 Jun 2024 15:26:17 -0700 Subject: [PATCH 04/12] Fix typo in `ESMFWeightGen` command. --- .../landice/tests/ismip6_GrIS_forcing/create_mapping_files.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py b/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py index 3ea737e08d..511983276e 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py @@ -149,7 +149,7 @@ def make_forcing_descriptor_and_weights(self, forcingDescriptor.to_scrip(forcing_scrip_fp) # generate mapping weights between forcing data and the MALI mesh - args = ['srun', '-n', '128', 'ESMF_RegridWeigh:tGen', + args = ['srun', '-n', '128', 'ESMF_RegridWeightGen', '-s', forcing_scrip_fp, '-d', mali_mesh_scrip_fp, '-w', weights_file_fp, From 078bd630c44d2976f03c64a36de681609669da6e Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Tue, 11 Jun 2024 15:35:37 -0700 Subject: [PATCH 05/12] Build out forcing processing step --- .../ismip6_GrIS_forcing/process_forcing.py | 108 +++++++++++++----- 1 file changed, 77 insertions(+), 31 deletions(-) diff --git a/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py b/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py index fc73f9f18c..b1b4ad346c 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py @@ -1,19 +1,22 @@ -import os - import xarray as xr +from mpas_tools.io import write_netcdf from mpas_tools.logging import check_call from compass.step import Step -def datetime_2_xtime(ds, var="time"): - return xr.apply_ufunc(lambda t: t.strftime("%Y-%m-%d_%H:%M:%S").ljust(64), - ds[var], vectorize=True, - output_dtypes=["S"]) +def datetime_2_xtime(ds, var="Time", xtime_fmt="%Y-%m-%d_%H:%M:%S"): + + return xr.apply_ufunc(lambda t: t.strftime(xtime_fmt).ljust(64), + ds[var], vectorize=True, output_dtypes=["S"]) -{"thermal_forcing": "ismip6_2dThermalForcing", - "basin_runoff": "ismip6Runoff"} +renaming_dict = {"thermal_forcing": "ismip6_2dThermalForcing", + "basin_runoff": "ismip6Runoff", + "aSMB": "sfcMassBal", + "dSMBdz": "sfcMassBal_lapseRate", + "aST": "surfaceAirTemperature", + "dSTdz": "surfaceAirTemperature_lapseRate"} class ProcessForcing(Step): @@ -52,40 +55,71 @@ def run(self): if expr_name in ["ctrl", "hist"]: continue + ocean_forcing_datastes = [] # loop over ocean variables seperately for var in ocean_vars: var_fp = self.test_case.findForcingFiles(GCM, scenario, var) - remapped_fp = f"gis_{var}_{GCM}_{scenario}_{start}-{end}.nc" + ocean_forcing_datastes.append( + self.process_variable(GCM, scenario, var, + var_fp, start, end)) + ocean_forcing_ds = xr.combined_by_coords( + ocean_forcing_datastes, combine_attrs="drop_conflicts") - self.remap_variable(var_fp, - remapped_fp, - self.test_case.ismip6_2_mali_weights) - - print(var_fp, os.path.exists(var_fp)) + write_netcdf(ocean_forcing_ds, "ocean_forcing_test.nc") # loop over atmosphere variables - for variable in atmosphere_vars: - pass + for var in atmosphere_vars: + # this pretty hacky with the tuple, fix the file finder + var_fp, _ = self.test_case.findForcingFiles(GCM, scenario, var) + + self.process_variable(GCM, scenario, var, var_fp, start, end) - def process_variable(self, variables): - # variable object be initialized in the setup? + def process_variable(self, GCM, scenario, var, forcing_fp, start, end): - # remap using the landice framework function + # create the remapped file name + remapped_fp = f"gis_{var}_{GCM}_{scenario}_{start}-{end}.nc" + + # remap the forcing file + self.remap_variable(forcing_fp, + remapped_fp, + self.test_case.ismip6_2_mali_weights) + + # Now that forcing has been regridded onto the MALI grid, we can drop + # ISMIP6 grid variables. Special note: the char field `mapping` + # particularily causes problem with `xtime` (40byte vs. 64byte chars) + vars_2_drop = ["time_bounds", "lat_vertices", "lon_vertices", "area", + "mapping", "lat", "lon"] # open the remapped file for post processing - # remapped_ds = xr.open_dataset(self.out_fn) + remapped_ds = xr.open_dataset(remapped_fp, + drop_variables=vars_2_drop, + use_cftime=True) + + # convert time to xtime + remapped_ds["xtime"] = datetime_2_xtime(remapped_ds, var="time") + + # rename the variable/dimensions to match MPAS/MALI conventions + remapped_ds = remapped_ds.rename({"ncol": "nCells", + "time": "Time", + var: renaming_dict[var]}) - # if var == "SMB": # SMB is not a var in ISMIP6 forcing files, need to use `SMB_ref` # and the processed `aSMB` to produce a `SMB` forcing + if var == "aSMB": + # open the reference climatology file + smb_ref = xr.open_dataset(self.test_case.smb_ref_climatology) + # squeeze the empty time dimension so that broadcasting works + smb_ref = smb_ref.squeeze("Time") + # add climatology to anomoly to make full forcing field + remapped_ds[renaming_dict[var]] += smb_ref["sfcMassBal"] - # rename varibales - # self.rename_variables(remapped_ds) + # need to process surfaceAirTemperature + if var == "surfaceAirTemperature": + pass - # convert time to xtime - # self.create_xtime(remapped_ds) - pass + # write_netcdf(remapped_ds, remapped_fp) + return remapped_ds def remap_variable(self, input_file, output_file, weights_file): """ @@ -99,12 +133,24 @@ def remap_variable(self, input_file, output_file, weights_file): check_call(args, logger=self.logger) - def rename_variable_and_trim_dataset(self, ds, var): + def rename_variable_and_trim_dataset(self, ds, var_src, var_dest): # drop unnecessary variables - ds = ds.drop_vars(["lat_vertices", "area", - "lon_vertices", "lat", "lon"]) - pass + ds = ds.drop_vars(["lat_vertices", "area", "lon_vertices", + "lat", "lon"]) + + ds = ds.rename({"ncol": "nCells", + "time": "Time", + var_src: var_dest}) + + # need to make `Time` the unlimited dimension, which also prevents + # `time` dimension for being added back to the dataset + # ds.encoding["unlimited_dims"] = {"Time"} + + return ds def create_xtime(self, ds): - pass + + ds["xtime"] = datetime_2_xtime(ds, var="Time") + + return ds From e1c97a40c5fec12482128d4aeee28b95385e36ae Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Wed, 12 Jun 2024 15:06:33 -0700 Subject: [PATCH 06/12] Build out the regridding and processing step. - Now packing atm/ocn forcing vars into their own respective forcing files - moved shared utilites to common file - fixed bug in returning conctenated atm files from file finder --- .../tests/ismip6_GrIS_forcing/file_finders.py | 2 +- .../ismip6_GrIS_forcing/process_forcing.py | 133 ++++++++---------- .../tests/ismip6_GrIS_forcing/utilities.py | 50 +++++++ 3 files changed, 112 insertions(+), 73 deletions(-) create mode 100644 compass/landice/tests/ismip6_GrIS_forcing/utilities.py diff --git a/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py b/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py index 7802874711..9f41ca0eb3 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py @@ -99,7 +99,7 @@ def get_filename(self, GCM, scenario, variable, start=2015, end=2100): # with output filename, combine the files to a single self.__combine_files(yearly_files, out_fp) - return out_fp, yearly_files + return out_fp def __find_yearly_files(self, GCM, scenario, variable, start, end): """ diff --git a/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py b/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py index b1b4ad346c..d8f8b0cf45 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py @@ -1,16 +1,14 @@ +import os + import xarray as xr from mpas_tools.io import write_netcdf -from mpas_tools.logging import check_call +from compass.landice.tests.ismip6_GrIS_forcing.utilities import ( + add_xtime, + remap_variables, +) from compass.step import Step - -def datetime_2_xtime(ds, var="Time", xtime_fmt="%Y-%m-%d_%H:%M:%S"): - - return xr.apply_ufunc(lambda t: t.strftime(xtime_fmt).ljust(64), - ds[var], vectorize=True, output_dtypes=["S"]) - - renaming_dict = {"thermal_forcing": "ismip6_2dThermalForcing", "basin_runoff": "ismip6Runoff", "aSMB": "sfcMassBal", @@ -55,102 +53,93 @@ def run(self): if expr_name in ["ctrl", "hist"]: continue - ocean_forcing_datastes = [] - # loop over ocean variables seperately - for var in ocean_vars: - var_fp = self.test_case.findForcingFiles(GCM, scenario, var) + atm_fn = f"gis_atm_forcing_{GCM}_{scenario}_{start}--{end}.nc" + self.process_variables(GCM, scenario, start, end, atmosphere_vars, + atm_fn) + + ocn_fn = f"gis_ocn_forcing_{GCM}_{scenario}_{start}--{end}.nc" + self.process_variables(GCM, scenario, start, end, ocean_vars, + ocn_fn) + + def process_variables(self, GCM, scenario, start, end, + variables, output_fn): - ocean_forcing_datastes.append( - self.process_variable(GCM, scenario, var, - var_fp, start, end)) - ocean_forcing_ds = xr.combined_by_coords( - ocean_forcing_datastes, combine_attrs="drop_conflicts") + forcing_datastes = [] - write_netcdf(ocean_forcing_ds, "ocean_forcing_test.nc") + for var in variables: + var_fp = self.test_case.findForcingFiles(GCM, scenario, var) - # loop over atmosphere variables - for var in atmosphere_vars: - # this pretty hacky with the tuple, fix the file finder - var_fp, _ = self.test_case.findForcingFiles(GCM, scenario, var) + ds = self.process_variable(GCM, scenario, var, var_fp, start, end) - self.process_variable(GCM, scenario, var, var_fp, start, end) + forcing_datastes.append(ds) + + forcing_ds = xr.merge(forcing_datastes) + + write_netcdf(forcing_ds, output_fn) def process_variable(self, GCM, scenario, var, forcing_fp, start, end): - # create the remapped file name - remapped_fp = f"gis_{var}_{GCM}_{scenario}_{start}-{end}.nc" + # + config = self.config + # + renamed_var = renaming_dict[var] + + # create the remapped file name, using original variable name + remapped_fn = f"remapped_{GCM}_{scenario}_{var}_{start}-{end}.nc" + # follow the directory structure of the concatenated source files + remapped_fp = os.path.join(f"{GCM}-{scenario}/{var}", remapped_fn) # remap the forcing file - self.remap_variable(forcing_fp, - remapped_fp, - self.test_case.ismip6_2_mali_weights) + remap_variables(forcing_fp, + remapped_fp, + self.test_case.ismip6_2_mali_weights, + variables=[var], + logger=self.logger) # Now that forcing has been regridded onto the MALI grid, we can drop # ISMIP6 grid variables. Special note: the char field `mapping` # particularily causes problem with `xtime` (40byte vs. 64byte chars) vars_2_drop = ["time_bounds", "lat_vertices", "lon_vertices", "area", - "mapping", "lat", "lon"] + "mapping", "lat", "lon", "polar_stereographic"] # open the remapped file for post processing remapped_ds = xr.open_dataset(remapped_fp, drop_variables=vars_2_drop, use_cftime=True) - # convert time to xtime - remapped_ds["xtime"] = datetime_2_xtime(remapped_ds, var="time") + # mask of desired time indices + mask = remapped_ds.time.dt.year.isin(range(start, end)) + # drop the non-desired time indices from remapped dataset + remapped_ds = remapped_ds.where(mask, drop=True) + # add xtime variable based on `time` + remapped_ds["xtime"] = add_xtime(remapped_ds, var="time") # rename the variable/dimensions to match MPAS/MALI conventions remapped_ds = remapped_ds.rename({"ncol": "nCells", "time": "Time", - var: renaming_dict[var]}) + var: renamed_var}) + + # drop the unneeded attributes from the dataset + for _var in remapped_ds: + remapped_ds[_var].attrs.pop("grid_mapping", None) + remapped_ds[_var].attrs.pop("cell_measures", None) # SMB is not a var in ISMIP6 forcing files, need to use `SMB_ref` # and the processed `aSMB` to produce a `SMB` forcing - if var == "aSMB": + if renamed_var == "sfcMassBal": # open the reference climatology file smb_ref = xr.open_dataset(self.test_case.smb_ref_climatology) # squeeze the empty time dimension so that broadcasting works smb_ref = smb_ref.squeeze("Time") # add climatology to anomoly to make full forcing field - remapped_ds[renaming_dict[var]] += smb_ref["sfcMassBal"] - - # need to process surfaceAirTemperature - if var == "surfaceAirTemperature": - pass - - # write_netcdf(remapped_ds, remapped_fp) - return remapped_ds - - def remap_variable(self, input_file, output_file, weights_file): - """ - """ - - # remap the forcing file onto the MALI mesh - args = ["ncremap", - "-i", input_file, - "-o", output_file, - "-m", weights_file] - - check_call(args, logger=self.logger) - - def rename_variable_and_trim_dataset(self, ds, var_src, var_dest): + remapped_ds[renamed_var] += smb_ref[renamed_var] - # drop unnecessary variables - ds = ds.drop_vars(["lat_vertices", "area", "lon_vertices", - "lat", "lon"]) - - ds = ds.rename({"ncol": "nCells", - "time": "Time", - var_src: var_dest}) - - # need to make `Time` the unlimited dimension, which also prevents - # `time` dimension for being added back to the dataset - # ds.encoding["unlimited_dims"] = {"Time"} - - return ds - - def create_xtime(self, ds): + if renamed_var == "surfaceAirTemperature": + # read the mesh path/name from the config file + mali_mesh_fp = config.get("ISMIP6_GrIS_Forcing", "MALI_mesh_fp") + # squeeze the empty time dimension so that broadcasting works + mali_ds = xr.open_dataset(mali_mesh_fp).squeeze("Time") - ds["xtime"] = datetime_2_xtime(ds, var="Time") + remapped_ds[renamed_var] += mali_ds[renamed_var] - return ds + return remapped_ds diff --git a/compass/landice/tests/ismip6_GrIS_forcing/utilities.py b/compass/landice/tests/ismip6_GrIS_forcing/utilities.py new file mode 100644 index 0000000000..b553ef8ebb --- /dev/null +++ b/compass/landice/tests/ismip6_GrIS_forcing/utilities.py @@ -0,0 +1,50 @@ +import xarray as xr +from mpas_tools.logging import check_call + + +def _datetime_2_xtime(ds, var="Time", xtime_fmt="%Y-%m-%d_%H:%M:%S"): + + return xr.apply_ufunc(lambda t: t.strftime(xtime_fmt).ljust(64), + ds[var], vectorize=True, output_dtypes=["S"]) + + +def add_xtime(ds, var="Time"): + """ + ds["xtime"] = add_xtime(ds) + """ + + # ensure time variable has been properly parsed as a datetime object + if not hasattr(ds[var], "dt"): + raise TypeError(f"The {var} variable passed has not been parsed as" + " a datetime object, so conversion to xtime string" + " will not work.\n\nTry using ds = xr.open_dataset" + "(..., use_cftime=True).") + + # xtime related attributes + attrs = {"units": "unitless", + "long_name": "model time, with format \'YYYY-MM-DD_HH:MM:SS\'"} + + # compute xtime dataarray from time variable passed + xtime = _datetime_2_xtime(ds, var=var) + + # update attributes of dataarray + xtime.attrs = attrs + + return xtime + + +def remap_variables(in_fp, out_fp, weights_fp, variables=None, logger=None): + """ + """ + + args = ["ncremap", + "-i", in_fp, + "-o", out_fp, + "-m", weights_fp] + + if variables and not isinstance(variables, list): + raise TypeError("`variables` kwarg must be a list of strings.") + elif variables: + args += ["-v"] + variables + + check_call(args, logger=logger) From b48c1e9d0bbe5e04fb30c391677bf4c302cb1781 Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Wed, 19 Jun 2024 10:15:14 -0700 Subject: [PATCH 07/12] Set nans in ocean forcing fields to zero. --- .../landice/tests/ismip6_GrIS_forcing/process_forcing.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py b/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py index d8f8b0cf45..4c991842ea 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py @@ -142,4 +142,13 @@ def process_variable(self, GCM, scenario, var, forcing_fp, start, end): remapped_ds[renamed_var] += mali_ds[renamed_var] + # ocean variable contain alot of nan's, replace with zeros + # ismip6Runoff : nan in catchments that are not marine terminating + # ismip6_2dThermalForcing : nan for all cells above sea level + if renamed_var in {"ismip6_2dThermalForcing", "ismip6Runoff"}: + # get the ocean dataarray of interest + da = remapped_ds[renamed_var] + # set nan values to zero in the parent dataset + remapped_ds[renamed_var] = xr.where(da.isnull(), 0, da) + return remapped_ds From ebe137e521b9a76bcbd3ee050a3502a7c572f3e3 Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Tue, 7 Oct 2025 14:26:12 -0700 Subject: [PATCH 08/12] uncommitted updates --- compass/landice/tests/ismip6_GrIS_forcing/file_finders.py | 2 +- .../tests/ismip6_GrIS_forcing/forcing_gen/forcing_gen.cfg | 5 +++++ compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py | 5 +++-- 3 files changed, 9 insertions(+), 3 deletions(-) diff --git a/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py b/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py index 9f41ca0eb3..92ead22d95 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py @@ -74,7 +74,7 @@ def __init__(self, archive_fp, version="v1", workdir="./"): # should do some checking here that all the fps work... self.dir2GCMs = dir2GCMs - def get_filename(self, GCM, scenario, variable, start=2015, end=2100): + def get_filename(self, GCM, scenario, variable, start=1950, end=2100): """ """ if GCM == "NorESM1-M": diff --git a/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/forcing_gen.cfg b/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/forcing_gen.cfg index becc05ce41..2dff24119d 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/forcing_gen.cfg +++ b/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/forcing_gen.cfg @@ -9,6 +9,11 @@ racmo_smb_fn = smb_rec.1958-2019.BN_RACMO2.3p2_FGRN055_GrIS.MM.nc climatology_start = 1960 climatology_end = 1989 +[TF_ref_climatology] +GCM = MIROC5 +climatology_start = 1990 +climatology_end = 2014 + [ISMIP6_GrIS_Forcing] # full path to MALI mesh forcing data will be regridded onto diff --git a/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py b/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py index 4c991842ea..deb62e6dfd 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py @@ -107,8 +107,9 @@ def process_variable(self, GCM, scenario, var, forcing_fp, start, end): drop_variables=vars_2_drop, use_cftime=True) - # mask of desired time indices - mask = remapped_ds.time.dt.year.isin(range(start, end)) + # create mask of desired time indices. Include forcing from year prior + # to requested start date since forcing in July but sims start in Jan + mask = remapped_ds.time.dt.year.isin(range(start - 1, end)) # drop the non-desired time indices from remapped dataset remapped_ds = remapped_ds.where(mask, drop=True) # add xtime variable based on `time` From 69ca57abd4e73dc44f57b81f6a9a2d5d46053e0f Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Wed, 8 Oct 2025 08:29:24 -0700 Subject: [PATCH 09/12] kwarg fix for pyremap and tmp fix for ESMF --- .../landice/tests/ismip6_GrIS_forcing/create_mapping_files.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py b/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py index 511983276e..b7d681d759 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py @@ -143,7 +143,7 @@ def make_forcing_descriptor_and_weights(self, # create pyremap `MeshDescriptor` obj for the forcing grid forcingDescriptor = ProjectionGridDescriptor.create( - projection=forcing_proj, x=x, y=y, meshName=forcing_file_fp) + projection=forcing_proj, x=x, y=y, mesh_name=forcing_file_fp) # write scrip file describing forcing grid forcingDescriptor.to_scrip(forcing_scrip_fp) @@ -153,7 +153,7 @@ def make_forcing_descriptor_and_weights(self, '-s', forcing_scrip_fp, '-d', mali_mesh_scrip_fp, '-w', weights_file_fp, - '--method', 'conserve', + '--method', 'bilinear', "--netcdf4", # "--no_log", "--dst_regional", "--src_regional", '--ignore_unmapped'] From 25f71db408f465c145b5891396d1c877588d5b05 Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Wed, 8 Oct 2025 08:30:50 -0700 Subject: [PATCH 10/12] Explicity do I/O with netcdf4 engine --- .../landice/tests/ismip6_GrIS_forcing/file_finders.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py b/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py index 92ead22d95..cd1bc629de 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py @@ -2,6 +2,7 @@ import os import xarray as xr +from xarray.coders import CFDatetimeCoder # create mapping dictionary of ISMIP6 variables to MALI variable names {"thermal_forcing": "ismip6_2dThermalForcing", @@ -82,6 +83,7 @@ def get_filename(self, GCM, scenario, variable, start=1950, end=2100): # get a list of all the yearly files within the period of intrest yearly_files = self.__find_yearly_files(GCM, scenario, variable, start, end) + # still need to make the output filename to write combined files to out_fn = f"MAR3.9_{GCM}_{scenario}_{variable}_{start}--{end}.nc" # relative to the workdir, which we've already checked if if existed @@ -129,9 +131,12 @@ def __find_yearly_files(self, GCM, scenario, variable, start, end): def __combine_files(self, files, out_fn): """ """ - ds = xr.open_mfdataset(files, concat_dim="time", combine="nested", + decoder = CFDatetimeCoder(use_cftime=True) + + ds = xr.open_mfdataset(files, decode_times=decoder, + concat_dim="time", combine="nested", data_vars='minimal', coords='minimal', compat="broadcast_equals", - combine_attrs="override") + combine_attrs="override", engine='netcdf4') - ds.to_netcdf(out_fn) + ds.to_netcdf(out_fn, engine='netcdf4') From 08c38df2b6f3135c619521d859c72ff3b69d8c2b Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Tue, 28 Oct 2025 14:42:59 -0700 Subject: [PATCH 11/12] Use shared functionailty for calling ncremap --- .../ref_smb_climatology.py | 30 ++++++------------- 1 file changed, 9 insertions(+), 21 deletions(-) diff --git a/compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py b/compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py index e9b6b718c1..e17c5db4f6 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py @@ -2,8 +2,8 @@ import xarray as xr from mpas_tools.io import write_netcdf -from mpas_tools.logging import check_call +from compass.landice.tests.ismip6_GrIS_forcing.utilities import remap_variables from compass.step import Step @@ -33,7 +33,6 @@ def setup(self): config = self.config smb_ref_section = config["smb_ref_climatology"] - # racmo_directory = smb_ref_section.get("racmo_directory") # this filename should probably just be hardcoded..... racmo_smb_fn = smb_ref_section.get("racmo_smb_fn") @@ -70,6 +69,9 @@ def setup(self): def run(self): """ """ + racmo_smb = self.racmo_smb + smb_ref_climatology = self.test_case.smb_ref_climatology + racmo_2_mali_weights = self.test_case.racmo_2_mali_weights # parse user specified parameters from the config config = self.config @@ -79,12 +81,11 @@ def run(self): clima_end = smb_ref_section.getint("climatology_end") # remap the gridded racmo data onto the mpas grid - self.remap_variable(self.racmo_smb, - self.test_case.smb_ref_climatology, - self.test_case.racmo_2_mali_weights) + remap_variables( + racmo_smb, smb_ref_climatology, racmo_2_mali_weights, ["SMB_rec"] + ) - ds = xr.open_dataset(self.test_case.smb_ref_climatology, - decode_times=False) + ds = xr.open_dataset(smb_ref_climatology, decode_times=False) # find indices of climatology start/end (TO DO: make more robust) s_idx = ((clima_start - 1958) * 12) - 1 @@ -105,17 +106,4 @@ def run(self): # expand sfcMassBal dimension to match what MALI expects ds["sfcMassBal"] = ds.sfcMassBal.expand_dims("Time") # write the file - write_netcdf(ds, self.test_case.smb_ref_climatology) - - def remap_variable(self, input_file, output_file, weights_file): - """ - """ - - # remap the forcing file onto the MALI mesh - args = ["ncremap", - "-i", input_file, - "-o", output_file, - "-m", weights_file, - "-v", "SMB_rec"] - - check_call(args, logger=self.logger) + write_netcdf(ds, smb_ref_climatology) From ef80161d423d560c52ed165dc3d609c6aae3cafa Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Fri, 7 Nov 2025 16:56:54 -0800 Subject: [PATCH 12/12] Formatting, docstrings, and minor fixes Made all variables snake case. Add docstings to all functions, methods, and classes. --- .../tests/ismip6_GrIS_forcing/__init__.py | 7 +- .../create_mapping_files.py | 84 ++++++-- .../tests/ismip6_GrIS_forcing/file_finders.py | 195 ++++++++++++++---- .../forcing_gen/__init__.py | 51 +++-- .../ismip6_GrIS_forcing/process_forcing.py | 94 +++++++-- .../ref_smb_climatology.py | 27 ++- .../tests/ismip6_GrIS_forcing/utilities.py | 60 ++++-- 7 files changed, 406 insertions(+), 112 deletions(-) diff --git a/compass/landice/tests/ismip6_GrIS_forcing/__init__.py b/compass/landice/tests/ismip6_GrIS_forcing/__init__.py index 232d8428fc..e978c60c91 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/__init__.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/__init__.py @@ -12,6 +12,8 @@ class Ismip6GrISForcing(TestGroup): """ def __init__(self, mpas_core): """ + Parameters + ---------- mpas_core : compass.landice.Landice the MPAS core that this test group belongs to """ @@ -30,7 +32,6 @@ def _read_and_validate_experiment_file(self): with open(f"{module_path}/experiments.yml", "r") as f: experiments = yaml.safe_load(f) - # should I filter through the experiments dictionary to make sure - # everything is valid.... - + # experiments dictionary is unverified... + # But the yaml file packages with compass shouldn't really be altered self.experiments = experiments diff --git a/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py b/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py index b7d681d759..62bedd1835 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/create_mapping_files.py @@ -13,10 +13,25 @@ class CreateMappingFiles(Step): """ + A step for creating mapping files for ISMIP6 GrIS forcing + + Attributes + ---------- + mali_mesh : str + File path of the MALI mesh ISMIP6 forcing will be remapped onto + + racmo_grid : str + File path to RACMO grid used to calculate the reference climatology """ def __init__(self, test_case): """ + Create the step + + Parameters + ---------- + test_case : compass.TestCase + The test case this step belongs to """ name = "create_mapping_files" subdir = "mapping_files" @@ -24,14 +39,18 @@ def __init__(self, test_case): super().__init__(test_case=test_case, name=name, subdir=subdir, cpus_per_task=128, min_cpus_per_task=1) - # initalize the attributes with empty values for now, attributes will - # be propely in `setup` method so user specified config options can - # be parsed + # attrs will set by the `setup` method, so user specified + # config options can be parsed self.mali_mesh = None self.racmo_grid = None def setup(self): """ + Parse config file for paths to input grids. + + Then add the paths to the scrip and weights files generated by this + step as attributes to the test case. This allows the mapping files + to be reused by other steps of the test case. """ # parse user specified parameters from the config @@ -95,6 +114,7 @@ def setup(self): def run(self): """ + Run this step of the test case """ # unpack the filepaths stored as attributes of the test_case @@ -112,29 +132,49 @@ def run(self): epsg_3413 = Proj("EPSG:3413") # make scrip file for racmo grid and mapping weights for racmo->mpas - self.make_forcing_descriptor_and_weights(self.racmo_grid, - racmo_gis_scrip, - mali_mesh_scrip, - epsg_3413, - racmo_2_mali_weights) + self.make_scrip_and_weights_files(self.racmo_grid, + racmo_gis_scrip, + mali_mesh_scrip, + epsg_3413, + racmo_2_mali_weights) # finding the right ismip6 forcing file is complicated, # so let's call a helper function to do that. ismip6_forcing_fp = self._find_ismip6_forcing_files() # make scrip file for ismip6 grid and mapping weights for ismip6->mpas - self.make_forcing_descriptor_and_weights(ismip6_forcing_fp, - ismip6_gis_scrip, - mali_mesh_scrip, - epsg_3413, - ismip6_2_mali_weights) - - def make_forcing_descriptor_and_weights(self, - forcing_file_fp, - forcing_scrip_fp, - mali_mesh_scrip_fp, - forcing_proj, - weights_file_fp): + self.make_scrip_and_weights_files(ismip6_forcing_fp, + ismip6_gis_scrip, + mali_mesh_scrip, + epsg_3413, + ismip6_2_mali_weights) + + def make_scrip_and_weights_files(self, + forcing_file_fp, + forcing_scrip_fp, + mali_mesh_scrip_fp, + forcing_proj, + weights_file_fp): + """ + Generate scrip and weight files needed for remapping + + Parameters + ---------- + forcing_file_fp: str + File path that constains mesh information of forcing data + + forcing_scrip_fp: str + File path to the scrip file describing the forcing data + + mali_mesh_scrip_fp: str + File path to scrip file describing the MPAS mesh + + forcing_proj: pyproj.Proj + Coordinate reference system of forcing data + + weights_file_fp: str + File path where the resulting weights file will be written + """ # open up forcing grid; get x and y coordinates as numpy arrays with xr.open_dataset(forcing_file_fp) as ds: @@ -184,6 +224,8 @@ def make_forcing_descriptor_and_weights(self, def _find_ismip6_forcing_files(self): """ + Helper function to find an ISMIP6 forcing file (not important which), + in order to create a mapping file to go from ISMIP6 --> MPAS grids """ # get a list of all projection experiments requesed. # (i.e. experiments names that use ISMIP6 forcing) @@ -202,6 +244,6 @@ def _find_ismip6_forcing_files(self): # get the filepath for any given forcing file (since all files are on # the same grid), in order to generate a scrip file for the grid - forcing_fp = self.test_case.findForcingFiles(GCM, scenario, variable) + forcing_fp = self.test_case.find_forcing_files(GCM, scenario, variable) return forcing_fp diff --git a/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py b/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py index cd1bc629de..c2f5263920 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/file_finders.py @@ -4,12 +4,22 @@ import xarray as xr from xarray.coders import CFDatetimeCoder -# create mapping dictionary of ISMIP6 variables to MALI variable names -{"thermal_forcing": "ismip6_2dThermalForcing", - "basin_runoff": "basin_runoff"} - def check_year(fn, start=2015, end=2100): + """ + Check is file is within desired period based on it's filename + + Parameters + ---------- + fn: str + Filename to check. + + start: int + Start of period + + end: str + End of period + """ # strip the file extension from filename and # convert the year at end of str to an int @@ -19,70 +29,170 @@ def check_year(fn, start=2015, end=2100): return start <= fn_year <= end -class oceanFileFinder: +class ISMIP6FileFinder: + """ + Base class for itterating over ISMIP6 forcing file archives - def __init__(self, archive_fp, version="v4"): + Attributes + ---------- + version: str + Version of the ISMIP6 ocean archive filename are from + + dir_w_GCMs: str + File path to directory containing the GCM data + """ + def __init__(self, version, dir_w_GCMs): + """ + Parameters + ---------- + version: str + Version of the ISMIP6 archive filename are from + + dir_w_GCMs: str + File path to directory containing the GCM data + """ self.version = version + self.dir_w_GCMs = self.check_file_exists(dir_w_GCMs) + + def get_filename(self): + """ + Return filepath to variable for requested GCM, scenario, and period + """ + raise NotImplementedError() + + def check_file_exists(self, fp): + """ + Ensure the filepath constructed actually exists + """ + + if os.path.exists(fp): + return fp + else: + msg = f"Cannot Find File: \n {fp} \n" + raise FileNotFoundError(msg) + + +class oceanFileFinder(ISMIP6FileFinder): + """ + Subclass for itterating ISMIP6 ocean archive + """ + + def __init__(self, archive_fp, version="v4"): + """ + Parameters + ---------- + version: str + Version of the ISMIP6 archive filename are from + + dir_w_GCMs: str + File path to directory containing the GCM data + """ + # file strucutre within Ocean_Forcing directory for navigating file_struct = f"Ocean_Forcing/Melt_Implementation/{version}" # file path to directory containing forcings files from each GCM - dir2GCMs = os.path.join(archive_fp, file_struct) + dir_w_GCMs = os.path.join(archive_fp, file_struct) - # should do some checking here that all the fps work... - self.dir2GCMs = dir2GCMs + super().__init__(version, dir_w_GCMs) def get_filename(self, GCM, scenario, variable): + """ + Return filepath to variable for requested GCM and scenario - # convert from variable name within NetCDF file, to variable name as it - # appears in the filename + Parameters + ---------- + GCM: str + General Circulation Model the forcing is derived from + + scenario: str + Emissions scenario + + var: str + Name of the variable to find file(s) for + """ + + # convert var name in NetCDF file, to var name in the filename if variable == "thermal_forcing": fn_var = "oceanThermalForcing" elif variable == "basin_runoff": fn_var = "basinRunoff" else: - raise ValueError(f"invalid varibale name: {variable}") + msg = f"invalid varibale name: {variable}" + raise ValueError(msg) # within filename scenario have for removed scenario_no_dot = "".join(scenario.split(".")) - fp = (f"{GCM.lower()}_{scenario}/" - f"MAR3.9_{GCM}_{scenario_no_dot}_{fn_var}_{self.version}.nc") + fp = ( + f"{GCM.lower()}_{scenario}/" + f"MAR3.9_{GCM}_{scenario_no_dot}_{fn_var}_{self.version}.nc" + ) - # check that file exists!!! - return os.path.join(self.dir2GCMs, fp) + out_fp = os.path.join(self.dir2GCMs, fp) + + return self.check_file_exists(out_fp) class atmosphereFileFinder: + """ + Subclass for itterating ISMIP6 atmosphere archive + """ def __init__(self, archive_fp, version="v1", workdir="./"): + """ + Parameters + ---------- + version: str + Version of the ISMIP6 archive filename are from - self.version = version - - if os.path.exists(workdir): - self.workdir = workdir - else: - print("gotta do error checking") + dir_w_GCMs: str + File path to directory containing the GCM data + """ # file strucutre within Atmosphere_Forcing directory for navigating file_struct = f"Atmosphere_Forcing/aSMB_observed/{version}" # file path to directory containing forcings files from each GCM - dir2GCMs = os.path.join(archive_fp, file_struct) + dir_w_GCMs = os.path.join(archive_fp, file_struct) - # should do some checking here that all the fps work... - self.dir2GCMs = dir2GCMs + super().__init__(version, dir_w_GCMs) + + if os.path.exists(workdir): + self.workdir = workdir + else: + print("gotta do error checking") def get_filename(self, GCM, scenario, variable, start=1950, end=2100): """ + Return filepath to variable for requested GCM, scenario and period + + Parameters + ---------- + GCM: str + General Circulation Model the forcing is derived from + + scenario: str + Emissions scenario + + var: str + Name of the variable to find file(s) for + + start: int + First year to process forcing for + + end: int + Final year to process forcing for """ if GCM == "NorESM1-M": GCM = "NorESM1" + # get a list of all the yearly files within the period of intrest - yearly_files = self.__find_yearly_files(GCM, scenario, variable, - start, end) + yearly_files = self._find_yearly_files( + GCM, scenario, variable, start, end + ) # still need to make the output filename to write combined files to out_fn = f"MAR3.9_{GCM}_{scenario}_{variable}_{start}--{end}.nc" @@ -99,19 +209,23 @@ def get_filename(self, GCM, scenario, variable, start=1950, end=2100): # only combine the files if they don't exist yet if not os.path.exists(out_fp): # with output filename, combine the files to a single - self.__combine_files(yearly_files, out_fp) + self._combine_files(yearly_files, out_fp) return out_fp - def __find_yearly_files(self, GCM, scenario, variable, start, end): + def _find_yearly_files(self, GCM, scenario, variable, start, end): """ + Each year of atmospheric forcing data is written to seperate file. + So, get a sorted list of files for the requested years """ # within filename scenario have for removed scenario_no_dot = "".join(scenario.split(".")) # create filename template to be globbed - fp = (f"{GCM}-{scenario_no_dot}/{variable}/" - f"{variable}_MARv3.9-yearly-{GCM}-{scenario_no_dot}-*.nc") + fp = ( + f"{GCM}-{scenario_no_dot}/{variable}/" + f"{variable}_MARv3.9-yearly-{GCM}-{scenario_no_dot}-*.nc" + ) # glob the files all_files = glob.glob(os.path.join(self.dir2GCMs, fp)) @@ -120,16 +234,27 @@ def __find_yearly_files(self, GCM, scenario, variable, start, end): # get the files within the desired period files = sorted(f for f in all_files if check_year(f, start, end)) + # make sure that each individual file exists + files = sorted(self.check_file_exists(f) for f in files) - # make sure the length of fns returned matches the length (in years) - # of the period + # number of fns returned must match length of the period (in years) period = (end - start) + 1 - assert len(files) == period + if len(files) != period: + msg = ( + f"Number of yearly files found for: \n" + f"\t GCM: {GCM}, Scenario: {scenario}, Period: {start}-{end}\n" + f"for variable: {variable} does not match the lenth of the" + f"period requested. Please investigate in: " + f"\t {self.dir_w_GCMs}" + ) + raise ValueError(msg) + return files - def __combine_files(self, files, out_fn): + def _combine_files(self, files, out_fn): """ + Concatenate multiple files along their time dimension """ decoder = CFDatetimeCoder(use_cftime=True) diff --git a/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py b/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py index 78a721df31..c4775f0135 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/forcing_gen/__init__.py @@ -34,14 +34,14 @@ class ForcingGen(TestCase): filepath to the `ESMF_RegirdWeightGen` generated mapping file experiments : dict - + Dictionary of ISMIP6 GrIS experiments to process data for """ def __init__(self, test_group): """ Parameters ---------- - test_group : compass.landice.tests... + test_group : compass.landice.tests.ismip6_GrIS_forcing The test group that this test case belongs to """ name = "forcing_gen" @@ -60,8 +60,8 @@ def __init__(self, test_group): # filename of the reference climatology self.smb_ref_climatology = None # place holder for file finders that will be initialized in `configure` - self.__atmFF = None - self.__ocnFF = None + self._atm_file_finder = None + self._ocn_file_finder = None # precusssory step that builds scrip file for mali mesh, # and generates a common weights file to be used in remapping @@ -75,10 +75,11 @@ def __init__(self, test_group): self.add_step(ProcessForcing(test_case=self)) def configure(self): - + """ + Set up the expriment dictionary, based on the requested experiments. + And initialize file finders based on archive filepath from config file. + """ config = self.config - # add ouputdir path to the remapping files - # get the list of requested experiments expr2run = config.getlist("ISMIP6_GrIS_Forcing", "experiments") # get the dictionary of experiments, as defined in the yaml file @@ -89,17 +90,41 @@ def configure(self): archive_fp = config.get("ISMIP6_GrIS_Forcing", "archive_fp") # initalize the oceanFileFinder - self.__ocnFF = oceanFileFinder(archive_fp) - self.__atmFF = atmosphereFileFinder(archive_fp) # , workdir=workdir) + self._atm_file_finder = oceanFileFinder(archive_fp) + self._ocn_file_finder = atmosphereFileFinder(archive_fp) - def findForcingFiles(self, GCM, scenario, variable): + def find_forcing_files( + self, GCM, scenario, variable, start=2015, end=2100 + ): """ + Parameters + ---------- + GCM: str + General Circulation Model the forcing is derived from + + scenario: str + Emissions scenario + + variable: str + Name of the variable to process + + forcing_fp: str + file path to read the forcing data from + + start: int + First year to process forcing for + + end: int + Final year to process forcing for """ if variable in ["basin_runoff", "thermal_forcing"]: - forcing_fp = self.__ocnFF.get_filename(GCM, scenario, variable) + file_finder = self._ocn_file_finder + # start and end have no effect for reading ocn forcing + kwargs = {} if variable in ["aSMB", "aST", "dSMBdz", "dSTdz"]: - forcing_fp = self.__atmFF.get_filename(GCM, scenario, variable) + file_finder = self._atm_file_finder + kwargs = {"start": start, "end": end} - return forcing_fp + return file_finder.get_filename(GCM, scenario, variable, **kwargs) diff --git a/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py b/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py index deb62e6dfd..2fd25af41a 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/process_forcing.py @@ -19,22 +19,28 @@ class ProcessForcing(Step): """ - + A step for remapping ISMIP6 GrIS forcing onto a MALI mesh. """ def __init__(self, test_case): """ + Create the step + + Parameters + ---------- + test_case : compass.TestCase + The test case this step belongs to """ name = "process_forcing" super().__init__(test_case=test_case, name=name, subdir=None, cpus_per_task=1, min_cpus_per_task=1) - # read and store the experiment dictionary - - # initalize the FileFinders are store them as dict? def run(self): + """ + Run this step of the test case + """ # list of variables to process atmosphere_vars = ["aSMB", "aST", "dSMBdz", "dSTdz"] @@ -42,7 +48,6 @@ def run(self): # loop over experiments passed via config file for expr_name, expr_dict in self.test_case.experiments.items(): - # print to logger which expriment we are on (i/n) GCM = expr_dict["GCM"] scenario = expr_dict["Scenario"] @@ -54,20 +59,46 @@ def run(self): continue atm_fn = f"gis_atm_forcing_{GCM}_{scenario}_{start}--{end}.nc" - self.process_variables(GCM, scenario, start, end, atmosphere_vars, - atm_fn) + self.process_variables( + GCM, scenario, start, end, atmosphere_vars, atm_fn + ) ocn_fn = f"gis_ocn_forcing_{GCM}_{scenario}_{start}--{end}.nc" - self.process_variables(GCM, scenario, start, end, ocean_vars, - ocn_fn) + self.process_variables( + GCM, scenario, start, end, ocean_vars, ocn_fn + ) + + def process_variables( + self, GCM, scenario, start, end, variables, output_fn + ): + """ + Parameters + ---------- + GCM: str + General Circulation Model the forcing is derived from + + scenario: str + Emissions scenario + + start: int + First year to process forcing for + + end: int + Final year to process forcing for - def process_variables(self, GCM, scenario, start, end, - variables, output_fn): + variables: list[str] + Variable names to process + + output_fn: str + File path where the processed forcing will be written + """ forcing_datastes = [] for var in variables: - var_fp = self.test_case.findForcingFiles(GCM, scenario, var) + var_fp = self.test_case.find_forcing_files( + GCM, scenario, var, start, end + ) ds = self.process_variable(GCM, scenario, var, var_fp, start, end) @@ -78,10 +109,29 @@ def process_variables(self, GCM, scenario, start, end, write_netcdf(forcing_ds, output_fn) def process_variable(self, GCM, scenario, var, forcing_fp, start, end): + """ + Parameters + ---------- + GCM: str + General Circulation Model the forcing is derived from + + scenario: str + Emissions scenario + + var: str + Name of the variable to process + + forcing_fp: str + file path to read the forcing data from + + start: int + First year to process forcing for + + end: int + Final year to process forcing for + """ - # config = self.config - # renamed_var = renaming_dict[var] # create the remapped file name, using original variable name @@ -103,9 +153,9 @@ def process_variable(self, GCM, scenario, var, forcing_fp, start, end): "mapping", "lat", "lon", "polar_stereographic"] # open the remapped file for post processing - remapped_ds = xr.open_dataset(remapped_fp, - drop_variables=vars_2_drop, - use_cftime=True) + remapped_ds = xr.open_dataset( + remapped_fp, drop_variables=vars_2_drop, use_cftime=True + ) # create mask of desired time indices. Include forcing from year prior # to requested start date since forcing in July but sims start in Jan @@ -116,17 +166,17 @@ def process_variable(self, GCM, scenario, var, forcing_fp, start, end): remapped_ds["xtime"] = add_xtime(remapped_ds, var="time") # rename the variable/dimensions to match MPAS/MALI conventions - remapped_ds = remapped_ds.rename({"ncol": "nCells", - "time": "Time", - var: renamed_var}) + remapped_ds = remapped_ds.rename( + {"ncol": "nCells", "time": "Time", var: renamed_var} + ) # drop the unneeded attributes from the dataset for _var in remapped_ds: remapped_ds[_var].attrs.pop("grid_mapping", None) remapped_ds[_var].attrs.pop("cell_measures", None) - # SMB is not a var in ISMIP6 forcing files, need to use `SMB_ref` - # and the processed `aSMB` to produce a `SMB` forcing + # surface mass balance and surface air temp. are provided as anomolies + # by ISMIP6. So remapped ISMIP6 fields must be added to a climatology if renamed_var == "sfcMassBal": # open the reference climatology file smb_ref = xr.open_dataset(self.test_case.smb_ref_climatology) diff --git a/compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py b/compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py index e17c5db4f6..d7269a9afd 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/ref_smb_climatology.py @@ -9,24 +9,38 @@ class SMBRefClimatology(Step): """ + A step to produce a surface mass balance climatology from RACMO data + + Attributes + ---------- + racmo_smb_fp: str + File path to RACMO SMB data proccessed by this step """ def __init__(self, test_case): """ + Create the step + + Parameters + ---------- + test_case : compass.TestCase + The test case this step belongs to """ name = "smb_ref_climatology" subdir = None super().__init__(test_case=test_case, name=name, subdir=subdir) - # initalize the attributes with empty values for now, attributes will - # be propely in `setup` method so user specified config options can - # be parsed + # attrs will set by the `setup` method, so user specified + # config options can be parsed self.racmo_smb_fp = None - self.smb_ref_climatology = None def setup(self): """ + Parse config file for path to RACMO data + + Then add path to the climatology produced by this step as an attribute + to the test case. Allowing climatology to be reused by other steps """ # parse user specified parameters from the config @@ -49,7 +63,7 @@ def setup(self): raise FileNotFoundError(f"{racmo_smb_fp} does not exist") # add the racmo smb as an attribute and as an input file - self.racmo_smb = racmo_smb_fp + self.racmo_smb_fp = racmo_smb_fp self.add_input_file(filename=racmo_smb_fp) # get the start and end dates for the climatological mean @@ -68,8 +82,9 @@ def setup(self): def run(self): """ + Run this step of the test case """ - racmo_smb = self.racmo_smb + racmo_smb = self.racmo_smb_fp smb_ref_climatology = self.test_case.smb_ref_climatology racmo_2_mali_weights = self.test_case.racmo_2_mali_weights diff --git a/compass/landice/tests/ismip6_GrIS_forcing/utilities.py b/compass/landice/tests/ismip6_GrIS_forcing/utilities.py index b553ef8ebb..6ecac007a7 100644 --- a/compass/landice/tests/ismip6_GrIS_forcing/utilities.py +++ b/compass/landice/tests/ismip6_GrIS_forcing/utilities.py @@ -2,27 +2,45 @@ from mpas_tools.logging import check_call -def _datetime_2_xtime(ds, var="Time", xtime_fmt="%Y-%m-%d_%H:%M:%S"): +def add_xtime(ds, var="Time"): + """ + Parse DateTime objects and convert to MPAS "xtime" format - return xr.apply_ufunc(lambda t: t.strftime(xtime_fmt).ljust(64), - ds[var], vectorize=True, output_dtypes=["S"]) + Parameters + ---------- + ds: xr.Dataset + Dataset containing the "var" data array of DateTime objects + var: str + Name of data array, composed of DateTime objects, to parse -def add_xtime(ds, var="Time"): - """ - ds["xtime"] = add_xtime(ds) + Returns + ------- + xtime: xr.DataArray + DateTime objects converted to 64 byte character strings """ + def _datetime_2_xtime(ds, var="Time", xtime_fmt="%Y-%m-%d_%H:%M:%S"): + + def f(t): + return t.strftime(xtime_fmt).ljust(64) + + return xr.apply_ufunc(f, ds[var], vectorize=True, output_dtypes=["S"]) + # ensure time variable has been properly parsed as a datetime object if not hasattr(ds[var], "dt"): - raise TypeError(f"The {var} variable passed has not been parsed as" - " a datetime object, so conversion to xtime string" - " will not work.\n\nTry using ds = xr.open_dataset" - "(..., use_cftime=True).") + msg = ( + f"The {var} variable passed has not been parsed as a datetime " + f"object, so conversion to xtime string will not work.\n\n" + f"Try using ds = xr.open_dataset(..., use_cftime=True)." + ) + raise TypeError(msg) # xtime related attributes - attrs = {"units": "unitless", - "long_name": "model time, with format \'YYYY-MM-DD_HH:MM:SS\'"} + attrs = { + "units": "unitless", + "long_name": "model time, with format \'YYYY-MM-DD_HH:MM:SS\'" + } # compute xtime dataarray from time variable passed xtime = _datetime_2_xtime(ds, var=var) @@ -35,6 +53,24 @@ def add_xtime(ds, var="Time"): def remap_variables(in_fp, out_fp, weights_fp, variables=None, logger=None): """ + Remap field using ncremp + + Parameters + ---------- + in_fp: str + File path to netcdf that has vars. on source grid + + out_fp: str + File path to netcdf where the vars. on destination grid will be written + + weight_fp: ste + File path to the weights file for remapping + + variables: list[str], optional + List of variable names (as strings) to be remapped. + + logger : logging.Logger + A logger for capturing the output from ncremap call """ args = ["ncremap",