Merge pull request #262 from correac/fix_stellar_abundance_plot

Updates in the stellar abundance plots!

colibre/config.yml

@@ -353,7 +353,7 @@ scripts:
       yvar: C_Fe
       dataset: GALAH
   - filename: scripts/stellar_abundances.py
-    caption: '[Fe/H] vs [C/Fe] using Asplund et al. (2009) values for [Fe/H]Sun = 7.5 and [C/H]Sun = 8.43. The median [C/Fe] vs median [Fe/H] is indicated by the solid curve(s). The scatter points show abundances of individual stellar particles. The observational data for MW compiles the data from the APOGEE survey (Holtzman et al. 2018). Contours use a log scale with 0.04 bin size and a minimum star count of 10.'
+    caption: '[Fe/H] vs [C/Fe] using Asplund et al. (2009) values for [Fe/H]Sun = 7.5 and [C/H]Sun = 8.43. The median [C/Fe] vs median [Fe/H] is indicated by the solid curve(s). The scatter points show abundances of individual stellar particles. The observational data for MW compiles the data from the APOGEE survey (Holtzman et al. 2018) and AstroNN added-value catalog (Leung, H.W. & Bovy, Jo 2019b). We create 6 stellar distributions by selecting stars from APOGEE based on galactocentric radial & azimuthal cuts, and combine them in order to derive a joint stellar abundance distribution that gives less weight to stars in the solar vicinity. The resulting contours use a log scale with 0.25 bin size and a minimum star count of 10.'
     output_file: stellar_abundances_FeH_CFe_APOGEE.png
     section: Stellar Metal Abundances
     title: '[Fe/H] vs [C/Fe]'
@@ -362,7 +362,7 @@ scripts:
       yvar: C_Fe
       dataset: APOGEE
   - filename: scripts/stellar_abundances.py
-    caption: '[Fe/H] vs [N/Fe] using Asplund et al. (2009) values for [Fe/H]Sun = 7.5 and [N/H]Sun = 7.83. The median [N/Fe] vs median [Fe/H] is indicated by the solid curve(s). The scatter points show abundances of individual stellar particles. The observational data for MW compiles the data from the APOGEE survey (Holtzman et al. 2018). Contours use a log scale with 0.04 bin size and a minimum star count of 10.'
+    caption: '[Fe/H] vs [N/Fe] using Asplund et al. (2009) values for [Fe/H]Sun = 7.5 and [N/H]Sun = 7.83. The median [N/Fe] vs median [Fe/H] is indicated by the solid curve(s). The scatter points show abundances of individual stellar particles. The observational data for MW compiles the data from the APOGEE survey (Holtzman et al. 2018) and AstroNN added-value catalog (Leung, H.W. & Bovy, Jo 2019b). We create 6 stellar distributions by selecting stars from APOGEE based on galactocentric radial & azimuthal cuts, and combine them in order to derive a joint stellar abundance distribution that gives less weight to stars in the solar vicinity. The resulting contours use a log scale with 0.25 bin size and a minimum star count of 10.'
     output_file: stellar_abundances_FeH_NFe_APOGEE.png
     section: Stellar Metal Abundances
     title: '[Fe/H] vs [N/Fe]'
@@ -380,7 +380,7 @@ scripts:
       yvar: O_Fe
       dataset: GALAH
   - filename: scripts/stellar_abundances.py
-    caption: '[Fe/H] vs [O/Fe] using Asplund et al. (2009) values for [Fe/H]Sun = 7.5 and [O/H]Sun = 8.69. The median [O/Fe] vs median [Fe/H] is indicated by the solid curve(s). The scatter points show abundances of individual stellar particles. The observational data for MW compiles the data from the APOGEE survey (Holtzman et al. 2018). Contours use a log scale with 0.04 bin size and a minimum star count of 10.'
+    caption: '[Fe/H] vs [O/Fe] using Asplund et al. (2009) values for [Fe/H]Sun = 7.5 and [O/H]Sun = 8.69. The median [O/Fe] vs median [Fe/H] is indicated by the solid curve(s). The scatter points show abundances of individual stellar particles. The observational data for MW compiles the data from the APOGEE survey (Holtzman et al. 2018) and AstroNN added-value catalog (Leung, H.W. & Bovy, Jo 2019b). We create 6 stellar distributions by selecting stars from APOGEE based on galactocentric radial & azimuthal cuts, and combine them in order to derive a joint stellar abundance distribution that gives less weight to stars in the solar vicinity. The resulting contours use a log scale with 0.25 bin size and a minimum star count of 10.'
     output_file: stellar_abundances_FeH_OFe_APOGEE.png
     section: Stellar Metal Abundances
     title: '[Fe/H] vs [O/Fe]'
