Nishikawa/io/#96

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

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