From 223ccf5b56fb4ef587c3e249e32ae794afa179d9 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sat, 18 Nov 2023 21:54:22 +0000 Subject: [PATCH 01/31] #135 Make include positional argument in select.by --- decode/select.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/decode/select.py b/decode/select.py index 53445da..8e816ad 100644 --- a/decode/select.py +++ b/decode/select.py @@ -20,11 +20,11 @@ def by( dems: xr.DataArray, coord_name: str, /, + include: Optional[Multiple[T]] = None, *, + exclude: Optional[Multiple[T]] = None, min: Optional[T] = None, max: Optional[T] = None, - include: Optional[Multiple[T]] = None, - exclude: Optional[Multiple[T]] = None, sort: bool = False, as_dim: bool = False, ) -> xr.DataArray: @@ -33,14 +33,14 @@ def by( Args: dems: DEMS DataArray to be selected. coord_name: Name of the coordinate for the selection. - min: Minimum selection bound (inclusive). - If not specified, no bound is set. - max: Maximum selection bound (exclusive). - If not specified, no bound is set. include: Coordinate values to be included. If not specified, all values are included. exclude: Coordinate values to be excluded. If not specified, any values are not excluded. + min: Minimum selection bound (inclusive). + If not specified, no bound is set. + max: Maximum selection bound (exclusive). + If not specified, no bound is set. sort: Whether to sort by the coordinate after selection. as_dim: Whether to use the coordinate as a dimension. @@ -58,18 +58,18 @@ def by( coord_dim = coord.dims[0] - if min is not None: - dems = dems.sel({coord_dim: dems[coord_name] >= min}) - - if max is not None: - dems = dems.sel({coord_dim: dems[coord_name] < max}) - if include is not None: dems = dems.sel({coord_dim: dems[coord_name].isin(include)}) if exclude is not None: dems = dems.sel({coord_dim: ~dems[coord_name].isin(exclude)}) + if min is not None: + dems = dems.sel({coord_dim: dems[coord_name] >= min}) + + if max is not None: + dems = dems.sel({coord_dim: dems[coord_name] < max}) + if sort: dems = dems.sortby(coord_name) From 2e642b9a117f04a52e380c40f58ded526b95debd Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sat, 18 Nov 2023 21:57:30 +0000 Subject: [PATCH 02/31] #135 Fix wrong (reversed) selections of src and sky --- decode/qlook.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 16c492c..bdb146e 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -535,13 +535,13 @@ def subtract_per_scan(dems: xr.DataArray) -> xr.DataArray: raise ValueError("State must be unique.") if (state := states[0]) == "ON": - src = select.by(dems, "beam", include="B") - sky = select.by(dems, "beam", include="A") + src = select.by(dems, "beam", include="A") + sky = select.by(dems, "beam", include="B") return src.mean("time") - sky.mean("time").data if state == "OFF": - src = select.by(dems, "beam", include="A") - sky = select.by(dems, "beam", include="B") + src = select.by(dems, "beam", include="B") + sky = select.by(dems, "beam", include="A") return src.mean("time") - sky.mean("time").data raise ValueError("State must be either ON or OFF.") From 154be123d5da27662eddc2c68d7e7aa17c01e459 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sat, 18 Nov 2023 22:28:47 +0000 Subject: [PATCH 03/31] #135 Fix assign.scan --- decode/assign.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/decode/assign.py b/decode/assign.py index ad73743..07d8b05 100644 --- a/decode/assign.py +++ b/decode/assign.py @@ -14,7 +14,7 @@ def scan( dems: xr.DataArray, /, *, - by: Literal["state"] = "state", + by: Literal["beam", "scan", "state"] = "state", dt: Optional[np.timedelta64] = None, inplace: bool = False, ) -> xr.DataArray: @@ -33,16 +33,16 @@ def scan( """ if not inplace: - dems = dems.copy() + # deepcopy except for data + dems = dems.copy(data=dems.data) - is_cut = np.zeros_like(dems.scan, dtype=bool) + is_div = xr.zeros_like(dems.scan, dtype=bool) - if by == "state": - state = dems.state.values - is_cut |= np.hstack([False, state[1:] != state[:-1]]) + ref = dems.coords[by].data + is_div[1:] |= ref[1:] != ref[:-1] if dt is not None: - is_cut |= np.hstack([False, np.diff(dems.time) >= dt]) + is_div[1:] |= np.diff(dems.time) >= dt - dems.scan.values[:] = np.cumsum(is_cut) - return dems + new_scan = is_div.cumsum().astype(dems.scan.dtype) + return dems.assign_coords(scan=new_scan) From ebdc85abe165878f67e721f7aab90d8dd2388d6f Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sun, 19 Nov 2023 01:17:09 +0000 Subject: [PATCH 04/31] #135 Remove inplace option from convert module --- decode/convert.py | 78 +++++++++++++++++++++++------------------------ 1 file changed, 38 insertions(+), 40 deletions(-) diff --git a/decode/convert.py b/decode/convert.py index 5399289..23b91b1 100644 --- a/decode/convert.py +++ b/decode/convert.py @@ -2,7 +2,7 @@ # standard library -from typing import Optional, Union +from typing import Optional, Sequence, TypeVar, Union # dependencies @@ -12,22 +12,23 @@ # type hints +T = TypeVar("T") +Multiple = Union[Sequence[T], T] UnitLike = Union[Unit, str] def coord_units( da: xr.DataArray, - coord_name: str, + coord_names: Multiple[str], new_units: UnitLike, /, - *, equivalencies: Optional[Equivalency] = None, ) -> xr.DataArray: - """Convert units of a coordinate of a DataArray. + """Convert the units of coordinate(s) of a DataArray. Args: da: Input DataArray. - coord_name: Name of the coordinate to be converted. + coord_names: Name(s) of the coordinate(s) to be converted. new_units: Units to be converted from the current ones. equivalencies: Optional Astropy equivalencies. @@ -35,63 +36,60 @@ def coord_units( DataArray with the units of the coordinate converted. """ - new_coord = units(da[coord_name], new_units, equivalencies=equivalencies) - return da.assign_coords({coord_name: new_coord}) + # deepcopy except for data + da = da.copy(data=da.data) + if isinstance(coord_names, str): + coord_names = [coord_names] + + for coord_name in coord_names: + coord = da.coords[coord_name] + new_coord = units(coord, new_units, equivalencies) + da = da.assign_coords({coord_name: new_coord}) + + return da -def frame( - dems: xr.DataArray, - new_frame: str, - /, - *, - inplace: bool = False, -) -> xr.DataArray: - """Convert skycoord frame of DEMS. + +def frame(da: xr.DataArray, new_frame: str, /) -> xr.DataArray: + """Convert the skycoord frame of a DataArray. Args: - dems: Target DEMS DataArray. - new_frame: Skycoord frame to be converted from the current ones. - inplace: Whether the skycoord frame are converted in-place. + da: Input DataArray. + new_frame: Frame to be converted from the current one. Returns: - DEMS DataArray with the skycoord frame converted. + DataArray with the skycoord frame converted. """ - if not inplace: - # deepcopy except for data - dems = dems.copy(data=dems.data) + # deepcopy except for data + da = da.copy(data=da.data) if not new_frame == "relative": raise ValueError("Relative is only available.") - lon, lon_origin = dems["lon"], dems["lon_origin"] - lat, lat_origin = dems["lat"], dems["lat_origin"] - cos = np.cos(Quantity(lat, lat.attrs["units"]).to("rad")) - - if lon.attrs["units"] != lon_origin.attrs["units"]: - raise ValueError("Units of lon and lon_origin must be same.") - - if lat.attrs["units"] != lat_origin.attrs["units"]: - raise ValueError("Units of lat and lat_origin must be same.") + lon = da.coords["lon"] + lat = da.coords["lat"] + lon_origin = da.coords["lon_origin"] + lat_origin = da.coords["lat_origin"] # do not change the order below! - lon -= lon_origin - lon *= cos - lat -= lat_origin - lon_origin[:] = 0.0 - lat_origin[:] = 0.0 + lon -= units(lon_origin, lon.attrs["units"]) + lon *= np.cos(units(lat, "rad")) + lat -= units(lat_origin, lat.attrs["units"]) + lon_origin *= 0.0 + lat_origin *= 0.0 - return dems.assign_coords(frame=dems.frame.copy(False, new_frame)) + new_frame = da.frame.copy(data=new_frame) + return da.assign_coords(frame=new_frame) def units( da: xr.DataArray, new_units: UnitLike, /, - *, equivalencies: Optional[Equivalency] = None, ) -> xr.DataArray: - """Convert units of a DataArray. + """Convert the units of a DataArray. Args: da: Input DataArray. @@ -106,4 +104,4 @@ def units( raise ValueError("Units must exist in DataArray attrs.") new_data = Quantity(da, units).to(new_units, equivalencies) - return da.copy(False, new_data).assign_attrs(units=str(new_units)) + return da.copy(data=new_data).assign_attrs(units=new_units) From 5814edc41570b594293333f7fd5d3b88a4b21f34 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sun, 19 Nov 2023 02:04:36 +0000 Subject: [PATCH 05/31] #135 Load DEMS with common loader --- decode/qlook.py | 210 ++++++++++++++++++++++++++---------------------- 1 file changed, 115 insertions(+), 95 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index bdb146e..9cf22de 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -15,13 +15,18 @@ # constants +DEFAULT_DATA_TYPE = "brightness" # fmt: off -BAD_MKID_IDS = ( +DEFAULT_EXCL_MKID_IDS = ( 18, 77, 117, 118, 140, 141, 161, 182, 183, 184, 201, 209, 211, 232, 233, 239, 258, 259, 278, 282, 283, 296, 297, 299, 301, 313, ) # fmt: on +DEFAULT_FREQUENCY_UNITS = "GHz" +DEFAULT_INCL_MKID_IDS = None +DEFAULT_SKYCOORD_GRID = "6 arcsec" +DEFAULT_SKYCOORD_UNITS = "arcsec" DFOF_TO_TSKY = (300 - 77) / 3e-5 TSKY_TO_DFOF = 3e-5 / (300 - 77) @@ -30,9 +35,9 @@ def still( dems: Path, /, *, - include_mkid_ids: Optional[Sequence[int]] = None, - exclude_mkid_ids: Optional[Sequence[int]] = BAD_MKID_IDS, - data_type: Literal["df/f", "brightness"] = "brightness", + include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, + exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, + data_type: Literal["df/f", "brightness"] = DEFAULT_DATA_TYPE, 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", cabin_temperature: float = 273.0, @@ -65,20 +70,11 @@ def still( out = Path(outdir) / dems.with_suffix(f".still.{format}").name # load DEMS - da = load.dems(dems, chunks=None) - da = assign.scan(da) - da = convert.frame(da, "relative") - - if data_type == "df/f": - da.attrs.update(long_name="df/f", units="dimensionless") - - # select DEMS - da = select.by(da, "d2_mkid_type", include="filter") - da = select.by( - da, - "d2_mkid_id", - include=include_mkid_ids, - exclude=exclude_mkid_ids, + da = load_dems( + dems, + include_mkid_ids=include_mkid_ids, + exclude_mkid_ids=exclude_mkid_ids, + data_type=data_type, ) da_off = select.by(da, "state", exclude=["ON", "SCAN"]) @@ -126,9 +122,9 @@ def pswsc( dems: Path, /, *, - include_mkid_ids: Optional[Sequence[int]] = None, - exclude_mkid_ids: Optional[Sequence[int]] = BAD_MKID_IDS, - data_type: Literal["df/f", "brightness"] = "brightness", + include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, + exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, + data_type: Literal["df/f", "brightness"] = DEFAULT_DATA_TYPE, frequency_units: str = "GHz", outdir: Path = Path(), format: str = "png", @@ -151,26 +147,15 @@ def pswsc( out = Path(outdir) / dems.with_suffix(f".pswsc.{format}").name # load DEMS - da = load.dems(dems, chunks=None) - da = assign.scan(da) - da = convert.frame(da, "relative") - da = convert.coord_units(da, "frequency", frequency_units) - da = convert.coord_units(da, "d2_mkid_frequency", frequency_units) - - if data_type == "df/f": - da = cast(xr.DataArray, np.abs(da)) - da.attrs.update(long_name="|df/f|", units="dimensionless") - - # select DEMS - da = select.by(da, "d2_mkid_type", include="filter") - da = select.by( - da, - "d2_mkid_id", - include=include_mkid_ids, - exclude=exclude_mkid_ids, + da = load_dems( + dems, + include_mkid_ids=include_mkid_ids, + exclude_mkid_ids=exclude_mkid_ids, + data_type=data_type, + frequency_units=frequency_units, ) - da = select.by(da, "state", include=["ON", "OFF"]) - da_sub = da.groupby("scan").map(subtract_per_scan) + da_scan = select.by(da, "state", ["ON", "OFF"]) + da_sub = da_scan.groupby("scan").map(subtract_per_scan) # export output spec = da_sub.mean("scan") @@ -216,13 +201,13 @@ def raster( dems: Path, /, *, - include_mkid_ids: Optional[Sequence[int]] = None, - exclude_mkid_ids: Optional[Sequence[int]] = BAD_MKID_IDS, - data_type: Literal["df/f", "brightness"] = "brightness", + include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, + exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, + data_type: Literal["df/f", "brightness"] = DEFAULT_DATA_TYPE, 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", - skycoord_grid: str = "6 arcsec", - skycoord_units: str = "arcsec", + skycoord_grid: str = DEFAULT_SKYCOORD_GRID, + skycoord_units: str = DEFAULT_SKYCOORD_UNITS, outdir: Path = Path(), format: str = "png", ) -> None: @@ -252,20 +237,12 @@ def raster( result = Path(outdir) / dems.with_suffix(f".raster.{format}").name # load DEMS - da = load.dems(dems, chunks=None) - da = assign.scan(da) - da = convert.frame(da, "relative") - - if data_type == "df/f": - da.attrs.update(long_name="df/f", units="dimensionless") - - # select DEMS - da = select.by(da, "d2_mkid_type", include="filter") - da = select.by( - da, - "d2_mkid_id", - include=include_mkid_ids, - exclude=exclude_mkid_ids, + da = load_dems( + dems, + include_mkid_ids=include_mkid_ids, + exclude_mkid_ids=exclude_mkid_ids, + data_type=data_type, + skycoord_units=skycoord_units, ) da_on = select.by(da, "state", include="SCAN") da_off = select.by(da, "state", exclude="SCAN") @@ -330,9 +307,9 @@ def skydip( dems: Path, /, *, - include_mkid_ids: Optional[Sequence[int]] = None, - exclude_mkid_ids: Optional[Sequence[int]] = BAD_MKID_IDS, - data_type: Literal["df/f", "brightness"] = "brightness", + include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, + exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, + data_type: Literal["df/f", "brightness"] = DEFAULT_DATA_TYPE, 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", outdir: Path = Path(), @@ -362,24 +339,19 @@ def skydip( result = Path(outdir) / dems.with_suffix(f".skydip.{format}").name # load DEMS - da = load.dems(dems, chunks=None) - - if data_type == "df/f": - da = cast(xr.DataArray, np.abs(da)) - da.attrs.update(long_name="|df/f|", units="dimensionless") - - # add sec(Z) coordinate - secz = 1 / np.cos(np.deg2rad(90.0 - da.lat)) - secz.attrs.update(long_name="sec(Z)", units="dimensionless") - da = da.assign_coords(secz=secz) - - # select DEMS - da = select.by(da, "d2_mkid_type", include="filter") - da = select.by( - da, - "d2_mkid_id", - include=include_mkid_ids, - exclude=exclude_mkid_ids, + da = load_dems( + dems, + include_mkid_ids=include_mkid_ids, + exclude_mkid_ids=exclude_mkid_ids, + data_type=data_type, + ) + z = np.pi / 2 - convert.units(da.lat, "rad") + secz = cast(xr.DataArray, 1 / np.cos(z)) + da = da.assign_coords( + secz=secz.assign_attrs( + long_name="sec(Z)", + units="dimensionless", + ) ) da_on = select.by(da, "state", include="SCAN") da_off = select.by(da, "state", exclude="SCAN") @@ -411,9 +383,9 @@ def zscan( dems: Path, /, *, - include_mkid_ids: Optional[Sequence[int]] = None, - exclude_mkid_ids: Optional[Sequence[int]] = BAD_MKID_IDS, - data_type: Literal["df/f", "brightness"] = "brightness", + include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, + exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, + data_type: Literal["df/f", "brightness"] = DEFAULT_DATA_TYPE, 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", outdir: Path = Path(), @@ -443,18 +415,11 @@ def zscan( result = Path(outdir) / dems.with_suffix(f".zscan.{format}").name # load DEMS - da = load.dems(dems, chunks=None) - - if data_type == "df/f": - da.attrs.update(long_name="df/f", units="dimensionless") - - # select DEMS - da = select.by(da, "d2_mkid_type", include="filter") - da = select.by( - da, - "d2_mkid_id", - include=include_mkid_ids, - exclude=exclude_mkid_ids, + da = load_dems( + dems, + include_mkid_ids=include_mkid_ids, + exclude_mkid_ids=exclude_mkid_ids, + data_type=data_type, ) da_on = select.by(da, "state", include="ON") da_off = select.by(da, "state", exclude="ON") @@ -547,6 +512,61 @@ def subtract_per_scan(dems: xr.DataArray) -> xr.DataArray: raise ValueError("State must be either ON or OFF.") +def load_dems( + dems: Path, + /, + *, + include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, + exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, + data_type: Literal["brightness", "df/f"] = DEFAULT_DATA_TYPE, + frequency_units: str = DEFAULT_FREQUENCY_UNITS, + skycoord_units: str = DEFAULT_SKYCOORD_UNITS, +) -> xr.DataArray: + """Load a DEMS with given conversions and selections. + + Args: + 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 bad MKID IDs found on 2023-11-07. + data_type: Data type of the input DEMS file. + frequency_units: Units of the frequency axis. + skycoord_units: Units of the sky coordinate axes. + + Returns: + DEMS as a DataArray with given conversion and selections. + + """ + da = load.dems(dems, chunks=None) + da = assign.scan(da, by="state") + da = convert.frame(da, "relative") + da = select.by(da, "d2_mkid_type", "filter") + da = select.by( + da, + "d2_mkid_id", + include=include_mkid_ids, + exclude=exclude_mkid_ids, + ) + da = convert.coord_units( + da, + ["d2_mkid_frequency", "frequency"], + frequency_units, + ) + da = convert.coord_units( + da, + ["lat", "lat_origin", "lon", "lon_origin"], + skycoord_units, + ) + + if data_type == "brightness": + return da.assign_attrs(long_name="brightness", units="K") + + if data_type == "df/f": + return -da.assign_attrs(long_name="-df/f", units="dimensionless") + + raise ValueError("Data type must be either df/f or brightness.") + + def main() -> None: """Entry point of the decode-qlook command.""" with xr.set_options(keep_attrs=True): From 152def1cd33b0f1178e62a5a5631671c1e70e4d0 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sun, 19 Nov 2023 15:38:57 +0000 Subject: [PATCH 06/31] #135 Use long_name in DEMS as data_type if it exists --- decode/qlook.py | 30 ++++++++++++++++++++---------- 1 file changed, 20 insertions(+), 10 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 9cf22de..00f001c 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -15,7 +15,7 @@ # constants -DEFAULT_DATA_TYPE = "brightness" +DEFAULT_DATA_TYPE = None # fmt: off DEFAULT_EXCL_MKID_IDS = ( 18, 77, 117, 118, 140, 141, 161, 182, 183, 184, @@ -37,7 +37,7 @@ def still( *, include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, - data_type: Literal["df/f", "brightness"] = DEFAULT_DATA_TYPE, + data_type: Literal["df/f", "brightness", None] = DEFAULT_DATA_TYPE, 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", cabin_temperature: float = 273.0, @@ -53,6 +53,7 @@ def still( exclude_mkid_ids: MKID IDs to be excluded in analysis. Defaults to bad MKID IDs found on 2023-11-07. data_type: Data type of the input DEMS file. + Defaults to the ``long_name`` attribute in it. chan_weight: Weighting method along the channel axis. uniform: Uniform weight (i.e. no channel dependence). std: Inverse square of temporal standard deviation of sky. @@ -124,7 +125,7 @@ def pswsc( *, include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, - data_type: Literal["df/f", "brightness"] = DEFAULT_DATA_TYPE, + data_type: Literal["df/f", "brightness", None] = DEFAULT_DATA_TYPE, frequency_units: str = "GHz", outdir: Path = Path(), format: str = "png", @@ -138,6 +139,7 @@ def pswsc( exclude_mkid_ids: MKID IDs to be excluded in analysis. Defaults to bad MKID IDs found on 2023-11-07. data_type: Data type of the input DEMS file. + Defaults to the ``long_name`` attribute in it. frequency_units: Units of the frequency axis. outdir: Output directory for the analysis result. format: Output data format of the analysis result. @@ -203,7 +205,7 @@ def raster( *, include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, - data_type: Literal["df/f", "brightness"] = DEFAULT_DATA_TYPE, + data_type: Literal["df/f", "brightness", None] = DEFAULT_DATA_TYPE, 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", skycoord_grid: str = DEFAULT_SKYCOORD_GRID, @@ -220,6 +222,7 @@ def raster( exclude_mkid_ids: MKID IDs to be excluded in analysis. Defaults to bad MKID IDs found on 2023-11-07. data_type: Data type of the input DEMS file. + Defaults to the ``long_name`` attribute in it. chan_weight: Weighting method along the channel axis. uniform: Uniform weight (i.e. no channel dependence). std: Inverse square of temporal standard deviation of sky. @@ -309,7 +312,7 @@ def skydip( *, include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, - data_type: Literal["df/f", "brightness"] = DEFAULT_DATA_TYPE, + data_type: Literal["df/f", "brightness", None] = DEFAULT_DATA_TYPE, 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", outdir: Path = Path(), @@ -324,6 +327,7 @@ def skydip( exclude_mkid_ids: MKID IDs to be excluded in analysis. Defaults to bad MKID IDs found on 2023-11-07. data_type: Data type of the input DEMS file. + Defaults to the ``long_name`` attribute in it. chan_weight: Weighting method along the channel axis. uniform: Uniform weight (i.e. no channel dependence). std: Inverse square of temporal standard deviation of sky. @@ -385,7 +389,7 @@ def zscan( *, include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, - data_type: Literal["df/f", "brightness"] = DEFAULT_DATA_TYPE, + data_type: Literal["df/f", "brightness", None] = DEFAULT_DATA_TYPE, 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", outdir: Path = Path(), @@ -400,6 +404,7 @@ def zscan( exclude_mkid_ids: MKID IDs to be excluded in analysis. Defaults to bad MKID IDs found on 2023-11-07. data_type: Data type of the input DEMS file. + Defaults to the ``long_name`` attribute in it. chan_weight: Weighting method along the channel axis. uniform: Uniform weight (i.e. no channel dependence). std: Inverse square of temporal standard deviation of sky. @@ -518,19 +523,21 @@ def load_dems( *, include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, - data_type: Literal["brightness", "df/f"] = DEFAULT_DATA_TYPE, + data_type: Literal["brightness", "df/f", None] = DEFAULT_DATA_TYPE, frequency_units: str = DEFAULT_FREQUENCY_UNITS, skycoord_units: str = DEFAULT_SKYCOORD_UNITS, ) -> xr.DataArray: """Load a DEMS with given conversions and selections. 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 bad MKID IDs found on 2023-11-07. data_type: Data type of the input DEMS file. - frequency_units: Units of the frequency axis. + Defaults to the ``long_name`` attribute in it. + frequency_units: Units of the frequency. skycoord_units: Units of the sky coordinate axes. Returns: @@ -558,13 +565,16 @@ def load_dems( skycoord_units, ) + if data_type is None and "units" in da.attrs: + return da + if data_type == "brightness": return da.assign_attrs(long_name="brightness", units="K") if data_type == "df/f": - return -da.assign_attrs(long_name="-df/f", units="dimensionless") + return da.assign_attrs(long_name="df/f", units="dimensionless") - raise ValueError("Data type must be either df/f or brightness.") + raise ValueError("Data type could not be inferred.") def main() -> None: From e0bbb642017a509e52d4a05d93bd35488232d0d6 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sun, 19 Nov 2023 15:44:55 +0000 Subject: [PATCH 07/31] #135 Update channel weight method --- decode/qlook.py | 30 ++++++++++++------------------ 1 file changed, 12 insertions(+), 18 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 00f001c..70a7698 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -80,8 +80,8 @@ def still( da_off = select.by(da, "state", exclude=["ON", "SCAN"]) # make continuum series - weight = get_weight(da_off, method=chan_weight, pwv=pwv) - series = (da * weight).sum("chan") / weight.sum("chan") + weight = calc_chan_weight(da_off, method=chan_weight, pwv=pwv) + series = da.weighted(weight).mean("chan") # export output if format == "csv": @@ -263,8 +263,8 @@ def raster( da_sub = da_on - da_base.values # make continuum series - weight = get_weight(da_off, method=chan_weight, pwv=pwv) - series = (da_sub * weight).sum("chan") / weight.sum("chan") + weight = calc_chan_weight(da_off, method=chan_weight, pwv=pwv) + series = da_sub.weighted(weight).mean("chan") # make continuum map cube = make.cube( @@ -361,8 +361,8 @@ def skydip( da_off = select.by(da, "state", exclude="SCAN") # plotting - weight = get_weight(da_off, method=chan_weight, pwv=pwv) - series = (da_on * weight).sum("chan") / weight.sum("chan") + weight = calc_chan_weight(da_off, method=chan_weight, pwv=pwv) + series = da_on.weighted(weight).mean("chan") fig, axes = plt.subplots(1, 2, figsize=(12, 4)) @@ -430,8 +430,8 @@ def zscan( da_off = select.by(da, "state", exclude="ON") # plotting - weight = get_weight(da_off, method=chan_weight, pwv=pwv) - series = (da_on * weight).sum("chan") / weight.sum("chan") + weight = calc_chan_weight(da_off, method=chan_weight, pwv=pwv) + series = da_on.weighted(weight).mean("chan") fig, axes = plt.subplots(1, 2, figsize=(12, 4)) @@ -450,7 +450,7 @@ def zscan( print(str(result)) -def get_weight( +def calc_chan_weight( dems: xr.DataArray, /, *, @@ -480,15 +480,9 @@ def get_weight( return dems.std("time") ** -2 if method == "std/tx": - tx = ( - load.atm(type="eta") - .sel(pwv=float(pwv)) - .interp( - freq=dems.d2_mkid_frequency, - method="linear", - ) - ) - return (dems.std("time") / tx) ** -2 + tx = load.atm(type="eta").sel(pwv=float(pwv)) + freq = convert.units(dems.d2_mkid_frequency, tx.freq.attrs["units"]) + return (dems.std("time") / tx.interp(freq=freq)) ** -2 raise ValueError("Method must be either uniform, std, or std/tx.") From 835f918586e75af923c04660aa4a26293c017d3d Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Thu, 23 Nov 2023 10:14:51 +0000 Subject: [PATCH 08/31] #135 Add common function for saving result --- decode/qlook.py | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/decode/qlook.py b/decode/qlook.py index 70a7698..8551c56 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -3,7 +3,7 @@ # standard library from pathlib import Path -from typing import Literal, Optional, Sequence, cast +from typing import Literal, Optional, Sequence, Union, cast # dependencies @@ -11,6 +11,7 @@ import xarray as xr import matplotlib.pyplot as plt from fire import Fire +from matplotlib.figure import Figure from . import assign, convert, load, make, plot, select, utils @@ -571,6 +572,30 @@ def load_dems( raise ValueError("Data type could not be inferred.") +def save_qlook(qlook: Union[Figure, xr.DataArray], filename: Path) -> Path: + """Save a quick look result to a file with given format. + + Args: + qlook: Matplotlib figure or DataArray to be saved. + filename: Path of the saved file. + + Returns: + Absolute path of the saved file. + + """ + if isinstance(qlook, Figure): + qlook.savefig(filename) + elif (ext := "".join(filename.suffixes)) == ".csv": + name = qlook.attrs["data_type"] + qlook.to_dataset(name=name).to_pandas().to_csv(filename) + elif ext == ".nc": + qlook.to_netcdf(filename) + elif ext == ".zarr" or format == ".zarr.zip": + qlook.to_zarr(filename, mode="w") + + return Path(filename).expanduser().resolve() + + def main() -> None: """Entry point of the decode-qlook command.""" with xr.set_options(keep_attrs=True): From 81feaee171aade1b385983b03daa701501a24877 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Thu, 23 Nov 2023 10:19:19 +0000 Subject: [PATCH 09/31] #135 Add common function for robust limit --- decode/qlook.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/decode/qlook.py b/decode/qlook.py index 8551c56..3372afa 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -28,6 +28,7 @@ DEFAULT_INCL_MKID_IDS = None DEFAULT_SKYCOORD_GRID = "6 arcsec" DEFAULT_SKYCOORD_UNITS = "arcsec" +SIGMA_OVER_MAD = 1.4826 DFOF_TO_TSKY = (300 - 77) / 3e-5 TSKY_TO_DFOF = 3e-5 / (300 - 77) @@ -512,6 +513,16 @@ def subtract_per_scan(dems: xr.DataArray) -> xr.DataArray: raise ValueError("State must be either ON or OFF.") +def get_robust_lim(da: xr.DataArray) -> tuple[float, float]: + """Calculate a robust limit for plotting.""" + sigma = SIGMA_OVER_MAD * utils.mad(da) + + return ( + float(np.percentile(da.data, 1) - sigma), + float(np.percentile(da.data, 99) + sigma), + ) + + def load_dems( dems: Path, /, From c79e5170c629b819c2aa572bf0de89fb819a771a Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Thu, 23 Nov 2023 10:20:04 +0000 Subject: [PATCH 10/31] #135 Assign sec(Z) when loading DEMS --- decode/qlook.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/decode/qlook.py b/decode/qlook.py index 3372afa..c46c901 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -551,6 +551,18 @@ def load_dems( """ da = load.dems(dems, chunks=None) + + if da.frame == "altaz": + z = np.pi / 2 - convert.units(da.lat, "rad") + secz = cast(xr.DataArray, 1 / np.cos(z)) + + da = da.assign_coords( + secz=secz.assign_attrs( + long_name="sec(Z)", + units="dimensionless", + ) + ) + da = assign.scan(da, by="state") da = convert.frame(da, "relative") da = select.by(da, "d2_mkid_type", "filter") From 9b229cd1fa54f666998a8f87e93cdefacaca16dd Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Thu, 23 Nov 2023 10:21:59 +0000 Subject: [PATCH 11/31] #135 Update still command using common functions --- decode/qlook.py | 67 +++++++++++++++++++------------------------------ 1 file changed, 26 insertions(+), 41 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index c46c901..55503ed 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -16,6 +16,7 @@ # constants +DATA_FORMATS = "csv", "nc", "zarr", "zarr.zip" DEFAULT_DATA_TYPE = None # fmt: off DEFAULT_EXCL_MKID_IDS = ( @@ -24,6 +25,8 @@ 283, 296, 297, 299, 301, 313, ) # fmt: on +DEFAULT_FIGSIZE = 12, 4 +DEFAULT_FORMAT = "png" DEFAULT_FREQUENCY_UNITS = "GHz" DEFAULT_INCL_MKID_IDS = None DEFAULT_SKYCOORD_GRID = "6 arcsec" @@ -42,10 +45,9 @@ def still( data_type: Literal["df/f", "brightness", None] = DEFAULT_DATA_TYPE, 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", - cabin_temperature: float = 273.0, + format: str = DEFAULT_FORMAT, outdir: Path = Path(), - format: str = "png", -) -> None: +) -> Path: """Quick-look at a still observation. Args: @@ -63,62 +65,45 @@ def still( transmission calculated by the ATM model. pwv: PWV in units of mm. Only used for the calculation of the atmospheric transmission when chan_weight is std/tx. - cabin_temperature: Temperature at the ASTE cabin. - Only used for the df/f-to-Tsky conversion. - outdir: Output directory for the analysis result. - format: Output data format of the analysis result. + format: Output data format of the quick-look result. + outdir: Output directory for the quick-look result. - """ - dems = Path(dems) - out = Path(outdir) / dems.with_suffix(f".still.{format}").name + Returns: + Absolute path of the saved file. - # load DEMS + """ da = load_dems( dems, include_mkid_ids=include_mkid_ids, exclude_mkid_ids=exclude_mkid_ids, data_type=data_type, ) - da_off = select.by(da, "state", exclude=["ON", "SCAN"]) # make continuum series + da_off = select.by(da, "state", exclude=["ON", "SCAN"]) weight = calc_chan_weight(da_off, method=chan_weight, pwv=pwv) series = da.weighted(weight).mean("chan") - # export output - if format == "csv": - series.to_dataset(name=data_type).to_pandas().to_csv(out) - elif format == "nc": - series.to_netcdf(out) - elif format.startswith("zarr"): - series.to_zarr(out) - else: - fig, axes = plt.subplots(1, 2, figsize=(12, 4)) + # save result + filename = Path(dems).with_suffix(f".still.{format}").name - ax = axes[0] - plot.state(da, add_colorbar=False, add_legend=False, ax=ax) - ax.set_title(Path(dems).name) - ax.grid(True) + if format in DATA_FORMATS: + return save_qlook(series, Path(outdir) / filename) - ax = axes[1] - plot.data(series, add_colorbar=False, ax=ax) - ax.set_title(Path(dems).name) - ax.grid(True) + fig, axes = plt.subplots(1, 2, figsize=DEFAULT_FIGSIZE) - if data_type == "df/f": - ax = ax.secondary_yaxis( - "right", - functions=( - lambda x: DFOF_TO_TSKY * x + cabin_temperature, - lambda x: TSKY_TO_DFOF * (x - cabin_temperature), - ), - ) - ax.set_ylabel("Approx. brightness [K]") + ax = axes[0] + plot.state(da, add_colorbar=False, add_legend=False, ax=ax) - fig.tight_layout() - fig.savefig(out) + ax = axes[1] + plot.data(series, add_colorbar=False, ax=ax) - print(str(out)) + for ax in axes: + ax.set_title(Path(dems).name) + ax.grid(True) + + fig.tight_layout() + return save_qlook(fig, Path(outdir) / filename) def pswsc( From a5abd2badb105b7028a0ea6078429507bad7f8a9 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Thu, 23 Nov 2023 10:22:39 +0000 Subject: [PATCH 12/31] #135 Update pswsc command using common functions --- decode/qlook.py | 71 +++++++++++++++++++------------------------------ 1 file changed, 28 insertions(+), 43 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 55503ed..47454bb 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -113,10 +113,10 @@ def pswsc( include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, data_type: Literal["df/f", "brightness", None] = DEFAULT_DATA_TYPE, - frequency_units: str = "GHz", + frequency_units: str = DEFAULT_FREQUENCY_UNITS, + format: str = DEFAULT_FORMAT, outdir: Path = Path(), - format: str = "png", -) -> None: +) -> Path: """Quick-look at a PSW observation with sky chopper. Args: @@ -128,14 +128,13 @@ def pswsc( data_type: Data type of the input DEMS file. Defaults to the ``long_name`` attribute in it. frequency_units: Units of the frequency axis. - outdir: Output directory for the analysis result. - format: Output data format of the analysis result. + format: Output data format of the quick-look result. + outdir: Output directory for the quick-look result. - """ - dems = Path(dems) - out = Path(outdir) / dems.with_suffix(f".pswsc.{format}").name + Returns: + Absolute path of the saved file. - # load DEMS + """ da = load_dems( dems, include_mkid_ids=include_mkid_ids, @@ -143,47 +142,33 @@ def pswsc( data_type=data_type, frequency_units=frequency_units, ) + + # make spectrum da_scan = select.by(da, "state", ["ON", "OFF"]) da_sub = da_scan.groupby("scan").map(subtract_per_scan) - - # export output spec = da_sub.mean("scan") - mad = utils.mad(spec) - - if format == "csv": - spec.to_dataset(name=data_type).to_pandas().to_csv(out) - elif format == "nc": - spec.to_netcdf(out) - elif format.startswith("zarr"): - spec.to_zarr(out) - else: - fig, axes = plt.subplots(1, 2, figsize=(12, 4)) - - ax = axes[0] - plot.data(da.scan, ax=ax) - ax.set_title(Path(dems).name) - ax.grid(True) - ax = axes[1] - plot.data(spec, x="frequency", s=5, hue=None, ax=ax) - ax.set_ylim(-mad, spec.max() + mad) - ax.set_title(Path(dems).name) - ax.grid(True) + # save result + filename = Path(dems).with_suffix(f".pswsc.{format}").name - if data_type == "df/f": - ax = ax.secondary_yaxis( - "right", - functions=( - lambda x: DFOF_TO_TSKY * x, - lambda x: TSKY_TO_DFOF * x, - ), - ) - ax.set_ylabel("Approx. brightness [K]") + if format in DATA_FORMATS: + return save_qlook(spec, Path(outdir) / filename) + + fig, axes = plt.subplots(1, 2, figsize=DEFAULT_FIGSIZE) + + ax = axes[0] + plot.data(da.scan, ax=ax) + + ax = axes[1] + plot.data(spec, x="frequency", s=5, hue=None, ax=ax) + ax.set_ylim(get_robust_lim(spec)) - fig.tight_layout() - fig.savefig(out) + for ax in axes: + ax.set_title(Path(dems).name) + ax.grid(True) - print(str(out)) + fig.tight_layout() + return save_qlook(fig, Path(outdir) / filename) def raster( From ba21e04d49ecd6eab008219b791319004d280ad3 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Thu, 23 Nov 2023 10:23:05 +0000 Subject: [PATCH 13/31] #135 Update raster command using common functions --- decode/qlook.py | 49 ++++++++++++++++++++++++------------------------- 1 file changed, 24 insertions(+), 25 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 47454bb..8803f5b 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -182,9 +182,9 @@ def raster( pwv: Literal["0.5", "1.0", "2.0", "3.0", "4.0", "5.0"] = "5.0", skycoord_grid: str = DEFAULT_SKYCOORD_GRID, skycoord_units: str = DEFAULT_SKYCOORD_UNITS, + format: str = DEFAULT_FORMAT, outdir: Path = Path(), - format: str = "png", -) -> None: +) -> Path: """Quick-look at a raster scan observation. Args: @@ -204,14 +204,13 @@ def raster( the atmospheric transmission when chan_weight is std/tx. skycoord_grid: Grid size of the sky coordinate axes. skycoord_units: Units of the sky coordinate axes. - outdir: Output directory for analysis results. format: Output image format of analysis results. + outdir: Output directory for analysis results. - """ - dems = Path(dems) - result = Path(outdir) / dems.with_suffix(f".raster.{format}").name + Returns: + Absolute path of the saved file. - # load DEMS + """ da = load_dems( dems, include_mkid_ids=include_mkid_ids, @@ -219,10 +218,10 @@ def raster( data_type=data_type, skycoord_units=skycoord_units, ) - da_on = select.by(da, "state", include="SCAN") - da_off = select.by(da, "state", exclude="SCAN") # subtract temporal baseline + da_on = select.by(da, "state", include="SCAN") + da_off = select.by(da, "state", exclude="SCAN") da_base = ( da_off.groupby("scan") .map(mean_in_time) @@ -244,38 +243,38 @@ def raster( skycoord_grid=skycoord_grid, skycoord_units=skycoord_units, ) - cont = (cube * weight).sum("chan") / weight.sum("chan") + cont = cube.weighted(weight).mean("chan") - if data_type == "df/f": - cont.attrs.update(long_name="df/f", units="dimensionless") + # save result + filename = Path(dems).with_suffix(f".pswsc.{format}").name - # plotting - map_lim = max(abs(cube.lon).max(), abs(cube.lat).max()) - max_pix = cont.where(cont == cont.max(), drop=True) - max_lon = float(max_pix.lon) - max_lat = float(max_pix.lat) + if format in DATA_FORMATS: + return save_qlook(cont, Path(outdir) / filename) fig, axes = plt.subplots(1, 2, figsize=(12, 5.5)) ax = axes[0] plot.data(series, ax=ax) ax.set_title(Path(dems).name) - ax.grid(True) ax = axes[1] + map_lim = max(abs(cube.lon).max(), abs(cube.lat).max()) + max_pix = cont.where(cont == cont.max(), drop=True) + cont.plot(ax=ax) # type: ignore + ax.set_xlim(-map_lim, map_lim) + ax.set_ylim(-map_lim, map_lim) ax.set_title( "Maxima: " - f"dAz = {max_lon:+.1f} {cont.lon.attrs['units']}, " - f"dEl = {max_lat:+.1f} {cont.lat.attrs['units']}" + f"dAz = {max_pix.lon:+.1f} {cont.lon.attrs['units']}, " + f"dEl = {max_pix.lat:+.1f} {cont.lat.attrs['units']}" ) - ax.set_xlim(-map_lim, map_lim) - ax.set_ylim(-map_lim, map_lim) - ax.grid() + + for ax in axes: + ax.grid(True) fig.tight_layout() - fig.savefig(result) - print(str(result)) + return save_qlook(fig, Path(outdir) / filename) def skydip( From 68773618f2b7761a6c6d645229140f7bcc4f2e31 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Thu, 23 Nov 2023 10:23:26 +0000 Subject: [PATCH 14/31] #135 Update skydip command using common functions --- decode/qlook.py | 44 ++++++++++++++++++++------------------------ 1 file changed, 20 insertions(+), 24 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 8803f5b..2095ef3 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -286,9 +286,9 @@ def skydip( data_type: Literal["df/f", "brightness", None] = DEFAULT_DATA_TYPE, 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", + format: str = DEFAULT_FORMAT, outdir: Path = Path(), - format: str = "png", -) -> None: +) -> Path: """Quick-look at a skydip observation. Args: @@ -306,52 +306,48 @@ def skydip( transmission calculated by the ATM model. pwv: PWV in units of mm. Only used for the calculation of the atmospheric transmission when chan_weight is std/tx. - outdir: Output directory for analysis results. format: Output image format of analysis results. + outdir: Output directory for analysis results. - """ - dems = Path(dems) - result = Path(outdir) / dems.with_suffix(f".skydip.{format}").name + Returns: + Absolute path of the saved file. - # load DEMS + """ da = load_dems( dems, include_mkid_ids=include_mkid_ids, exclude_mkid_ids=exclude_mkid_ids, data_type=data_type, ) - z = np.pi / 2 - convert.units(da.lat, "rad") - secz = cast(xr.DataArray, 1 / np.cos(z)) - da = da.assign_coords( - secz=secz.assign_attrs( - long_name="sec(Z)", - units="dimensionless", - ) - ) + + # make continuum series da_on = select.by(da, "state", include="SCAN") da_off = select.by(da, "state", exclude="SCAN") - - # plotting weight = calc_chan_weight(da_off, method=chan_weight, pwv=pwv) series = da_on.weighted(weight).mean("chan") - fig, axes = plt.subplots(1, 2, figsize=(12, 4)) + # save result + filename = Path(dems).with_suffix(f".skydip.{format}").name + + if format in DATA_FORMATS: + return save_qlook(series, Path(outdir) / filename) + + fig, axes = plt.subplots(1, 2, figsize=DEFAULT_FIGSIZE) ax = axes[0] plot.data(series, hue="secz", ax=ax) - ax.set_title(Path(dems).name) - ax.grid(True) ax = axes[1] plot.data(series, x="secz", ax=ax) - ax.set_title(Path(dems).name) ax.set_xscale("log") ax.set_yscale("log") - ax.grid(True) + + for ax in axes: + ax.set_title(Path(dems).name) + ax.grid(True) fig.tight_layout() - fig.savefig(result) - print(str(result)) + return save_qlook(fig, Path(outdir) / filename) def zscan( From db53cf20b64b3aced2c84ff9f4a5e2f4d9418a73 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Thu, 23 Nov 2023 10:23:47 +0000 Subject: [PATCH 15/31] #135 Update zscan command using common functions --- decode/qlook.py | 36 ++++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 2095ef3..f1f50ea 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -359,9 +359,9 @@ def zscan( data_type: Literal["df/f", "brightness", None] = DEFAULT_DATA_TYPE, 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", + format: str = DEFAULT_FORMAT, outdir: Path = Path(), - format: str = "png", -) -> None: +) -> Path: """Quick-look at an observation of subref axial focus scan. Args: @@ -379,42 +379,46 @@ def zscan( transmission calculated by the ATM model. pwv: PWV in units of mm. Only used for the calculation of the atmospheric transmission when chan_weight is std/tx. - outdir: Output directory for analysis results. format: Output image format of analysis results. + outdir: Output directory for analysis results. - """ - dems = Path(dems) - result = Path(outdir) / dems.with_suffix(f".zscan.{format}").name + Returns: + Absolute path of the saved file. - # load DEMS + """ da = load_dems( dems, include_mkid_ids=include_mkid_ids, exclude_mkid_ids=exclude_mkid_ids, data_type=data_type, ) + + # make continuum series da_on = select.by(da, "state", include="ON") da_off = select.by(da, "state", exclude="ON") - - # plotting weight = calc_chan_weight(da_off, method=chan_weight, pwv=pwv) series = da_on.weighted(weight).mean("chan") - fig, axes = plt.subplots(1, 2, figsize=(12, 4)) + # save output + filename = Path(dems).with_suffix(f".zscan.{format}").name + + if format in DATA_FORMATS: + return save_qlook(series, Path(outdir) / filename) + + fig, axes = plt.subplots(1, 2, figsize=DEFAULT_FIGSIZE) ax = axes[0] plot.data(series, hue="aste_subref_z", ax=ax) - ax.set_title(Path(dems).name) - ax.grid(True) ax = axes[1] plot.data(series, x="aste_subref_z", ax=ax) - ax.set_title(Path(dems).name) - ax.grid(True) + + for ax in axes: + ax.set_title(Path(dems).name) + ax.grid(True) fig.tight_layout() - fig.savefig(result) - print(str(result)) + return save_qlook(fig, Path(outdir) / filename) def calc_chan_weight( From b4d34e11dd6e85f2c323662d6ebd735a8c0aded8 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Thu, 23 Nov 2023 10:25:47 +0000 Subject: [PATCH 16/31] #135 Reorder commands alphabetically --- decode/qlook.py | 144 ++++++++++++++++++++++++------------------------ 1 file changed, 72 insertions(+), 72 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index f1f50ea..55be607 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -1,4 +1,4 @@ -__all__ = ["still", "pswsc", "raster", "skydip", "zscan"] +__all__ = ["pswsc", "raster", "skydip", "still", "zscan"] # standard library @@ -36,76 +36,6 @@ TSKY_TO_DFOF = 3e-5 / (300 - 77) -def still( - dems: Path, - /, - *, - include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, - exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, - data_type: Literal["df/f", "brightness", None] = DEFAULT_DATA_TYPE, - 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", - format: str = DEFAULT_FORMAT, - outdir: Path = Path(), -) -> Path: - """Quick-look at a still observation. - - 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 bad MKID IDs found on 2023-11-07. - data_type: Data type of the input DEMS file. - Defaults to the ``long_name`` attribute in it. - chan_weight: Weighting method along the channel axis. - uniform: Uniform weight (i.e. no channel dependence). - std: Inverse square of temporal standard deviation of sky. - std/tx: Same as std but std is divided by the atmospheric - transmission calculated by the ATM model. - pwv: PWV in units of mm. Only used for the calculation of - the atmospheric transmission when chan_weight is std/tx. - format: Output data format of the quick-look result. - outdir: Output directory for the quick-look result. - - Returns: - Absolute path of the saved file. - - """ - da = load_dems( - dems, - include_mkid_ids=include_mkid_ids, - exclude_mkid_ids=exclude_mkid_ids, - data_type=data_type, - ) - - # make continuum series - da_off = select.by(da, "state", exclude=["ON", "SCAN"]) - weight = calc_chan_weight(da_off, method=chan_weight, pwv=pwv) - series = da.weighted(weight).mean("chan") - - # save result - filename = Path(dems).with_suffix(f".still.{format}").name - - if format in DATA_FORMATS: - return save_qlook(series, Path(outdir) / filename) - - fig, axes = plt.subplots(1, 2, figsize=DEFAULT_FIGSIZE) - - ax = axes[0] - plot.state(da, add_colorbar=False, add_legend=False, ax=ax) - - ax = axes[1] - plot.data(series, add_colorbar=False, ax=ax) - - for ax in axes: - ax.set_title(Path(dems).name) - ax.grid(True) - - fig.tight_layout() - return save_qlook(fig, Path(outdir) / filename) - - def pswsc( dems: Path, /, @@ -350,6 +280,76 @@ def skydip( return save_qlook(fig, Path(outdir) / filename) +def still( + dems: Path, + /, + *, + include_mkid_ids: Optional[Sequence[int]] = DEFAULT_INCL_MKID_IDS, + exclude_mkid_ids: Optional[Sequence[int]] = DEFAULT_EXCL_MKID_IDS, + data_type: Literal["df/f", "brightness", None] = DEFAULT_DATA_TYPE, + 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", + format: str = DEFAULT_FORMAT, + outdir: Path = Path(), +) -> Path: + """Quick-look at a still observation. + + 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 bad MKID IDs found on 2023-11-07. + data_type: Data type of the input DEMS file. + Defaults to the ``long_name`` attribute in it. + chan_weight: Weighting method along the channel axis. + uniform: Uniform weight (i.e. no channel dependence). + std: Inverse square of temporal standard deviation of sky. + std/tx: Same as std but std is divided by the atmospheric + transmission calculated by the ATM model. + pwv: PWV in units of mm. Only used for the calculation of + the atmospheric transmission when chan_weight is std/tx. + format: Output data format of the quick-look result. + outdir: Output directory for the quick-look result. + + Returns: + Absolute path of the saved file. + + """ + da = load_dems( + dems, + include_mkid_ids=include_mkid_ids, + exclude_mkid_ids=exclude_mkid_ids, + data_type=data_type, + ) + + # make continuum series + da_off = select.by(da, "state", exclude=["ON", "SCAN"]) + weight = calc_chan_weight(da_off, method=chan_weight, pwv=pwv) + series = da.weighted(weight).mean("chan") + + # save result + filename = Path(dems).with_suffix(f".still.{format}").name + + if format in DATA_FORMATS: + return save_qlook(series, Path(outdir) / filename) + + fig, axes = plt.subplots(1, 2, figsize=DEFAULT_FIGSIZE) + + ax = axes[0] + plot.state(da, add_colorbar=False, add_legend=False, ax=ax) + + ax = axes[1] + plot.data(series, add_colorbar=False, ax=ax) + + for ax in axes: + ax.set_title(Path(dems).name) + ax.grid(True) + + fig.tight_layout() + return save_qlook(fig, Path(outdir) / filename) + + def zscan( dems: Path, /, @@ -594,10 +594,10 @@ def main() -> None: Fire( { "default": still, - "still": still, "pswsc": pswsc, "raster": raster, "skydip": skydip, + "still": still, "zscan": zscan, } ) From fdfd685dbf8b3a9391ea412dc90fcfb78621197d Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Thu, 23 Nov 2023 10:35:04 +0000 Subject: [PATCH 17/31] #135 Pin Poetry version (1.7.1) --- .devcontainer/devcontainer.json | 2 +- .github/workflows/pypi.yaml | 2 +- .github/workflows/tests.yaml | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json index 95dee1c..fcbc0f4 100644 --- a/.devcontainer/devcontainer.json +++ b/.devcontainer/devcontainer.json @@ -1,7 +1,7 @@ { "name": "decode", "image":"python:3.11", - "onCreateCommand": "pip install poetry==1.6.1", + "onCreateCommand": "pip install poetry==1.7.1", "postCreateCommand": "poetry install", "containerEnv": { "POETRY_VIRTUALENVS_CREATE": "false" diff --git a/.github/workflows/pypi.yaml b/.github/workflows/pypi.yaml index 9d90eae..ce39d69 100644 --- a/.github/workflows/pypi.yaml +++ b/.github/workflows/pypi.yaml @@ -16,4 +16,4 @@ jobs: - uses: actions/setup-python@v4 with: python-version: "3.11" - - run: pip install poetry && poetry publish --build + - run: pip install poetry==1.7.1 && poetry publish --build diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml index 5f00058..3a1ffd7 100644 --- a/.github/workflows/tests.yaml +++ b/.github/workflows/tests.yaml @@ -27,7 +27,7 @@ jobs: - uses: actions/setup-python@v4 with: python-version: ${{ matrix.python }} - - run: pip install poetry && poetry install + - run: pip install poetry==1.7.1 && poetry install - run: black --check decode docs tests - run: pyright decode docs tests - run: pytest -v tests From 3218e4bbfb84940b7f7c7adf2de6def7cc97b487 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Fri, 24 Nov 2023 10:09:16 +0000 Subject: [PATCH 18/31] #135 Fix common function for robust limit --- decode/qlook.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 55be607..6752ac4 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -487,8 +487,8 @@ def get_robust_lim(da: xr.DataArray) -> tuple[float, float]: sigma = SIGMA_OVER_MAD * utils.mad(da) return ( - float(np.percentile(da.data, 1) - sigma), - float(np.percentile(da.data, 99) + sigma), + float(np.nanpercentile(da.data, 1) - sigma), + float(np.nanpercentile(da.data, 99) + sigma), ) From c68da6e6c5803937fbdada949bb4c9a9510b9767 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Fri, 24 Nov 2023 10:10:01 +0000 Subject: [PATCH 19/31] #135 Fix channel weight method --- decode/qlook.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 6752ac4..c83b9ea 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -165,7 +165,7 @@ def raster( # make continuum series weight = calc_chan_weight(da_off, method=chan_weight, pwv=pwv) - series = da_sub.weighted(weight).mean("chan") + series = da_sub.weighted(weight.fillna(0)).mean("chan") # make continuum map cube = make.cube( @@ -173,7 +173,7 @@ def raster( skycoord_grid=skycoord_grid, skycoord_units=skycoord_units, ) - cont = cube.weighted(weight).mean("chan") + cont = cube.weighted(weight.fillna(0)).mean("chan") # save result filename = Path(dems).with_suffix(f".pswsc.{format}").name @@ -254,7 +254,7 @@ def skydip( da_on = select.by(da, "state", include="SCAN") da_off = select.by(da, "state", exclude="SCAN") weight = calc_chan_weight(da_off, method=chan_weight, pwv=pwv) - series = da_on.weighted(weight).mean("chan") + series = da_on.weighted(weight.fillna(0)).mean("chan") # save result filename = Path(dems).with_suffix(f".skydip.{format}").name @@ -326,7 +326,7 @@ def still( # make continuum series da_off = select.by(da, "state", exclude=["ON", "SCAN"]) weight = calc_chan_weight(da_off, method=chan_weight, pwv=pwv) - series = da.weighted(weight).mean("chan") + series = da.weighted(weight.fillna(0)).mean("chan") # save result filename = Path(dems).with_suffix(f".still.{format}").name @@ -397,7 +397,7 @@ def zscan( da_on = select.by(da, "state", include="ON") da_off = select.by(da, "state", exclude="ON") weight = calc_chan_weight(da_off, method=chan_weight, pwv=pwv) - series = da_on.weighted(weight).mean("chan") + series = da_on.weighted(weight.fillna(0)).mean("chan") # save output filename = Path(dems).with_suffix(f".zscan.{format}").name From 3b0d9442be8bc162899bf7b150022ca0d2d3f29c Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Fri, 24 Nov 2023 10:39:15 +0000 Subject: [PATCH 20/31] #135 Ignore warnings related to standard deviation --- decode/qlook.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index c83b9ea..aef206b 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -4,6 +4,7 @@ # standard library from pathlib import Path from typing import Literal, Optional, Sequence, Union, cast +from warnings import catch_warnings, simplefilter # dependencies @@ -448,12 +449,17 @@ def calc_chan_weight( return xr.ones_like(dems.mean("time")) if method == "std": - return dems.std("time") ** -2 + with catch_warnings(): + simplefilter("ignore") + return dems.std("time") ** -2 if method == "std/tx": tx = load.atm(type="eta").sel(pwv=float(pwv)) freq = convert.units(dems.d2_mkid_frequency, tx.freq.attrs["units"]) - return (dems.std("time") / tx.interp(freq=freq)) ** -2 + + with catch_warnings(): + simplefilter("ignore") + return (dems.std("time") / tx.interp(freq=freq)) ** -2 raise ValueError("Method must be either uniform, std, or std/tx.") From 8cec3f8dd7d0c966f269e7711bd052127ce38a75 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Fri, 24 Nov 2023 10:40:40 +0000 Subject: [PATCH 21/31] #135 Fix figure title of raster command --- decode/qlook.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index aef206b..7a1263b 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -197,8 +197,8 @@ def raster( ax.set_ylim(-map_lim, map_lim) ax.set_title( "Maxima: " - f"dAz = {max_pix.lon:+.1f} {cont.lon.attrs['units']}, " - f"dEl = {max_pix.lat:+.1f} {cont.lat.attrs['units']}" + 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: From f8b9d41821ffb2fac117e370daa9d15bff870a10 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sat, 25 Nov 2023 00:16:41 +0000 Subject: [PATCH 22/31] #135 Fix assigned long_name in load_dems --- decode/qlook.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/decode/qlook.py b/decode/qlook.py index 7a1263b..cf66d2f 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -562,7 +562,7 @@ def load_dems( return da if data_type == "brightness": - return da.assign_attrs(long_name="brightness", units="K") + return da.assign_attrs(long_name="Brightness", units="K") if data_type == "df/f": return da.assign_attrs(long_name="df/f", units="dimensionless") From 35ce579159663af78a177f5bd26b8e000e508f2d Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sat, 25 Nov 2023 00:17:05 +0000 Subject: [PATCH 23/31] #135 Fix skydip command --- decode/qlook.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index cf66d2f..7a44ffc 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -270,8 +270,6 @@ def skydip( ax = axes[1] plot.data(series, x="secz", ax=ax) - ax.set_xscale("log") - ax.set_yscale("log") for ax in axes: ax.set_title(Path(dems).name) From 7d919ff00543e0399ddfc6a9e43a65c45fa79f05 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sat, 25 Nov 2023 00:17:31 +0000 Subject: [PATCH 24/31] #135 Update docstrings of qlook commands --- decode/qlook.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 7a44ffc..dacf93e 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -135,8 +135,8 @@ def raster( the atmospheric transmission when chan_weight is std/tx. skycoord_grid: Grid size of the sky coordinate axes. skycoord_units: Units of the sky coordinate axes. - format: Output image format of analysis results. - outdir: Output directory for analysis results. + format: Output image format of quick-look result. + outdir: Output directory for the quick-look result. Returns: Absolute path of the saved file. @@ -237,8 +237,8 @@ def skydip( transmission calculated by the ATM model. pwv: PWV in units of mm. Only used for the calculation of the atmospheric transmission when chan_weight is std/tx. - format: Output image format of analysis results. - outdir: Output directory for analysis results. + format: Output image format of quick-look result. + outdir: Output directory for the quick-look result. Returns: Absolute path of the saved file. @@ -378,8 +378,8 @@ def zscan( transmission calculated by the ATM model. pwv: PWV in units of mm. Only used for the calculation of the atmospheric transmission when chan_weight is std/tx. - format: Output image format of analysis results. - outdir: Output directory for analysis results. + format: Output image format of quick-look result. + outdir: Output directory for the quick-look result. Returns: Absolute path of the saved file. From 150fddc63d23b075807da43317109dbcc17a974c Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sat, 25 Nov 2023 00:21:05 +0000 Subject: [PATCH 25/31] #135 Update default MKID IDs to be excluded --- decode/qlook.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index dacf93e..3cf79da 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -21,9 +21,8 @@ DEFAULT_DATA_TYPE = None # fmt: off DEFAULT_EXCL_MKID_IDS = ( - 18, 77, 117, 118, 140, 141, 161, 182, 183, 184, - 201, 209, 211, 232, 233, 239, 258, 259, 278, 282, - 283, 296, 297, 299, 301, 313, + 0, 18, 26, 73, 130, 184, 118, 119, 201, 202, + 208, 214, 261, 266, 280, 283, 299, 304, 308, 321, ) # fmt: on DEFAULT_FIGSIZE = 12, 4 @@ -55,7 +54,7 @@ def pswsc( 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 bad MKID IDs found on 2023-11-07. + Defaults to bad MKID IDs found on 2023-11-19. data_type: Data type of the input DEMS file. Defaults to the ``long_name`` attribute in it. frequency_units: Units of the frequency axis. @@ -123,7 +122,7 @@ def raster( 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 bad MKID IDs found on 2023-11-07. + Defaults to bad MKID IDs found on 2023-11-19. data_type: Data type of the input DEMS file. Defaults to the ``long_name`` attribute in it. chan_weight: Weighting method along the channel axis. @@ -227,7 +226,7 @@ def skydip( 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 bad MKID IDs found on 2023-11-07. + Defaults to bad MKID IDs found on 2023-11-19. data_type: Data type of the input DEMS file. Defaults to the ``long_name`` attribute in it. chan_weight: Weighting method along the channel axis. @@ -298,7 +297,7 @@ def still( 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 bad MKID IDs found on 2023-11-07. + Defaults to bad MKID IDs found on 2023-11-19. data_type: Data type of the input DEMS file. Defaults to the ``long_name`` attribute in it. chan_weight: Weighting method along the channel axis. @@ -368,7 +367,7 @@ def zscan( 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 bad MKID IDs found on 2023-11-07. + Defaults to bad MKID IDs found on 2023-11-19. data_type: Data type of the input DEMS file. Defaults to the ``long_name`` attribute in it. chan_weight: Weighting method along the channel axis. @@ -513,7 +512,7 @@ def load_dems( 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 bad MKID IDs found on 2023-11-07. + Defaults to bad MKID IDs found on 2023-11-19. data_type: Data type of the input DEMS file. Defaults to the ``long_name`` attribute in it. frequency_units: Units of the frequency. From 39c2c41b7345f0062e3747a40fafdb5dcec81cb7 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sat, 25 Nov 2023 00:21:22 +0000 Subject: [PATCH 26/31] #135 Remove unused constants --- decode/qlook.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/decode/qlook.py b/decode/qlook.py index 3cf79da..3dfd2fe 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -32,8 +32,6 @@ DEFAULT_SKYCOORD_GRID = "6 arcsec" DEFAULT_SKYCOORD_UNITS = "arcsec" SIGMA_OVER_MAD = 1.4826 -DFOF_TO_TSKY = (300 - 77) / 3e-5 -TSKY_TO_DFOF = 3e-5 / (300 - 77) def pswsc( From 12feca014bdeff12f889847ef8357b437a4c9745 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sat, 25 Nov 2023 00:22:09 +0000 Subject: [PATCH 27/31] #135 Fix figure title of raster command --- decode/qlook.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/decode/qlook.py b/decode/qlook.py index 3dfd2fe..00bf60c 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -193,7 +193,7 @@ def raster( ax.set_xlim(-map_lim, map_lim) ax.set_ylim(-map_lim, map_lim) ax.set_title( - "Maxima: " + "Maximum: " f"dAz = {float(max_pix.lon):+.1f} {cont.lon.attrs['units']}, " f"dEl = {float(max_pix.lat):+.1f} {cont.lat.attrs['units']}" ) From 9558e0c9f9fe30d83e0e5650185f27352195bf6a Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sat, 25 Nov 2023 00:32:23 +0000 Subject: [PATCH 28/31] #135 Accept DataArray with units as new_units --- decode/convert.py | 11 ++++++++--- decode/qlook.py | 2 +- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/decode/convert.py b/decode/convert.py index 23b91b1..37f399e 100644 --- a/decode/convert.py +++ b/decode/convert.py @@ -14,7 +14,7 @@ # type hints T = TypeVar("T") Multiple = Union[Sequence[T], T] -UnitLike = Union[Unit, str] +UnitLike = Union[xr.DataArray, Unit, str] def coord_units( @@ -30,6 +30,7 @@ def coord_units( da: Input DataArray. coord_names: Name(s) of the coordinate(s) to be converted. new_units: Units to be converted from the current ones. + A DataArray that has units attribute is also accepted. equivalencies: Optional Astropy equivalencies. Returns: @@ -73,9 +74,9 @@ def frame(da: xr.DataArray, new_frame: str, /) -> xr.DataArray: lat_origin = da.coords["lat_origin"] # do not change the order below! - lon -= units(lon_origin, lon.attrs["units"]) + lon -= units(lon_origin, lon) lon *= np.cos(units(lat, "rad")) - lat -= units(lat_origin, lat.attrs["units"]) + lat -= units(lat_origin, lat) lon_origin *= 0.0 lat_origin *= 0.0 @@ -94,6 +95,7 @@ def units( Args: da: Input DataArray. new_units: Units to be converted from the current ones. + A DataArray that has units attribute is also accepted. equivalencies: Optional Astropy equivalencies. Returns: @@ -103,5 +105,8 @@ def units( if (units := da.attrs.get("units")) is None: raise ValueError("Units must exist in DataArray attrs.") + if isinstance(new_units, xr.DataArray): + new_units = new_units.attrs["units"] + new_data = Quantity(da, units).to(new_units, equivalencies) return da.copy(data=new_data).assign_attrs(units=new_units) diff --git a/decode/qlook.py b/decode/qlook.py index 00bf60c..e5cec6c 100644 --- a/decode/qlook.py +++ b/decode/qlook.py @@ -450,7 +450,7 @@ def calc_chan_weight( if method == "std/tx": tx = load.atm(type="eta").sel(pwv=float(pwv)) - freq = convert.units(dems.d2_mkid_frequency, tx.freq.attrs["units"]) + freq = convert.units(dems.d2_mkid_frequency, tx.freq) with catch_warnings(): simplefilter("ignore") From c1e9b873413799fce4fd228ec8d2cbe277ed56e4 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sat, 25 Nov 2023 00:33:14 +0000 Subject: [PATCH 29/31] #135 Update docstrings of utils.mad --- decode/utils.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/decode/utils.py b/decode/utils.py index 53e4f74..11df1a6 100644 --- a/decode/utils.py +++ b/decode/utils.py @@ -15,17 +15,17 @@ def mad( keep_attrs: Optional[bool] = None, **kwargs: Any, ) -> xr.DataArray: - """Calculate median absolute deviation (MAD) of a DataArray. + """Calculate the median absolute deviation (MAD) of a DataArray. Args: da: Input DataArray. - dim: Name of dimension(s) along which MAD is calculated. + dim: Name of dimension(s) along which the MAD is calculated. skipna: Same-name option to be passed to ``DataArray.median``. keep_attrs: Same-name option to be passed to ``DataArray.median``. kwargs: Same-name option(s) to be passed to ``DataArray.median``. Returns: - MAD of the input DataArray. + The MAD of the input DataArray. """ From 848b5fd5c9b17e915dd214925cda5f1a3568915c Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sat, 25 Nov 2023 00:45:14 +0000 Subject: [PATCH 30/31] =?UTF-8?q?#135=20Update=20package=20version=20(2.8.?= =?UTF-8?q?0=20=E2=86=92=202.9.0)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CITATION.cff | 2 +- README.md | 2 +- decode/__init__.py | 2 +- pyproject.toml | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/CITATION.cff b/CITATION.cff index 0ebb54f..6baab48 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -3,7 +3,7 @@ message: "If you use this software, please cite it as below." title: "de:code" abstract: "DESHIMA code for data analysis" -version: 2.8.0 +version: 2.9.0 date-released: 2023-11-12 license: "MIT" doi: "10.5281/zenodo.3384216" diff --git a/README.md b/README.md index ce962ce..fc2b423 100644 --- a/README.md +++ b/README.md @@ -11,5 +11,5 @@ DESHIMA code for data analysis ## Installation ```shell -pip install decode==2.8.0 +pip install decode==2.9.0 ``` diff --git a/decode/__init__.py b/decode/__init__.py index 3c4771d..3277d48 100644 --- a/decode/__init__.py +++ b/decode/__init__.py @@ -10,7 +10,7 @@ "select", "utils", ] -__version__ = "2.8.0" +__version__ = "2.9.0" # submodules diff --git a/pyproject.toml b/pyproject.toml index 23128bc..edde958 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "decode" -version = "2.8.0" +version = "2.9.0" description = "DESHIMA code for data analysis" authors = ["Akio Taniguchi "] keywords = [ From ad9ee4df45f157489fc00f3955fb0065526a20f5 Mon Sep 17 00:00:00 2001 From: Akio Taniguchi Date: Sat, 25 Nov 2023 00:45:34 +0000 Subject: [PATCH 31/31] #135 Update release date --- CITATION.cff | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CITATION.cff b/CITATION.cff index 6baab48..3d24c7b 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -4,7 +4,7 @@ message: "If you use this software, please cite it as below." title: "de:code" abstract: "DESHIMA code for data analysis" version: 2.9.0 -date-released: 2023-11-12 +date-released: 2023-11-25 license: "MIT" doi: "10.5281/zenodo.3384216" url: "https://github.com/deshima-dev/decode"