@@ -414,7 +414,7 @@ scripts:
       yvar: Mg_Fe
       dataset: GALAH
   - filename: scripts/stellar_abundances.py
-    caption: '[Fe/H] vs [Mg/Fe] using Asplund et al. (2009) values for [Fe/H]Sun = 7.5 and [Mg/H]Sun = 7.6. The median [Mg/Fe] vs median [Fe/H] is indicated by the solid curve(s). The scatter points show abundances of individual stellar particles. The observational data for MW compiles the data from the APOGEE survey (Holtzman et al. 2018). Contours use a log scale with 0.04 bin size and a minimum star count of 10.'
+    caption: '[Fe/H] vs [Mg/Fe] using Asplund et al. (2009) values for [Fe/H]Sun = 7.5 and [Mg/H]Sun = 7.6. The median [Mg/Fe] vs median [Fe/H] is indicated by the solid curve(s). The scatter points show abundances of individual stellar particles. The observational data for MW compiles the data from the APOGEE survey (Holtzman et al. 2018) and AstroNN added-value catalog (Leung, H.W. & Bovy, Jo 2019b). We create 6 stellar distributions by selecting stars from APOGEE based on galactocentric radial & azimuthal cuts, and combine them in order to derive a joint stellar abundance distribution that gives less weight to stars in the solar vicinity. The resulting contours use a log scale with 0.25 bin size and a minimum star count of 10.'
     output_file: stellar_abundances_FeH_MgFe_APOGEE.png
     section: Stellar Metal Abundances
     title: '[Fe/H] vs [Mg/Fe]'
@@ -466,7 +466,7 @@ scripts:
       yvar: Eu_Fe
       dataset: GALAH
   - filename: scripts/stellar_abundances.py
-    caption: '[O/H] vs [O/Fe] using Asplund et al. (2009) values for [Fe/H]Sun = 7.5 and [O/H]Sun = 8.69. The median [O/Fe] vs median [O/H] is indicated by the solid curve(s). The scatter points show abundances of individual stellar particles. The observational data for MW compiles the data from the APOGEE survey (Holtzman et al. 2018). Contours use a log scale with 0.04 bin size and a minimum star count of 10.'
+    caption: '[O/H] vs [O/Fe] using Asplund et al. (2009) values for [Fe/H]Sun = 7.5 and [O/H]Sun = 8.69. The median [O/Fe] vs median [O/H] is indicated by the solid curve(s). The scatter points show abundances of individual stellar particles. The observational data for MW compiles the data from the APOGEE survey (Holtzman et al. 2018) and AstroNN added-value catalog (Leung, H.W. & Bovy, Jo 2019b). We create 6 stellar distributions by selecting stars from APOGEE based on galactocentric radial & azimuthal cuts, and combine them in order to derive a joint stellar abundance distribution that gives less weight to stars in the solar vicinity. The resulting contours use a log scale with 0.25 bin size and a minimum star count of 10.'
     output_file: stellar_abundances_OH_OFe_APOGEE.png
     section: Stellar Metal Abundances
     title: '[O/H] vs [O/Fe]'
@@ -475,7 +475,7 @@ scripts:
       yvar: O_Fe
       dataset: APOGEE
   - filename: scripts/stellar_abundances.py
