debug and updates

.pre-commit-config.yaml

@@ -25,7 +25,7 @@ repos:
  - id: black
 
  - repo: https://github.com/timothycrosley/isort
- rev: 5.12.0
+ rev: 5.13.2
  hooks:
  - id: isort
  args: ["--honor-noqa"]
@@ -34,7 +34,7 @@ repos:
  # Python linting
  # =======================
  - repo: https://github.com/pycqa/flake8
- rev: 6.0.0
+ rev: 7.0.0
  hooks:
  - id: flake8
  args: ["--config=setup.cfg"]

pcmdi_metrics/io/__init__.py

@@ -7,6 +7,7 @@
  da_to_ds,
  get_axis_list,
  get_data_list,
+ get_grid,
  get_latitude_bounds_key,
  get_latitude_key,
  get_latitude,
@@ -21,4 +22,4 @@
  get_time_key,
  select_subset,
 )
-from .default_regions_define import load_regions_specs, region_subset # noqa
+from .regions import load_regions_specs, region_subset # noqa

pcmdi_metrics/io/default_regions_define.py → pcmdi_metrics/io/regions.py

@@ -3,7 +3,7 @@
 import xarray as xr
 import xcdat as xc
 
-from pcmdi_metrics.io import da_to_ds
+from pcmdi_metrics.io import da_to_ds, get_longitude, select_subset
 
 
 def load_regions_specs() -> dict:
@@ -76,63 +76,67 @@ def load_regions_specs() -> dict:
 
 
 def region_subset(
- ds: Union[xr.Dataset, xr.DataArray], region: str, regions_specs: dict = None
+ ds: Union[xr.Dataset, xr.DataArray],
+ region: str,
+ data_var: str = "variable",
+ regions_specs: dict = None,
 ) -> Union[xr.Dataset, xr.DataArray]:
