Skip to content

Commit

Permalink
Update qlook.py (issue187)
Browse files Browse the repository at this point in the history
  • Loading branch information
@ShinjiFujita
ShinjiFujita authored Jul 29, 2024
1 parent da3a29d commit ce46b92
Showing 1 changed file with 153 additions and 13 deletions.
166 changes: 153 additions & 13 deletions decode/qlook.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -185,7 +188,6 @@ def daisy(
data_type=data_type,
skycoord_units=skycoord_units,
)
da = select.by(da, "state", exclude="GRAD")

# fmt: off
is_source = (
Expand Down Expand Up @@ -222,6 +224,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
Expand All @@ -240,13 +255,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: {min_frequency/1e9:+.1f} GHz, max_frequency: {max_frequency/1e9:+.1f} GHz"
)

for ax in axes: # type: ignore
ax.grid(True)
Expand Down Expand Up @@ -430,8 +472,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)
Expand All @@ -445,6 +486,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
Expand All @@ -463,13 +517,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: {min_frequency/1e9:+.1f} GHz, max_frequency: {max_frequency/1e9:+.1f} GHz"
)

for ax in axes: # type: ignore
ax.grid(True)
Expand Down Expand Up @@ -1232,6 +1308,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."""

Expand Down

0 comments on commit ce46b92

Please sign in to comment.