update pointing (yaml, function, .py) and make test set

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()

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)

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
 

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

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
+  theta: 0.00          # position angle
+  floor: -0.02