- """
- ds: xarray.Dataset
- regions_specs: dict
- region: string
+ """_summary_
+
+ Parameters
+ ----------
+ ds : Union[xr.Dataset, xr.DataArray]
+ _description_
+ region : str
+ _description_
+ data_var : str, optional
+ _description_, by default None
+ regions_specs : dict, optional
+ _description_, by default None
+
+ Returns
+ -------
+ Union[xr.Dataset, xr.DataArray]
+ _description_
  """
  if isinstance(ds, xr.DataArray):
  is_dataArray = True
- varname = "variable"
- ds = da_to_ds(ds, varname)
+ ds = da_to_ds(ds, data_var)
  else:
  is_dataArray = False
 
  if regions_specs is None:
  regions_specs = load_regions_specs()
 
- if "domain" in list(regions_specs[region].keys()):
- if "latitude" in list(regions_specs[region]["domain"].keys()):
+ if "domain" in regions_specs[region]:
+ if "latitude" in regions_specs[region]["domain"]:
  lat0 = regions_specs[region]["domain"]["latitude"][0]
  lat1 = regions_specs[region]["domain"]["latitude"][1]
  # proceed subset
- if "latitude" in (ds.coords.dims):
- ds = ds.sel(latitude=slice(lat0, lat1))
- elif "lat" in (ds.coords.dims):
- ds = ds.sel(lat=slice(lat0, lat1))
+ ds = select_subset(ds, lat=(lat0, lat1))
 
- if "longitude" in list(regions_specs[region]["domain"].keys()):
+ if "longitude" in regions_specs[region]["domain"]:
  lon0 = regions_specs[region]["domain"]["longitude"][0]
  lon1 = regions_specs[region]["domain"]["longitude"][1]
 
  # check original dataset longitude range
- if "longitude" in (ds.coords.dims):
- lon_min = ds.longitude.min()
- lon_max = ds.longitude.max()
- elif "lon" in (ds.coords.dims):
- lon_min = ds.lon.min()
- lon_max = ds.lon.max()
-
- # longitude range swap if needed
- if (
- min(lon0, lon1) < 0
- ): # when subset region lon is defined in (-180, 180) range
- if (
- min(lon_min, lon_max) < 0
- ): # if original data lon range is (-180, 180) no treatment needed
+ lon_min = get_longitude(ds).min().values.item()
+ lon_max = get_longitude(ds).max().values.item()
+
+ # Check if longitude range swap is needed
+ if min(lon0, lon1) < 0:
+ # when subset region lon is defined in (-180, 180) range
+ if min(lon_min, lon_max) < 0:
+ # if original data lon range is (-180, 180), no treatment needed
  pass
- else: # if original data lon range is (0, 360), convert swap lon
+ else:
+ # if original data lon range is (0, 360), convert and swap lon
  ds = xc.swap_lon_axis(ds, to=(-180, 180))
 
  # proceed subset
- if "longitude" in (ds.coords.dims):
- ds = ds.sel(longitude=slice(lon0, lon1))
- elif "lon" in (ds.coords.dims):
- ds = ds.sel(lon=slice(lon0, lon1))
+ ds = select_subset(ds, lon=(lon0, lon1))
 
  if is_dataArray:
- return ds["variable"]
+ return ds[data_var]
  else:
  return ds

pcmdi_metrics/io/xcdat_dataset_io.py

@@ -197,3 +197,25 @@ def da_to_ds(d: Union[xr.Dataset, xr.DataArray], var: str = "variable") -> xr.Da
  raise TypeError(
  "Input must be an instance of either xarrary.DataArray or xarrary.Dataset"
  )
+
+
+def get_grid(
+ ds: xr.Dataset,
+) -> xr.Dataset:
+ """Get grid information
+
+ Parameters
+ ----------
+ ds : xr.Dataset
+ xarray dataset to extract grid information that has latitude, longitude, and their bounds included
+
+ Returns
+ -------
+ xr.Dataset
+ xarray dataset with grid information
+ """
+ lat_key = get_latitude_key(ds)
+ lon_key = get_longitude_key(ds)
+ lat_bnds_key = get_latitude_bounds_key(ds)
+ lon_bnds_key = get_longitude_bounds_key(ds)
+ return ds[[lat_key, lon_key, lat_bnds_key, lon_bnds_key]]

pcmdi_metrics/stats/compute_statistics_dataset.py

@@ -6,7 +6,7 @@
 import xcdat as xc
 
 
-def _check_data_convert_to_ds_if_needed(
+def da_to_ds(
  d: Union[xr.Dataset, xr.DataArray], var: str = "variable"
 ):
  if isinstance(d, xr.Dataset):
@@ -29,8 +29,8 @@ def annual_mean(dm, do, var="variable"):
  "Comments": "Assumes input are 12 months climatology",
  }
 
- dm = _check_data_convert_to_ds_if_needed(dm, var)
- do = _check_data_convert_to_ds_if_needed(do, var)
+ dm = da_to_ds(dm, var)
+ do = da_to_ds(do, var)
 
  dm_am = dm.temporal.average(var)
  do_am = do.temporal.average(var)
@@ -84,8 +84,8 @@ def bias_xy(dm, do, var="variable", weights=None):
  "Contact": "[email protected]",
  }
 
- dm = _check_data_convert_to_ds_if_needed(dm, var)
- do = _check_data_convert_to_ds_if_needed(do, var)
+ dm = da_to_ds(dm, var)
+ do = da_to_ds(do, var)
 
  dif = dm[var] - do[var]
  if weights is None:
@@ -104,8 +104,8 @@ def bias_xyt(dm, do, var="variable"):
  "Contact": "[email protected]",
  }
 
- dm = _check_data_convert_to_ds_if_needed(dm, var)
- do = _check_data_convert_to_ds_if_needed(do, var)
+ dm = da_to_ds(dm, var)
+ do = da_to_ds(do, var)
 
  ds = dm.copy(deep=True)
  ds["dif"] = dm[var] - do[var]
@@ -124,8 +124,8 @@ def cor_xy(dm, do, var="variable", weights=None):
  "Contact": "[email protected]",
  }
 
- dm = _check_data_convert_to_ds_if_needed(dm, var)
- do = _check_data_convert_to_ds_if_needed(do, var)
+ dm = da_to_ds(dm, var)
+ do = da_to_ds(do, var)
 
  if weights is None:
  weights = dm.spatial.get_weights(axis=["X", "Y"])
@@ -155,7 +155,7 @@ def mean_xy(d, var="variable", weights=None):
  "Contact": "[email protected]",
  }
 
- d = _check_data_convert_to_ds_if_needed(d, var)
+ d = da_to_ds(d, var)
 
  lat_key = xc.axis.get_dim_keys(d, axis="Y")
  lon_key = xc.axis.get_dim_keys(d, axis="X")
@@ -176,8 +176,8 @@ def meanabs_xy(dm, do, var="variable", weights=None):
  "Contact": "[email protected]",
  }
 
- dm = _check_data_convert_to_ds_if_needed(dm, var)
- do = _check_data_convert_to_ds_if_needed(do, var)
+ dm = da_to_ds(dm, var)
+ do = da_to_ds(do, var)
 
  if weights is None:
  weights = dm.spatial.get_weights(axis=["X", "Y"])
@@ -197,8 +197,8 @@ def meanabs_xyt(dm, do, var="variable"):
  "Contact": "[email protected]",
  }
 
- dm = _check_data_convert_to_ds_if_needed(dm, var)
- do = _check_data_convert_to_ds_if_needed(do, var)
+ dm = da_to_ds(dm, var)
+ do = da_to_ds(do, var)
 
  ds = dm.copy(deep=True)
  ds["absdif"] = abs(dm[var] - do[var])
@@ -219,8 +219,8 @@ def rms_0(dm, do, var="variable", weighted=True):
  "Contact": "[email protected]",
  }
 
- dm = _check_data_convert_to_ds_if_needed(dm, var)
- do = _check_data_convert_to_ds_if_needed(do, var)
+ dm = da_to_ds(dm, var)
+ do = da_to_ds(do, var)
 
  dif_square = (dm[var] - do[var]) ** 2
  if weighted:
@@ -240,8 +240,8 @@ def rms_xy(dm, do, var="variable", weights=None):
  "Contact": "[email protected]",
  }
 
- dm = _check_data_convert_to_ds_if_needed(dm, var)
- do = _check_data_convert_to_ds_if_needed(do, var)
+ dm = da_to_ds(dm, var)
+ do = da_to_ds(do, var)
 
  dif_square = (dm[var] - do[var]) ** 2
  if weights is None:
@@ -259,8 +259,8 @@ def rms_xyt(dm, do, var="variable"):
  "Contact": "[email protected]",
  }
 
- dm = _check_data_convert_to_ds_if_needed(dm, var)
- do = _check_data_convert_to_ds_if_needed(do, var)
+ dm = da_to_ds(dm, var)
+ do = da_to_ds(do, var)
 
  ds = dm.copy(deep=True)
  ds["diff_square"] = (dm[var] - do[var]) ** 2
@@ -280,8 +280,8 @@ def rmsc_xy(dm, do, var="variable", weights=None, NormalizeByOwnSTDV=False):
  "Contact": "[email protected]",
  }
 
- dm = _check_data_convert_to_ds_if_needed(dm, var)
- do = _check_data_convert_to_ds_if_needed(do, var)
+ dm = da_to_ds(dm, var)
+ do = da_to_ds(do, var)
 
  if weights is None:
  weights = dm.spatial.get_weights(axis=["X", "Y"])
@@ -310,7 +310,7 @@ def std_xy(ds, var="variable", weights=None):
  "Contact": "[email protected]",
  }
 
- ds = _check_data_convert_to_ds_if_needed(ds, var)
+ ds = da_to_ds(ds, var)
 
  if weights is None:
  weights = ds.spatial.get_weights(axis=["X", "Y"])
@@ -334,7 +334,7 @@ def std_xyt(d, var="variable"):
  "Contact": "[email protected]",
  }
  ds = d.copy(deep=True)
- ds = _check_data_convert_to_ds_if_needed(ds, var)
+ ds = da_to_ds(ds, var)
  average = d.spatial.average(var, axis=["X", "Y"]).temporal.average(var)[var]
  ds["anomaly"] = (d[var] - average) ** 2
  variance = (
@@ -353,8 +353,8 @@ def zonal_mean(dm, do, var="variable"):
  "Contact": "[email protected]",
  "Comments": "",
  }
- dm = _check_data_convert_to_ds_if_needed(dm, var)
- do = _check_data_convert_to_ds_if_needed(do, var)
+ dm = da_to_ds(dm, var)
+ do = da_to_ds(do, var)
 
  dm_zm = dm.spatial.average(var, axis=["X"])
  do_zm = do.spatial.average(var, axis=["X"])

pcmdi_metrics/variability_mode/lib/__init__.py

@@ -21,6 +21,7 @@
 )
 from .landmask import data_land_mask_out, estimate_landmask # noqa
 from .lib_variability_mode import ( # noqa
+ check_start_end_year,
  debug_print,
  get_domain_range,
  read_data_in,

pcmdi_metrics/variability_mode/lib/calc_stat.py

@@ -1,11 +1,12 @@
 from time import gmtime, strftime
 
-from pcmdi_metrics.io import region_subset
+from pcmdi_metrics.io import get_grid, region_subset
 from pcmdi_metrics.stats import bias_xy as calcBias
 from pcmdi_metrics.stats import cor_xy as calcSCOR
 from pcmdi_metrics.stats import mean_xy
 from pcmdi_metrics.stats import rms_xy as calcRMS
 from pcmdi_metrics.stats import rmsc_xy as calcRMSc
+from pcmdi_metrics.utils import regrid
 
 
 def calc_stats_save_dict(
@@ -60,12 +61,13 @@ def calc_stats_save_dict(
  # . . . . . . . . . . . . . . . . . . . . . . . . .
  if obs_compare:
  if method in ["eof", "cbf"]:
- ref_grid_global = eof_lr_obs.getGrid()
+ ref_grid_global = get_grid(eof_lr_obs)
  # Regrid (interpolation, model grid to ref grid)
  debug_print("regrid (global) start", debug)
- eof_model_global = eof_lr.regrid(
- ref_grid_global, regridTool="regrid2", mkCyclic=True
- )
+ # eof_model_global = eof_lr.regrid(eof_lr,
+ # ref_grid_global, regridTool="regrid2", mkCyclic=True
+ # )
+ eof_model_global = regrid(eof_lr, ref_grid_global)
  debug_print("regrid end", debug)
  # Extract subdomain
  # eof_model = eof_model_global(region_subdomain)