From d3915b370672702fb55608a9d4529093cca6e9b6 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 21 May 2024 17:11:09 +0200 Subject: [PATCH] replace src/analysis with execution from src/lib/plots and move unused scripts to utils --- src/lib/plots.py | 37 +++++++++++++++----- {src => utils}/analysis.py | 4 +-- {src => utils}/comparison_Kishimoto.py | 10 +++--- utils/get_cdelt.py | 47 ++++++++++++++++++++++++++ 4 files changed, 82 insertions(+), 16 deletions(-) rename {src => utils}/analysis.py (95%) rename {src => utils}/comparison_Kishimoto.py (98%) create mode 100755 utils/get_cdelt.py diff --git a/src/lib/plots.py b/src/lib/plots.py index e357348..3386c41 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -1,3 +1,4 @@ +#!/usr/bin/python """ Library functions for displaying informations using matplotlib @@ -54,7 +55,7 @@ from astropy.wcs import WCS from astropy.io import fits from astropy.coordinates import SkyCoord from scipy.ndimage import zoom as sc_zoom -from lib.utils import rot2D, princ_angle, sci_not +from utils import rot2D, princ_angle, sci_not def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="", **kwargs): @@ -423,7 +424,7 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c X, Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) U, V = poldata*np.cos(np.pi/2.+pangdata*np.pi/180.), poldata*np.sin(np.pi/2.+pangdata*np.pi/180.) ax.quiver(X[::step_vec, ::step_vec], Y[::step_vec, ::step_vec], U[::step_vec, ::step_vec], V[::step_vec, ::step_vec], units='xy', angles='uv', - scale=1./vec_scale, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.2, linewidth=0.3, color='w', edgecolor='k') + scale=1./vec_scale, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='w', edgecolor='k') pol_sc = AnchoredSizeBar(ax.transData, vec_scale, r"$P$= 100 %", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w') ax.add_artist(pol_sc) @@ -734,7 +735,7 @@ class overplot_radio(align_maps): self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) self.U, self.V = pol*np.cos(np.pi/2.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.) self.Q = self.ax_overplot.quiver(self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], self.U[::step_vec, ::step_vec], self.V[::step_vec, ::step_vec], units='xy', angles='uv', scale=1./self.vec_scale, - scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.1, linewidth=0.5, color='white', edgecolor='black', label="{0:s} polarisation map".format(self.map_observer)) + scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='white', edgecolor='black', label="{0:s} polarisation map".format(self.map_observer)) self.ax_overplot.autoscale(False) # Display other map as contours @@ -853,7 +854,7 @@ class overplot_chandra(align_maps): self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) self.U, self.V = pol*np.cos(np.pi/2.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.) self.Q = self.ax_overplot.quiver(self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], self.U[::step_vec, ::step_vec], self.V[::step_vec, ::step_vec], units='xy', angles='uv', scale=1./self.vec_scale, - scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.1, linewidth=0.5, color='white', edgecolor='black', label="{0:s} polarisation map".format(self.map_observer)) + scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='white', edgecolor='black', label="{0:s} polarisation map".format(self.map_observer)) self.ax_overplot.autoscale(False) # Display other map as contours @@ -971,8 +972,8 @@ class overplot_pol(align_maps): px_scale = self.other_wcs.wcs.get_cdelt()[0]/self.wcs_UV.wcs.get_cdelt()[0] self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) self.U, self.V = pol*np.cos(np.pi/2.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.) - self.Q = self.ax_overplot.quiver(self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], self.U[::step_vec, ::step_vec], self.V[::step_vec, ::step_vec], units='xy', angles='uv', scale=px_scale/self.vec_scale, scale_units='xy', pivot='mid', - headwidth=0., headlength=0., headaxislength=0., width=0.1/px_scale, linewidth=0.5, color='white', edgecolor='black', transform=self.ax_overplot.get_transform(self.wcs_UV), label="{0:s} polarisation map".format(self.map_observer)) + self.Q = self.ax_overplot.quiver(self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], self.U[::step_vec, ::step_vec], self.V[::step_vec, ::step_vec], units='xy', angles='uv', scale=1./self.vec_scale, scale_units='xy', pivot='mid', + headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='white', edgecolor='black', transform=self.ax_overplot.get_transform(self.wcs_UV), label="{0:s} polarisation map".format(self.map_observer)) # Display Stokes I as contours if levels is None: @@ -1128,7 +1129,7 @@ class align_pol(object): X, Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) U, V = pol*np.cos(np.pi/2.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.) ax.quiver(X[::step_vec, ::step_vec], Y[::step_vec, ::step_vec], U[::step_vec, ::step_vec], V[::step_vec, ::step_vec], units='xy', - angles='uv', scale=0.5, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.1, color='w') + angles='uv', scale=0.5, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='w') pol_sc = AnchoredSizeBar(ax.transData, 2., r"$P$= 100 %", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w') ax.add_artist(pol_sc) @@ -2349,12 +2350,12 @@ class pol_map(object): if hasattr(self, 'quiver'): self.quiver.remove() self.quiver = ax.quiver(X, Y, XY_U, XY_V, units='xy', scale=1./self.vec_scale, scale_units='xy', pivot='mid', headwidth=0., - headlength=0., headaxislength=0., width=0.2, linewidth=0.3, color='white', edgecolor='black') + headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='white', edgecolor='black') fig.canvas.draw_idle() return self.quiver else: ax.quiver(X, Y, XY_U, XY_V, units='xy', scale=1./self.vec_scale, scale_units='xy', pivot='mid', headwidth=0., - headlength=0., headaxislength=0., width=0.2, linewidth=0.3, color='white', edgecolor='black') + headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='white', edgecolor='black') fig.canvas.draw_idle() def pol_int(self, fig=None, ax=None): @@ -2463,3 +2464,21 @@ class pol_map(object): if self.region is not None: ax.contour(self.region.astype(float), levels=[0.5], colors='white', linewidths=0.8) fig.canvas.draw_idle() + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser(description='Interactively plot the pipeline products') + parser.add_argument('-f', '--file', metavar='path', required=False, help='the full or relative path to the data product', type=str, default=None) + parser.add_argument('-p', '--snrp', metavar='snrp_cut', required=False, help='the cut in signal-to-noise for the polarisation degree', type=float, default=3.) + parser.add_argument('-i', '--snri', metavar='snri_cut', required=False, help='the cut in signal-to-noise for the intensity', type=float, default=3.) + parser.add_argument('-l', '--lim', metavar='flux_lim', nargs=2, required=False, help='limits for the intensity map', default=None) + args = parser.parse_args() + + if args.file is not None: + Stokes_UV = fits.open(args.file, mode='readonly') + p = pol_map(Stokes_UV, SNRp_cut=args.snrp, SNRi_cut=args.snri, flux_lim=args.lim) + + else: + print("python3 plots.py -f -p -i -l ") diff --git a/src/analysis.py b/utils/analysis.py similarity index 95% rename from src/analysis.py rename to utils/analysis.py index e2d2fea..6daecb8 100755 --- a/src/analysis.py +++ b/utils/analysis.py @@ -1,4 +1,4 @@ -#!/usr/bin/python3 +#!/usr/bin/python from getopt import getopt, error as get_error from sys import argv @@ -30,7 +30,7 @@ except get_error as err: if fits_path is not None: from astropy.io import fits - from lib.plots import pol_map + from src.lib.plots import pol_map Stokes_UV = fits.open(fits_path) p = pol_map(Stokes_UV, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim) diff --git a/src/comparison_Kishimoto.py b/utils/comparison_Kishimoto.py similarity index 98% rename from src/comparison_Kishimoto.py rename to utils/comparison_Kishimoto.py index 56b67ae..75b7073 100755 --- a/src/comparison_Kishimoto.py +++ b/utils/comparison_Kishimoto.py @@ -1,8 +1,8 @@ -# !/usr/bin/env python -from lib.background import gauss, bin_centers -from lib.deconvolve import zeropad -from lib.reduction import align_data -from lib.plots import princ_angle +#!/usr/bin/python +from src.lib.background import gauss, bin_centers +from src.lib.deconvolve import zeropad +from src.lib.reduction import align_data +from src.lib.plots import princ_angle from matplotlib.colors import LogNorm from os.path import join as path_join from astropy.io import fits diff --git a/utils/get_cdelt.py b/utils/get_cdelt.py new file mode 100755 index 0000000..45e526b --- /dev/null +++ b/utils/get_cdelt.py @@ -0,0 +1,47 @@ +#!/usr/bin/python + +def main(infiles=None): + """ + Retrieve native spatial resolution from given observation. + """ + from os.path import join as path_join + from warnings import catch_warnings, filterwarnings + from astropy.io.fits import getheader + from astropy.wcs import WCS, FITSFixedWarning + from numpy.linalg import eig + + if infiles is None: + print("Usage: \"python get_cdelt.py -f infiles\"") + return 1 + prod = [["/".join(filepath.split('/')[:-1]), filepath.split('/')[-1]] for filepath in infiles] + data_folder = prod[0][0] + infiles = [p[1] for p in prod] + + cdelt = {} + size = {} + for currfile in infiles: + with catch_warnings(): + filterwarnings('ignore', message="'datfix' made the change", category=FITSFixedWarning) + wcs = WCS(getheader(path_join(data_folder, currfile))).celestial + key = currfile[:-5] + size[key] = wcs.array_shape + if wcs.wcs.has_cd(): + cdelt[key] = eig(wcs.wcs.cd)[0]*3600. + else: + cdelt[key] = wcs.wcs.cdelt*3600. + + print("Image name, native resolution in arcsec and shape") + for currfile in infiles: + key = currfile[:-5] + print(key, cdelt[key], size[key]) + + return 0 + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser(description='Query MAST for target products') + parser.add_argument('-f', '--files', metavar='path', required=False, nargs='*', help='the full or relative path to the data products', default=None) + args = parser.parse_args() + exitcode = main(infiles=args.files)