From 5c0add7e14a9de83a8d27625cd9d085fc0778727 Mon Sep 17 00:00:00 2001 From: EvgeniiChaikin Date: Thu, 2 May 2024 12:06:41 +0200 Subject: [PATCH 1/4] Make pipeline work for SOAP output instead of VR output --- morphology-pipeline | 18 ++++- morpholopy/KS.py | 61 ++++++++++++++- morpholopy/filtered_catalogue.py | 11 ++- morpholopy/galaxy_data.py | 130 ++++++++++++++++++++----------- morpholopy/galaxy_plots.py | 63 +++++++++------ morpholopy/morphology.py | 24 +++--- morpholopy/orientation.py | 13 ++-- requirements.txt | 2 +- velociraptor-comparison-data | 2 +- 9 files changed, 220 insertions(+), 104 deletions(-) diff --git a/morphology-pipeline b/morphology-pipeline index bc476bd..2150003 100755 --- a/morphology-pipeline +++ b/morphology-pipeline @@ -12,6 +12,7 @@ Usage: -C/--config \ -i/--input \ -s/--snapshots \ + -g/--groups \ -c/--catalogues \ -o/--output \ [-n/--run-names ] \ @@ -41,6 +42,7 @@ parser = ap.ArgumentParser( " -C/--config \ \n" " -i/--input \ \n" " -s/--snapshots \ \n" + " -g/--groups \ \n" " -c/--catalogues \ \n" " -o/--output \ \n" " [-n/--run-names ] \ \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, diff --git a/morpholopy/KS.py b/morpholopy/KS.py index f2d1bc5..c07b773 100644 --- a/morpholopy/KS.py +++ b/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) diff --git a/morpholopy/filtered_catalogue.py b/morpholopy/filtered_catalogue.py index e2acb91..d8f481a 100644 --- a/morpholopy/filtered_catalogue.py +++ b/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 diff --git a/morpholopy/galaxy_data.py b/morpholopy/galaxy_data.py index f97c912..2faa515 100644 --- a/morpholopy/galaxy_data.py +++ b/morpholopy/galaxy_data.py @@ -14,11 +14,12 @@ """ import numpy as np -from swiftsimio import load as load_snapshot +import swiftsimio from swiftsimio.objects import cosmo_array from velociraptor import load as load_catalogue -from velociraptor.particles import load_groups -from velociraptor.swift.swift import to_swiftsimio_dataset +import h5py as h5 +from collections import namedtuple + import unyt import yaml @@ -98,9 +99,9 @@ # for brevity, we append individual plots in logical chunks: # - spatially resolved and azimuthally averaged KS plots for type, label in [ - ("neutral", "$\\Sigma{}_{{\\rm{}HI}+{\\rm{}H2}}$"), - ("HI", "$\\Sigma{}_{\\rm{}HI}$"), - ("H2", "$\\Sigma{}_{\\rm{}H2}$"), + ("neutral", r"$\Sigma_{\rm HI+H2}$"), + ("HI", r"$\Sigma_{\rm HI}$"), + ("H2", r"$\Sigma_{\rm H2}$"), ]: for Zmask in ["all", "Zm1", "Z0", "Zp1"]: medians[f"sigma_{type}_tgas_spatial_{Zmask}"] = { @@ -113,7 +114,7 @@ "x units": "Msun/pc**2", "y units": "yr", "x label": label, - "y label": label + "$/\\Sigma{}_{{\\rm{}SFR}}}$", + "y label": label + r"$/ \Sigma_{\rm SFR}$", } for method in ["spatial", "azimuthal"]: medians[f"sigma_{type}_SFR_{method}_{Zmask}"] = { @@ -126,7 +127,7 @@ "x units": "Msun/pc**2", "y units": "Msun/yr/kpc**2", "x label": label, - "y label": "$\\Sigma{}_{{\\rm{}SFR}}}$", + "y label": r"$\Sigma_{\rm SFR}$", } # - Schruba like plots for Zmask in ["all", "Zm1", "Z0", "Zp1"]: @@ -139,8 +140,8 @@ "range in y": [-8.0, 1.0], "x units": "Msun/pc**2", "y units": "dimensionless", - "x label": "$\\Sigma{}_{{\\rm{}HI}+{\\rm{}H2}}$", - "y label": "$\\Sigma{}_{\\rm{}H2}/\\Sigma{}_{{\\rm{}HI}+{\\rm{}H2}}$", + "x label": r"$\Sigma_{\rm HI+H2}$", + "y label": r"$\Sigma_{\rm H2}/\Sigma_{\rm HI+H2}$", } medians[f"H2_to_HI_vs_neutral_spatial_{Zmask}"] = { "number of bins x": 20, @@ -151,8 +152,8 @@ "range in y": [-2.0, 3.0], "x units": "Msun/pc**2", "y units": "dimensionless", - "x label": "$\\Sigma{}_{{\\rm{}HI}+{\\rm{}H2}}$", - "y label": "$\\Sigma{}_{\\rm{}H2}/\\Sigma{}_{\\rm{}HI}$", + "x label": r"$\Sigma_{\rm HI+H2}$", + "y label": r"$\Sigma_{\rm H2}/\Sigma_{\rm HI}$", } # - Sanchez like plots @@ -165,8 +166,8 @@ "range in y": [-3.0, 1.0], "x units": "Msun/pc**2", "y units": "dimensionless", - "x label": "$\\Sigma{}_{\\bigstar{}}$", - "y label": "$\\Sigma{}_{\\rm{}H2}/\\Sigma{}_{\\bigstar{}}$", + "x label": r"$\Sigma_{\bigstar}$", + "y label": r"$\Sigma_{\rm H2}/\Sigma_{\bigstar}$", } medians["SFR_to_H2_vs_H2_spatial"] = { "number of bins x": 20, @@ -177,8 +178,8 @@ "range in y": [-11.0, -7.0], "x units": "Msun/pc**2", "y units": "1/yr", - "x label": "$\\Sigma{}_{\\rm{}H2}$", - "y label": "$\\Sigma{}_{\\rm{}SFR}/\\Sigma{}_{\\rm{}H2}$", + "x label": r"$\Sigma_{\rm H2}$", + "y label": r"$\Sigma_{\rm SFR}/\Sigma_{\rm H2}$", } medians["SFR_to_star_vs_star_spatial"] = { "number of bins x": 20, @@ -189,8 +190,8 @@ "range in y": [-13.0, -7.0], "x units": "Msun/pc**2", "y units": "1/yr", - "x label": "$\\Sigma{}_{\\bigstar{}}$", - "y label": "$\\Sigma{}_{\\rm{}SFR}/\\Sigma{}_{\\bigstar{}}$", + "x label": r"$\Sigma{}_{\bigstar}$", + "y label": r"$\Sigma{}_{\rm SFR}/\Sigma_{\bigstar}$", } # construct a structured array data type for the median arrays @@ -413,6 +414,7 @@ def process_galaxy( index, galaxy_index, catalogue_filename, + halo_membership_filename, snapshot_filename, output_path, observational_data_path, @@ -433,18 +435,55 @@ def process_galaxy( # (re)load the catalogue to read the properties of this particular galaxy # we need to disregard the units to avoid problems with the SFR unit catalogue = load_catalogue(catalogue_filename, disregard_units=True) - groups = load_groups( - catalogue_filename.replace(".properties", ".catalog_groups"), - catalogue=catalogue, - ) + membership_data = h5.File(halo_membership_filename, "r") # get the particles belonging to this galaxy - particles, _ = groups.extract_halo(halo_index=galaxy_index) + # particles, _ = groups.extract_halo(halo_index=galaxy_index) # turn this information into a swiftsimio mask and read the data # using swiftsimio - data, mask = to_swiftsimio_dataset( - particles, snapshot_filename, generate_extra_mask=True + + ##data, mask = to_swiftsimio_dataset( + ## particles, snapshot_filename, generate_extra_mask=True + ##) + + swift_mask = swiftsimio.mask(snapshot_filename, spatial_only=True) + # SWIFT data is stored in comoving units, so we need to un-correct + # the velociraptor data if it is stored in physical. + length_factor = catalogue.a + + halo_x = catalogue.get_quantity("positions.xc")[galaxy_index] / length_factor + halo_y = catalogue.get_quantity("positions.yc")[galaxy_index] / length_factor + halo_z = catalogue.get_quantity("positions.zc")[galaxy_index] / length_factor + r_size = ( + catalogue.get_quantity("searchradius.search_radius")[galaxy_index] + / length_factor ) + + spatial_mask = [ + [halo_x - r_size, halo_x + r_size], + [halo_y - r_size, halo_y + r_size], + [halo_z - r_size, halo_z + r_size], + ] + + swift_mask.constrain_spatial(spatial_mask) + data = swiftsimio.load(snapshot_filename) + + particle_name_masks = {} + for particle_name, particle_type in zip( + data.metadata.present_particle_names, data.metadata.present_particle_types + ): + mask_per_type = ( + membership_data[f"PartType{particle_type}"]["GroupNr_bound"][:] + == galaxy_index + ) + particle_name_masks[particle_name] = mask_per_type + MaskTuple = namedtuple("MaskCollection", data.metadata.present_particle_names) + mask = MaskTuple(**particle_name_masks) + membership_data.close() + + for key, val in particle_name_masks.items(): + print(key, np.sum(val)) + # mask out all particles that are actually bound to the galaxy # while at it: convert everything to physical coordinates for parttype in ["gas", "dark_matter", "stars", "black_holes"]: @@ -470,9 +509,9 @@ def process_galaxy( galaxy_center = cosmo_array( unyt.unyt_array( [ - catalogue.positions.xcmbp[galaxy_index], - catalogue.positions.ycmbp[galaxy_index], - catalogue.positions.zcmbp[galaxy_index], + catalogue.get_quantity("positions.xcmbp")[galaxy_index], + catalogue.get_quantity("positions.ycmbp")[galaxy_index], + catalogue.get_quantity("positions.zcmbp")[galaxy_index], ] ), comoving=False, @@ -481,9 +520,9 @@ def process_galaxy( galaxy_velocity = cosmo_array( unyt.unyt_array( [ - catalogue.velocities.vxc[galaxy_index], - catalogue.velocities.vyc[galaxy_index], - catalogue.velocities.vzc[galaxy_index], + catalogue.get_quantity("velocities.vxc")[galaxy_index], + catalogue.get_quantity("velocities.vyc")[galaxy_index], + catalogue.get_quantity("velocities.vzc")[galaxy_index], ] ), comoving=False, @@ -492,11 +531,13 @@ def process_galaxy( # store some of the catalogue properties directly into the galaxy data array galaxy_data["stellar_mass"] = cosmo_array( - catalogue.apertures.mass_star_50_kpc[galaxy_index].to("Msun"), + catalogue.get_quantity("apertures.mass_star_50_kpc")[galaxy_index].to("Msun"), comoving=False, cosmo_factor=data.gas.masses.cosmo_factor, ) - r_halfmass_star = catalogue.radii.r_halfmass_star[galaxy_index] + r_halfmass_star = catalogue.get_quantity("apertures.rhalfmass_star_50_kpc")[ + galaxy_index + ] galaxy_data["half_mass_radius_star"] = cosmo_array( r_halfmass_star.to("kpc"), comoving=False, @@ -510,8 +551,12 @@ def process_galaxy( # Lowest sSFR below which the galaxy is considered passive marginal_ssfr = unyt.unyt_quantity(1e-11, units=1 / unyt.year) - stellar_mass = catalogue.apertures.mass_star_50_kpc[galaxy_index].to("Msun") - star_formation_rate = catalogue.apertures.sfr_gas_50_kpc[galaxy_index] + stellar_mass = catalogue.get_quantity("apertures.mass_star_50_kpc")[ + galaxy_index + ].to("Msun") + star_formation_rate = catalogue.get_quantity("apertures.sfr_gas_50_kpc")[ + galaxy_index + ] if stellar_mass == unyt.unyt_quantity(0.0, units=unyt.msun): ssfr = unyt.unyt_quantity(0.0, units=1 / unyt.year) @@ -530,17 +575,12 @@ def process_galaxy( # get other radius quantities from the catalogue that the orientation # calculation might use Rhalf = cosmo_array( - catalogue.apertures.rhalfmass_star_50_kpc[galaxy_index], + catalogue.get_quantity("apertures.rhalfmass_star_50_kpc")[galaxy_index], comoving=False, cosmo_factor=data.gas.coordinates.cosmo_factor, ) R200crit = cosmo_array( - catalogue.spherical_overdensities.r_200_rhocrit[galaxy_index], - comoving=False, - cosmo_factor=data.gas.coordinates.cosmo_factor, - ) - Rvir = cosmo_array( - catalogue.radii.rvir[galaxy_index], + catalogue.get_quantity("spherical_overdensities.r_200_rhocrit")[galaxy_index], comoving=False, cosmo_factor=data.gas.coordinates.cosmo_factor, ) @@ -600,7 +640,7 @@ def process_galaxy( # determine the angular momentum vector and corresponding face-on and # edge-on rotation matrices orientation_vector, face_on_rmatrix, edge_on_rmatrix = get_orientation_matrices( - data, Rhalf, R200crit, Rvir, orientation_type + data, Rhalf, R200crit, orientation_type ) # done with the pre-processing: @@ -608,6 +648,7 @@ def process_galaxy( # - axis ratios and angular momenta (a, b, c), z_axis = get_axis_lengths_reduced_tensor(galaxy_log, data.stars, Rhalf) + print(a, b, c, z_axis) if (a > 0.0) and (b > 0.0) and (c > 0.0): galaxy_data["stars_axis_ca_reduced"] = c / a galaxy_data["stars_axis_cb_reduced"] = c / b @@ -632,7 +673,7 @@ def process_galaxy( galaxy_data["stars_z_axis_normal"][:] = np.nan stars_momentum, kappa_corot = get_kappa_corot( - data.stars, Rhalf, R200crit, Rvir, orientation_type, orientation_vector + data.stars, Rhalf, R200crit, orientation_type, orientation_vector ) galaxy_data["stars_kappa_co"] = kappa_corot galaxy_data["orientation_vector"] = orientation_vector @@ -670,7 +711,6 @@ def process_galaxy( data.gas, Rhalf, R200crit, - Rvir, orientation_type, orientation_vector, mass_variable="H_neutral_mass", diff --git a/morpholopy/galaxy_plots.py b/morpholopy/galaxy_plots.py index 4fc93c6..ea868d5 100644 --- a/morpholopy/galaxy_plots.py +++ b/morpholopy/galaxy_plots.py @@ -12,16 +12,13 @@ matplotlib.use("Agg") import matplotlib.pyplot as pl -from matplotlib import gridspec import unyt -from swiftsimio.visualisation.rotation import rotation_matrix_from_vector from swiftsimio.visualisation.projection import project_gas from swiftsimio.visualisation.projection import project_pixel_grid from swiftsimio.visualisation.slice import kernel_gamma from swiftsimio.visualisation.smoothing_length_generation import ( generate_smoothing_lengths, ) -from swiftsimio import swift_cosmology_to_astropy from astropy.visualization import make_lupton_rgb from scipy.optimize import curve_fit @@ -87,7 +84,6 @@ def calculate_scaleheight_pointmasses(zcoords, weights, zrange, resolution): ) z_1D = (bin_edges[1:] + bin_edges[:-1]) / 2.0 - p0 = (hist.max(), 1.0, 0.0) try: popt, pcov = curve_fit( @@ -107,8 +103,6 @@ def calculate_scaleheight_fit(mass_map, r_img_kpc, r_abs_max_kpc): z_1D = np.ravel(z[:, (np.abs(xx) < r_abs_max_kpc)]) S_1D = np.ravel(mass_map[:, (np.abs(xx) < r_abs_max_kpc)]) - p0 = (mass_map.max(), 1.0, 0.0) - try: popt, pcov = curve_fit( exponential, @@ -512,42 +506,60 @@ def plot_galaxy( image_info += " resolution: (%i x %i) pixel." % (npix, npix) galaxy_coords_text = ( - f"[ {catalogue.positions.xcmbp[halo_id].to('Mpc').value:.02f}, " - + f"{catalogue.positions.ycmbp[halo_id].to('Mpc').value:.02f}, " - + f"{catalogue.positions.zcmbp[halo_id].to('Mpc').value:.02f} ] cMpc" + f"[ {catalogue.get_quantity('positions.xcmbp')[halo_id].to('Mpc').value:.02f}, " + + f"{catalogue.get_quantity('positions.ycmbp')[halo_id].to('Mpc').value:.02f}, " + + f"{catalogue.get_quantity('positions.zcmbp')[halo_id].to('Mpc').value:.02f} ] cMpc" ) galaxy_info_title = f"Galaxy {halo_id:08d} " + galaxy_coords_text galaxy_info_short = ( r"M$_{\mathrm{200,crit}}$ = " - + sci_notation(catalogue.masses.mass_200crit[halo_id].to("Msun").value) + + sci_notation( + catalogue.get_quantity("masses.mass_200crit")[halo_id].to("Msun").value + ) + r" M$_{\odot}$, " + r"M$_{\mathrm{*,30kpc}}$ = " - + sci_notation(catalogue.masses.mass_star_30kpc[halo_id].to("Msun").value) + + sci_notation( + catalogue.get_quantity("apertures.mass_star_30_kpc")[halo_id] + .to("Msun") + .value + ) + r" M$_{\odot}$, " ) galaxy_info = " Coordinates (x,y,z) = " + galaxy_coords_text + ", " galaxy_info += ( r"M$_{\mathrm{200,crit}}$ = " - + sci_notation(catalogue.masses.mass_200crit[halo_id].to("Msun").value) + + sci_notation( + catalogue.get_quantity("masses.mass_200crit")[halo_id].to("Msun").value + ) + r" M$_{\odot}$, " ) galaxy_info += ( r"M$_{\mathrm{*,30kpc}}$ = " - + sci_notation(catalogue.masses.mass_star_30kpc[halo_id].to("Msun").value) + + sci_notation( + catalogue.get_quantity("apertures.mass_star_30_kpc")[halo_id] + .to("Msun") + .value + ) + r" M$_{\odot}$, " ) galaxy_info += ( r"M$_{\mathrm{gas,30kpc}}$ = " - + sci_notation(catalogue.masses.mass_gas_30kpc[halo_id].to("Msun").value) + + sci_notation( + catalogue.get_quantity("apertures.mass_gas_30_kpc")[halo_id] + .to("Msun") + .value + ) + r" M$_{\odot}$, " ) galaxy_info += ( r"M$_{\mathrm{HI,30kpc}}$ = " + sci_notation( - catalogue.gas_hydrogen_species_masses.HI_mass_30_kpc[halo_id] + catalogue.get_quantity("gas_hydrogen_species_masses.HI_mass_30_kpc")[ + halo_id + ] .to("Msun") .value ) @@ -556,15 +568,17 @@ def plot_galaxy( galaxy_info += ( r"M$_{\mathrm{H2,30kpc}}$ = " + sci_notation( - catalogue.gas_hydrogen_species_masses.H2_mass_30_kpc[halo_id] + catalogue.get_quantity("gas_hydrogen_species_masses.H2_mass_30_kpc")[ + halo_id + ] .to("Msun") .value ) + r" M$_{\odot}$, " ) sSFR = ( - catalogue.apertures.sfr_gas_100_kpc[halo_id] - / catalogue.apertures.mass_star_100_kpc[halo_id] + catalogue.get_quantity("apertures.sfr_gas_100_kpc")[halo_id] + / catalogue.get_quantity("apertures.mass_star_100_kpc")[halo_id] ) if np.isfinite(sSFR): galaxy_info += ( @@ -573,12 +587,13 @@ def plot_galaxy( + r" Gyr$^{-1}$" ) - stars_faceon_filename = "galaxy_%3.3i_map_stars_faceon.png" % (index) - stars_edgeon_filename = "galaxy_%3.3i_map_stars_edgeon.png" % (index) - HI_faceon_filename = "galaxy_%3.3i_map_HI_faceon.png" % (index) - HI_edgeon_filename = "galaxy_%3.3i_map_HI_edgeon.png" % (index) - H2_faceon_filename = "galaxy_%3.3i_map_H2_faceon.png" % (index) - H2_edgeon_filename = "galaxy_%3.3i_map_H2_edgeon.png" % (index) + base_file_name = f"galaxy_{index:03d}_map" + stars_faceon_filename = f"{base_file_name}_stars_faceon.png" + stars_edgeon_filename = f"{base_file_name}_stars_edgeon.png" + HI_faceon_filename = f"{base_file_name}_HI_faceon.png" + HI_edgeon_filename = f"{base_file_name}_HI_edgeon.png" + H2_faceon_filename = f"{base_file_name}_H2_faceon.png" + H2_edgeon_filename = f"{base_file_name}_H2_edgeon.png" plots = {} diff --git a/morpholopy/morphology.py b/morpholopy/morphology.py index 4224fa7..524d5a1 100644 --- a/morpholopy/morphology.py +++ b/morpholopy/morphology.py @@ -9,7 +9,6 @@ import numpy as np import unyt import swiftsimio as sw -from swiftsimio.visualisation.rotation import rotation_matrix_from_vector from swiftsimio import SWIFTDataset from .orientation import get_orientation_mask @@ -22,7 +21,6 @@ from typing import Dict, List, Tuple, Union from numpy.typing import NDArray from .logging import GalaxyLog -from .medians import plot_median_on_axis_as_line from scipy.optimize import curve_fit from scipy import stats @@ -176,8 +174,6 @@ def calculate_scaleheight_from_map(mass_map, img_size): z_1D = np.ravel(z) S_1D = np.ravel(mass_map) - p0 = (mass_map.max(), 1.0, 0.0) - try: popt, pcov = curve_fit( exponential, @@ -358,7 +354,9 @@ def get_axis_lengths_tensor( galaxy_log.debug("Total weight of 0, so not calculating axis lengths.") return unyt.unyt_array(np.zeros(3), position.units), np.zeros(3) - Itensor = (weight[:, None, None] / weight.sum()) * np.ones((weight.shape[0], 3, 3)) + Itensor = ( + (weight[:, None, None] / weight.sum()) * np.ones((weight.shape[0], 3, 3)) + ).value # Note: unyt currently ignores the position units in the *= # i.e. Itensor is dimensionless throughout (even though it should not be) for i in range(3): @@ -441,9 +439,9 @@ def get_axis_lengths_tensor( ) break - Itensor = (weight[:, None, None] / weight.sum()) * np.ones( - (weight.shape[0], 3, 3) - ) + Itensor = ( + (weight[:, None, None] / weight.sum()) * np.ones((weight.shape[0], 3, 3)) + ).value # Note: unyt currently ignores the position units in the *= # i.e. Itensor is dimensionless throughout (even though it should not be) for i in range(3): @@ -512,7 +510,6 @@ def get_kappa_corot( partdata: "SWIFTParticleDataset", half_mass_radius: unyt.unyt_quantity, R200crit: unyt.unyt_quantity, - Rvir: unyt.unyt_quantity, orientation_type: str, orientation_vector: NDArray[float], mass_variable: str = "masses", @@ -536,7 +533,7 @@ def get_kappa_corot( radius = np.sqrt((position ** 2).sum(axis=1)) mask = get_orientation_mask( - radius, half_mass_radius, R200crit, Rvir, inner_aperture, outer_aperture + radius, half_mass_radius, R200crit, inner_aperture, outer_aperture ) position = position[mask] @@ -544,14 +541,15 @@ def get_kappa_corot( mass = mass[mask] radius = radius[mask] + print(f"{position=}") + print(f"{velocity=}") + K = 0.5 * (mass[:, None] * velocity ** 2).sum() if K == 0.0: return 0.0, 0.0 # np.cross does not preserve units, so we need to multiply them back in - angular_momentum = ( - (mass[:, None] * np.cross(position, velocity)) * position.units * velocity.units - ) + angular_momentum = mass[:, None] * np.cross(position, velocity) Lz = (angular_momentum * orientation_vector[None, :]).sum(axis=1) rdotL = (position * orientation_vector[None, :]).sum(axis=1) R2 = radius ** 2 - rdotL ** 2 diff --git a/morpholopy/orientation.py b/morpholopy/orientation.py index 1638b16..dc7098a 100644 --- a/morpholopy/orientation.py +++ b/morpholopy/orientation.py @@ -11,7 +11,7 @@ where: - component type can be stars, gas, ISM, HI, baryons - inner mask radius can be 0xR0.5, 0.5R0.5 - - outer mask radius can be R0.5, 2xR0.5, 4xR0.5, 50kpc, R200crit, 0.1Rvir + - outer mask radius can be R0.5, 2xR0.5, 4xR0.5, 50kpc, R200crit, 0.1R200crit - sigma clipping can be 0sigma (no clipping), 1sigma, 2sigma This file contains a number of functions that can be called in other @@ -36,7 +36,6 @@ def get_orientation_mask( radius: unyt.unyt_array, half_mass_radius: unyt.unyt_quantity, R200crit: unyt.unyt_quantity, - Rvir: unyt.unyt_quantity, inner_aperture: str, outer_aperture: str, ) -> NDArray[bool]: @@ -67,8 +66,8 @@ def get_orientation_mask( ) elif outer_aperture == "R200crit": mask &= radius < R200crit - elif outer_aperture == "0.1Rvir": - mask &= radius < 0.1 * Rvir + elif outer_aperture == "0.1R200crit": + mask &= radius < 0.1 * R200crit else: raise RuntimeError(f"Unknown outer aperture: {outer_aperture}!") @@ -145,7 +144,6 @@ def get_mass_position_velocity( data: SWIFTDataset, half_mass_radius: unyt.unyt_quantity, R200crit: unyt.unyt_quantity, - Rvir: unyt.unyt_quantity, orientation_type: str, ) -> Tuple[unyt.unyt_array, unyt.unyt_array, unyt.unyt_array]: """ @@ -161,7 +159,7 @@ def get_mass_position_velocity( radius = np.sqrt((position ** 2).sum(axis=1)) mask = get_orientation_mask( - radius, half_mass_radius, R200crit, Rvir, inner_aperture, outer_aperture + radius, half_mass_radius, R200crit, inner_aperture, outer_aperture ) position = position[mask] @@ -196,7 +194,6 @@ def get_orientation_matrices( data: SWIFTDataset, half_mass_radius: unyt.unyt_quantity, R200crit: unyt.unyt_quantity, - Rvir: unyt.unyt_quantity, orientation_type: str, ) -> Tuple[NDArray[float], NDArray[float], NDArray[float]]: """ @@ -210,7 +207,7 @@ def get_orientation_matrices( """ mass, position, velocity = get_mass_position_velocity( - data, half_mass_radius, R200crit, Rvir, orientation_type + data, half_mass_radius, R200crit, orientation_type ) angular_momentum = (mass[:, None] * np.cross(position, velocity)).sum(axis=0) diff --git a/requirements.txt b/requirements.txt index 34b0036..a303a7d 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ -velociraptor>=0.15.8 +velociraptor>=0.16.0 swiftsimio>=6.0.0 scipy>=1.9.3 matplotlib>=3.6.2 diff --git a/velociraptor-comparison-data b/velociraptor-comparison-data index 9dc4cb3..ec6e35f 160000 --- a/velociraptor-comparison-data +++ b/velociraptor-comparison-data @@ -1 +1 @@ -Subproject commit 9dc4cb3a2b30736a213d116f7dd6fc91065c96f4 +Subproject commit ec6e35f4ff5b163b67b6718ef2906cc435afeafc From f935c2fd3cf9b73d6b7b0ba6104c4f992c51a1ed Mon Sep 17 00:00:00 2001 From: EvgeniiChaikin Date: Thu, 2 May 2024 13:06:10 +0200 Subject: [PATCH 2/4] Remove debug prints and add try-except to prevent possible unyt errors --- morpholopy/galaxy_data.py | 4 ---- morpholopy/morphology.py | 25 ++++++++++++++++--------- 2 files changed, 16 insertions(+), 13 deletions(-) diff --git a/morpholopy/galaxy_data.py b/morpholopy/galaxy_data.py index 2faa515..03ac19f 100644 --- a/morpholopy/galaxy_data.py +++ b/morpholopy/galaxy_data.py @@ -481,9 +481,6 @@ def process_galaxy( mask = MaskTuple(**particle_name_masks) membership_data.close() - for key, val in particle_name_masks.items(): - print(key, np.sum(val)) - # mask out all particles that are actually bound to the galaxy # while at it: convert everything to physical coordinates for parttype in ["gas", "dark_matter", "stars", "black_holes"]: @@ -648,7 +645,6 @@ def process_galaxy( # - axis ratios and angular momenta (a, b, c), z_axis = get_axis_lengths_reduced_tensor(galaxy_log, data.stars, Rhalf) - print(a, b, c, z_axis) if (a > 0.0) and (b > 0.0) and (c > 0.0): galaxy_data["stars_axis_ca_reduced"] = c / a galaxy_data["stars_axis_cb_reduced"] = c / b diff --git a/morpholopy/morphology.py b/morpholopy/morphology.py index 524d5a1..0c70473 100644 --- a/morpholopy/morphology.py +++ b/morpholopy/morphology.py @@ -541,25 +541,32 @@ def get_kappa_corot( mass = mass[mask] radius = radius[mask] - print(f"{position=}") - print(f"{velocity=}") - K = 0.5 * (mass[:, None] * velocity ** 2).sum() if K == 0.0: return 0.0, 0.0 - # np.cross does not preserve units, so we need to multiply them back in - angular_momentum = mass[:, None] * np.cross(position, velocity) + # np.cross may not preserve units, so we need to multiply them back in + try: + angular_momentum = mass[:, None] * np.cross(position, velocity) + j = angular_momentum.sum(axis=0) / mass.sum() + j.convert_to_units("kpc*km/s") + except unyt.exceptions.UnitConversionError: + angular_momentum = ( + mass[:, None] + * np.cross(position, velocity) + * position.units + * velocity.units + ) + j = angular_momentum.sum(axis=0) / mass.sum() + j.convert_to_units("kpc*km/s") + j = np.sqrt((j ** 2).sum()) + Lz = (angular_momentum * orientation_vector[None, :]).sum(axis=1) rdotL = (position * orientation_vector[None, :]).sum(axis=1) R2 = radius ** 2 - rdotL ** 2 mask = (Lz > 0.0) & (R2 > 0.0) Kcorot = 0.5 * (Lz[mask] ** 2 / (mass[mask] * R2[mask])).sum() - j = angular_momentum.sum(axis=0) / mass.sum() - j.convert_to_units("kpc*km/s") - j = np.sqrt((j ** 2).sum()) - return j, Kcorot / K From 9026cf4c4c02dbbbe7fca4649d3759edeebd541a Mon Sep 17 00:00:00 2001 From: EvgeniiChaikin Date: Thu, 2 May 2024 13:08:42 +0200 Subject: [PATCH 3/4] Remove commeted-out, outdated code --- morpholopy/galaxy_data.py | 9 --------- 1 file changed, 9 deletions(-) diff --git a/morpholopy/galaxy_data.py b/morpholopy/galaxy_data.py index 03ac19f..a9b6b9b 100644 --- a/morpholopy/galaxy_data.py +++ b/morpholopy/galaxy_data.py @@ -437,15 +437,6 @@ def process_galaxy( catalogue = load_catalogue(catalogue_filename, disregard_units=True) membership_data = h5.File(halo_membership_filename, "r") - # get the particles belonging to this galaxy - # particles, _ = groups.extract_halo(halo_index=galaxy_index) - # turn this information into a swiftsimio mask and read the data - # using swiftsimio - - ##data, mask = to_swiftsimio_dataset( - ## particles, snapshot_filename, generate_extra_mask=True - ##) - swift_mask = swiftsimio.mask(snapshot_filename, spatial_only=True) # SWIFT data is stored in comoving units, so we need to un-correct # the velociraptor data if it is stored in physical. From 04046f5f408500cec14ecd9b7e95f5095490d73c Mon Sep 17 00:00:00 2001 From: EvgeniiChaikin Date: Mon, 6 May 2024 14:15:03 +0200 Subject: [PATCH 4/4] Changes following review --- morpholopy/morphology.py | 20 ++++---------------- requirements.txt | 2 ++ 2 files changed, 6 insertions(+), 16 deletions(-) diff --git a/morpholopy/morphology.py b/morpholopy/morphology.py index 0c70473..b56be8d 100644 --- a/morpholopy/morphology.py +++ b/morpholopy/morphology.py @@ -357,8 +357,6 @@ def get_axis_lengths_tensor( Itensor = ( (weight[:, None, None] / weight.sum()) * np.ones((weight.shape[0], 3, 3)) ).value - # Note: unyt currently ignores the position units in the *= - # i.e. Itensor is dimensionless throughout (even though it should not be) for i in range(3): for j in range(3): Itensor[:, i, j] *= position[:, i] * position[:, j] @@ -545,20 +543,10 @@ def get_kappa_corot( if K == 0.0: return 0.0, 0.0 - # np.cross may not preserve units, so we need to multiply them back in - try: - angular_momentum = mass[:, None] * np.cross(position, velocity) - j = angular_momentum.sum(axis=0) / mass.sum() - j.convert_to_units("kpc*km/s") - except unyt.exceptions.UnitConversionError: - angular_momentum = ( - mass[:, None] - * np.cross(position, velocity) - * position.units - * velocity.units - ) - j = angular_momentum.sum(axis=0) / mass.sum() - j.convert_to_units("kpc*km/s") + angular_momentum = mass[:, None] * np.cross(position, velocity) + j = angular_momentum.sum(axis=0) / mass.sum() + j.convert_to_units("kpc*km/s") + j = np.sqrt((j ** 2).sum()) Lz = (angular_momentum * orientation_vector[None, :]).sum(axis=1) diff --git a/requirements.txt b/requirements.txt index a303a7d..d7affa9 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,3 +3,5 @@ swiftsimio>=6.0.0 scipy>=1.9.3 matplotlib>=3.6.2 swiftpipeline>=0.3.5 +h5py>=3.9.0 +unyt>=3.0.0 \ No newline at end of file