From e28b288b4a8c489d4d53574fde2eb9a800f0d621 Mon Sep 17 00:00:00 2001 From: Fujita Shinji <109708402+ShinjiFujita@users.noreply.github.com> Date: Fri, 2 Aug 2024 22:20:06 +0900 Subject: [PATCH 1/6] Update qlook.py --- decode/qlook.py | 76 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 75 insertions(+), 1 deletion(-) diff --git a/decode/qlook.py b/decode/qlook.py index 83c2cb9..fe49067 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -27,7 +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 DATA_FORMATS = "csv", "nc", "zarr", "zarr.zip" @@ -223,6 +225,10 @@ def daisy( ) cont = cube.weighted(weight.fillna(0)).mean("chan") + ### 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 @@ -446,6 +452,10 @@ def raster( ) cont = cube.weighted(weight.fillna(0)).mean("chan") + ### 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 @@ -1233,6 +1243,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 953748fa4f27ef896162412d73fce00fefc7c459 Mon Sep 17 00:00:00 2001 From: Fujita Shinji <109708402+ShinjiFujita@users.noreply.github.com> Date: Tue, 6 Aug 2024 17:31:21 +0900 Subject: [PATCH 2/6] Update fit.py --- decode/fit.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/decode/fit.py b/decode/fit.py index 01986bb..6eda07c 100644 --- a/decode/fit.py +++ b/decode/fit.py @@ -93,3 +93,32 @@ def dtau_dpwv(freq: NDArray[np.float_]) -> xr.DataArray: tau = load.atm(type="tau").interp(freq=freq, method="linear") fit = tau.curvefit("pwv", lambda x, a, b: a * x + b) return fit["curvefit_coefficients"].sel(param="a", drop=True) + + +def cube( + cube: xr.DataArray, + /, + *, + init_amp: float = 1.0, + init_x0: float = 0.0, + init_y0: float = 0.0, + init_sigma_x: float = 20.0, + init_sigma_y: float = 20.0, + init_theta: float = 0.0, + init_offset: float = 0.0, +) -> xr.DataArray: + """Apply 2D Gaussian fit to each channel of a 3D spectral cube.""" + return cube.curvefit( + coords=("lon", "lat"), + func=gaussian_2d, + p0={ + "amp": init_amp, + "x0": init_x0, + "y0": init_y0, + "sigma_x": init_sigma_x, + "sigma_y": init_sigma_y, + "theta": init_theta, + "offset": init_offset, + }, + errors="ignore", + ) From c7f700744d3245f0714e800734e953d6e00dc4c8 Mon Sep 17 00:00:00 2001 From: Fujita Shinji <109708402+ShinjiFujita@users.noreply.github.com> Date: Tue, 6 Aug 2024 17:31:49 +0900 Subject: [PATCH 3/6] Update qlook.py --- decode/qlook.py | 62 ++++++------------------------------------------- 1 file changed, 7 insertions(+), 55 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index fe49067..9d4ad7d 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -17,7 +17,7 @@ from pathlib import Path from typing import Any, Literal, Optional, Sequence, Union, cast from warnings import catch_warnings, simplefilter - +import copy # dependencies import numpy as np @@ -26,9 +26,8 @@ from astropy.units import Quantity 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 +from . import assign, convert, load, make, plot, select, utils import pandas as pd # constants @@ -226,7 +225,7 @@ def daisy( cont = cube.weighted(weight.fillna(0)).mean("chan") ### GaussFit (all chan) - df_gauss_fit = fit_cube(cube) + fitted_cube = fit.cube(cube) # to toml here # save result @@ -1243,10 +1242,10 @@ def save_qlook( raise ValueError("Extension of filename is not valid.") -def gaussian_2d(xy, amp, xo, yo, sigma_x, sigma_y, theta, offset): +def gaussian_2d(xy, amp, x0, y0, sigma_x, sigma_y, theta, offset): x, y = xy - xo = float(xo) - yo = float(yo) + x0 = float(x0) + y0 = float(y0) a = (np.cos(theta) ** 2) / (2 * sigma_x**2) + (np.sin(theta) ** 2) / ( 2 * sigma_y**2 ) @@ -1255,58 +1254,11 @@ def gaussian_2d(xy, amp, xo, yo, sigma_x, sigma_y, theta, offset): 2 * sigma_y**2 ) g = offset + amp * np.exp( - -(a * ((x - xo) ** 2) + 2 * b * (x - xo) * (y - yo) + c * ((y - yo) ** 2)) + -(a * ((x - x0) ** 2) + 2 * b * (x - x0) * (y - y0) + c * ((y - y0) ** 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 be5472b6a9ec562f3eb9552cd5bbcbde4d81643d Mon Sep 17 00:00:00 2001 From: Fujita Shinji <109708402+ShinjiFujita@users.noreply.github.com> Date: Tue, 6 Aug 2024 17:43:13 +0900 Subject: [PATCH 4/6] Update qlook.py --- decode/qlook.py | 21 ++------------------- 1 file changed, 2 insertions(+), 19 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 9d4ad7d..6fd32ed 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -27,7 +27,7 @@ from fire import Fire from matplotlib.figure import Figure from scipy.optimize import curve_fit -from . import assign, convert, load, make, plot, select, utils +from . import assign, convert, load, make, plot, select, utils, fit import pandas as pd # constants @@ -452,7 +452,7 @@ def raster( cont = cube.weighted(weight.fillna(0)).mean("chan") ### GaussFit (all chan) - df_gauss_fit = fit_cube(cube) + fitted_cube = fit.cube(cube) # to toml here # save result @@ -1242,23 +1242,6 @@ def save_qlook( raise ValueError("Extension of filename is not valid.") -def gaussian_2d(xy, amp, x0, y0, sigma_x, sigma_y, theta, offset): - x, y = xy - x0 = float(x0) - y0 = float(y0) - 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 - x0) ** 2) + 2 * b * (x - x0) * (y - y0) + c * ((y - y0) ** 2)) - ) - return g.ravel() - - def main() -> None: """Entry point of the decode-qlook command.""" From e33041eb7ad1786d1eaa182c86e9fa0013dd0dd4 Mon Sep 17 00:00:00 2001 From: Fujita Shinji <109708402+ShinjiFujita@users.noreply.github.com> Date: Tue, 6 Aug 2024 17:43:36 +0900 Subject: [PATCH 5/6] Update fit.py --- decode/fit.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/decode/fit.py b/decode/fit.py index 6eda07c..b66ed67 100644 --- a/decode/fit.py +++ b/decode/fit.py @@ -122,3 +122,20 @@ def cube( }, errors="ignore", ) + + +def gaussian_2d(xy, amp, x0, y0, sigma_x, sigma_y, theta, offset): + x, y = xy + x0 = float(x0) + y0 = float(y0) + 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 - x0) ** 2) + 2 * b * (x - x0) * (y - y0) + c * ((y - y0) ** 2)) + ) + return g.ravel() From a52f664e40fc890d509f281eb14097dabe29681e Mon Sep 17 00:00:00 2001 From: Fujita Shinji <109708402+ShinjiFujita@users.noreply.github.com> Date: Tue, 6 Aug 2024 17:46:53 +0900 Subject: [PATCH 6/6] Update fit.py --- decode/fit.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/decode/fit.py b/decode/fit.py index b66ed67..db2b90a 100644 --- a/decode/fit.py +++ b/decode/fit.py @@ -106,7 +106,7 @@ def cube( init_sigma_y: float = 20.0, init_theta: float = 0.0, init_offset: float = 0.0, -) -> xr.DataArray: +) -> xr.Dataset: """Apply 2D Gaussian fit to each channel of a 3D spectral cube.""" return cube.curvefit( coords=("lon", "lat"),