Merge pull request #10 from SWIFTSIM/pipeline_with_soap_catalogues

Working version of the morphology pipeline for SOAP output

morphology-pipeline

@@ -12,6 +12,7 @@ Usage:
     -C/--config      <configuration directory> \
     -i/--input       <input folder that contains the snapshots and catalogues> \
     -s/--snapshots   <input snapshot> \
+    -g/--groups      <input particle membership> \
     -c/--catalogues  <input catalogue> \
     -o/--output      <directory where output images and web pages are stored> \
    [-n/--run-names   <name for the run>] \
@@ -41,6 +42,7 @@ parser = ap.ArgumentParser(
         "    -C/--config      <configuration directory> \ \n"
         "    -i/--input       <input folder that contains the snapshots and catalogues> \ \n"
         "    -s/--snapshots   <input snapshot> \ \n"
+        "    -g/--groups      <input particle membership> \ \n"
         "    -c/--catalogues  <input catalogue> \ \n"
         "    -o/--output      <directory where output images and web pages are stored> \ \n"
         "   [-n/--run-names   <name for the run>] \ \n"
@@ -66,7 +68,16 @@ parser.add_argument(
     "--catalogues",
     type=str,
     required=True,
-    help="Name of the VELOCIraptor HDF5 .properties file(s). Required.",
+    help="Name of the SOAP HDF5 properties file(s). Required.",
+    nargs="*",
+)
+
+parser.add_argument(
+    "-g",
+    "--groups",
+    type=str,
+    required=True,
+    help="Name of the SOAP HDF5 memberships file(s). Required.",
     nargs="*",
 )
 
@@ -302,6 +313,7 @@ if __name__ == "__main__":
             # compute the galaxy properties
             # start by getting the catalogue and snapshot file name
             halo_catalogue_filename = f"{args.input[0]}/{args.catalogues[0]}"
+            halo_membership_filename = f"{args.input[0]}/{args.groups[0]}"
             snapshot_filename = f"{args.input[0]}/{args.snapshots[0]}"
 
             # load the catalogue and create the filtered catalogue that only
@@ -329,8 +341,9 @@ if __name__ == "__main__":
             # these should correspond to the argument needed by
             # process_galaxy():
             #  - the index in the global galaxy list
-            #  - the catalogue index of the galaxy (its VR/SOAP index)
+            #  - the catalogue index of the galaxy (its SOAP index)
             #  - the catalogue file name
+            #  - the halo membership file name
             #  - the snapshot file name
             #  - the output directory name
             #  - the observational data path
@@ -342,6 +355,7 @@ if __name__ == "__main__":
                     index,
                     galaxy_index,
                     halo_catalogue_filename,
+                    halo_membership_filename,
                     snapshot_filename,
                     args.output,
                     observational_data_path,

morpholopy/KS.py

@@ -362,12 +362,12 @@ def plot_KS_relations(
        Data for all the simulations that need to be plotted (can be a single galaxy for
        individual galaxy plots).
      - prefix: str
-       Unique prefix that is prepended to images file names. Useful to distinguish
+       Unique prefix that is prepended to image file names. Useful to distinguish
        individual galaxy images.
      - always_plot_scatter: bool
-       Wheter or not to always plot individual median histrogram bins as well as median lines.
+       Whether to always plot individual median histogram bins as well as median lines.
      - plot_integrated_quantities: bool
-       Whether or not to plot integrated KS plots. Set to False for individual galaxy plots,
+       Whether to plot integrated KS plots. Set to False for individual galaxy plots,
        since a galaxy is only a single point on those plots.
 
     Returns an image data dictionary compatible with MorphologyConfig.add_images().
@@ -444,7 +444,6 @@ def plot_KS_relations(
                     ax.plot([], [], marker=marker, color="k", linestyle="None")[0]
                 )
                 sim_labels.append(name)
-
             if dataname is not None:
                 observational_data = load_observations(
                     sorted(glob.glob(f"{observational_data_path}/{dataname}/*.hdf5"))
@@ -578,14 +577,32 @@ def plot_KS_relations(
         ax.add_artist(sim_legend)
 
     neut_filename = f"{prefix}azimuthal_KS_neutral.png"
+    ax_neut.set_xlabel(
+        r"Neutral Surface Density $\Sigma_{\rm HI+H2}$ $\rm \left[\frac{M_\odot}{pc^{2}}\right]$"
+    )
+    ax_neut.set_ylabel(
+        r"Star Formation Rate Surface Density $\rm \left[\frac{M_\odot}{yr \cdot kpc^{2}}\right]$"
+    )
     fig_neut.savefig(f"{output_path}/{neut_filename}", dpi=300)
     pl.close(fig_neut)
 
     at_filename = f"{prefix}azimuthal_KS_atomic.png"
+    ax_at.set_xlabel(
+        r"Atomic Surface Density $\Sigma_{\rm HI}$ $\rm \left[\frac{M_\odot}{pc^{2}}\right]$"
+    )
+    ax_at.set_ylabel(
+        r"Star Formation Rate Surface Density $\rm \left[\frac{M_\odot}{yr \cdot kpc^{2}}\right]$"
+    )
     fig_at.savefig(f"{output_path}/{at_filename}", dpi=300)
     pl.close(fig_at)
 
     mol_filename = f"{prefix}azimuthal_KS_molecular.png"
+    ax_mol.set_xlabel(
+        r"Molecular Surface Density $\Sigma_{\rm H2}$ $\rm \left[\frac{M_\odot}{pc^{2}}\right]$"
+    )
+    ax_mol.set_ylabel(
+        r"Star Formation Rate Surface Density $\rm \left[\frac{M_\odot}{yr \cdot kpc^{2}}\right]$"
+    )
     fig_mol.savefig(f"{output_path}/{mol_filename}", dpi=300)
     pl.close(fig_mol)
 
@@ -706,14 +723,32 @@ def plot_KS_relations(
         ax.add_artist(sim_legend)
 
     neut_filename = f"{prefix}spatial_KS_neutral.png"
+    ax_neut.set_xlabel(
+        r"Neutral Surface Density $\Sigma_{\rm HI+H2}$ $\rm \left[\frac{M_\odot}{pc^{2}}\right]$"
+    )
+    ax_neut.set_ylabel(
+        r"Star Formation Rate Surface Density $\rm \left[\frac{M_\odot}{yr \cdot kpc^{2}}\right]$"
+    )
     fig_neut.savefig(f"{output_path}/{neut_filename}", dpi=300)
     pl.close(fig_neut)
 
     at_filename = f"{prefix}spatial_KS_atomic.png"
+    ax_at.set_xlabel(
+        r"Atomic Surface Density $\Sigma_{\rm HI}$ $\rm \left[\frac{M_\odot}{pc^{2}}\right]$"
+    )
+    ax_at.set_ylabel(
+        r"Star Formation Rate Surface Density $\rm \left[\frac{M_\odot}{yr \cdot kpc^{2}}\right]$"
+    )
     fig_at.savefig(f"{output_path}/{at_filename}", dpi=300)
     pl.close(fig_at)
 
     mol_filename = f"{prefix}spatial_KS_molecular.png"
+    ax_mol.set_xlabel(
+        r"Molecular Surface Density $\Sigma_{\rm H2}$ $\rm \left[\frac{M_\odot}{pc^{2}}\right]$"
+    )
+    ax_mol.set_ylabel(
+        r"Star Formation Rate Surface Density $\rm \left[\frac{M_\odot}{yr \cdot kpc^{2}}\right]$"
+    )
     fig_mol.savefig(f"{output_path}/{mol_filename}", dpi=300)
     pl.close(fig_mol)
 
@@ -821,14 +856,32 @@ def plot_KS_relations(
         ax.add_artist(sim_legend)
 
     neut_filename = f"{prefix}spatial_tgas_neutral.png"
+    ax_neut.set_xlabel(
+        r"Neutral Surface Density $\Sigma_{\rm HI+H2}$ $\rm \left[\frac{M_\odot}{pc^{2}}\right]$"
+    )
+    ax_neut.set_ylabel(
+        r"Neutral gas depletion time $\Sigma_{\rm HI+H2} / \Sigma_{\rm SFR}$ [yr]"
+    )
     fig_neut.savefig(f"{output_path}/{neut_filename}", dpi=300)
     pl.close(fig_neut)
 
     at_filename = f"{prefix}spatial_tgas_atomic.png"
+    ax_at.set_xlabel(
+        r"Atomic Surface Density $\Sigma_{\rm HI}$ $\rm \left[\frac{M_\odot}{pc^{2}}\right]$"
+    )
+    ax_at.set_ylabel(
+        r"Atomic gas depletion time $\Sigma_{\rm HI} / \Sigma_{\rm SFR}$ [yr]"
+    )
     fig_at.savefig(f"{output_path}/{at_filename}", dpi=300)
     pl.close(fig_at)
 
     mol_filename = f"{prefix}spatial_tgas_molecular.png"
+    ax_mol.set_xlabel(
+        r"Molecular Surface Density $\Sigma_{\rm H2}$ $\rm \left[\frac{M_\odot}{pc^{2}}\right]$"
+    )
+    ax_mol.set_ylabel(
+        r"Molecular gas depletion time $\Sigma_{\rm H2} / \Sigma_{\rm SFR}$ [yr]"
+    )
     fig_mol.savefig(f"{output_path}/{mol_filename}", dpi=300)
     pl.close(fig_mol)
 

morpholopy/filtered_catalogue.py

@@ -12,10 +12,9 @@
 
 import numpy as np
 import unyt
-from functools import reduce
 
 from numpy.typing import NDArray
-from velociraptor.catalogue.catalogue import VelociraptorCatalogue
+from velociraptor.catalogue.catalogue import Catalogue
 
 
 class FilteredCatalogue:
@@ -30,7 +29,7 @@ class FilteredCatalogue:
 
     def __init__(
         self,
-        catalogue: VelociraptorCatalogue,
+        catalogue: Catalogue,
         minimum_mass_stars: unyt.unyt_quantity,
         mass_variable_stars: str,
         minimum_mass_gas: unyt.unyt_quantity,
@@ -74,13 +73,13 @@ def __init__(
         """
 
         # get the masses from the catalogue
-        Mstar = reduce(getattr, mass_variable_stars.split("."), catalogue)
-        Mgas = reduce(getattr, mass_variable_gas.split("."), catalogue)
+        Mstar = catalogue.get_quantity(mass_variable_stars)
+        Mgas = catalogue.get_quantity(mass_variable_gas)
 
         # compute the mask
         mask = (
             (Mstar >= minimum_mass_stars)
-            & (catalogue.structure_type.structuretype == 10)
+            & (catalogue.get_quantity("structure_type.structuretype") == 10)
             & (Mgas > minimum_mass_gas)
         )
         # turn the mask into a list of indices