From 16c5968cdb55e5bc8fcd191cf3a0272109d03bfe Mon Sep 17 00:00:00 2001 From: hapipapijuris <“yamanaka.juri.h0@s.mail.nagoya-u.ac.jp”> Date: Sun, 17 Nov 2024 16:50:37 +0000 Subject: [PATCH 1/9] #220 Add ABBA method code --- decode/qlook.py | 110 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 110 insertions(+) diff --git a/decode/qlook.py b/decode/qlook.py index 1da71fa..78c9a23 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -403,6 +403,106 @@ def pswsc( return save_qlook(fig, file, overwrite=overwrite, **options) +def abba( + dems: Path, + /, + *, + # options for loading + include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, + exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, + min_frequency: Optional[str] = DEFAULT_MIN_FREQUENCY, + max_frequency: Optional[str] = DEFAULT_MAX_FREQUENCY, + data_type: Literal["auto", "brightness", "df/f"] = DEFAULT_DATA_TYPE, + frequency_units: str = DEFAULT_FREQUENCY_UNITS, + # options for saving + format: str = DEFAULT_FORMAT, + outdir: Path = DEFAULT_OUTDIR, + overwrite: bool = DEFAULT_OVERWRITE, + suffix: str = "abba_pswsc", + # other options + debug: bool = DEFAULT_DEBUG, + **options: Any, +) -> Path: + """Quick-look at a ABBA PSW observation with sky chopper. + + Args: + dems: Input DEMS file (netCDF or Zarr). + include_mkid_ids: MKID IDs to be included in analysis. + Defaults to all MKID IDs. + exclude_mkid_ids: MKID IDs to be excluded in analysis. + Defaults to no MKID IDs. + min_frequency: Minimum frequency to be included in analysis. + Defaults to no minimum frequency bound. + max_frequency: Maximum frequency to be included in analysis. + Defaults to no maximum frequency bound. + data_type: Data type of the input DEMS file. + Defaults to the ``long_name`` attribute in it. + frequency_units: Units of the frequency axis. + format: Output data format of the quick-look result. + outdir: Output directory for the quick-look result. + overwrite: Whether to overwrite the output if it exists. + suffix: Suffix that precedes the file extension. + debug: Whether to print detailed logs for debugging. + **options: Other options for saving the output (e.g. dpi). + + Returns: + Absolute path of the saved file. + """ + + with set_logger(debug): + for key, val in locals().items(): + LOGGER.debug(f"{key}: {val!r}") + + with xr.set_options(keep_attrs=True): + da = load_dems( + dems, + include_mkid_ids=include_mkid_ids, + exclude_mkid_ids=exclude_mkid_ids, + min_frequency=min_frequency, + max_frequency=max_frequency, + data_type=data_type, + frequency_units=frequency_units, + ) + + da_despike = despike(da) + + # extracting nodAB dataset + da_onoff = select.by(da_despike, "state", ["ON", "OFF"]) + # make cuntinus ID label of nod (Equal to scan label,but the latter is discrete) + scan_onoff = utils.phaseof(da_onoff.state) + + chop_per_scan = da_onoff.beam.groupby(scan_onoff).apply(utils.phaseof) + is_second_half = chop_per_scan.groupby(scan_onoff).apply( + lambda group: (group >= group.mean()) + ) + + # make is_second_half coord usinf the func each_half_cal + # make abba_cycle label and abba_phase using quotient and remainder devided by 4. + abba_cycle = (scan_onoff.data * 2 + is_second_half - 1) // 4 + abba_phase = (scan_onoff * 2 + is_second_half - 1) % 4 + da_abba = da_onoff.assign_coords(abba_cycle=abba_cycle, abba_phase=abba_phase) + da_abba_select = ( + da_abba.groupby("abba_cycle") + .map(subtract_per_abba_cycle) + .mean("abba_cycle") + ) + + # save result + suffixes = f".{suffix}.{format}" + file = Path(outdir) / Path(dems).with_suffix(suffixes).name + if format in DATA_FORMATS: + return save_qlook(da_abba_select, file, overwrite=overwrite, **options) + + fig, ax = plt.subplots(figsize=(6, 4)) + plot.data(da_abba_select, x="frequency", s=5, hue=None) + ax.set_ylim(get_robust_lim(da_abba_select)) + ax.grid(True) + ax.set_title(f"{Path(dems).name}\n({da.observation})") + + fig.tight_layout() + return save_qlook(fig, file, overwrite=overwrite, **options) + + def raster( dems: Path, /, @@ -1160,6 +1260,16 @@ def subtract_per_scan(dems: xr.DataArray) -> xr.DataArray: raise ValueError("State must be either ON or OFF.") +def subtract_per_abba_cycle(dems: xr.DataArray) -> xr.DataArray: + """Apply abba atm-corrction to a single-scan DEMS.""" + required_elements = {0, 1, 2, 3} + if required_elements == set(np.unique(dems.abba_phase.values)): + per_abba = dems.groupby("abba_phase").map(subtract_per_scan) + return per_abba.mean("abba_phase") + else: + return xr.DataArray(np.nan) + + def get_chan_weight( dems: xr.DataArray, /, From 8add9d44626bad1713c21137a36e4824c95c1e0d Mon Sep 17 00:00:00 2001 From: hapipapijuris <“yamanaka.juri.h0@s.mail.nagoya-u.ac.jp”> Date: Sun, 17 Nov 2024 17:31:19 +0000 Subject: [PATCH 2/9] #220 Clean up ABBA method code --- decode/qlook.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 78c9a23..6b653c8 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -2,6 +2,7 @@ "auto", "daisy", "pswsc", + "abba", "raster", "skydip", "still", @@ -89,6 +90,9 @@ def auto(dems: Path, /, **options: Any) -> Path: if "pswsc" in obs: return pswsc(dems, **options) + if "abba" in obs: + return abba(dems, **options) + if "raster" in obs: return raster(dems, **options) @@ -418,7 +422,7 @@ def abba( format: str = DEFAULT_FORMAT, outdir: Path = DEFAULT_OUTDIR, overwrite: bool = DEFAULT_OVERWRITE, - suffix: str = "abba_pswsc", + suffix: str = "abba", # other options debug: bool = DEFAULT_DEBUG, **options: Any, @@ -466,18 +470,13 @@ def abba( da_despike = despike(da) - # extracting nodAB dataset + # make spectrum da_onoff = select.by(da_despike, "state", ["ON", "OFF"]) - # make cuntinus ID label of nod (Equal to scan label,but the latter is discrete) scan_onoff = utils.phaseof(da_onoff.state) - chop_per_scan = da_onoff.beam.groupby(scan_onoff).apply(utils.phaseof) is_second_half = chop_per_scan.groupby(scan_onoff).apply( lambda group: (group >= group.mean()) ) - - # make is_second_half coord usinf the func each_half_cal - # make abba_cycle label and abba_phase using quotient and remainder devided by 4. abba_cycle = (scan_onoff.data * 2 + is_second_half - 1) // 4 abba_phase = (scan_onoff * 2 + is_second_half - 1) % 4 da_abba = da_onoff.assign_coords(abba_cycle=abba_cycle, abba_phase=abba_phase) @@ -490,10 +489,12 @@ def abba( # save result suffixes = f".{suffix}.{format}" file = Path(outdir) / Path(dems).with_suffix(suffixes).name + if format in DATA_FORMATS: return save_qlook(da_abba_select, file, overwrite=overwrite, **options) - fig, ax = plt.subplots(figsize=(6, 4)) + fig, ax = plt.subplots(figsize=(6, 4)) # type: ignore + plot.data(da_abba_select, x="frequency", s=5, hue=None) ax.set_ylim(get_robust_lim(da_abba_select)) ax.grid(True) @@ -1492,6 +1493,7 @@ def main() -> None: "auto": auto, "daisy": daisy, "pswsc": pswsc, + "abba": abba, "raster": raster, "skydip": skydip, "still": still, From df8154da454d3a68b5edb5c8017716335ef95cbf Mon Sep 17 00:00:00 2001 From: hapipapijuris <“yamanaka.juri.h0@s.mail.nagoya-u.ac.jp”> Date: Mon, 18 Nov 2024 11:54:34 +0000 Subject: [PATCH 3/9] #220 incorporate ABBA-based analysis code into pswsc func --- decode/qlook.py | 123 +++++++----------------------------------------- 1 file changed, 16 insertions(+), 107 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 6b653c8..349cb2b 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -2,7 +2,6 @@ "auto", "daisy", "pswsc", - "abba", "raster", "skydip", "still", @@ -50,6 +49,7 @@ DEFAULT_SKYCOORD_UNITS = "arcsec" SIGMA_OVER_MAD = 1.4826 LOGGER = getLogger(__name__) +ABBA_PHASE = {0, 1, 2, 3} @contextmanager @@ -90,9 +90,6 @@ def auto(dems: Path, /, **options: Any) -> Path: if "pswsc" in obs: return pswsc(dems, **options) - if "abba" in obs: - return abba(dems, **options) - if "raster" in obs: return raster(dems, **options) @@ -379,9 +376,20 @@ def pswsc( da_despiked = despike(da) # make spectrum - da_scan = select.by(da_despiked, "state", ["ON", "OFF"]) - da_sub = da_scan.groupby("scan").map(subtract_per_scan) - spec = da_sub.mean("scan") + da_onoff = select.by(da_despiked, "state", ["ON", "OFF"]) + scan_onoff = utils.phaseof(da_onoff.state) + chop_per_scan = da_onoff.beam.groupby(scan_onoff).apply(utils.phaseof) + is_second_half = chop_per_scan.groupby(scan_onoff).apply( + lambda group: (group >= group.mean()) + ) + abba_cycle = (scan_onoff.data * 2 + is_second_half - 1) // 4 + abba_phase = (scan_onoff.data * 2 + is_second_half - 1) % 4 + da_abba = da_onoff.assign_coords(abba_cycle=abba_cycle, abba_phase=abba_phase) + spec = ( + da_abba.groupby("abba_cycle") + .map(subtract_per_abba_cycle) + .mean("abba_cycle") + ) # save result suffixes = f".{suffix}.{format}" @@ -407,103 +415,6 @@ def pswsc( return save_qlook(fig, file, overwrite=overwrite, **options) -def abba( - dems: Path, - /, - *, - # options for loading - include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, - exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, - min_frequency: Optional[str] = DEFAULT_MIN_FREQUENCY, - max_frequency: Optional[str] = DEFAULT_MAX_FREQUENCY, - data_type: Literal["auto", "brightness", "df/f"] = DEFAULT_DATA_TYPE, - frequency_units: str = DEFAULT_FREQUENCY_UNITS, - # options for saving - format: str = DEFAULT_FORMAT, - outdir: Path = DEFAULT_OUTDIR, - overwrite: bool = DEFAULT_OVERWRITE, - suffix: str = "abba", - # other options - debug: bool = DEFAULT_DEBUG, - **options: Any, -) -> Path: - """Quick-look at a ABBA PSW observation with sky chopper. - - Args: - dems: Input DEMS file (netCDF or Zarr). - include_mkid_ids: MKID IDs to be included in analysis. - Defaults to all MKID IDs. - exclude_mkid_ids: MKID IDs to be excluded in analysis. - Defaults to no MKID IDs. - min_frequency: Minimum frequency to be included in analysis. - Defaults to no minimum frequency bound. - max_frequency: Maximum frequency to be included in analysis. - Defaults to no maximum frequency bound. - data_type: Data type of the input DEMS file. - Defaults to the ``long_name`` attribute in it. - frequency_units: Units of the frequency axis. - format: Output data format of the quick-look result. - outdir: Output directory for the quick-look result. - overwrite: Whether to overwrite the output if it exists. - suffix: Suffix that precedes the file extension. - debug: Whether to print detailed logs for debugging. - **options: Other options for saving the output (e.g. dpi). - - Returns: - Absolute path of the saved file. - """ - - with set_logger(debug): - for key, val in locals().items(): - LOGGER.debug(f"{key}: {val!r}") - - with xr.set_options(keep_attrs=True): - da = load_dems( - dems, - include_mkid_ids=include_mkid_ids, - exclude_mkid_ids=exclude_mkid_ids, - min_frequency=min_frequency, - max_frequency=max_frequency, - data_type=data_type, - frequency_units=frequency_units, - ) - - da_despike = despike(da) - - # make spectrum - da_onoff = select.by(da_despike, "state", ["ON", "OFF"]) - scan_onoff = utils.phaseof(da_onoff.state) - chop_per_scan = da_onoff.beam.groupby(scan_onoff).apply(utils.phaseof) - is_second_half = chop_per_scan.groupby(scan_onoff).apply( - lambda group: (group >= group.mean()) - ) - abba_cycle = (scan_onoff.data * 2 + is_second_half - 1) // 4 - abba_phase = (scan_onoff * 2 + is_second_half - 1) % 4 - da_abba = da_onoff.assign_coords(abba_cycle=abba_cycle, abba_phase=abba_phase) - da_abba_select = ( - da_abba.groupby("abba_cycle") - .map(subtract_per_abba_cycle) - .mean("abba_cycle") - ) - - # save result - suffixes = f".{suffix}.{format}" - file = Path(outdir) / Path(dems).with_suffix(suffixes).name - - if format in DATA_FORMATS: - return save_qlook(da_abba_select, file, overwrite=overwrite, **options) - - fig, ax = plt.subplots(figsize=(6, 4)) # type: ignore - - plot.data(da_abba_select, x="frequency", s=5, hue=None) - ax.set_ylim(get_robust_lim(da_abba_select)) - ax.grid(True) - ax.set_title(f"{Path(dems).name}\n({da.observation})") - - fig.tight_layout() - return save_qlook(fig, file, overwrite=overwrite, **options) - - def raster( dems: Path, /, @@ -1263,8 +1174,7 @@ def subtract_per_scan(dems: xr.DataArray) -> xr.DataArray: def subtract_per_abba_cycle(dems: xr.DataArray) -> xr.DataArray: """Apply abba atm-corrction to a single-scan DEMS.""" - required_elements = {0, 1, 2, 3} - if required_elements == set(np.unique(dems.abba_phase.values)): + if ABBA_PHASE == set(np.unique(dems.abba_phase.values)): per_abba = dems.groupby("abba_phase").map(subtract_per_scan) return per_abba.mean("abba_phase") else: @@ -1493,7 +1403,6 @@ def main() -> None: "auto": auto, "daisy": daisy, "pswsc": pswsc, - "abba": abba, "raster": raster, "skydip": skydip, "still": still, From 8009ddde7ac701f21c9f0d39f86fe652331b3571 Mon Sep 17 00:00:00 2001 From: hapipapijuris <“yamanaka.juri.h0@s.mail.nagoya-u.ac.jp”> Date: Tue, 19 Nov 2024 01:21:41 +0000 Subject: [PATCH 4/9] #220 minor change of ABBA code --- decode/qlook.py | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 349cb2b..3806c62 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -32,6 +32,7 @@ # constants +ABBA_PHASES = {0, 1, 2, 3} DATA_FORMATS = "csv", "nc", "zarr", "zarr.zip" DEFAULT_DATA_TYPE = "auto" DEFAULT_DEBUG = False @@ -49,7 +50,6 @@ DEFAULT_SKYCOORD_UNITS = "arcsec" SIGMA_OVER_MAD = 1.4826 LOGGER = getLogger(__name__) -ABBA_PHASE = {0, 1, 2, 3} @contextmanager @@ -1145,8 +1145,21 @@ def mean_in_time(dems: xr.DataArray) -> xr.DataArray: return xr.zeros_like(middle) + dems.mean("time") -def subtract_per_scan(dems: xr.DataArray) -> xr.DataArray: - """Apply source-sky subtraction to a single-scan DEMS.""" +def subtract_per_abba_cycle(dems: xr.DataArray, /) -> xr.DataArray: + """Take the whole abba-cycle average which are applied atm corecction""" + if not set(np.unique(dems.abba_phase)) == ABBA_PHASES: + return xr.DataArray(np.nan) + + return ( + dems + .groupby("abba_phase") + .map(subtract_per_abba_phase) + .mean("abba_phase") + ) + + +def subtract_per_abba_phase(dems: xr.DataArray) -> xr.DataArray: + """Apply atm corecction to each abba-phase DEMS.""" t_amb = 273.15 if len(states := np.unique(dems.state)) != 1: raise ValueError("State must be unique.") @@ -1172,15 +1185,6 @@ def subtract_per_scan(dems: xr.DataArray) -> xr.DataArray: raise ValueError("State must be either ON or OFF.") -def subtract_per_abba_cycle(dems: xr.DataArray) -> xr.DataArray: - """Apply abba atm-corrction to a single-scan DEMS.""" - if ABBA_PHASE == set(np.unique(dems.abba_phase.values)): - per_abba = dems.groupby("abba_phase").map(subtract_per_scan) - return per_abba.mean("abba_phase") - else: - return xr.DataArray(np.nan) - - def get_chan_weight( dems: xr.DataArray, /, From df80d7686e8982b9c241864cbfc98d944cfdcfe3 Mon Sep 17 00:00:00 2001 From: hapipapijuris <“yamanaka.juri.h0@s.mail.nagoya-u.ac.jp”> Date: Tue, 19 Nov 2024 01:38:55 +0000 Subject: [PATCH 5/9] #220 reformated 1 file --- decode/qlook.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 3806c62..70d934e 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -1150,12 +1150,7 @@ def subtract_per_abba_cycle(dems: xr.DataArray, /) -> xr.DataArray: if not set(np.unique(dems.abba_phase)) == ABBA_PHASES: return xr.DataArray(np.nan) - return ( - dems - .groupby("abba_phase") - .map(subtract_per_abba_phase) - .mean("abba_phase") - ) + return dems.groupby("abba_phase").map(subtract_per_abba_phase).mean("abba_phase") def subtract_per_abba_phase(dems: xr.DataArray) -> xr.DataArray: From 7c23dca12d623d564fcdfb3ba48199982027269e Mon Sep 17 00:00:00 2001 From: hapipapijuris <“yamanaka.juri.h0@s.mail.nagoya-u.ac.jp”> Date: Tue, 19 Nov 2024 21:36:57 +0000 Subject: [PATCH 6/9] #220 change explanation and data type --- decode/qlook.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 70d934e..b7fd270 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -382,8 +382,8 @@ def pswsc( is_second_half = chop_per_scan.groupby(scan_onoff).apply( lambda group: (group >= group.mean()) ) - abba_cycle = (scan_onoff.data * 2 + is_second_half - 1) // 4 - abba_phase = (scan_onoff.data * 2 + is_second_half - 1) % 4 + abba_cycle = (scan_onoff * 2 + is_second_half - 1) // 4 + abba_phase = (scan_onoff * 2 + is_second_half - 1) % 4 da_abba = da_onoff.assign_coords(abba_cycle=abba_cycle, abba_phase=abba_phase) spec = ( da_abba.groupby("abba_cycle") @@ -1146,7 +1146,7 @@ def mean_in_time(dems: xr.DataArray) -> xr.DataArray: def subtract_per_abba_cycle(dems: xr.DataArray, /) -> xr.DataArray: - """Take the whole abba-cycle average which are applied atm corecction""" + """Take each abba-cycle's average which are applied atm corecction""" if not set(np.unique(dems.abba_phase)) == ABBA_PHASES: return xr.DataArray(np.nan) From fcb194c2b372d49d0fb252df186deea600cb5bdf Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Thu, 21 Nov 2024 07:28:20 +0000 Subject: [PATCH 7/9] #220 Refactor subtract_per_abba_cycle --- decode/qlook.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index b7fd270..f9311dc 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -1146,9 +1146,19 @@ def mean_in_time(dems: xr.DataArray) -> xr.DataArray: def subtract_per_abba_cycle(dems: xr.DataArray, /) -> xr.DataArray: - """Take each abba-cycle's average which are applied atm corecction""" + """Subtract sky from source with atmospheric correction for each ABBA cycle. + + Args: + dems: 2D DataArray (time x chan) of DEMS per ABBA cycle. + + Returns: + 1D DataArray (chan) of the mean spectrum after subtraction and correction. + If ABBA phases per cycle are incomplete, i.e., some phases are missing, + a spectrum filled with NaN will be returned instead. + + """ if not set(np.unique(dems.abba_phase)) == ABBA_PHASES: - return xr.DataArray(np.nan) + return dems.mean("time") * np.nan return dems.groupby("abba_phase").map(subtract_per_abba_phase).mean("abba_phase") From 2621436a28fcff13a9d701d94212f466e552f6e5 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Thu, 21 Nov 2024 07:28:48 +0000 Subject: [PATCH 8/9] #220 Refactor subtract_per_abba_phase --- decode/qlook.py | 43 +++++++++++++++++++++++-------------------- 1 file changed, 23 insertions(+), 20 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index f9311dc..dffa6a6 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -1163,31 +1163,34 @@ def subtract_per_abba_cycle(dems: xr.DataArray, /) -> xr.DataArray: return dems.groupby("abba_phase").map(subtract_per_abba_phase).mean("abba_phase") -def subtract_per_abba_phase(dems: xr.DataArray) -> xr.DataArray: - """Apply atm corecction to each abba-phase DEMS.""" +def subtract_per_abba_phase(dems: xr.DataArray, /) -> xr.DataArray: + """Subtract sky from source with atmospheric correction for each ABBA phase. + + Args: + dems: 2D DataArray (time x chan) of DEMS per ABBA phase. + + Returns: + 1D DataArray (chan) of the mean spectrum after subtraction and correction. + + Raises: + ValueError: Raised if ``dems.state`` is not ON-only nor OFF-only. + + """ t_amb = 273.15 + if len(states := np.unique(dems.state)) != 1: raise ValueError("State must be unique.") - if (state := states[0]) == "ON": - src = select.by(dems, "beam", include="A") - sky = select.by(dems, "beam", include="B") - return ( - t_amb - * (src.mean("time") - sky.mean("time").data) - / ((t_amb - sky.mean("time"))) - ) - - if state == "OFF": - src = select.by(dems, "beam", include="B") - sky = select.by(dems, "beam", include="A") - return ( - t_amb - * (src.mean("time") - sky.mean("time").data) - / ((t_amb - sky.mean("time"))) - ) + if states[0] == "ON": + src = select.by(dems, "beam", include="A").mean("time") + sky = select.by(dems, "beam", include="B").mean("time") + elif states[0] == "OFF": + src = select.by(dems, "beam", include="B").mean("time") + sky = select.by(dems, "beam", include="A").mean("time") + else: + raise ValueError("State must be either ON or OFF.") - raise ValueError("State must be either ON or OFF.") + return t_amb * (src - sky.data) / (t_amb - sky) def get_chan_weight( From a27e853dad05c682792c4337be4e0515019c6971 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Thu, 21 Nov 2024 07:51:40 +0000 Subject: [PATCH 9/9] #220 Refactor pswsc --- decode/qlook.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index dffa6a6..328549d 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -375,7 +375,7 @@ def pswsc( da_despiked = despike(da) - # make spectrum + # calculate ABBA cycles and phases da_onoff = select.by(da_despiked, "state", ["ON", "OFF"]) scan_onoff = utils.phaseof(da_onoff.state) chop_per_scan = da_onoff.beam.groupby(scan_onoff).apply(utils.phaseof) @@ -384,9 +384,11 @@ def pswsc( ) abba_cycle = (scan_onoff * 2 + is_second_half - 1) // 4 abba_phase = (scan_onoff * 2 + is_second_half - 1) % 4 - da_abba = da_onoff.assign_coords(abba_cycle=abba_cycle, abba_phase=abba_phase) + + # make spectrum spec = ( - da_abba.groupby("abba_cycle") + da_onoff.assign_coords(abba_cycle=abba_cycle, abba_phase=abba_phase) + .groupby("abba_cycle") .map(subtract_per_abba_cycle) .mean("abba_cycle") )