From 0690d953970f96f84d0bc691ca4ba1fad2296d9e Mon Sep 17 00:00:00 2001 From: Miles Lucas Date: Wed, 29 May 2024 20:40:48 -1000 Subject: [PATCH] add m3 diat to polarimetry model footprint --- .gitignore | 4 +++- pyproject.toml | 3 ++- src/vampires_dpp/__init__.py | 2 +- src/vampires_dpp/image_processing.py | 20 +++++++++++++++----- src/vampires_dpp/indexing.py | 2 +- src/vampires_dpp/pdi/models.py | 15 +++++++++------ src/vampires_dpp/pdi/processing.py | 5 +++-- src/vampires_dpp/pipeline/pipeline.py | 8 ++++++-- src/vampires_dpp/wcs.py | 12 +++++++++--- 9 files changed, 49 insertions(+), 22 deletions(-) diff --git a/.gitignore b/.gitignore index b4c139f..d44c96e 100644 --- a/.gitignore +++ b/.gitignore @@ -166,4 +166,6 @@ docs/**/data docs/.jupyter_cache docs/jupyter_execute -.DS_store \ No newline at end of file +.DS_store + +tests/data \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 85a0029..d7dec4e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,10 +25,11 @@ dependencies = [ "photutils", "pydantic", "reproject", - "synphot", + "result", "scikit_image", "scipy", "sep", + "synphot", "tomli_w", "tomli", "tqdm", diff --git a/src/vampires_dpp/__init__.py b/src/vampires_dpp/__init__.py index 8f52a6e..b665bfc 100644 --- a/src/vampires_dpp/__init__.py +++ b/src/vampires_dpp/__init__.py @@ -1 +1 @@ -__version__ = "0.12.9" +__version__ = "0.12.10" diff --git a/src/vampires_dpp/image_processing.py b/src/vampires_dpp/image_processing.py index 251f468..7988830 100644 --- a/src/vampires_dpp/image_processing.py +++ b/src/vampires_dpp/image_processing.py @@ -164,7 +164,7 @@ def weighted_collapse(data: ArrayLike, angles: ArrayLike, **kwargs) -> NDArray: # if the variance is zero, return the mean if np.allclose(variance_frame, 0): derotated = derotate_cube(data, angles, **kwargs) - return np.nanmean(derotated, 0) + return bn.nanmean(derotated, 0) # expand the variance frame into a cube variance_cube = np.repeat(variance_frame, data.shape[0], axis=0) @@ -173,8 +173,8 @@ def weighted_collapse(data: ArrayLike, angles: ArrayLike, **kwargs) -> NDArray: derotated_variance = derotate_cube(variance_cube, angles, **kwargs) derotated_variance[derotated_variance == 0] = np.inf # calculate weighted sum - numer = np.nansum(derotated_data / derotated_variance, axis=0) - denom = np.nansum(1 / derotated_variance, axis=0) + numer = bn.nansum(derotated_data / derotated_variance, axis=0) + denom = bn.nansum(1 / derotated_variance, axis=0) weighted_frame = numer / denom return weighted_frame @@ -217,10 +217,14 @@ def collapse_cube( weights = 1 / bn.nanvar(cube, axis=(1, 2), keepdims=True) frame = bn.nansum(cube * weights, axis=0) / bn.nansum(weights) case "biweight": - frame = biweight_location(cube, axis=0, c=6, skip_nans=True) + # suppress all-nan axis warnings + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + frame = biweight_location(cube, axis=0, c=6, ignore_nan=True) if header is not None: - header["COL_METH"] = method, "DPP cube collapse method" + header["hierarch DPP COL METH"] = method, "DPP cube collapse method" + header["hierarch DPP NCOADD"] = cube.shape[0], "Num. of coadded frames in cube" return frame, header @@ -505,3 +509,9 @@ def adaptive_sigma_clip_mask(data, sigma=10, boxsize=8): output_mask[inds] = np.abs(cutout - med) > sigma * std return output_mask + + +def create_footprint(cube, angles): + mask = np.isfinite(cube) + derot = derotate_cube(mask.astype(float), angles) + return bn.nanmean(derot, axis=0) diff --git a/src/vampires_dpp/indexing.py b/src/vampires_dpp/indexing.py index c74896b..c98fd78 100644 --- a/src/vampires_dpp/indexing.py +++ b/src/vampires_dpp/indexing.py @@ -80,7 +80,7 @@ def frame_angles_image(frame, center): def frame_angles_astro(frame, center): Ys, Xs = np.ogrid[0 : frame.shape[-2], 0 : frame.shape[-1]] # degrees East of North: phi = arctan(-x, y) - thetas = np.arctan2(center[-2] - Xs, Ys - center[-1]) + thetas = np.arctan2(center[-1] - Xs, Ys - center[-2]) return thetas diff --git a/src/vampires_dpp/pdi/models.py b/src/vampires_dpp/pdi/models.py index 3a2dcd3..0d7bfc9 100644 --- a/src/vampires_dpp/pdi/models.py +++ b/src/vampires_dpp/pdi/models.py @@ -35,7 +35,8 @@ def load_calibration_file(header): class VAMPIRESMuellerMatrix(BaseModel): name: str = "ideal" - pbs_ratio: float = 1 # cam1 / cam2 ratio + m3_diat: float = 0 + m3_offset: float = 0 # deg hwp_offset: float = 0 # deg hwp_phi: float = 0.5 # wave imr_offset: float = 0 # deg @@ -49,7 +50,8 @@ class VAMPIRESMuellerMatrix(BaseModel): def common_path_mm(self, pa, alt, az, hwp, imr, hwp_adi_sync=True): # telescope pa_theta = np.deg2rad(pa) - m3 = mm.generic(epsilon=0, delta=np.pi) + m3_theta = np.deg2rad(self.m3_offset) + m3 = mm.generic(epsilon=self.m3_diat, theta=m3_theta, delta=np.pi) alt_theta = np.deg2rad(alt) tel_mm = mm.rotator(-alt_theta) @ m3 @ mm.rotator(pa_theta) @@ -82,11 +84,9 @@ def __call__(self, flc_state: str, camera: int, *args, **kwargs) -> NDArray: flc_theta = np.deg2rad(self.flc_theta[flc_state]) flc_mm = mm.waveplate(flc_theta, self.flc_phi * 2 * np.pi) - # beamsplitter + # beamsplitter - horiztonal to camera 1 is_ordinary = camera == 1 pbs_mm = mm.wollaston(is_ordinary) - if is_ordinary: - pbs_mm *= self.pbs_ratio M = pbs_mm @ flc_mm @ cp_mm return M.astype("f4") @@ -111,7 +111,8 @@ def load_calib_file(cls, header): flc_theta = {k: t + table["flc_theta"] for k, t in zip(("A", "B"), (0, 45), strict=True)} return cls( name=table.name, - pbs_ratio=table["emgain"], + m3_diat=table["m3_diat"], + m3_offset=table["m3_theta"], hwp_offset=table["hwp_delta"], hwp_phi=table["hwp_phi"], imr_offset=table["imr_delta"], @@ -171,6 +172,8 @@ def load_calib_file(cls, header): table = load_calibration_file(header) return cls( name=table.name, + m3_diat=table["m3_diat"], + m3_offset=table["m3_theta"], hwp_offset=table["hwp_delta"], hwp_phi=table["hwp_phi"], imr_offset=table["imr_delta"], diff --git a/src/vampires_dpp/pdi/processing.py b/src/vampires_dpp/pdi/processing.py index 5efb623..1f0a0f2 100644 --- a/src/vampires_dpp/pdi/processing.py +++ b/src/vampires_dpp/pdi/processing.py @@ -343,7 +343,7 @@ def get_triplediff_set(table, row) -> dict | None: (table["RET-ANG1"] == key[0]) & (table["U_FLC"] == key[1]) & (table["U_CAMERA"] == key[2]) - & (delta_angs < 2) # 1 l/D rotation @ 45 l/d sep + & (delta_angs < 4) # 2 l/D rotation @ 45 l/d sep ) if np.any(mask): idx = delta_time[mask].argmin() @@ -351,7 +351,6 @@ def get_triplediff_set(table, row) -> dict | None: return None output_set[key] = table.loc[mask, "path"].iloc[idx] - return output_set @@ -507,6 +506,8 @@ def polarization_ip_correct(stokes_data, phot_rad, method, header=None): if header is not None: header["IP_PQ"] = pQ, "I -> Q IP correction value" header["IP_PU"] = pU, "I -> U IP correction value" + header["IP_POL"] = np.hypot(pU, pQ) * 100, "[%] Residual IP" + header["IP_ANG"] = 0.5 * np.rad2deg(np.arctan2(pU, pQ)), "[deg] Residual IP angle" header["IP_METH"] = method, "IP measurement method" return stokes_data, header diff --git a/src/vampires_dpp/pipeline/pipeline.py b/src/vampires_dpp/pipeline/pipeline.py index dcd6232..120fe8e 100644 --- a/src/vampires_dpp/pipeline/pipeline.py +++ b/src/vampires_dpp/pipeline/pipeline.py @@ -460,7 +460,6 @@ def polarimetry_difference(self, table, method, num_proc=None, force=False): full_paths.append(tuple(sorted(stokes_set.values()))) full_path_set = list(set(full_paths)) - stokes_files = [] # stokes_sets_path = self.paths.pdi_dir / f"{self.config.name}_stokes_sets.txt" # with stokes_sets_path.open("w") as fh: @@ -523,7 +522,12 @@ def polarimetry_difference(self, table, method, num_proc=None, force=False): ## Collapse outputs stokes_data = np.array(stokes_data) stokes_err = np.array(stokes_err) - coll_frame, _ = collapse_frames(stokes_data) + coll_frame, _ = collapse_frames(np.nan_to_num(stokes_data)) + footprint = np.mean(np.isfinite(stokes_data).astype("f4"), axis=0) + fits.writeto( + self.paths.pdi_dir / f"{self.config.name}_footprint.fits", footprint, overwrite=True + ) + # coll_frame *= footprint with warnings.catch_warnings(): warnings.simplefilter("ignore") coll_err = np.sqrt(np.nansum(stokes_err**2, axis=0)) / stokes_err.shape[0] diff --git a/src/vampires_dpp/wcs.py b/src/vampires_dpp/wcs.py index a2cb15a..eca57a8 100644 --- a/src/vampires_dpp/wcs.py +++ b/src/vampires_dpp/wcs.py @@ -47,10 +47,15 @@ def get_coord_header(header, time=None): GAIA_CATALOGS = {"dr1": "I/337/gaia", "dr2": "I/345/gaia2", "dr3": "I/355/gaiadr3"} -def get_gaia_astrometry(target, catalog="dr3", radius=1): +def get_gaia_astrometry(target: str, catalog="dr3", radius=1): + """ + Get coordinate from GAIA catalogue with proper motions and parallax + + Radius is in arcminute. + """ # get precise RA and DEC gaia_catalog_list = Vizier.query_object( - target, radius=radius * u.arcsec, catalog=GAIA_CATALOGS[catalog.lower()] + target, radius=radius * u.arcminute, catalog=GAIA_CATALOGS[catalog.lower()] ) if len(gaia_catalog_list) == 0: return None @@ -68,6 +73,7 @@ def get_gaia_astrometry(target, catalog="dr3", radius=1): return coord -def get_precise_coord(coord, time, scale="utc"): +def get_precise_coord(coord: SkyCoord, time: str, scale="utc"): + """Use astropy to get proper-motion corrected coordinate""" _time = Time(time, scale=scale, location=SUBARU_LOC) return coord.apply_space_motion(_time)