PyGDSM 1.5.0

pygdsm/gsm08.py

@@ -30,12 +30,13 @@
 from .base_observer import BaseObserver
 from .base_skymodel import BaseSkyModel
 
+T_CMB = 2.725
 
 class GlobalSkyModel(BaseSkyModel):
     """ Global sky model (GSM) class for generating sky models.
     """
 
-    def __init__(self, freq_unit='MHz', basemap='haslam', interpolation='pchip'):
+    def __init__(self, freq_unit='MHz', basemap='haslam', interpolation='pchip', include_cmb=False):
         """ Global sky model (GSM) class for generating sky models.
 
         Upon initialization, the map PCA data are loaded into memory and interpolation
@@ -58,6 +59,9 @@ def __init__(self, freq_unit='MHz', basemap='haslam', interpolation='pchip'):
             piecewise cubic hermitian interpolating polynomial (PCHIP).
             PCHIP is designed to never locally overshoot data, whereas
             splines are designed to have smooth first and second derivatives.
+        include_cmb: bool (default False)
+            Choose whether to include the CMB. Defaults to False. A value of
+            T_CMB = 2.725 K is used.
 
         Notes
         -----
@@ -82,6 +86,7 @@ def __init__(self, freq_unit='MHz', basemap='haslam', interpolation='pchip'):
         super(GlobalSkyModel, self).__init__('GSM2008', GSM_FILEPATH, freq_unit, data_unit, basemap)
 
         self.interpolation_method = interpolation
+        self.include_cmb = include_cmb
 
 
         self.pca_map_data = None
@@ -160,11 +165,17 @@ def generate(self, freqs):
         # c=comp, f=freq, p=pixel. We want to dot product over c for each freq
         #print comps.shape, self.pca_map_data.shape, scaling.shape
         map_out = np.einsum('cf,pc,f->fp', comps, self.pca_map_data, scaling)
+
+        if self.include_cmb:
+            map_out += T_CMB
 
         if map_out.shape[0] == 1:
             map_out = map_out[0]
         self.generated_map_data = map_out
         self.generated_map_freqs = freqs
+
+
+
         return map_out
 
     def set_basemap(self, new_basemap):

pygdsm/gsm16.py

@@ -49,7 +49,7 @@ class GlobalSkyModel16(BaseSkyModel):
     """ Global sky model (GSM) class for generating sky models.
     """
 
-    def __init__(self, freq_unit='MHz', data_unit='TCMB', resolution='hi', theta_rot=0, phi_rot=0, interpolation='pchip'):
+    def __init__(self, freq_unit='MHz', data_unit='TCMB', resolution='hi', theta_rot=0, phi_rot=0, interpolation='pchip', include_cmb=False):
         """ Global sky model (GSM) class for generating sky models.
 
         Upon initialization, the map PCA data are loaded into memory and interpolation
