Merge pull request #1100 from PCMDI/979_cloudfeedback_xcdat

Cloud feedbacks: updates to code backend

pcmdi_metrics/cloud_feedback/cloud_feedback_driver.py

@@ -113,7 +113,7 @@
         json.dump(filenames, f, sort_keys=True, indent=4)
 
 # calculate all feedback components and Klein et al (2013) error metrics:
-fbk_dict, obsc_fbk_dict, err_dict = CloudRadKernel(filenames)
+fbk_dict, obsc_fbk_dict, err_dict = CloudRadKernel(filenames, data_path)
 
 print("calc done")
 

pcmdi_metrics/cloud_feedback/lib/cal_CloudRadKernel_xr.py

@@ -7,20 +7,20 @@
 # to expert-assessed values from Sherwood et al (2020)
 # =============================================
 
+import os
 from datetime import date
 
 # IMPORT STUFF:
 import numpy as np
 import xarray as xr
 import xcdat as xc
-from global_land_mask import globe
+
+from pcmdi_metrics.utils import create_land_sea_mask
 
 # =============================================
 # define necessary information
 # =============================================
 
-datadir = "./data/"
-
 # Define a python dictionary containing the sections of the histogram to consider
 # These are the same as in Zelinka et al, GRL, 2016
 sections = ["ALL", "HI680", "LO680"]
@@ -35,7 +35,15 @@
 cutoff = int(len(binmids) / 2)
 # [:cutoff] = ascent; [cutoff:-1] = descent; [-1] = land
 
+# Define the grid
+LAT = np.arange(-89, 91, 2.0)
+LON = np.arange(1.25, 360, 2.5)
+lat_axis = xc.create_axis("lat", LAT)
+lon_axis = xc.create_axis("lon", LON)
+output_grid = xc.create_grid(x=lon_axis, y=lat_axis)
+
 # load land sea mask (90x144):
+"""
 lat = np.arange(-89, 90, 2.0)
 lon = np.arange(1.25, 360, 2.5) - 180
 lon_grid, lat_grid = np.meshgrid(lon, lat)
@@ -52,6 +60,9 @@
 land_mask = xr.where(land_mask, 1.0, 0.0)
 ocean_mask = xr.where(land_mask == 1.0, 0.0, 1.0)
 del E, W
+"""
+land_mask = create_land_sea_mask(output_grid)
+ocean_mask = xr.where(land_mask == 1.0, 0.0, 1.0)
 
 
 ###########################################################################
@@ -68,14 +79,14 @@ def get_amip_data(filename, var, lev=None):
     ).load()
     if lev:
         f = f.sel(time=tslice, plev=lev)
-        f = f.drop_vars(["plev", "plev_bnds"])
+        # f = f.drop_vars(["plev", "plev_bnds"])
     else:
         f = f.sel(time=tslice)
 
     # Compute climatological monthly means
     avg = f.temporal.climatology(var, freq="month", weighted=True)
     # Regrid to cloud kernel grid
-    output_grid = xc.create_grid(land_mask.lat.values, land_mask.lon.values)
+    # output_grid = xc.create_grid(land_mask.lat.values, land_mask.lon.values)
     output = avg.regridder.horizontal(
         var, output_grid, tool="xesmf", method="bilinear", extrap_method="inverse_dist"
     )
@@ -131,14 +142,14 @@ def get_model_data(filepath):
 
 
 ###########################################################################
-def get_CRK_data(filepath):
+def get_CRK_data(filepath, datadir):
     # Read in data, regrid
 
     # Load in regridded monthly mean climatologies from control and perturbed simulation
     print("get data")
     ctl, fut = get_model_data(filepath)
     print("get LW and SW kernel")
-    LWK, SWK = get_kernel_regrid(ctl)
+    LWK, SWK = get_kernel_regrid(ctl, datadir)
 
     # global mean and annual average delta tas
     avgdtas0 = fut["tas"] - ctl["tas"]
@@ -188,10 +199,12 @@ def get_CRK_data(filepath):
 
 
 ###########################################################################
-def get_kernel_regrid(ctl):
+def get_kernel_regrid(ctl, datadir):
     # Read in data and map kernels to lat/lon
 
-    f = xc.open_mfdataset(datadir + "cloud_kernels2.nc", decode_times=False)
+    f = xc.open_mfdataset(
+        os.path.join(datadir, "cloud_kernels2.nc"), decode_times=False
+    )
     f = f.rename({"mo": "time", "tau_midpt": "tau", "p_midpt": "plev"})
     f["time"] = ctl["time"].copy()
     f["tau"] = np.arange(7)
@@ -641,7 +654,8 @@ def xc_to_dataset(idata):
     idata = idata.to_dataset(name="data")
     if "height" in idata.coords:
         idata = idata.drop("height")