-    caption: '[O/H] vs [Mg/Fe] using Asplund et al. (2009) values for [O/H]Sun = 8.69 and [Mg/H]Sun = 7.6. The median [Mg/Fe] vs median [O/H] is indicated by the solid curve(s). The scatter points show abundances of individual stellar particles. The observational data for MW compiles the data from the APOGEE survey (Holtzman et al. 2018). Contours use a log scale with 0.04 bin size and a minimum star count of 10.'
+    caption: '[O/H] vs [Mg/Fe] using Asplund et al. (2009) values for [O/H]Sun = 8.69 and [Mg/H]Sun = 7.6. The median [Mg/Fe] vs median [O/H] is indicated by the solid curve(s). The scatter points show abundances of individual stellar particles. The observational data for MW compiles the data from the APOGEE survey (Holtzman et al. 2018) and AstroNN added-value catalog (Leung, H.W. & Bovy, Jo 2019b). We create 6 stellar distributions by selecting stars from APOGEE based on galactocentric radial & azimuthal cuts, and combine them in order to derive a joint stellar abundance distribution that gives less weight to stars in the solar vicinity. The resulting contours use a log scale with 0.25 bin size and a minimum star count of 10.'
     output_file: stellar_abundances_OH_MgFe_APOGEE.png
     section: Stellar Metal Abundances
     title: '[O/H] vs [Mg/Fe]'

colibre/scripts/stellar_abundances.py

@@ -14,6 +14,41 @@
 import h5py
 
 
+def make_hist(x, y, cut, xi, yi):
+
+    selection = np.where(cut)[0]
+
+    # Create a histogram
+    h, xedges, yedges = np.histogram2d(
+        x[selection], y[selection], bins=(xi, yi), normed=True
+    )
+
+    return h, xedges, yedges
+
+
+def make_stellar_abundance_distribution(x, y, R, z, xi, yi):
+
+    # Galactocentric radius (R) in kpc units
+    # Galactocentric azimuthal distance (z) in kpc units
+    h = np.zeros((len(xi) - 1, len(yi) - 1))
+
+    # We apply masks to select stars in R & z bins
+    for Ri in range(0, 9, 3):
+        for zi in range(0, 2, 1):
+            mask_R = (R >= Ri) & (R < Ri + 3)
+            mask_z = (np.abs(z) >= zi) & (np.abs(z) < zi + 1)
+            distance_cut = np.logical_and(mask_R, mask_z)
+
+            hist, xedges, yedges = make_hist(x, y, distance_cut, xi, yi)
+
+            # We combine (add) the 6 histograms to give less weight to stars in the solar vicinity.
+            # In this manner all stars in the radial/azimuthal bins contribute with equal weight to
+            # the final stellar distribution
+            h += hist
+
+    return h, xedges, yedges
+
+
 def read_data(data, xvar, yvar):
     """
     Grabs the data
@@ -274,9 +309,13 @@ def read_data(data, xvar, yvar):
         ymax = 2
     else:
         raise AttributeError(f"No APOGEE dataset for x variable {xvar}!")
-    obs_data = load_observations([observational_data])[0]
-    x = obs_data.x
-    y = obs_data.y
+
+    # Reading APOGEE data
+    with h5py.File(observational_data, "r") as obs_data:
+        x = obs_data["x"][:]
+        y = obs_data["y"][:]
+        GalR = obs_data["GalR"][:]  # in kpc units
+        Galz = obs_data["Galz"][:]  # in kpc units
 
     ngridx = 100
     ngridy = 50
@@ -285,8 +324,11 @@ def read_data(data, xvar, yvar):
     xi = np.linspace(xmin, xmax, ngridx)
     yi = np.linspace(ymin, ymax, ngridy)
 
-    # Create a histogram
-    h, xedges, yedges = np.histogram2d(x.value, y.value, bins=(xi, yi))
+    # We apply radial & azimuthal cuts, and combine the stellar distributions
+    # to give less weight to stars in the solar vicinity. We create a histogram
+    # for each distance cut and then combine them
+    h, xedges, yedges = make_stellar_abundance_distribution(x, y, GalR, Galz, xi, yi)
+
     xbins = 0.5 * (xedges[1:] + xedges[:-1])
     ybins = 0.5 * (yedges[1:] + yedges[:-1])
 
@@ -303,6 +345,7 @@ def read_data(data, xvar, yvar):
     )
 
     ax.annotate("APOGEE data", (-3.8, -1.3))
+
 elif dataset == "GALAH":
     observational_data = (
         f"{path_to_obs_data}/data/StellarAbundances/raw/Buder21_data.hdf5"