Merge pull request #17 from LSSTDESC/issue/16/copy_obscond

Copying observing_condition_degrader from rail_attic

Author: aimalz

src/rail/creation/degradation/observing_condition_degrader.py

@@ -45,7 +45,7 @@ class ObsCondition(Degrader):
         uses the same arguments as LSSTErrorModel (from PhotErr).
         The following arguments, if supplied, may contain either
         a single number (as in the case of LSSTErrorModel), or a path:
-        [m5, nVisYr, airmass, gamma, msky, theta, km, tvis]
+        [m5, nVisYr, airmass, gamma, msky, theta, km, tvis, EBV]
         For the following keys:
         [m5, nVisYr, gamma, msky, theta, km]
         numbers/paths for specific bands must be passed.
@@ -121,6 +121,7 @@ def __init__(self, args, comm=None):
             "theta",
             "km",
             "tvis",
+            "EBV",
         ]
 
         # validate input parameters
@@ -165,10 +166,12 @@ def _validate_obs_config(self):
 
         # Check if extra keys are passed
         # get lsst_error_model keys
-        lsst_error_model_keys = [field.name for field in fields(LsstErrorParams)]
+        lsst_error_model_keys = list(LsstErrorParams.__dataclass_fields__.keys())
         if len(set(self.config["map_dict"].keys()) - set(lsst_error_model_keys)) != 0:
             extra_keys = set(self.config["map_dict"].keys()) - set(lsst_error_model_keys)
-            raise ValueError("Extra keywords are passed to the configuration: \n" + str(extra_keys))
+            # now we added EBV, which is not in LsstErrorParams:
+            if extra_keys != {"EBV"}:
+                raise ValueError("Extra keywords are passed to the configuration: \n" + str(extra_keys))
 
         # Check data type for the keys:
         # Note that LSSTErrorModel checks
@@ -189,7 +192,7 @@ def _validate_obs_config(self):
                 elif key in self.obs_cond_keys:
 
                     # band-independent keys:
-                    if key in ["airmass", "tvis"]:
+                    if key in ["airmass", "tvis", "EBV"]:
 
                         # check if the input is a string or number
                         if not (
@@ -261,7 +264,7 @@ def _get_maps(self):
             for key in self.config["map_dict"]:
                 if key in self.obs_cond_keys:
                     # band-independent keys:
-                    if key in ["airmass", "tvis"]:
+                    if key in ["airmass", "tvis", "EBV"]:
                         if isinstance(self.config["map_dict"][key], str):
                             maps[key] = hp.read_map(self.config["map_dict"][key])[pixels]
                         elif isinstance(self.config["map_dict"][key], float):
@@ -304,9 +307,10 @@ def get_pixel_conditions(self, pixel: int) -> dict:
             # For keys that may contain the survey condition maps
             if key in self.obs_cond_keys:
                 # band-independent keys:
-                if key in ["airmass", "tvis"]:
-                    obs_conditions[key] = float(self.maps[key][ind])
-                # band-dependent keys:
+                if key in ["airmass", "tvis", "EBV"]:
+                    if key != "EBV":# exclude EBV because it is not in LsstErrorModel
+                        obs_conditions[key] = float(self.maps[key][ind])
+                # band-dependent keys
                 else:
                     obs_conditions[key] = {}
                     for band in (self.maps[key]).keys():
@@ -333,6 +337,45 @@ def assign_pixels(self, catalog: pd.DataFrame) -> pd.DataFrame:
         catalog = pd.concat([catalog, assigned_pix], axis=1)
 
         return catalog
+    
+    # this is milky way extinction, should be added before other observing conditions is applied
+    def apply_galactic_extinction(self, pixel: int, pixel_cat: pd.DataFrame) -> pd.DataFrame:
+        """
+        MW extinction reddening of the magnitudes
+        """
+        # set the A_lamba/E(B-V) values for the six LSST filters 
+        band_a_ebv = {
+            "u":4.81,
+            "g":3.64,
+            "r":2.70,
+            "i":2.06,
+            "z":1.58,
+            "y":1.31,
+        }
+
+        # find the corresponding ebv for the pixel
+        ind = self.maps["pixels"]==pixel
+        ebvvec = self.maps["EBV"][ind]
+
+        if "renameDict" in self.maps.keys():
+            bands = self.maps["renameDict"]
+        elif "renameDict" not in self.maps.keys():
+            bands = {
+                "u":"u",
+                "g":"g",
+                "r":"r",
+                "i":"i",
+                "z":"z",
+                "y":"y",
+            }
+        
+        for b in bands.keys():
+            key=bands[b]
+            # update pixel_cat to the reddened magnitudes
+            pixel_cat[key] = (pixel_cat[key].copy())+ebvvec*band_a_ebv[b]
+
+        return pixel_cat
+
 
     def run(self):
         """
@@ -356,12 +399,17 @@ def run(self):
             # assign each galaxy to a pixel
             print("Assigning pixels.")
             catalog = self.assign_pixels(catalog)
-
+            
             # loop over each pixel
             pixel_cat_list = []
             for pixel, pixel_cat in catalog.groupby("pixel"):
                 # get the observing conditions for this pixel
                 obs_conditions = self.get_pixel_conditions(pixel)
+                
+                # apply MW extinction if supplied, 
+                # replace the Mag column with reddened magnitudes:
+                if "EBV" in self.maps.keys():
+                    pixel_cat = self.apply_galactic_extinction(pixel, pixel_cat)
 
                 # creating the error model for this pixel
                 errorModel = LsstErrorModel(**obs_conditions)
@@ -402,4 +450,4 @@ def __repr__(self):
 
         printMsg += str(self.config["map_dict"])
 
-        return printMsg
+        return printMsg

tests/astro_tools/test_degraders.py

@@ -213,6 +213,7 @@ def test_ObsCondition_extended(data):
     weight = ""
     map_dict = {
         "airmass": find_rail_file('examples_data/creation_data/data/survey_conditions/minion_1016_dc2_Median_airmass_i_and_nightlt1825_HEAL.fits'),
+        "EBV": 0.0,
         "nVisYr": {"u": 50.0},
         "tvis": 30.0,
     }
@@ -242,6 +243,16 @@ def test_ObsCondition_empty_map_dict(data):
     assert degraded_data1.equals(degraded_data2)
 
     os.remove(degrader1.get_output(degrader1.get_aliased_tag("output"), final_name=True))
+    
+
+def test_ObsCondition_renameDict(data):
+    """Test with no renameDict included"""
+    degrader1 = ObsCondition.make_stage(random_seed=0, map_dict={"EBV": 0.0,"renameDict": {"u": "u"},})
+
+    # make sure setting the same seeds yields the same output
+    degraded_data1 = degrader1(data).data
+
+    os.remove(degrader1.get_output(degrader1.get_aliased_tag("output"), final_name=True))
 
 
 def test_LSSTErrorModel_returns_correct_columns(data):