-    idata = idata.bounds.add_missing_bounds()
+    idata = idata.bounds.add_missing_bounds(axes=["X", "Y", "T"])
+    # idata = idata.bounds.add_missing_bounds()
     return idata
 
 
@@ -652,8 +666,9 @@ def monthly_anomalies(idata):
     usage:
     anom,avg = monthly_anomalies(data)
     """
-    idata = xc_to_dataset(idata)
     idata["time"].encoding["calendar"] = "noleap"
+    idata = xc_to_dataset(idata)
+    # idata["time"].encoding["calendar"] = "noleap"
     clim = idata.temporal.climatology("data", freq="month", weighted=True)
     anom = idata.temporal.departures("data", freq="month", weighted=True)
 
@@ -995,13 +1010,13 @@ def klein_metrics(obs_clisccp, gcm_clisccp, LWkern, SWkern, WTS):
 
 ###########################################################################
 def cal_Klein_metrics(
-    ctl_clisccp, LWK, SWK, area_wts, ctl_clisccp_wap, LWK_wap, SWK_wap
+    ctl_clisccp, LWK, SWK, area_wts, ctl_clisccp_wap, LWK_wap, SWK_wap, datadir
 ):
     ##########################################################
     ##### Load in ISCCP HGG clisccp climo annual cycle  ######
     ##########################################################
 
-    f = xc.open_dataset(datadir + "AC_clisccp_ISCCP_HGG_198301-200812.nc")
+    f = xc.open_dataset(os.path.join(datadir, "AC_clisccp_ISCCP_HGG_198301-200812.nc"))
     # f.history='Written by /work/zelinka1/scripts/load_ISCCP_HGG.py on feedback.llnl.gov'
     # f.comment='Monthly-resolved climatological annual cycle over 198301-200812'
     obs_clisccp_AC = f["AC_clisccp"]
@@ -1011,7 +1026,9 @@ def cal_Klein_metrics(
     obs_clisccp_AC["plev"] = np.arange(7)
     del f
 
-    f = xc.open_dataset(datadir + "AC_clisccp_wap_ISCCP_HGG_198301-200812.nc")
+    f = xc.open_dataset(
+        os.path.join(datadir, "AC_clisccp_wap_ISCCP_HGG_198301-200812.nc")
+    )
     # f.history='Written by /work/zelinka1/scripts/load_ISCCP_HGG.py on feedback.llnl.gov'
     # f.comment='Monthly-resolved climatological annual cycle over 198301-200812, in omega500 space'
     obs_clisccp_AC_wap = f["AC_clisccp_wap"]
@@ -1112,7 +1129,7 @@ def regional_breakdown(data, OCN, LND, area_wts, norm=False):
 
 
 ###########################################################################
-def CloudRadKernel(filepath):
+def CloudRadKernel(filepath, datadir):
     (
         ctl_clisccp,
         fut_clisccp,
@@ -1126,7 +1143,7 @@ def CloudRadKernel(filepath):
         SWK_wap,
         ctl_N,
         fut_N,
-    ) = get_CRK_data(filepath)
+    ) = get_CRK_data(filepath, datadir)
 
     OCN = ocean_mask.copy()
     LND = land_mask.copy()
@@ -1136,7 +1153,7 @@ def CloudRadKernel(filepath):
     ##############################################################################
     print("Compute Klein et al error metrics")
     KEM_dict = cal_Klein_metrics(
-        ctl_clisccp, LWK, SWK, area_wts, ctl_clisccp_wap, LWK_wap, SWK_wap
+        ctl_clisccp, LWK, SWK, area_wts, ctl_clisccp_wap, LWK_wap, SWK_wap, datadir
     )
     # [sec][region][sfc][ID], ID=E_TCA,E_ctpt,E_LW,E_SW,E_NET
 

pcmdi_metrics/cloud_feedback/param/my_param.py

@@ -5,6 +5,8 @@
 
 input_files_json = "./param/input_files.json"
 
+data_path = "./data"
+
 # Flag to compute ECS
 # True: compute ECS using abrupt-4xCO2 run
 # False: do not compute, instead rely on ECS value present in the json file (if it exists)

pcmdi_metrics/cloud_feedback/param/my_param_pcmdi_internal.py

@@ -11,7 +11,7 @@
 # True: compute ECS using abrupt-4xCO2 run
 # False: do not compute, instead rely on ECS value present in the json file (if it exists)
 get_ecs = True
-data_pathh = "/home/lee1043/git/pcmdi_metrics_20221013/pcmdi_metrics/pcmdi_metrics/cloud_feedback/data"
+data_path = "/home/lee1043/git/pcmdi_metrics_20221013/pcmdi_metrics/pcmdi_metrics/cloud_feedback/data"
 
 # Output directory path (directory will be generated if it does not exist yet.)
 xml_path = "/work/lee1043/cdat/pmp/cloud_feedback/xmls/"