@@ -69,6 +69,9 @@ def __init__(self, freq_unit='MHz', data_unit='TCMB', resolution='hi', theta_rot
             piecewise cubic hermitian interpolating polynomial (PCHIP).
             PCHIP is designed to never locally overshoot data, whereas
             splines are designed to have smooth first and second derivatives.
+        include_cmb: bool (default False)
+            Choose whether to include the CMB. Defaults to False. A value of
+            T_CMB = 2.725 K is used.
 
         Notes
         -----
@@ -94,6 +97,7 @@ def __init__(self, freq_unit='MHz', data_unit='TCMB', resolution='hi', theta_rot
 
         self.interpolation_method = interpolation
         self.resolution = resolution
+        self.include_cmb = include_cmb
 
         # Map data to load
         labels = ['Synchrotron', 'CMB', 'HI', 'Dust1', 'Dust2', 'Free-Free']
@@ -187,10 +191,15 @@ def generate(self, freqs):
 
         output = np.single(np.einsum('cf,pc,f->fp', comps, map_ni.T, scaling))
 
-        for ifreq, freq in enumerate(freqs_ghz):
 
+        for ifreq, freq in enumerate(freqs_ghz):
+
             output[ifreq] = hp.pixelfunc.reorder(output[ifreq], n2r=True)
 
+            # DCP 2023.04.28 - Remove CMB as default
+            if self.include_cmb is False:
+                output -=  K_CMB2MJysr(T, 1e9 * freq)
+
             if self.data_unit == 'TCMB':
                 conversion = 1. / K_CMB2MJysr(1., 1e9 * freq)
             elif self.data_unit == 'TRJ':
@@ -199,7 +208,6 @@ def generate(self, freqs):
                 conversion = 1.
             output[ifreq] *= conversion
 
-#            output.append(result)
 
         if len(output) == 1:
             output = output[0]

pygdsm/haslam.py

@@ -7,18 +7,20 @@
 from .base_skymodel import BaseSkyModel
 from .base_observer import BaseObserver
 
+T_CMB = 2.725
 
 class HaslamSkyModel(BaseSkyModel):
     """ Haslam destriped, desourced sky model """
 
-    def __init__(self,  freq_unit='MHz', spectral_index=-2.6):
+    def __init__(self,  freq_unit='MHz', spectral_index=-2.6, include_cmb=False):
         """ Generate a sky model at a given frequency based off Haslam 408 MHz map
 
         Parameters
         ----------
         freq_unit (str): Frequency unit to use, defaults to MHz
         spectral_index (float): Spectral index to use for calculations.
-
+        include_cmb (bool):  Choose whether to include the CMB. Defaults to False. A value of
+                             T_CMB = 2.725 K is used if True.
         Notes
         -----
         The basemap is the destriped and desoruced Haslam map from Remazeilles (2015)
@@ -47,12 +49,14 @@ def __init__(self,  freq_unit='MHz', spectral_index=-2.6):
         """
         data_unit = 'K'
         basemap = 'Haslam'
-        T_cmb = 2.725
+
         super(HaslamSkyModel, self).__init__('Haslam', HASLAM_FILEPATH, freq_unit, data_unit, basemap)
         self.spectral_index = spectral_index
-        self.data = hp.read_map(self.fits, verbose=False, dtype=np.float64) - T_cmb
+        self.data = hp.read_map(self.fits, verbose=False, dtype=np.float64) - T_CMB
         self.nside = 512
 
+        self.include_cmb = include_cmb
+
     def generate(self, freqs):
         """ Generate a global sky model at a given frequency or frequencies
 
@@ -77,6 +81,8 @@ def generate(self, freqs):
 
         map_out = np.outer((freqs_mhz / 408.0) ** (self.spectral_index), self.data).squeeze()
 
+        if self.include_cmb:
+            map_out += T_CMB
         self.generated_map_data = map_out
         self.generated_map_freqs = freqs
         return map_out

pygdsm/lfsm.py

@@ -7,6 +7,7 @@
 from .base_skymodel import BaseSkyModel
 from .base_observer import BaseObserver
 
+T_CMB = 2.725
 
 def rotate_equatorial_to_galactic(map):
     """
@@ -24,8 +25,14 @@ class LowFrequencySkyModel(BaseSkyModel):
     """ LWA1 Low Frequency Sky Model
     """
 
-    def __init__(self,  freq_unit='MHz'):
+    def __init__(self,  freq_unit='MHz', include_cmb=False):
         """ Global sky model (GSM) class for generating sky models.
+
+        Parameters
+        ----------
+        freq_unit (str): Frequency unit to use, defaults to MHz
+        include_cmb (bool):  Choose whether to include the CMB. Defaults to False. A value of
+                             T_CMB = 2.725 K is used if True.
         """
         data_unit = 'K'
         basemap = 'LFSS'
@@ -35,6 +42,8 @@ def __init__(self,  freq_unit='MHz'):
         self.pca_components =  self.h5['lfsm_components.dat'][:]
         self.nside = 256
 
+        self.include_cmb = include_cmb
+
         freqs  = self.pca_components[:, 0]
         sigmas = self.pca_components[:, 1]
         comps  = self.pca_components[:, 2:]
@@ -90,6 +99,10 @@ def generate(self, freqs):
             map_out[ff] =  rotate_equatorial_to_galactic(map_out[ff])
 
         map_out = map_out.squeeze()
+
+        if self.include_cmb == False:
+            map_out -= T_CMB
+
         self.generated_map_data = map_out
         self.generated_map_freqs = freqs
         return map_out

tests/test_gsm.py

@@ -104,10 +104,21 @@ def test_write_fits():
 
     os.remove("test_write_fits.fits")
 
+def test_cmb_removal():
+    g = GlobalSkyModel(freq_unit='MHz', include_cmb=False)
+    sky_no_cmb = g.generate(400)
+    g = GlobalSkyModel(freq_unit='MHz',  include_cmb=True)
+    sky_with_cmb = g.generate(400)    
+
+    T_cmb = (sky_with_cmb - sky_no_cmb).mean()
+    print(T_cmb)
+    assert np.isclose(T_cmb, 2.725)
+
 
 if __name__ == "__main__":
     test_gsm_generate()
     compare_to_gsm()
     test_speed()
     test_write_fits()
     test_set_methods()
+    test_cmb_removal()

tests/test_gsm2016.py

@@ -82,6 +82,16 @@ def test_gsm_opts():
     with pytest.raises(RuntimeError):
         g = GlobalSkyModel16(data_unit='furlongs/fortnight')
 
+def test_cmb_removal():
+    g = GlobalSkyModel16(freq_unit='GHz', resolution='lo', data_unit='TCMB', include_cmb=False)
+    sky_no_cmb = g.generate(400)
+    g = GlobalSkyModel16(freq_unit='GHz', resolution='lo', data_unit='TCMB', include_cmb=True)
+    sky_with_cmb = g.generate(400)    
+
+    T_cmb = (sky_with_cmb - sky_no_cmb).mean()
+    print(T_cmb)
+    assert np.isclose(T_cmb, 2.725)
+
 
 if __name__ == "__main__":
     test_compare_gsm_to_old()

tests/test_haslam.py

@@ -41,6 +41,17 @@ def test_observer_test():
     d = ov.view(logged=True)
     plt.show()
 
+def test_cmb_removal():
+    g = HaslamSkyModel(freq_unit='MHz', include_cmb=False)
+    sky_no_cmb = g.generate(400)
+    g = HaslamSkyModel(freq_unit='MHz',  include_cmb=True)
+    sky_with_cmb = g.generate(400)    
+
+    T_cmb = (sky_with_cmb - sky_no_cmb).mean()
+    print(T_cmb)
+    assert np.isclose(T_cmb, 2.725)
+
+
 if __name__ == "__main__":
     test_compare_gsm_to_old()
     test_observer_test()

tests/test_lfsm.py

@@ -43,6 +43,16 @@ def test_observer_test():
     d = ov.view(logged=True)
     plt.show()
 
+def test_cmb_removal():
+    g = LowFrequencySkyModel(freq_unit='MHz', include_cmb=False)
+    sky_no_cmb = g.generate(400)
+    g = LowFrequencySkyModel(freq_unit='MHz',  include_cmb=True)
+    sky_with_cmb = g.generate(400)    
+
+    T_cmb = (sky_with_cmb - sky_no_cmb).mean()
+    print(T_cmb)
+    assert np.isclose(T_cmb, 2.725)
+
 if __name__ == "__main__":
     test_compare_gsm_to_old()
     # test_observer_test()