From ced27d0146745619ba01147b23171d2f1b125239 Mon Sep 17 00:00:00 2001 From: Fujita Shinji <109708402+ShinjiFujita@users.noreply.github.com> Date: Wed, 31 Jul 2024 19:13:51 +0900 Subject: [PATCH 1/3] Update qlook.py (issue186) --- decode/qlook.py | 177 ++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 162 insertions(+), 15 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 83c2cb9..e97d776 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -27,6 +27,9 @@ from fire import Fire from matplotlib.figure import Figure from . import assign, convert, load, make, plot, select, utils +import copy +from scipy.optimize import curve_fit +import pandas as pd # constants @@ -40,6 +43,7 @@ DEFAULT_INCL_MKID_IDS = None DEFAULT_MIN_FREQUENCY = None DEFAULT_MAX_FREQUENCY = None +DEFAULT_ROLLING_TIME = 200 DEFAULT_OUTDIR = Path() DEFAULT_OVERWRITE = False DEFAULT_SKYCOORD_GRID = "6 arcsec" @@ -119,6 +123,7 @@ def daisy( exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, min_frequency: Optional[str] = DEFAULT_MIN_FREQUENCY, max_frequency: Optional[str] = DEFAULT_MAX_FREQUENCY, + rolling_time: int = DEFAULT_ROLLING_TIME, data_type: Literal["auto", "brightness", "df/f"] = DEFAULT_DATA_TYPE, # options for analysis source_radius: str = "60 arcsec", @@ -185,7 +190,6 @@ def daisy( data_type=data_type, skycoord_units=skycoord_units, ) - da = select.by(da, "state", exclude="GRAD") # fmt: off is_source = ( @@ -208,8 +212,13 @@ def daisy( kwargs={"fill_value": "extrapolate"}, ) ) - t_atm = da_on.temperature - da_sub = t_atm * (da_on - da_base) / (t_atm - da_base) + da_sub = da_on - da_base.values + + ### Rolling + da_sub_rolled = ( + copy.deepcopy(da_sub).rolling(time=int(rolling_time), center=True).mean() + ) + da_sub = da_sub - da_sub_rolled # make continuum series weight = get_chan_weight(da_off, method=chan_weight, pwv=pwv) @@ -223,6 +232,19 @@ def daisy( ) cont = cube.weighted(weight.fillna(0)).mean("chan") + ### GaussFit (cont) + data = np.array(copy.deepcopy(cont).data) + data[data != data] = 0.0 + x, y = np.meshgrid(np.array(cube["lon"]), np.array(cube["lat"])) + initial_guess = (1, 0, 0, 30, 30, 0, 0) + popt, pcov = curve_fit(gaussian_2d, (x, y), data.ravel(), p0=initial_guess) + perr = np.sqrt(np.diag(pcov)) + data_fitted = gaussian_2d((x, y), *popt).reshape(x.shape) + + ### GaussFit (all chan) + df_gauss_fit = fit_cube(cube) + # to toml here + # save result suffixes = f".{suffix}.{format}" file = Path(outdir) / Path(dems).with_suffix(suffixes).name @@ -241,13 +263,40 @@ def daisy( max_pix = cont.where(cont == cont.max(), drop=True) cont.plot(ax=ax) # type: ignore + cont_fit = ax.contour( + data_fitted, + extent=(x.min(), x.max(), y.min(), y.max()), + origin="lower", + levels=np.linspace(0, np.nanmax(data), 8), + colors="k", + ) ax.set_xlim(-map_lim, map_lim) ax.set_ylim(-map_lim, map_lim) - ax.set_title( - f"Maximum {cont.long_name.lower()} = {cont.max():.2e} [{cont.units}]\n" - f"(dAz = {float(max_pix.lon):+.1f} [{cont.lon.attrs['units']}], " - f"dEl = {float(max_pix.lat):+.1f} [{cont.lat.attrs['units']}])" - ) + # ax.set_title( + # f"Maximum {cont.long_name.lower()} = {cont.max():.2e} [{cont.units}]\n" + # f"(dAz = {float(max_pix.lon):+.1f} [{cont.lon.attrs['units']}], " + # f"dEl = {float(max_pix.lat):+.1f} [{cont.lat.attrs['units']}])" + # ) + if min_frequency == None or max_frequency == None: + ax.set_title( + f"Maximum {cont.long_name.lower()} = {popt[0]:+.2f} [{cont.units}]\n" + f"(dAz = {popt[1]:+.2f} [{cont.lon.attrs['units']}], " + f"dEl = {popt[2]:+.2f} [{cont.lat.attrs['units']}]), \n" + f"(sigma_x = {popt[3]:+.2f}, " + f"sigma_y = {popt[4]:+.2f}," + f"theta = {np.rad2deg(popt[5]):+.1f}, \n" + f"min_frequency: None, max_frequency: None" + ) + else: + ax.set_title( + f"Maximum {cont.long_name.lower()} = {popt[0]:+.2f} [{cont.units}]\n" + f"(dAz = {popt[1]:+.2f} [{cont.lon.attrs['units']}], " + f"dEl = {popt[2]:+.2f} [{cont.lat.attrs['units']}]), \n" + f"(sigma_x = {popt[3]:+.2f}, " + f"sigma_y = {popt[4]:+.2f}," + f"theta = {np.rad2deg(popt[5]):+.1f}, \n" + f"min_frequency: {float(min_frequency):+.1f} GHz, max_frequency: {float(max_frequency):+.1f} GHz" + ) for ax in axes: # type: ignore ax.grid(True) @@ -431,8 +480,7 @@ def raster( kwargs={"fill_value": "extrapolate"}, ) ) - t_atm = da_on.temperature - da_sub = t_atm * (da_on - da_base) / (t_atm - da_base) + da_sub = da_on - da_base.values # make continuum series weight = get_chan_weight(da_off, method=chan_weight, pwv=pwv) @@ -446,6 +494,19 @@ def raster( ) cont = cube.weighted(weight.fillna(0)).mean("chan") + ### GaussFit (cont) + data = np.array(copy.deepcopy(cont).data) + data[data != data] = 0.0 + x, y = np.meshgrid(np.array(cube["lon"]), np.array(cube["lat"])) + initial_guess = (1, 0, 0, 30, 30, 0, 0) + popt, pcov = curve_fit(gaussian_2d, (x, y), data.ravel(), p0=initial_guess) + perr = np.sqrt(np.diag(pcov)) + data_fitted = gaussian_2d((x, y), *popt).reshape(x.shape) + + ### GaussFit (all chan) + df_gauss_fit = fit_cube(cube) + # to toml here + # save result suffixes = f".{suffix}.{format}" file = Path(outdir) / Path(dems).with_suffix(suffixes).name @@ -464,13 +525,35 @@ def raster( max_pix = cont.where(cont == cont.max(), drop=True) cont.plot(ax=ax) # type: ignore + cont_fit = ax.contour( + data_fitted, + extent=(x.min(), x.max(), y.min(), y.max()), + origin="lower", + levels=np.linspace(0, np.nanmax(data), 8), + colors="k", + ) ax.set_xlim(-map_lim, map_lim) ax.set_ylim(-map_lim, map_lim) - ax.set_title( - f"Maximum {cont.long_name.lower()} = {cont.max():.2e} [{cont.units}]\n" - f"(dAz = {float(max_pix.lon):+.1f} [{cont.lon.attrs['units']}], " - f"dEl = {float(max_pix.lat):+.1f} [{cont.lat.attrs['units']}])" - ) + if min_frequency == None or max_frequency == None: + ax.set_title( + f"Maximum {cont.long_name.lower()} = {popt[0]:+.2f} [{cont.units}]\n" + f"(dAz = {popt[1]:+.2f} [{cont.lon.attrs['units']}], " + f"dEl = {popt[2]:+.2f} [{cont.lat.attrs['units']}]), \n" + f"(sigma_x = {popt[3]:+.2f}, " + f"sigma_y = {popt[4]:+.2f}," + f"theta = {np.rad2deg(popt[5]):+.1f}, \n" + f"min_frequency: None, max_frequency: None" + ) + else: + ax.set_title( + f"Maximum {cont.long_name.lower()} = {popt[0]:+.2f} [{cont.units}]\n" + f"(dAz = {popt[1]:+.2f} [{cont.lon.attrs['units']}], " + f"dEl = {popt[2]:+.2f} [{cont.lat.attrs['units']}]), \n" + f"(sigma_x = {popt[3]:+.2f}, " + f"sigma_y = {popt[4]:+.2f}," + f"theta = {np.rad2deg(popt[5]):+.1f}, \n" + f"min_frequency: {float(min_frequency):+.1f} GHz, max_frequency: {float(max_frequency):+.1f} GHz" + ) for ax in axes: # type: ignore ax.grid(True) @@ -1233,6 +1316,70 @@ def save_qlook( raise ValueError("Extension of filename is not valid.") +def gaussian_2d(xy, amp, xo, yo, sigma_x, sigma_y, theta, offset): + x, y = xy + xo = float(xo) + yo = float(yo) + a = (np.cos(theta) ** 2) / (2 * sigma_x**2) + (np.sin(theta) ** 2) / ( + 2 * sigma_y**2 + ) + b = -(np.sin(2 * theta)) / (4 * sigma_x**2) + (np.sin(2 * theta)) / (4 * sigma_y**2) + c = (np.sin(theta) ** 2) / (2 * sigma_x**2) + (np.cos(theta) ** 2) / ( + 2 * sigma_y**2 + ) + g = offset + amp * np.exp( + -(a * ((x - xo) ** 2) + 2 * b * (x - xo) * (y - yo) + c * ((y - yo) ** 2)) + ) + return g.ravel() + + +def fit_cube(cube): + res_list = [] + header = [ + "mkid_id", + "amp", + "center_lon", + "center_lat", + "sigma_x", + "sigma_y", + "theta_deg", + "floor", + "err_amp", + "err_center_lon", + "err_center_lat", + "err_sigma_x", + "err_sigma_y", + "err_theta_deg", + "err_floor", + ] + for i in range(len(cube)): + res = [] + xr_tempo = cube[i] + mkid_id = int(xr_tempo["d2_mkid_id"]) + res.append(mkid_id) + try: + data_tempo = np.array(xr_tempo.data) + data_tempo[data_tempo != data_tempo] = 0.0 + x, y = np.meshgrid(np.array(cube["lon"]), np.array(cube["lat"])) + initial_guess = (1, 0, 0, 30, 30, 0, 0) + popt, pcov = curve_fit( + gaussian_2d, (x, y), data_tempo.ravel(), p0=initial_guess + ) + perr = np.sqrt(np.diag(pcov)) + for j in range(7): + res.append(popt[j]) + for j in range(7): + res.append(perr[j]) + except: + for j in range(7): + res.append(np.nan) + for j in range(7): + res.append(np.nan) + res_list.append(res) + res_df = pd.DataFrame(np.array(res_list), columns=np.array(header)) + return res_df + + def main() -> None: """Entry point of the decode-qlook command.""" From c63cf06c8d7ceccfd6696975745b1412ce650e49 Mon Sep 17 00:00:00 2001 From: Fujita Shinji <109708402+ShinjiFujita@users.noreply.github.com> Date: Fri, 2 Aug 2024 22:20:44 +0900 Subject: [PATCH 2/3] Update qlook.py --- decode/qlook.py | 173 +++++------------------------------------------- 1 file changed, 16 insertions(+), 157 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index e97d776..ddfc5ae 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -27,9 +27,6 @@ from fire import Fire from matplotlib.figure import Figure from . import assign, convert, load, make, plot, select, utils -import copy -from scipy.optimize import curve_fit -import pandas as pd # constants @@ -190,6 +187,7 @@ def daisy( data_type=data_type, skycoord_units=skycoord_units, ) + da = select.by(da, "state", exclude="GRAD") # fmt: off is_source = ( @@ -212,12 +210,11 @@ def daisy( kwargs={"fill_value": "extrapolate"}, ) ) - da_sub = da_on - da_base.values + t_atm = da_on.temperature + da_sub = t_atm * (da_on - da_base) / (t_atm - da_base) ### Rolling - da_sub_rolled = ( - copy.deepcopy(da_sub).rolling(time=int(rolling_time), center=True).mean() - ) + da_sub_rolled = da_sub.rolling(time=int(rolling_time), center=True).mean() da_sub = da_sub - da_sub_rolled # make continuum series @@ -232,19 +229,6 @@ def daisy( ) cont = cube.weighted(weight.fillna(0)).mean("chan") - ### GaussFit (cont) - data = np.array(copy.deepcopy(cont).data) - data[data != data] = 0.0 - x, y = np.meshgrid(np.array(cube["lon"]), np.array(cube["lat"])) - initial_guess = (1, 0, 0, 30, 30, 0, 0) - popt, pcov = curve_fit(gaussian_2d, (x, y), data.ravel(), p0=initial_guess) - perr = np.sqrt(np.diag(pcov)) - data_fitted = gaussian_2d((x, y), *popt).reshape(x.shape) - - ### GaussFit (all chan) - df_gauss_fit = fit_cube(cube) - # to toml here - # save result suffixes = f".{suffix}.{format}" file = Path(outdir) / Path(dems).with_suffix(suffixes).name @@ -263,40 +247,13 @@ def daisy( max_pix = cont.where(cont == cont.max(), drop=True) cont.plot(ax=ax) # type: ignore - cont_fit = ax.contour( - data_fitted, - extent=(x.min(), x.max(), y.min(), y.max()), - origin="lower", - levels=np.linspace(0, np.nanmax(data), 8), - colors="k", - ) ax.set_xlim(-map_lim, map_lim) ax.set_ylim(-map_lim, map_lim) - # ax.set_title( - # f"Maximum {cont.long_name.lower()} = {cont.max():.2e} [{cont.units}]\n" - # f"(dAz = {float(max_pix.lon):+.1f} [{cont.lon.attrs['units']}], " - # f"dEl = {float(max_pix.lat):+.1f} [{cont.lat.attrs['units']}])" - # ) - if min_frequency == None or max_frequency == None: - ax.set_title( - f"Maximum {cont.long_name.lower()} = {popt[0]:+.2f} [{cont.units}]\n" - f"(dAz = {popt[1]:+.2f} [{cont.lon.attrs['units']}], " - f"dEl = {popt[2]:+.2f} [{cont.lat.attrs['units']}]), \n" - f"(sigma_x = {popt[3]:+.2f}, " - f"sigma_y = {popt[4]:+.2f}," - f"theta = {np.rad2deg(popt[5]):+.1f}, \n" - f"min_frequency: None, max_frequency: None" - ) - else: - ax.set_title( - f"Maximum {cont.long_name.lower()} = {popt[0]:+.2f} [{cont.units}]\n" - f"(dAz = {popt[1]:+.2f} [{cont.lon.attrs['units']}], " - f"dEl = {popt[2]:+.2f} [{cont.lat.attrs['units']}]), \n" - f"(sigma_x = {popt[3]:+.2f}, " - f"sigma_y = {popt[4]:+.2f}," - f"theta = {np.rad2deg(popt[5]):+.1f}, \n" - f"min_frequency: {float(min_frequency):+.1f} GHz, max_frequency: {float(max_frequency):+.1f} GHz" - ) + ax.set_title( + f"Maximum {cont.long_name.lower()} = {cont.max():.2e} [{cont.units}]\n" + f"(dAz = {float(max_pix.lon):+.1f} [{cont.lon.attrs['units']}], " + f"dEl = {float(max_pix.lat):+.1f} [{cont.lat.attrs['units']}])" + ) for ax in axes: # type: ignore ax.grid(True) @@ -480,7 +437,8 @@ def raster( kwargs={"fill_value": "extrapolate"}, ) ) - da_sub = da_on - da_base.values + t_atm = da_on.temperature + da_sub = t_atm * (da_on - da_base) / (t_atm - da_base) # make continuum series weight = get_chan_weight(da_off, method=chan_weight, pwv=pwv) @@ -494,19 +452,6 @@ def raster( ) cont = cube.weighted(weight.fillna(0)).mean("chan") - ### GaussFit (cont) - data = np.array(copy.deepcopy(cont).data) - data[data != data] = 0.0 - x, y = np.meshgrid(np.array(cube["lon"]), np.array(cube["lat"])) - initial_guess = (1, 0, 0, 30, 30, 0, 0) - popt, pcov = curve_fit(gaussian_2d, (x, y), data.ravel(), p0=initial_guess) - perr = np.sqrt(np.diag(pcov)) - data_fitted = gaussian_2d((x, y), *popt).reshape(x.shape) - - ### GaussFit (all chan) - df_gauss_fit = fit_cube(cube) - # to toml here - # save result suffixes = f".{suffix}.{format}" file = Path(outdir) / Path(dems).with_suffix(suffixes).name @@ -525,35 +470,13 @@ def raster( max_pix = cont.where(cont == cont.max(), drop=True) cont.plot(ax=ax) # type: ignore - cont_fit = ax.contour( - data_fitted, - extent=(x.min(), x.max(), y.min(), y.max()), - origin="lower", - levels=np.linspace(0, np.nanmax(data), 8), - colors="k", - ) ax.set_xlim(-map_lim, map_lim) ax.set_ylim(-map_lim, map_lim) - if min_frequency == None or max_frequency == None: - ax.set_title( - f"Maximum {cont.long_name.lower()} = {popt[0]:+.2f} [{cont.units}]\n" - f"(dAz = {popt[1]:+.2f} [{cont.lon.attrs['units']}], " - f"dEl = {popt[2]:+.2f} [{cont.lat.attrs['units']}]), \n" - f"(sigma_x = {popt[3]:+.2f}, " - f"sigma_y = {popt[4]:+.2f}," - f"theta = {np.rad2deg(popt[5]):+.1f}, \n" - f"min_frequency: None, max_frequency: None" - ) - else: - ax.set_title( - f"Maximum {cont.long_name.lower()} = {popt[0]:+.2f} [{cont.units}]\n" - f"(dAz = {popt[1]:+.2f} [{cont.lon.attrs['units']}], " - f"dEl = {popt[2]:+.2f} [{cont.lat.attrs['units']}]), \n" - f"(sigma_x = {popt[3]:+.2f}, " - f"sigma_y = {popt[4]:+.2f}," - f"theta = {np.rad2deg(popt[5]):+.1f}, \n" - f"min_frequency: {float(min_frequency):+.1f} GHz, max_frequency: {float(max_frequency):+.1f} GHz" - ) + ax.set_title( + f"Maximum {cont.long_name.lower()} = {cont.max():.2e} [{cont.units}]\n" + f"(dAz = {float(max_pix.lon):+.1f} [{cont.lon.attrs['units']}], " + f"dEl = {float(max_pix.lat):+.1f} [{cont.lat.attrs['units']}])" + ) for ax in axes: # type: ignore ax.grid(True) @@ -1316,70 +1239,6 @@ def save_qlook( raise ValueError("Extension of filename is not valid.") -def gaussian_2d(xy, amp, xo, yo, sigma_x, sigma_y, theta, offset): - x, y = xy - xo = float(xo) - yo = float(yo) - a = (np.cos(theta) ** 2) / (2 * sigma_x**2) + (np.sin(theta) ** 2) / ( - 2 * sigma_y**2 - ) - b = -(np.sin(2 * theta)) / (4 * sigma_x**2) + (np.sin(2 * theta)) / (4 * sigma_y**2) - c = (np.sin(theta) ** 2) / (2 * sigma_x**2) + (np.cos(theta) ** 2) / ( - 2 * sigma_y**2 - ) - g = offset + amp * np.exp( - -(a * ((x - xo) ** 2) + 2 * b * (x - xo) * (y - yo) + c * ((y - yo) ** 2)) - ) - return g.ravel() - - -def fit_cube(cube): - res_list = [] - header = [ - "mkid_id", - "amp", - "center_lon", - "center_lat", - "sigma_x", - "sigma_y", - "theta_deg", - "floor", - "err_amp", - "err_center_lon", - "err_center_lat", - "err_sigma_x", - "err_sigma_y", - "err_theta_deg", - "err_floor", - ] - for i in range(len(cube)): - res = [] - xr_tempo = cube[i] - mkid_id = int(xr_tempo["d2_mkid_id"]) - res.append(mkid_id) - try: - data_tempo = np.array(xr_tempo.data) - data_tempo[data_tempo != data_tempo] = 0.0 - x, y = np.meshgrid(np.array(cube["lon"]), np.array(cube["lat"])) - initial_guess = (1, 0, 0, 30, 30, 0, 0) - popt, pcov = curve_fit( - gaussian_2d, (x, y), data_tempo.ravel(), p0=initial_guess - ) - perr = np.sqrt(np.diag(pcov)) - for j in range(7): - res.append(popt[j]) - for j in range(7): - res.append(perr[j]) - except: - for j in range(7): - res.append(np.nan) - for j in range(7): - res.append(np.nan) - res_list.append(res) - res_df = pd.DataFrame(np.array(res_list), columns=np.array(header)) - return res_df - - def main() -> None: """Entry point of the decode-qlook command.""" From f24d25eaa1e7aa5e15936691f880bb6c29c56420 Mon Sep 17 00:00:00 2001 From: Fujita Shinji <109708402+ShinjiFujita@users.noreply.github.com> Date: Tue, 6 Aug 2024 17:32:29 +0900 Subject: [PATCH 3/3] Update qlook.py --- decode/qlook.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index ddfc5ae..2abe74f 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -120,9 +120,9 @@ def daisy( exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, min_frequency: Optional[str] = DEFAULT_MIN_FREQUENCY, max_frequency: Optional[str] = DEFAULT_MAX_FREQUENCY, - rolling_time: int = DEFAULT_ROLLING_TIME, data_type: Literal["auto", "brightness", "df/f"] = DEFAULT_DATA_TYPE, # options for analysis + rolling_time: int = DEFAULT_ROLLING_TIME, source_radius: str = "60 arcsec", chan_weight: Literal["uniform", "std", "std/tx"] = "std/tx", pwv: Literal["0.5", "1.0", "2.0", "3.0", "4.0", "5.0"] = "5.0", @@ -151,6 +151,7 @@ def daisy( Defaults to no maximum frequency bound. data_type: Data type of the input DEMS file. Defaults to the ``long_name`` attribute in it. + rolling_time: Moving window size. source_radius: Radius of the on-source area. Other areas are considered off-source in sky subtraction. chan_weight: Weighting method along the channel axis. @@ -189,6 +190,10 @@ def daisy( ) da = select.by(da, "state", exclude="GRAD") + ### Rolling + da_rolled = da.rolling(time=int(rolling_time), center=True).mean() + da = da - da_rolled + # fmt: off is_source = ( (da.lon**2 + da.lat**2) @@ -213,10 +218,6 @@ def daisy( t_atm = da_on.temperature da_sub = t_atm * (da_on - da_base) / (t_atm - da_base) - ### Rolling - da_sub_rolled = da_sub.rolling(time=int(rolling_time), center=True).mean() - da_sub = da_sub - da_sub_rolled - # make continuum series weight = get_chan_weight(da_off, method=chan_weight, pwv=pwv) series = da_sub.weighted(weight.fillna(0)).mean("chan")