Added NH2 vs integrated intensities plots and datasets. (b3296900) · Commits · Andrea Giannetti / swiss_army_knife_stable

etl/prs/prs_compute_integrated_fluxes_and_ratios.py

+6 −2

Original line number	Diff line number	Diff line
		@@ -3,6 +3,7 @@ import sqlalchemy
		import numpy as np
		from datetime import datetime
		from astropy.io import fits
		from astropy import constants
		from typing import Union, List
		from assets.commons import (validate_parameter,
		upsert,
		@@ -121,9 +122,12 @@ def compute_moment_zero(cube: str,
		fitsfile = open_fits_file_duck_typing(fitsfile=cube, fits_path=_cube_path)
		header = fitsfile[hdu_idx].header.copy()
		data = fitsfile[hdu_idx].data
		mom0 = data.sum(axis=0) * abs(header['CDELT3'])
		conversion = header['CUNIT3'].strip() == 'HZ'
		spectral_unit = "km/s" if conversion is True else header["CUNIT3"].strip()
		conversion_factor = 1 if conversion is False else (-header['CDELT3'] / header['CRVAL3'] * constants.c).to('km/s').value
		mom0 = (data.sum(axis=0) * abs(conversion_factor))
		keywords_to_delete = ['NAXIS3', 'CRPIX3', 'CDELT3', 'CRVAL3', 'CUNIT3', 'CTYPE3']
		header['BUNIT'] += f' {header["CUNIT3"].replace(" ", "")}'
		header['BUNIT'] += f' {spectral_unit}'
		header['BTYPE'] = 'MOMENT0'
		header['NAXIS'] = 2
		for key in keywords_to_delete:

etl/prs/prs_inspect_results.py

+83 −18

Original line number	Diff line number	Diff line
		@@ -4,10 +4,11 @@ import os
		import seaborn as sns
		import pandas as pd
		import xarray as xr
		from typing import Union, Tuple
		from typing import Union, Tuple, List
		from contextlib import closing
		from itertools import product
		from astropy.io import fits
		from astropy import units as u
		from sqlalchemy.orm import Session
		from sqlalchemy import and_, or_
		from sqlalchemy import engine as sqla_engine
		@@ -64,18 +65,18 @@ def get_results(
		lines: Union[list, tuple],
		run_id: str,
		session: Session,
		is_isothermal: bool = False) -> Tuple[float, str, str, tuple, np.array]:
		is_isothermal: bool = False) -> Tuple[float, str, str, float, str, str, tuple, np.array]:
		"""
		Get results from the database, given the parameters of the model
		:param dust_temperature: the characteristic dust temperature
		:param gas_density: the characteristic number density of molecualar hydrogen
		:param gas_density: the characteristic number density of molecular hydrogen
		:param lines: the lines used to compute the ratio
		:param run_id: the run_id of the model
		:param session: a session to query the database
		:param is_isothermal: whether the model is isothermal
		:return: a tuple containing the aggregated line ratios, the ratio map name, the fits grid name (for the number
		density), the scaling factor between the grid and image pixel (must be integer!), and the pixel-by-pixel
		measured line ratios
		:return: a tuple containing the aggregated line ratios, the moment zero maps names, the pixel size, the ratio map
		name, the fits grid name (for the number density), the scaling factor between the grid and image pixel (must be
		integer!), and the pixel-by-pixel measured line ratios
		"""
		# Select correct column for filtering
		if is_isothermal is True:
		@@ -104,6 +105,9 @@ def get_results(
		ratio_values = ratio_fitsfile[0].data
		assert (results[3].npix % results[0].grid_shape_1 == 0) and (results[3].npix % results[0].grid_shape_2 == 0)
		return results[1].aggregated_ratio, \
		results[1].mom_zero_name_1, \
		results[1].mom_zero_name_2, \
		(results[0].grid_size_3 * u.pc).to(u.cm).value / results[0].grid_shape_3, \
		results[1].ratio_map_name, \
		results[2].fits_grid_name, \
		(int(results[3].npix / results[0].grid_shape_1), int(results[3].npix / results[0].grid_shape_2)), \
		@@ -146,6 +150,25 @@ def get_grid_values(
		return data


		def get_column_density_vs_mom0(volume_density_grid: np.array,
		mom_zero_name_1: str,
		mom_zero_name_2: str,
		grid_pixel_size: float,
		lines: Tuple[str, str],
		gas_density: float,
		dust_temperature: float,
		zoom_ratios: Tuple[int, int]) -> pd.DataFrame:
		column_density_map = np.nansum(volume_density_grid * grid_pixel_size, axis=2).flatten()
		df_results = pd.DataFrame(data=column_density_map, columns=['H2_column_density'])
		df_results['nh2'] = gas_density
		df_results['tdust'] = dust_temperature
		with fits.open(os.path.join('prs', 'fits', 'moments', mom_zero_name_1)) as fitsfile:
		df_results[f'mom_zero_{lines[0]}'] = fitsfile[0].data[::zoom_ratios[0], ::zoom_ratios[1]].flatten()
		with fits.open(os.path.join('prs', 'fits', 'moments', mom_zero_name_2)) as fitsfile:
		df_results[f'mom_zero_{lines[1]}'] = fitsfile[0].data[::zoom_ratios[0], ::zoom_ratios[1]].flatten()
		return df_results


		def main(run_id: str,
		is_isothermal: bool,
		engine: sqla_engine = None):
		@@ -167,9 +190,9 @@ def main(run_id: str,
		'gas_density': central_densities
		})
		results_dict = {}

		with closing(Session(engine)) as session:
		for lines in line_pairs:
		df_coldens_mom0_list = []
		results_dict[f'{"-".join(lines)}'] = {}
		for (tdust, nh2) in product(dust_temperatures, central_densities):
		aggregated_ratio = get_aggregated_ratio_from_db(dust_temperature=tdust,
		@@ -179,7 +202,14 @@ def main(run_id: str,
		run_id=run_id)
		results.loc[tdust, nh2] = aggregated_ratio
		logger.debug(f'The aggregated ratio for lines {lines}, using {nh2}, {tdust} is: {aggregated_ratio}')
		ratio, ratio_fitsname, grid_fitsname, zoom_ratios, ratio_values = \
		ratio, \
		mom0_name_1, \
		mom0_name_2, \
		grid_pixel_size, \
		ratio_fitsname, \
		grid_fitsname, \
		zoom_ratios, \
		ratio_values = \
		get_results(
		dust_temperature=tdust,
		gas_density=nh2,
		@@ -214,20 +244,27 @@ def main(run_id: str,
		ratio_values[::zoom_ratios[0], ::zoom_ratios[1]].flatten()])
		results_dict[f'{"-".join(lines)}'][f'{"_".join([str(tdust), str(nh2)])}'] = correlation_data
		plt.scatter(correlation_data[0], correlation_data[2])
		df_coldens_mom0_list.append(get_column_density_vs_mom0(volume_density_grid=density_grid,
		mom_zero_name_1=mom0_name_1,
		mom_zero_name_2=mom0_name_2,
		lines=lines,
		gas_density=nh2,
		dust_temperature=tdust,
		zoom_ratios=zoom_ratios,
		grid_pixel_size=grid_pixel_size))
		plot_ratio_vs_density(lines=lines,
		results=results)

		df_coldens_mom0 = plot_coldens_vs_integrated_intensity(df_coldens_mom0_list=df_coldens_mom0_list,
		lines=lines)
		df_coldens_mom0.to_csv(f'full_dataset_NH2_lines_{"-".join(lines)}.csv')

		plt.semilogx()
		plt.xlabel('avg(gas_density)')
		plt.ylabel('simulated ratio')
		plt.savefig(os.path.join('prs', 'output', f'ratio_vs_avg_density_los_{"-".join(lines)}.png'))
		plt.clf()
		results.plot(x='dust_temperature', y='gas_density', yscale='log')
		plt.savefig(os.path.join('prs', 'output', f'ratio_grid_lines_{"-".join(lines)}.png'))
		plt.clf()
		df_list_concat = []
		for (tdust, nh2) in product(dust_temperatures, central_densities):
		df_list = []
		for lines in line_pairs:
		df = pd.DataFrame(data=results_dict[f'{"-".join(lines)}'][f'{"_".join([str(tdust), str(nh2)])}'].transpose(),
		df = pd.DataFrame(
		data=results_dict[f'{"-".join(lines)}'][f'{"_".join([str(tdust), str(nh2)])}'].transpose(),
		columns=['avg_nh2', 'avg_tdust', f'ratio_{"-".join(lines)}'])
		df_list.append(df.dropna())
		df_tmp = pd.concat(df_list, axis=1)
		@@ -243,6 +280,34 @@ def main(run_id: str,
		df_all.to_csv('full_dataset.csv')


		def plot_coldens_vs_integrated_intensity(df_coldens_mom0_list: List[pd.DataFrame],
		lines=Tuple[str, str]) -> pd.DataFrame:
		df_coldens_mom0 = pd.concat(df_coldens_mom0_list)
		plt.xlabel('H_2 column density)')
		plt.ylabel('Moment 0')
		plt.loglog()
		plt.scatter(df_coldens_mom0['H2_column_density'], df_coldens_mom0[f'mom_zero_{lines[0]}'],
		label=f'mom_zero_{lines[0]}')
		plt.scatter(df_coldens_mom0['H2_column_density'], df_coldens_mom0[f'mom_zero_{lines[1]}'],
		label=f'mom_zero_{lines[1]}')
		plt.legend()
		plt.savefig(os.path.join('prs', 'output', f'coldens_moments_lines_{"-".join(lines)}.png'))
		plt.clf()
		return df_coldens_mom0


		def plot_ratio_vs_density(lines: Tuple[str, str],
		results: xr.DataArray):
		plt.semilogx()
		plt.xlabel('avg(gas_density)')
		plt.ylabel('simulated ratio')
		plt.savefig(os.path.join('prs', 'output', f'ratio_vs_avg_density_los_{"-".join(lines)}.png'))
		plt.clf()
		results.plot(x='dust_temperature', y='gas_density', yscale='log')
		plt.savefig(os.path.join('prs', 'output', f'ratio_grid_lines_{"-".join(lines)}.png'))
		plt.clf()


		if __name__ == '__main__':
		main(run_id='0d0761f8-fcf3-4293-9d9d-ba1a37f79a19', is_isothermal=False)
		# main(run_id='ba7fd3ed-5947-4dc6-bcef-38f151f19b77', is_isothermal=False)

etl/prs/prs_prepare_backup.py

+6 −2

Original line number	Diff line number	Diff line
		import glob
		import os
		import shutil
		import pandas as pd
		@@ -103,8 +104,11 @@ def main(run_id: str):
		for object_category in objects_to_save:
		os.rmdir(os.path.join(backup_dir, run_id, object_category))

		shutil.copyfile('full_dataset.csv',
		os.path.join(backup_dir, run_id, 'full_dataset.csv'))
		full_datasets_list = glob.glob('full_dataset_NH2*.csv')
		full_datasets_list.append('full_dataset.csv')
		for dataset_file in full_datasets_list:
		shutil.copyfile(dataset_file,
		os.path.join(backup_dir, run_id, dataset_file))
		save_output_images(backup_dir=backup_dir,
		run_id=run_id)