bugfix: incorrect notation of the rotation vector

floenavi/__init__.py

@@ -140,24 +140,29 @@ class FloeNavi(object):
         :return: icedrift.GeoReferenceStation (+ FloeNaviGrid if return_grid_solution is True)
         """
 
-        # Use default base station setup
+        # Use default base station setup if not specifically supplied
         if basestation_setup is None:
             basestation_setup = self.cfg.get_floenavi_cfg(dataset_time_bounds)
 
-        # Load the collection of FloeNavi basestation files
-        logger.info("--- Establish FloeNavi grid ---")
-        floenavi_grid = FloeNaviGrid(self.cfg, basestation_setup)
-        floenavi_grid.set_time_bounds(dataset_time_bounds)
-
-        # Get Reference station data from two positions
-        pos_ref, pos_hdg, fngrid = floenavi_grid.get_virtual_reference_station_positions(
-            dataset_time_bounds, **refstat_kwargs)
-        if pos_ref is not None:
-            refstat_ais = GeoReferenceStation.from_two_basestations(pos_ref, pos_hdg)
-            refstat_ais.set_uncertainty(fngrid.origin_uncertainty_m, fngrid.heading_uncertainty_deg)
-        else:
+        # Check if base station setup exists
+        if basestation_setup is None:
             refstat_ais, fngrid = None, None
 
+        else:
+            # Load the collection of FloeNavi basestation files
+            logger.info("--- Establish FloeNavi grid ---")
+            floenavi_grid = FloeNaviGrid(self.cfg, basestation_setup)
+            floenavi_grid.set_time_bounds(dataset_time_bounds)
+
+            # Get Reference station data from two positions
+            pos_ref, pos_hdg, fngrid = floenavi_grid.get_virtual_reference_station_positions(
+                dataset_time_bounds, **refstat_kwargs)
+            if pos_ref is not None:
+                refstat_ais = GeoReferenceStation.from_two_basestations(pos_ref, pos_hdg)
+                refstat_ais.set_uncertainty(fngrid.origin_uncertainty_m, fngrid.heading_uncertainty_deg)
+            else:
+                refstat_ais, fngrid = None, None
+
         # Return value
         if not return_grid_solution:
             return refstat_ais
@@ -1125,12 +1130,20 @@ class FloeNaviGridSolution(object):
         origins_xc = np.full((2, self.bs1.n_records), np.nan)
         origins_yc = origins_xc.copy()
         for i, bs in enumerate([self.bs1, self.bs2]):
-            xc, yc = self.get_origin(bs.geo.xc, bs.geo.yc, bs.range_m, bs.ref_angle, self.rotation)
+            xc, yc = self.get_origin(bs.geo.longitude, bs.geo.latitude, bs.range_m, bs.ref_angle, self.rotation)
             invalid = np.where(np.logical_or(xc > 1.0e10, yc > 1.0e10))[0]
             xc[invalid], yc[invalid] = np.nan, np.nan
             origins_xc[i, :] = xc
             origins_yc[i, :] = yc
 
+        # import matplotlib.pyplot as plt
+        # plt.figure()
+        # for i, bs in enumerate([self.bs1, self.bs2]):
+        #     # plt.scatter(bs.geo.xc, bs.geo.yc, s=40, label="geo {}".format(bs.mmsi))
+        #     plt.scatter(origins_xc[i, :], origins_yc[i, :], s=5, label="origin {}".format(bs.mmsi))
+        # plt.legend()
+        # plt.show()
+
         # Compute the average deviation (uncertainty) of the origin
         origin_xc_delta = origins_xc[1]-origins_xc[0]
         origin_yc_delta = origins_yc[1]-origins_yc[0]
@@ -1168,8 +1181,7 @@ class FloeNaviGridSolution(object):
             angle = np.rad2deg(angle)
         return angle
 
-    @staticmethod
-    def get_origin(xc, yc, range_m, ref_heading, rotation):
+    def get_origin(self, lon, lat, range_m, ref_heading, rotation):
         """
         Compute the origin from basestation data and the known rotation of the
         ice reference system
@@ -1180,11 +1192,17 @@ class FloeNaviGridSolution(object):
         :param rotation:
         :return:
         """
-        direction_to_origin_deg = 180. - (rotation - ref_heading)
-        direction_to_origin_rad = np.deg2rad(direction_to_origin_deg)
-        vector_x = range_m * np.cos(direction_to_origin_rad)
-        vector_y = range_m * np.sin(direction_to_origin_rad)
-        return xc + vector_x, yc + vector_y
+
+        # Step 1: Project the origin positions in geographic space
+        direction_to_origin_deg = rotation - ref_heading + 180.
+        g = pyproj.Geod(ellps='WGS84')
+        origin_lon, origin_lat, _ = g.fwd(lon, lat, direction_to_origin_deg, np.full(lon.shape, range_m))
+
+        # Step 2: Compute x, y in local projection coordinates
+        p = pyproj.Proj(**self.proj_dict)
+        xc_orig, yc_orig = p(origin_lon, origin_lat)
+
+        return xc_orig, yc_orig
 
     @property
     def proj_dict(self):