Merge pull request #3259 from samsrabin/subset_data-lon-fixes

subset_data point: Fix --create-datm and Longitude TypeErrors

Author: samsrabin

.git-blame-ignore-revs

@@ -67,3 +67,4 @@ cdf40d265cc82775607a1bf25f5f527bacc97405
 3b7a2876933263f8986e4069f5d23bd45635756f
 3dd489af7ebe06566e2c6a1c7ade18550f1eb4ba
 742cfa606039ab89602fde5fef46458516f56fd4
+4ad46f46de7dde753b4653c15f05326f55116b73

doc/source/users_guide/running-single-points/generic-single-point-regional.rst

@@ -41,7 +41,7 @@ You can also have the script subset land-use data. See the help (``tools/site_an
 .. note::
    This script defaults to subsetting specific surface data, land-use timeseries, and the CRUJRA2024 DATM data. It can currently only be run as-is on Derecho. If you're not on Derecho, use ``--inputdata-dir`` to specify where the top level of your CESM input data is. 
 
-   Also, to subset GSWP3 instead of CRUJRA2024 DATM data, you currently need to hardwire ``datm_type = "datm_gswp3"`` (instead of the default ``"datm_crujra"``) in ``python/ctsm/subset_data.py``.
+   Using ``--create-datm`` with GSWP3 data is no longer supported; see `CTSM issue #3269 <https://github.com/ESCOMP/CTSM/issues/3269>`_.
 
 
 

python/ctsm/longitude.py

@@ -58,7 +58,7 @@ def _convert_lon_type_180_to_360(lon_in):
     return lon_out
 
 
-def _detect_lon_type(lon_in):
+def detect_lon_type(lon_in):
     """
     Detect longitude type of a given numeric. If lon_in contains more than one number (as in a list
     or Numpy array), this function will assume all members are of the same type if (a) there is at

python/ctsm/site_and_regional/regional_case.py

@@ -19,7 +19,7 @@
 from ctsm.utils import add_tag_to_filename
 from ctsm.utils import abort
 from ctsm.config_utils import check_lon1_lt_lon2
-from ctsm.longitude import Longitude, _detect_lon_type
+from ctsm.longitude import Longitude, detect_lon_type
 
 logger = logging.getLogger(__name__)
 
@@ -142,7 +142,7 @@ def _subset_lon_lat(self, x_dim, y_dim, f_in):
 
         # Detect longitude type (180 or 360) of input file, throwing a helpful error if it can't be
         # determined.
-        f_lon_type = _detect_lon_type(lon)
+        f_lon_type = detect_lon_type(lon)
         lon1_type = self.lon1.lon_type()
         lon2_type = self.lon2.lon_type()
         if lon1_type != lon2_type:

python/ctsm/site_and_regional/single_point_case.py

@@ -15,17 +15,14 @@
 # -- import local classes for this script
 from ctsm.site_and_regional.base_case import BaseCase, USRDAT_DIR, DatmFiles
 from ctsm.utils import add_tag_to_filename, ensure_iterable
+from ctsm.longitude import detect_lon_type
 
 logger = logging.getLogger(__name__)
 
 NAT_PFT = 15  # natural pfts
 NUM_PFT = 17  # for runs with generic crops
 MAX_PFT = 78  # for runs with explicit crops
 
-# -- constants to represent months of year
-FIRST_MONTH = 1
-LAST_MONTH = 12
-
 
 class SinglePointCase(BaseCase):
     """
@@ -151,6 +148,26 @@ def __init__(
         # self.check_nonveg()
         self.check_pct_pft()
 
+    def convert_plon_to_filetype_if_needed(self, lon_da):
+        """
+        Check that point and input file longitude types are equal. If not, convert point to match
+        file.
+        """
+        plon_in = self.plon
+        f_lon_type = detect_lon_type(lon_da)
+        plon_type = plon_in.lon_type()
+        if f_lon_type == plon_type:
+            plon_out = plon_in.get(plon_type)
+        else:
+            plon_orig = plon_in.get(plon_type)
+            plon_out = plon_in.get(f_lon_type)
+            if plon_orig != plon_out:
+                print(
+                    f"Converted plon from type {plon_type} (value {plon_orig}) "
+                    f"to type {f_lon_type} (value {plon_out})"
+                )
+        return plon_out
+
     def create_tag(self):
         """
         Create a tag for single point which is the site name
@@ -363,8 +380,11 @@ def create_landuse_at_point(self, indir, file, user_mods_dir):
         # create 1d coordinate variables to enable sel() method
         f_in = self.create_1d_coord(fluse_in, "LONGXY", "LATIXY", "lsmlon", "lsmlat")
 
+        # get point longitude, converting to match file type if needed
+        plon_float = self.convert_plon_to_filetype_if_needed(f_in["lsmlon"])
+
         # extract gridcell closest to plon/plat
-        f_out = f_in.sel(lsmlon=self.plon, lsmlat=self.plat, method="nearest")
+        f_out = f_in.sel(lsmlon=plon_float, lsmlat=self.plat, method="nearest")
 
         # expand dimensions
         f_out = f_out.expand_dims(["lsmlat", "lsmlon"])
@@ -498,8 +518,11 @@ def create_surfdata_at_point(self, indir, file, user_mods_dir, specify_fsurf_out
         # create 1d coordinate variables to enable sel() method
         f_in = self.create_1d_coord(fsurf_in, "LONGXY", "LATIXY", "lsmlon", "lsmlat")
 
+        # get point longitude, converting to match file type if needed
+        plon_float = self.convert_plon_to_filetype_if_needed(f_in["lsmlon"])
+
         # extract gridcell closest to plon/plat
-        f_tmp = f_in.sel(lsmlon=self.plon, lsmlat=self.plat, method="nearest")
+        f_tmp = f_in.sel(lsmlon=plon_float, lsmlat=self.plat, method="nearest")
 
         # expand dimensions
         f_tmp = f_tmp.expand_dims(["lsmlat", "lsmlon"]).copy(deep=True)
@@ -525,10 +548,10 @@ def create_surfdata_at_point(self, indir, file, user_mods_dir, specify_fsurf_out
         # update lsmlat and lsmlon to match site specific instead of the nearest point
         # we do this so that if we create user_mods the PTS_LON and PTS_LAT in CIME match
         # the surface data coordinates - which is required
-        f_out["lsmlon"] = np.atleast_1d(self.plon)
+        f_out["lsmlon"] = np.atleast_1d(plon_float)
         f_out["lsmlat"] = np.atleast_1d(self.plat)
         f_out["LATIXY"][:, :] = self.plat
-        f_out["LONGXY"][:, :] = self.plon
+        f_out["LONGXY"][:, :] = plon_float
 
         # update attributes
         self.update_metadata(f_out)
@@ -568,8 +591,11 @@ def create_datmdomain_at_point(self, datm_tuple: DatmFiles):
         # create 1d coordinate variables to enable sel() method
         f_in = self.create_1d_coord(fdatmdomain_in, "xc", "yc", "ni", "nj")
 
+        # get point longitude, converting to match file type if needed
+        plon_float = self.convert_plon_to_filetype_if_needed(f_in["lon"])
+
         # extract gridcell closest to plon/plat
-        f_out = f_in.sel(ni=self.plon, nj=self.plat, method="nearest")
+        f_out = f_in.sel(ni=plon_float, nj=self.plat, method="nearest")
 
         # expand dimensions
         f_out = f_out.expand_dims(["nj", "ni"])
@@ -591,14 +617,17 @@ def extract_datm_at(self, file_in, file_out):
         # create 1d coordinate variables to enable sel() method
         f_in = self.create_1d_coord(file_in, "LONGXY", "LATIXY", "lon", "lat")
 
+        # get point longitude, converting to match file type if needed
+        plon_float = self.convert_plon_to_filetype_if_needed(f_in["lon"])
+
         # extract gridcell closest to plon/plat
-        f_out = f_in.sel(lon=self.plon, lat=self.plat, method="nearest")
+        f_out = f_in.sel(lon=plon_float, lat=self.plat, method="nearest")
 
         # expand dimensions
         f_out = f_out.expand_dims(["lat", "lon"])
 
         # specify dimension order
-        f_out = f_out.transpose("scalar", "time", "lat", "lon")
+        f_out = f_out.transpose("time", "lat", "lon")
 
         # update attributes
         self.update_metadata(f_out)
@@ -653,46 +682,36 @@ def create_datm_at_point(self, datm_tuple: DatmFiles, datm_syr, datm_eyr, datm_s
         tpqwfiles = []
         for year in range(datm_syr, datm_eyr + 1):
             ystr = str(year)
-            for month in range(FIRST_MONTH, LAST_MONTH + 1):
-                mstr = str(month)
-                if month < 10:
-                    mstr = "0" + mstr
-
-                dtag = ystr + "-" + mstr
 
-                fsolar = os.path.join(
-                    datm_tuple.indir,
-                    datm_tuple.dir_solar,
-                    "{}{}.nc".format(datm_tuple.tag_solar, dtag),
-                )
-                fsolar2 = "{}{}.{}.nc".format(datm_tuple.tag_solar, self.tag, dtag)
-                fprecip = os.path.join(
-                    datm_tuple.indir,
-                    datm_tuple.dir_prec,
-                    "{}{}.nc".format(datm_tuple.tag_prec, dtag),
-                )
-                fprecip2 = "{}{}.{}.nc".format(datm_tuple.tag_prec, self.tag, dtag)
-                ftpqw = os.path.join(
-                    datm_tuple.indir,
-                    datm_tuple.dir_tpqw,
-                    "{}{}.nc".format(datm_tuple.tag_tpqw, dtag),
-                )
-                ftpqw2 = "{}{}.{}.nc".format(datm_tuple.tag_tpqw, self.tag, dtag)
-
-                outdir = os.path.join(self.out_dir, datm_tuple.outdir)
-                infile += [fsolar, fprecip, ftpqw]
-                outfile += [
-                    os.path.join(outdir, fsolar2),
-                    os.path.join(outdir, fprecip2),
-                    os.path.join(outdir, ftpqw2),
-                ]
-                solarfiles.append(
-                    os.path.join("${}".format(USRDAT_DIR), datm_tuple.outdir, fsolar2)
-                )
-                precfiles.append(
-                    os.path.join("${}".format(USRDAT_DIR), datm_tuple.outdir, fprecip2)
-                )
-                tpqwfiles.append(os.path.join("${}".format(USRDAT_DIR), datm_tuple.outdir, ftpqw2))
+            fsolar = os.path.join(
+                datm_tuple.indir,
+                datm_tuple.dir_solar,
+                "{}{}.nc".format(datm_tuple.tag_solar, ystr),
+            )
+            fsolar2 = "{}{}.{}.nc".format(datm_tuple.tag_solar, self.tag, ystr)
+            fprecip = os.path.join(
+                datm_tuple.indir,
+                datm_tuple.dir_prec,
+                "{}{}.nc".format(datm_tuple.tag_prec, ystr),
+            )
+            fprecip2 = "{}{}.{}.nc".format(datm_tuple.tag_prec, self.tag, ystr)
+            ftpqw = os.path.join(
+                datm_tuple.indir,
+                datm_tuple.dir_tpqw,
+                "{}{}.nc".format(datm_tuple.tag_tpqw, ystr),
+            )
+            ftpqw2 = "{}{}.{}.nc".format(datm_tuple.tag_tpqw, self.tag, ystr)
+
+            outdir = os.path.join(self.out_dir, datm_tuple.outdir)
+            infile += [fsolar, fprecip, ftpqw]
+            outfile += [
+                os.path.join(outdir, fsolar2),
+                os.path.join(outdir, fprecip2),
+                os.path.join(outdir, ftpqw2),
+            ]
+            solarfiles.append(os.path.join("${}".format(USRDAT_DIR), datm_tuple.outdir, fsolar2))
+            precfiles.append(os.path.join("${}".format(USRDAT_DIR), datm_tuple.outdir, fprecip2))
+            tpqwfiles.append(os.path.join("${}".format(USRDAT_DIR), datm_tuple.outdir, ftpqw2))
 
         for idx, out_f in enumerate(outfile):
             logger.debug(out_f)

python/ctsm/subset_data.py

@@ -69,7 +69,7 @@
 from ctsm.path_utils import path_to_ctsm_root
 from ctsm.utils import abort
 from ctsm.config_utils import check_lon1_lt_lon2
-from ctsm.longitude import Longitude, _detect_lon_type
+from ctsm.longitude import Longitude, detect_lon_type
 
 # -- import ctsm logging flags
 from ctsm.ctsm_logging import (
@@ -629,32 +629,41 @@ def setup_files(args, defaults, cesmroot):
     file_dict = {"main_dir": clmforcingindir}
 
     # DATM data
-    # TODO Issue #2960: Make datm_type a user option at the command
-    # line. For reference, this option affects three .cfg files:
-    #      tools/site_and_regional/default_data_1850.cfg
-    #      tools/site_and_regional/default_data_2000.cfg
-    #      python/ctsm/test/testinputs/default_data.cfg
+    # To find the affected files, from the top level of ctsm, do:
+    #     grep "\[datm\]" $(find . -type f -name "*cfg")
     if args.create_datm:
-        datm_type = "datm_crujra"  # also available: datm_type = "datm_gswp3"
+        datm_cfg_section = "datm"
+
+        # Issue #3269: Changes in PR #3259 mean that --create-datm won't work with GSWP3
+        settings_to_check_for_gswp3 = ["solartag", "prectag", "tpqwtag"]
+        for setting in settings_to_check_for_gswp3:
+            value = defaults.get(datm_cfg_section, setting)
+            if "gswp3" in value.lower():
+                msg = (
+                    "--create-datm is no longer supported for GSWP3 data; "
+                    "see https://github.com/ESCOMP/CTSM/issues/3269"
+                )
+                raise NotImplementedError(msg)
+
         dir_output_datm = "datmdata"
-        dir_input_datm = os.path.join(clmforcingindir, defaults.get(datm_type, "dir"))
+        dir_input_datm = os.path.join(clmforcingindir, defaults.get(datm_cfg_section, "dir"))
         if not os.path.isdir(os.path.join(args.out_dir, dir_output_datm)):
             os.mkdir(os.path.join(args.out_dir, dir_output_datm))
         logger.info("dir_input_datm : %s", dir_input_datm)
         logger.info("dir_output_datm: %s", os.path.join(args.out_dir, dir_output_datm))
         file_dict["datm_tuple"] = DatmFiles(
             dir_input_datm,
             dir_output_datm,
-            defaults.get(datm_type, "domain"),
-            defaults.get(datm_type, "solardir"),
-            defaults.get(datm_type, "precdir"),
-            defaults.get(datm_type, "tpqwdir"),
-            defaults.get(datm_type, "solartag"),
-            defaults.get(datm_type, "prectag"),
-            defaults.get(datm_type, "tpqwtag"),
-            defaults.get(datm_type, "solarname"),
-            defaults.get(datm_type, "precname"),
-            defaults.get(datm_type, "tpqwname"),
+            defaults.get(datm_cfg_section, "domain"),
+            defaults.get(datm_cfg_section, "solardir"),
+            defaults.get(datm_cfg_section, "precdir"),
+            defaults.get(datm_cfg_section, "tpqwdir"),
+            defaults.get(datm_cfg_section, "solartag"),
+            defaults.get(datm_cfg_section, "prectag"),
+            defaults.get(datm_cfg_section, "tpqwtag"),
+            defaults.get(datm_cfg_section, "solarname"),
+            defaults.get(datm_cfg_section, "precname"),
+            defaults.get(datm_cfg_section, "tpqwname"),
         )
 
     # if the crop flag is on - we need to use a different land use and surface data file
@@ -833,10 +842,10 @@ def process_args(args):
     if any(lon_arg_values):
         if args.lon_type is None:
             if hasattr(args, "plon"):
-                args.lon_type = _detect_lon_type(args.plon)
+                args.lon_type = detect_lon_type(args.plon)
             else:
-                lon1_type = _detect_lon_type(args.lon1)
-                lon2_type = _detect_lon_type(args.lon2)
+                lon1_type = detect_lon_type(args.lon1)
+                lon2_type = detect_lon_type(args.lon2)
                 if lon1_type != lon2_type:
                     raise argparse.ArgumentTypeError(
                         "--lon1 and --lon2 seem to be of different types"