diff --git a/nasco_analysis/doppler.py b/nasco_analysis/doppler.py index 9ef75d1..bd07835 100644 --- a/nasco_analysis/doppler.py +++ b/nasco_analysis/doppler.py @@ -243,11 +243,9 @@ def ch_speed(self): """ spec_freq = self.heterodyne() freq_resolution = self.band_width / self.ch_num - speed_resolution = -1 * self.band_width / self.rest_freq / self.ch_num * const.c - v_0GHz = ( - -1 * speed_resolution * spec_freq / freq_resolution * self.sideband_factor - ) - v_base = np.arange(self.ch_num) * speed_resolution * self.sideband_factor + speed_resolution = self.band_width / self.rest_freq / self.ch_num * const.c + v_0GHz = speed_resolution * spec_freq / freq_resolution * self.sideband_factor + v_base = -1 * np.arange(self.ch_num) * speed_resolution * self.sideband_factor v_apparent = v_base + v_0GHz v_obs = self.v_obs() v_lsr = v_apparent + np.atleast_2d(v_obs).T diff --git a/nasco_analysis/io.py b/nasco_analysis/io.py index 808b6c0..2e4ba5a 100644 --- a/nasco_analysis/io.py +++ b/nasco_analysis/io.py @@ -83,7 +83,7 @@ def convert( data["timestamp"] = timestamp2datetime(data["timestamp"].astype(float)) obsmode["received_time"] = timestamp2datetime( obsmode["received_time"].astype(float) - ) + ) # no timestamp recorded encoder["timestamp"] = timestamp2datetime(encoder["timestamp"].astype(float)) weather["timestamp"] = timestamp2datetime(weather["timestamp"].astype(float)) @@ -126,12 +126,12 @@ def correct_kisa(self) -> Tuple[xr.DataArray]: enc_az_array = self.encoder_set.enc_az enc_el_array = self.encoder_set.enc_el - d_az, d_el = apply_kisa_test( + dAz, dEl = apply_kisa_test( azel=(enc_az_array, enc_el_array), hosei=self.kisa_path ) - az = enc_az_array + d_az / 3600 - el = enc_el_array + d_el / 3600 + az = enc_az_array + dAz / 3600 + el = enc_el_array + dEl / 3600 self.encoder_set = self.encoder_set.assign( az=("t", az), el=("t", el) @@ -140,12 +140,55 @@ def correct_kisa(self) -> Tuple[xr.DataArray]: return az, el def correct_collimation_error( - self, collimation_params: Dict[str, float] + self, collimation_params: n2const.Constants = None ) -> Tuple[xr.DataArray]: - return NotImplemented - - def combine_metadata(self, int_obsmode: bool = False) -> xr.DataArray: + """ + Examples + -------- + >>> from n_const.constants import Constants + >>> params = Constants.set_values(r=300, theta=80, d1=-10, d2=35) + >>> InitialArray(...).correct_collimation_error(params) + """ self.correct_kisa() + if collimation_params is None: + self.encoder_set = self.encoder_set.assign_attrs( + {"collimation_applied": False} + ) + return + + # specify the beam by self.topic_name using newly defined Constants? + + r = collimation_params.r # arcsec + theta = np.deg2rad(collimation_params.theta) + d1 = collimation_params.d1 # arcsec + d2 = collimation_params.d2 # arcsec + el = np.deg2rad(self.data.el) + + dx = (r * np.cos(theta - el) + d1) / 3600 # deg + dAz = dx / np.cos(el) # deg + dEl = (r * np.sin(theta - el) + d2) / 3600 # deg + + az = self.encoder_set.az + dAz + el = self.encoder_set.el + dEl + + self.encoder_set = self.encoder_set.assign( + az=("t", az), el=("t", el) + ).assign_attrs({"collimation_applied": True}) + return az, el + + def combine_metadata( + self, + collimation_params: Optional[n2const.Constants] = None, + int_obsmode: bool = False, + ) -> xr.DataArray: + """ + Notes + ----- + TODO: remove collimation_params and add self.corrim_path, + determine beam using n2const.BOARD2BEAM[self.topic_name] + """ + + self.correct_collimation_error(collimation_params) # kisa correction also done if int_obsmode: @@ -171,7 +214,7 @@ def convert(mode: bytes) -> float: reindexed_obsmode_set = reindexed_obsmode_set.where( reindexed_obsmode_set == self.obsmode_set.reindex_like(self.data_set, method="bfill"), - b"None", + # b"None", ) reindexed_encoder_set = self.encoder_set.interp_like( self.data_set, method="linear" @@ -197,10 +240,25 @@ def convert(mode: bytes) -> float: for key in reindexed_weather_set.keys(): coords[key] = ("t", reindexed_weather_set[key]) + unit = { + "in_humi": "percent", + "out_humi": "percent", + "wind_sp": "m/s", + "wind_dir": "deg", + "press": "hPa", + "rain": "mm", + "v_lsr": "km/s", + } + unit.update(dict.fromkeys(["in_temp", "out_temp"], "K")) + unit.update(dict.fromkeys(["cabin_temp1", "cabin_temp2"], "K")) + unit.update(dict.fromkeys(["dome_temp1", "dome_temp2"], "K")) + unit.update(dict.fromkeys(["gen_temp1", "gen_temp2"], "K")) + unit.update(dict.fromkeys(["az", "el", "ra", "dec", "l", "b"], "deg")) + self.data = ( self.data_set[main_data_field] .assign_coords(coords) - .assign_attrs(self.encoder_set.attrs) + .assign_attrs(self.encoder_set.attrs, unit=unit) ) return self.data @@ -217,7 +275,7 @@ def channel2velocity(self, args: Optional[dict] = None, **kwargs) -> xr.DataArra if "spec" not in self.data_set.keys(): raise ValueError( "Spectral data not found." - "Cannot calculate doppler velocity to total power data." + "Cannot calculate doppler velocity on total power data." ) kwargs.update( @@ -228,21 +286,22 @@ def channel2velocity(self, args: Optional[dict] = None, **kwargs) -> xr.DataArra } ) dp = Doppler(args=args, **kwargs) - self.velocity = dp.ch_speed() + velocity = dp.ch_speed() - self.data = self.data.assign_coords(v_lsr=(["t", "spectra"], self.velocity)) + self.data = self.data.assign_coords(v_lsr=(["t", "spectra"], velocity)) return self.data def convert_coordinates(self) -> xr.DataArray: - location = n2const.LOC_NANTEN2 - horizontal_coord = SkyCoord( - az=self.data.az, - alt=self.data.el, + az=self.data.az.data * u.deg, + alt=self.data.el.data * u.deg, frame=AltAz, obstime=self.data.t, - location=location, - unit="deg", + location=n2const.LOC_NANTEN2, + pressure=self.data.press.data * u.hPa, + temperature=self.data.out_temp.data * u.deg_C, + relative_humidity=self.data.out_humi.data * u.percent, + # refraction have almost no dependency on obs_wl at RADIO frequency ) equatorial_coord = horizontal_coord.transform_to(FK5) galactic_coord = equatorial_coord.transform_to(Galactic) @@ -251,7 +310,7 @@ def convert_coordinates(self) -> xr.DataArray: ra=("t", equatorial_coord.ra), dec=("t", equatorial_coord.dec), l=("t", galactic_coord.l), - b=("t", galactic_coord.l), + b=("t", galactic_coord.b), ) return self.data @@ -260,11 +319,22 @@ def get_spectral_data( data_path: PathLike, topic_name: str, kisa_path: Optional[PathLike] = None, + collimation_params: Optional[n2const.Constants] = None, trans_coord: Optional[bool] = False, args: Optional[dict] = None, - **kwargs: Optional[Any] + **kwargs: Optional[Any], ) -> xr.DataArray: """ + Parameters + ---------- + data_path: PathLike + topic_name: str + kisa_path: PathLike + trans_coord: bool + args: dict + **kwargs: Any + kwargs for doppler.Doppler class + Notes ----- Wall time: @@ -273,7 +343,9 @@ def get_spectral_data( """ ia = InitialArray(data_path, topic_name, kisa_path) ia.create_data_array() - ret = ia.combine_metadata() # kisa correction implicitly done here + ret = ia.combine_metadata( + collimation_params=collimation_params + ) # kisa correction implicitly done here if not trans_coord: return ret ret = ia.convert_coordinates() @@ -286,6 +358,7 @@ def get_totalpower_data( data_path: PathLike, topic_name: str, kisa_path: Optional[PathLike] = None, + collimation_params: Optional[n2const.Constants] = None, trans_coord: Optional[bool] = False, ) -> xr.DataArray: """ @@ -297,7 +370,9 @@ def get_totalpower_data( """ ia = InitialArray(data_path, topic_name, kisa_path) ia.create_tp_array() - ret = ia.combine_metadata() # kisa correction implicitly done here + ret = ia.combine_metadata( + collimation_params=collimation_params + ) # kisa correction implicitly done here if not trans_coord: return ret return ia.convert_coordinates() @@ -350,8 +425,8 @@ def concatenate(self): def ch2velo(self, arg_dict): self.channel2velocity(arg_dict) - # variable mapping - self.velo_list = self.velocity + # variable mappings + self.velo_list = self.data.v_lsr self.concatenated_array = self.data return self.velo_list