From 0e65ff9a46d7b1e2a528c74ad0a346948f2275ea Mon Sep 17 00:00:00 2001 From: ttakekoshi Date: Wed, 23 Aug 2023 16:53:45 +0900 Subject: [PATCH] update pointing (yaml, function, .py) and make test set --- pipelines/beammap.py | 31 ++++------ pipelines/pointing.py | 112 ++++++++++++++++++++++++----------- pipelines/utils/functions.py | 32 ++++++---- yaml/beammap_params.yaml | 5 +- yaml/pointing_params.yaml | 8 +-- 5 files changed, 117 insertions(+), 71 deletions(-) diff --git a/pipelines/beammap.py b/pipelines/beammap.py index a931a7a..9549bc6 100644 --- a/pipelines/beammap.py +++ b/pipelines/beammap.py @@ -15,7 +15,6 @@ import astropy.units as u import aplpy from tqdm import tqdm - # original package from utils import functions as fc @@ -27,13 +26,13 @@ # command line arguments parser = argparse.ArgumentParser() parser.add_argument("dfits_file", help="DFITS name") -parser.add_argument("flux_file", help="flux text") parser.add_argument("yaml_file", help="parameter file") +parser.add_argument("flux_file", help="flux list file made by planetFlux") args = parser.parse_args() dfits_file = pathlib.Path(args.dfits_file) -flux_file = pathlib.Path(args.flux_file) yaml_file = pathlib.Path(args.yaml_file) +flux_file = pathlib.Path(args.flux_file) with open(yaml_file) as file: params = yaml.load(file, Loader=yaml.SafeLoader) @@ -337,26 +336,22 @@ print("#6: planet properties") plt.style.use("seaborn-darkgrid") - -planet_data = np.loadtxt(flux_file, comments="#", usecols=(1, 2, 3, 4, 5, 6)) -# dtable = table.Table(planet_data) -istart = params["planet"]["istart"] -iend = params["planet"]["iend"] - +planet_data = np.loadtxt(flux_file, comments="#", usecols=(1, 2, 3, 4)) fig, ax = plt.subplots(2, 1, figsize=(10, 5), dpi=dpi) -x = np.linspace(300, 399, 100) # frequency -ll = interp1d(planet_data[istart:iend, 0], planet_data[istart:iend, 1], kind="cubic") -theta_s = planet_data[istart:iend, 1].mean() +#x = np.linspace(300, 399, 100) # frequency +x = planet_data[:,0] # frequency +ll = interp1d(planet_data[:, 0], planet_data[:, 1], kind="cubic") +theta_s = planet_data[:, 1].mean() -ax[0].plot(planet_data[istart:iend, 0], planet_data[istart:iend, 1], ".") -ax[0].plot(planet_data[istart:iend, 0], ll(x), "k-") +ax[0].plot(planet_data[:, 0], planet_data[:, 1], ".") +ax[0].plot(planet_data[:, 0], ll(x), "k-") ax[0].set_xlabel("frequency [GHz]") ax[0].set_ylabel("angular diameter [arcsec]") -mm = interp1d(planet_data[istart:iend, 0], planet_data[istart:iend, 2], kind="cubic") +mm = interp1d(planet_data[:, 0], planet_data[:, 2], kind="cubic") -ax[1].plot(planet_data[istart:iend, 0], planet_data[istart:iend, 2], ".") -ax[1].plot(planet_data[istart:iend, 0], mm(x), "k-") +ax[1].plot(planet_data[:, 0], planet_data[:, 2], ".") +ax[1].plot(planet_data[:, 0], mm(x), "k-") ax[1].set_xlabel("frequency [GHz]") ax[1].set_ylabel("Tb [K]") @@ -471,7 +466,7 @@ label="beam efficiency", ) ax.axhline(eta_med, color="C1", label=f"median: {eta_med:.2f}") -ax.set_ylim(0, 1) +ax.set_ylim(0, 1.2) ax.set_xlabel("frequency [GHz]") ax.set_ylabel(r"$\eta_{mb}$") ax.legend() diff --git a/pipelines/pointing.py b/pipelines/pointing.py index 01dab0b..2198c12 100644 --- a/pipelines/pointing.py +++ b/pipelines/pointing.py @@ -148,8 +148,8 @@ # 3rd step: make cube/continuum print("#3: make cube/continuum") +scanarray_cal.kidtp[list(params["imaging"]["exchs"])] = -1 -scanarray_cal.kidtp[[16, 18, 44, 46]] = -1 scanarray_cal = scanarray_cal.where(scanarray_cal.kidtp == 1, drop=True) scanarray_cal = scanarray_cal[:-20000] @@ -160,16 +160,13 @@ ymin = scanarray_cal.y.min().values ymax = scanarray_cal.y.max().values -cube_array = dc.tocube( - scanarray_cal, gx=gx, gy=gy, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax -) +cube_array = dc.tocube( scanarray_cal, gx=gx, gy=gy, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) dc.io.savefits(cube_array, cube_obs_fits, overwrite=True) -# exchs = params["imaging"]["exchs"] -# mask = np.full_like(scanarray_cal.kidid.values, True, dtype=np.bool) -# mask[exchs] = False -# mask[np.where(scanarray_cal.kidtp != 1)] = False -# masked_cube_array = cube_array[:, :, mask] +# exchs = params["imaging"]["exchs"] mask = +# np.full_like(scanarray_cal.kidid.values, True, dtype=np.bool) mask[exchs] = +# False mask[np.where(scanarray_cal.kidtp != 1)] = False masked_cube_array = +# cube_array[:, :, mask] # weight = dc.ones_like(masked_cube_array) weight = dc.ones_like(cube_array) @@ -178,19 +175,21 @@ dc.io.savefits(cont_array, cont_obs_fits, dropdeg=True, overwrite=True) # fits -fig = plt.figure(figsize=(5, 5)) -ax = aplpy.FITSFigure(str(cont_obs_fits), figure=fig, subplot=(1, 1, 1)) - -ax.show_colorscale(cmap="viridis", stretch="linear") -ax.add_colorbar(width=0.15) - -fig.tight_layout() -fig.savefig(output_dir / f"continuum_image.{image_format}") -if do_plot: - plt.show() -else: - plt.clf() - plt.close() +#fig = plt.figure(figsize=(5, 5)) +#ax = aplpy.FITSFigure(str(cont_obs_fits), +# figure=fig, subplot=(1, 1, 1)) +# +#ax.show_colorscale(cmap="viridis", stretch="linear") +#ax.add_colorbar(width=0.15) +# +#fig.tight_layout() +#fig.savefig(output_dir / f"continuum_image.{image_format}") +# +#if do_plot: +# plt.show() +#else: +# plt.clf() +# plt.close() # 4th step: 2D-Gauss fit on the continuum map print("#4: 2D-Gauss fit on the continuum map") @@ -203,7 +202,24 @@ "x_stddev", "y_stddev", "theta", - ) + "e_peak", + "e_x_mean", + "e_y_mean", + "e_x_stddev", + "e_y_stddev", + "e_theta", + "e_floor", + "x_mean_arcsec", + "e_x_mean_arcsec", + "y_mean_arcsec", + "e_y_mean_arcsec", + "fwhm_major_arcsec", + "e_fwhm_major_arcsec", + "fwhm_minor_arcsec", + "e_fwhm_minor_arcsec", + "pa_deg", + "e_pa_deg" + ) ) amplitude = float(cont_array.max().values) @@ -212,7 +228,7 @@ x_stddev = params["fitting"]["x_stddev"] y_stddev = params["fitting"]["y_stddev"] theta = params["fitting"]["theta"] -noise = params["fitting"]["noise"] +floor = params["fitting"]["floor"] f = fc.gauss_fit( cont_array, @@ -224,16 +240,29 @@ x_stddev=x_stddev, y_stddev=y_stddev, theta=theta, - noise=noise + floor=floor ) -sigma2hpbw = 2 * np.sqrt(2 * np.log(2)) -hpbw_major_arcsec = float(f.x_stddev * 3600 * sigma2hpbw) -hpbw_major_rad = float(f.x_stddev * sigma2hpbw * np.pi / 180) -hpbw_minor_arcsec = float(f.y_stddev * 3600 * sigma2hpbw) -hpbw_minor_rad = float(f.y_stddev * sigma2hpbw * np.pi / 180) -print(f"hpbw_major: {hpbw_major_arcsec:.1f} [arcsec], {hpbw_major_rad:.1e} [rad]") -print(f"hpbw_minor: {hpbw_minor_arcsec:.1f} [arcsec], {hpbw_minor_rad:.1e} [rad]") +sigma2fwhm = 2 * np.sqrt(2 * np.log(2)) +fwhm_major_arcsec = float(f.x_stddev * 3600 * sigma2fwhm) +fwhm_minor_arcsec = float(f.y_stddev * 3600 * sigma2fwhm) +e_fwhm_major_arcsec = float(f.e_x_stddev * 3600 * sigma2fwhm) +e_fwhm_minor_arcsec = float(f.e_y_stddev * 3600 * sigma2fwhm) +fwhm_major_rad = float(f.x_stddev * sigma2fwhm * np.pi / 180) +fwhm_minor_rad = float(f.y_stddev * sigma2fwhm * np.pi / 180) +x_mean_arcsec = float(f.x_mean * 3600) +e_x_mean_arcsec = float(f.e_x_mean * 3600) +y_mean_arcsec = float(f.y_mean * 3600) +e_y_mean_arcsec = float(f.e_y_mean * 3600) +pa_deg = float(f.theta * 180/np.pi) +e_pa_deg = float(f.e_theta * 180/np.pi) +print(f"Peak: {float(f.peak):.3f}+/-{float(f.e_peak):.3f} [K]") +print(f"X mean {x_mean_arcsec:.1f}+/-{e_x_mean_arcsec:.1f} [arcsec]") +print(f"Y mean: {y_mean_arcsec:.1f}+/-{e_y_mean_arcsec:.1f} [arcsec]") +print(f"fwhm_major: {fwhm_major_arcsec:.1f}+/-{e_fwhm_major_arcsec:.1f} [arcsec]") +print(f"fwhm_minor: {fwhm_minor_arcsec:.1f}+/-{e_fwhm_minor_arcsec:.1f} [arcsec]") +print(f"PA: {pa_deg:.1f}+/-{e_pa_deg:.1f} [deg]") + header = fits.getheader(cont_obs_fits) fits.writeto(cont_mod_fits, f[:, :, 0].values.T, header, overwrite=True) @@ -278,6 +307,23 @@ f.x_stddev, f.y_stddev, f.theta, - ] + f.e_peak, + f.e_x_mean, + f.e_y_mean, + f.e_x_stddev, + f.e_y_stddev, + f.e_theta, + f.e_floor, + x_mean_arcsec, + e_x_mean_arcsec, + y_mean_arcsec, + e_y_mean_arcsec, + fwhm_major_arcsec, + e_fwhm_major_arcsec, + fwhm_minor_arcsec, + e_fwhm_minor_arcsec, + pa_deg, + e_pa_deg + ] ) alldata.write(result_file, format="ascii", overwrite=True) diff --git a/pipelines/utils/functions.py b/pipelines/utils/functions.py index e20b84d..0e90b11 100644 --- a/pipelines/utils/functions.py +++ b/pipelines/utils/functions.py @@ -320,7 +320,7 @@ def gauss_fit( y_stddev=None, theta=None, cov_matrix=None, - noise=0, + floor=0, **kwargs ): if chs is None: @@ -347,7 +347,7 @@ def gauss_fit( theta=theta, cov_matrix=cov_matrix, **kwargs - ) + models.Const2D(noise) + ) + models.Const2D(floor) fit_g = fitting.LevMarLSQFitter() g = fit_g(g_init, mX, mY, subdata) @@ -436,7 +436,7 @@ def gauss_fit( theta=theta, cov_matrix=cov_matrix, **kwargs - ) + models.Const2D(noise) + ) + models.Const2D(floor) fit_g = fitting.LevMarLSQFitter() g = fit_g(g_init, mX, mY, subdata) @@ -472,11 +472,15 @@ def gauss_fit( x_stddevs = np.array([g.x_stddev_0.value]) y_stddevs = np.array([g.y_stddev_0.value]) thetas = np.array([g.theta_0.value]) - noises = np.array([g.amplitude_1.value]) + floors = np.array([g.amplitude_1.value]) error = np.diag(fit_g.fit_info["param_cov"]) ** 0.5 - uncerts0 = np.array(error[0]) # <- added - uncerts3 = np.array(error[3]) # <- added - uncerts4 = np.array(error[4]) # <- added + uncerts0 = np.array(error[0]) + uncerts1 = np.array(error[1]) + uncerts2 = np.array(error[2]) + uncerts3 = np.array(error[3]) + uncerts4 = np.array(error[4]) + uncerts5 = np.array(error[5]) + uncerts6 = np.array(error[6]) result = map_data.copy() result.values = np.transpose(results) @@ -488,12 +492,16 @@ def gauss_fit( "x_stddev": x_stddevs, "y_stddev": y_stddevs, "theta": thetas, - "noise": noises, - "uncert0": uncerts0, # <- added - "uncert3": uncerts3, # <- added - "uncert4": uncerts4, + "floor": floors, + "e_peak": uncerts0, + "e_x_mean": uncerts1, + "e_y_mean": uncerts2, + "e_x_stddev": uncerts3, + "e_y_stddev": uncerts4, + "e_theta": uncerts5, + "e_floor": uncerts6, } - ) # <- added + ) return result diff --git a/yaml/beammap_params.yaml b/yaml/beammap_params.yaml index 60c1848..bda26cf 100644 --- a/yaml/beammap_params.yaml +++ b/yaml/beammap_params.yaml @@ -1,5 +1,5 @@ file: - output_dir: ~/Documents/deshima-dev/qlook-pipeline/results # output directory (absolute path) + output_dir: /home/deshima/raid1/desql/git/qlook-pipeline/test/result # output directory (absolute path) result_file: results.txt # file name of result image_format: png # image format (pdf, png, svg, etc...) do_plot: False # plot flag (True or False) @@ -37,6 +37,3 @@ fitting: # initial values for the continuum fitting y_stddev: 0.003 # standard deviation of y theta: 1.25 # position angle -planet: - istart: 1000 # line index in the beginning of fluxtxt - iend: 1100 # line index in the end of fluxtxt \ No newline at end of file diff --git a/yaml/pointing_params.yaml b/yaml/pointing_params.yaml index a12fdf9..37177f7 100644 --- a/yaml/pointing_params.yaml +++ b/yaml/pointing_params.yaml @@ -1,8 +1,8 @@ file: - output_dir: ~/Documents/deshima-dev/qlook-pipeline/results # output directory + output_dir: /home/deshima/raid1/desql/git/qlook-pipeline/test/result #otput directory result_file: results.txt # file name of result image_format: png # image format (pdf, png, svg, etc...) - do_plot: True # plot flag (True or False) + do_plot: False # plot flag (True or False) dpi: 100 # dpi loaddfits: # parameters for dc.io.loaddfits() @@ -30,5 +30,5 @@ imaging: fitting: # initial values for the continuum fitting x_stddev: 0.006 # standard deviation of x y_stddev: 0.006 # standard deviation of y - theta: 0.60357481 # position angle - noise: -0.02 \ No newline at end of file + theta: 0.00 # position angle + floor: -0.02