diff --git a/include/gnss_utils/observations.h b/include/gnss_utils/observations.h
index 852b35d646caed9d0d7967d2181158bb68f3ee8f..db2efc751470a5d37c781c243747a4193058968c 100644
--- a/include/gnss_utils/observations.h
+++ b/include/gnss_utils/observations.h
@@ -5,7 +5,6 @@
 
 namespace GnssUtils
 {
-
 class Observations
 {
 public:
@@ -15,34 +14,26 @@ public:
   ~Observations();
 
   void clearObservations();
-
   void addObservation(const obsd_t& obs);
-  void
-  loadFromRinex(const std::string& rnx_file, gtime_t t_start, gtime_t t_end, double dt = 0.0, const char* opt = "");
-
   void removeObservationByIdx(const int& _idx);
   void removeObservationBySat(const int& _sat);
-
   std::vector<obsd_t>&       getObservations();
   const std::vector<obsd_t>& getObservations() const;
-
   obsd_t&       getObservationBySat(const unsigned char& sat_number);
   const obsd_t& getObservationBySat(const unsigned char& sat_number) const;
-
   obsd_t&       getObservationByIdx(const int& idx);
   const obsd_t& getObservationByIdx(const int& idx) const;
 
   obsd_t*       data();
   const obsd_t* data() const;
-
   size_t size() const;
 
   bool hasSatellite(const unsigned char& i) const;
 
   static void print(obsd_t& _obs);
-  void printBySat(const int& _sat);
-  void printByIdx(const int& _idx);
-  void print();
+  void        printBySat(const int& _sat);
+  void        printByIdx(const int& _idx);
+  void        print();
 
   // Filter observations
   std::set<int> filterByEphemeris(const SatellitesPositions& sats_pos);
@@ -50,38 +41,35 @@ public:
   std::set<int> filterByCode();
   std::set<int> filterByCarrierPhase();
   std::set<int> filterByConstellations(const int& navsys);
-  std::set<int> filterByElevationSnr(const Eigen::Vector3d& x_r,
+  std::set<int> filterByElevationSnr(const Eigen::Vector3d&     x_r,
                                      const SatellitesPositions& sats_pos,
-                                     const snrmask_t& snrmask,
-                                     const double& elmin);
-  std::set<int> filter(const SatellitesPositions&   sats_pos,
-                       const std::set<int>&         discarded_sats,
-                       const Eigen::Vector3d&       x_r,
-                       const bool&                  check_code,
-                       const bool&                  check_carrier_phase,
-                       const prcopt_t&              opt);
-  std::set<int> filter(const SatellitesPositions&   sats_pos,
-                       const std::set<int>&         discarded_sats,
-                       const Eigen::Vector3d&       x_r,
-                       const bool&                  check_code,
-                       const bool&                  check_carrier_phase,
-                       const int&                   navsys,
-                       const snrmask_t&             snrmask,
-                       const double&                elmin);
-
-  static std::set<int> findCommonObservations(const Observations& obs_1,
-                                              const Observations& obs_2);
+                                     const snrmask_t&           snrmask,
+                                     const double&              elmin);
+  std::set<int> filter(const SatellitesPositions& sats_pos,
+                       const std::set<int>&       discarded_sats,
+                       const Eigen::Vector3d&     x_r,
+                       const bool&                check_code,
+                       const bool&                check_carrier_phase,
+                       const prcopt_t&            opt);
+  std::set<int> filter(const SatellitesPositions& sats_pos,
+                       const std::set<int>&       discarded_sats,
+                       const Eigen::Vector3d&     x_r,
+                       const bool&                check_code,
+                       const bool&                check_carrier_phase,
+                       const int&                 navsys,
+                       const snrmask_t&           snrmask,
+                       const double&              elmin);
+
+  static std::set<int> findCommonObservations(const Observations& obs_1, const Observations& obs_2);
 
   bool operator==(const Observations& other_obs) const;
-  bool operator !=(const Observations &other_obs) const;
+  bool operator!=(const Observations& other_obs) const;
 
 private:
   // Private objects
   std::map<unsigned char, int> sat_2_idx_;  //< key: corresponding sat number, value: idx in obs_ vector
   std::vector<unsigned char>   idx_2_sat_;  //< key: idx in obs_ vector, value: corresponding sat number
   std::vector<obsd_t>          obs_;        //< vector of RTKLIB observations for a given epoch
-
-  // Private methods
 };
 
 }  // namespace GnssUtils
@@ -92,7 +80,6 @@ private:
 
 namespace GnssUtils
 {
-
 inline const std::vector<obsd_t>& Observations::getObservations() const
 {
   return this->obs_;
@@ -147,9 +134,9 @@ inline bool Observations::hasSatellite(const unsigned char& i) const
   return sat_2_idx_.count(i) != 0;
 }
 
-inline bool Observations::operator !=(const Observations &other_obs) const
+inline bool Observations::operator!=(const Observations& other_obs) const
 {
-    return !(*this == other_obs);
+  return !(*this == other_obs);
 }
 
 }  // namespace GnssUtils
diff --git a/include/gnss_utils/utils/utils.h b/include/gnss_utils/utils/utils.h
index a107051641c6063ebfd6383e0682a741d8d6b605..72a1d77739e4a5574f07a779592122b9af9a11d8 100644
--- a/include/gnss_utils/utils/utils.h
+++ b/include/gnss_utils/utils/utils.h
@@ -8,6 +8,7 @@
 
 #include "gnss_utils/internal/config.h"
 #include "gnss_utils/gnss_utils.h"
+#include "gnss_utils/observations.h"
 
 #define ARRAY_SIZE(arr) sizeof(arr) / sizeof(arr[0])
 #define ARRAY2D_NROWS(arr) sizeof(arr) / sizeof(arr[0])
@@ -165,6 +166,9 @@ bool equalArray3d(const T*   array_1,
   return true;
 }
 
+std::vector<Observations>
+loadObservationsFromRinex(const std::string& rnx_file, gtime_t t_start, gtime_t t_end, double dt, const char* opt);
+
 }  // namespace GnssUtils
 
 bool operator==(const gtime_t& time1, const gtime_t& time2);
diff --git a/src/examples/gnss_utils_test.cpp b/src/examples/gnss_utils_test.cpp
index fba3b7f10045cd17d2238b30834b8f0c8f53e6f6..ce7a617b78a147ef258b70064d82823cac6c2e93 100644
--- a/src/examples/gnss_utils_test.cpp
+++ b/src/examples/gnss_utils_test.cpp
@@ -1,6 +1,7 @@
 #include "gnss_utils/observations.h"
 #include "gnss_utils/navigation.h"
 #include "gnss_utils/utils/transformations.h"
+#include "gnss_utils/utils/utils.h"
 
 #include <typeinfo>
 
@@ -33,14 +34,16 @@ int main(int argc, char* argv[])
   const char* opt = "";         // only GPS | GPS+GAL: "-SYS=G,L" | ALL: ""
 
   // load observations from RINEX file
-  observations.loadFromRinex("../src/examples/sample_data.obs", t_start, t_end, dt, opt);
+  std::vector<Observations> obs_rinex =
+      loadObservationsFromRinex("../src/examples/sample_data.obs", t_start, t_end, dt, opt);
+  observations = obs_rinex[0];
   observations.print();
 
   // RTKLIB trace
-  // char str_file[80];
-  // snprintf(str_file, sizeof str_file, "../src/examples/trace");
-  // tracelevel(5);
-  // traceopen(str_file);
+  char str_file[80];
+  snprintf(str_file, sizeof str_file, "../src/examples/trace");
+  tracelevel(5);
+  traceopen(str_file);
 
   // load navigation from RINEX file
   navigation.loadFromRinex("../src/examples/sample_data.nav", t_start, t_end, dt, opt);
@@ -89,15 +92,17 @@ int main(int argc, char* argv[])
 
   char msg[128] = "";
 
-  stat = pntpos(observations.data(), observations.size(), &navigation.getNavigation(), &prcopt, &solb, NULL, NULL, msg);
+  stat = pntpos(observations.data(), observations.size(), &navigation.getNavigation(), &prcopt, &solb, NULL, NULL,
+  msg);
 
   std::cout << "msg: " << msg << std::endl;
   std::cout << "sol.stat: " << int(solb.stat) << std::endl;
   std::cout << "Position:      " << solb.rr[0] << ", " << solb.rr[1] << ", " << solb.rr[2] << std::endl;
   std::cout << "Position (GT): " << latLonAltToEcef(lla_gt).transpose() << std::endl;
-  std::cout << "Position LLA:      " << ecefToLatLonAlt(Eigen::Vector3d(solb.rr[0], solb.rr[1], solb.rr[2])).transpose()
+  std::cout << "Position LLA:      " << ecefToLatLonAlt(Eigen::Vector3d(solb.rr[0], solb.rr[1],
+  solb.rr[2])).transpose()
             << std::endl;
   std::cout << "Position LLA (GT): " << lla_gt.transpose() << std::endl;
 
-  traceclose();
+  // traceclose();
 }
diff --git a/src/observations.cpp b/src/observations.cpp
index 0ca80fd9e9612ec569a9e295375a40cbb3e76748..a90ea7084d576746ccf7fb1ba746447fcc6af09a 100644
--- a/src/observations.cpp
+++ b/src/observations.cpp
@@ -74,41 +74,6 @@ void Observations::addObservation(const obsd_t& obs)
   assert(sat_2_idx_.size() == obs_.size());
 }
 
-void Observations::loadFromRinex(const std::string& rnx_file,
-                                 gtime_t            t_start,
-                                 gtime_t            t_end,
-                                 double             dt,
-                                 const char*        opt)
-{
-  obs_t obs;
-  obs.data = (obsd_t*)malloc(sizeof(obsd_t) * MAXSAT);
-  obs.n    = 0;
-  obs.nmax = MAXSAT;
-
-  // const char* opt = "";
-  auto stat = readrnxt(rnx_file.c_str(), 1, t_start, t_end, dt, opt, &obs, NULL, NULL);
-  if (stat == 1)
-    sortobs(&obs);
-  else
-  {
-    std::cout << "Observation: couldn't load provided observation file, reason: " << (stat == 0 ? "no data" : "error") << stat << std::endl;
-    return;
-  }
-
-  for (int i = 0; i < obs.n; i++)
-  {
-    // std::cout << "time: " << time_str(obs.data[i].time, 3) <<  " | sat: " << int(obs.data[i].sat) << " | rcv: " <<
-    // int(obs.data[i].rcv) <<
-    //           " | SNR: " << int(obs.data[i].SNR[0]) << " | LLI: " << int(obs.data[i].LLI[0]) << " | code: " <<
-    //           int(obs.data[i].code[0]) <<
-    //             " | L: " << obs.data[i].L[0] <<  " | P: " << obs.data[i].P[0] << " | D: " << obs.data[i].D[0] <<
-    //             std::endl;
-    addObservation(obs.data[i]);
-  }
-
-  free(obs.data);
-}
-
 void Observations::removeObservationByIdx(const int& _idx)
 {
   // std::cout << "removing observation of idx " << _idx << std::endl;
diff --git a/src/utils/utils.cpp b/src/utils/utils.cpp
index ce565b86ceed99855aefb54358ec006b38a41c90..2cb59ab57b524b097f0e756542bb64ea235aff1b 100644
--- a/src/utils/utils.cpp
+++ b/src/utils/utils.cpp
@@ -9,6 +9,55 @@ void print(std::string& _msg)
   std::cout << msg << "\n";
 }
 
+std::vector<Observations>
+loadObservationsFromRinex(const std::string& rnx_file, gtime_t t_start, gtime_t t_end, double dt, const char* opt)
+{
+  obs_t obs;
+  obs.data = (obsd_t*)malloc(sizeof(obsd_t) * MAXSAT);
+  obs.n    = 0;
+  obs.nmax = MAXSAT;
+
+  // const char* opt = "";
+  auto stat = readrnxt(rnx_file.c_str(), 1, t_start, t_end, dt, opt, &obs, NULL, NULL);
+  int  n_epochs;
+  if (stat == 1)
+    n_epochs = sortobs(&obs);
+  else
+  {
+    std::cout << "Observation: couldn't load provided observation file, reason: " << (stat == 0 ? "no data" : "error")
+              << stat << std::endl;
+    return;
+  }
+
+  std::vector<Observations> obs_vector;
+  obs_vector.reserve(n_epochs);
+
+  Observations obs_epoch;
+  for (int i = 0; i < obs.n; i++)
+  {
+    obs_epoch.addObservation(obs.data[i]);
+
+    bool last_obs_epoch = false;
+    if (i < obs.n - 1)
+    {
+      double tt = timediff(obs.data[i].time, obs.data[i + 1].time);
+      if (fabs(tt) > DTTOL)
+      {
+        obs_vector.push_back(obs_epoch);
+        obs_epoch.clearObservations();
+      }
+    }
+    else
+    {
+      obs_vector.push_back(obs_epoch);
+    }
+  }
+
+  free(obs.data);
+
+  return obs_vector;
+}
+
 }  // namespace GnssUtils
 
 bool operator==(const gtime_t& time1, const gtime_t& time2)
diff --git a/test/gtest_observations.cpp b/test/gtest_observations.cpp
index cbb25d72d1863ea094dd2bd83cc28c3a8f60a5d6..4b0699d64a63e1993756412e1970209b1a3ebe1d 100644
--- a/test/gtest_observations.cpp
+++ b/test/gtest_observations.cpp
@@ -55,7 +55,9 @@ TEST(ObservationsTest, LoadFromRinex)
 
   // GnssUtils utilities
   Observations observations;
-  observations.loadFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+
+  std::vector<Observations> obs_rinex = loadObservationsFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+  observations                        = obs_rinex[0];
 
   ASSERT_EQ(obs.n, 18);
 
@@ -75,7 +77,9 @@ TEST(ObservationsTest, GetObservationBySat)
   loadRinex(obs);
 
   Observations observations;
-  observations.loadFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+
+  std::vector<Observations> obs_rinex = loadObservationsFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+  observations                        = obs_rinex[0];
 
   for (int ii = 0; ii < obs.n; ++ii)
   {
@@ -90,7 +94,9 @@ TEST(ObservationsTest, GetObservationByIdx)
   loadRinex(obs);
 
   Observations observations;
-  observations.loadFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+
+  std::vector<Observations> obs_rinex = loadObservationsFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+  observations                        = obs_rinex[0];
 
   for (int ii = 0; ii < obs.n; ++ii)
   {
@@ -104,7 +110,9 @@ TEST(ObservationsTest, data)
   loadRinex(obs);
 
   Observations observations;
-  observations.loadFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+
+  std::vector<Observations> obs_rinex = loadObservationsFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+  observations                        = obs_rinex[0];
 
   for (int ii = 0; ii < obs.n; ++ii)
   {
@@ -118,7 +126,9 @@ TEST(ObservationsTest, HasSatellite)
   loadRinex(obs);
 
   Observations observations;
-  observations.loadFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+
+  std::vector<Observations> obs_rinex = loadObservationsFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+  observations                        = obs_rinex[0];
 
   for (int ii = 0; ii < obs.n; ++ii)
   {
@@ -129,7 +139,10 @@ TEST(ObservationsTest, HasSatellite)
 TEST(ObservationsTest, FindCommonObservations)
 {
   Observations observations1;
-  observations1.loadFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+
+  std::vector<Observations> obs_rinex = loadObservationsFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+  observations1                       = obs_rinex[0];
+
   Observations observations2(observations1);
 
   ASSERT_TRUE(observations1 == observations2);
@@ -140,15 +153,17 @@ TEST(ObservationsTest, FindCommonObservations)
   ASSERT_TRUE(common_sats.size() == observations2.size());
   for (auto sat : common_sats)
   {
-      ASSERT_TRUE(observations1.hasSatellite(sat));
-      ASSERT_TRUE(observations2.hasSatellite(sat));
+    ASSERT_TRUE(observations1.hasSatellite(sat));
+    ASSERT_TRUE(observations2.hasSatellite(sat));
   }
 }
 
 TEST(ObservationsTest, FindCommonObservationsRemoved)
 {
   Observations observations1;
-  observations1.loadFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+
+  std::vector<Observations> obs_rinex = loadObservationsFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+  observations1                       = obs_rinex[0];
   Observations observations2(observations1);
 
   // Remove first observation of observations2
@@ -161,15 +176,18 @@ TEST(ObservationsTest, FindCommonObservationsRemoved)
   ASSERT_TRUE(common_sats.size() == observations2.size());
   for (auto sat : common_sats)
   {
-      ASSERT_TRUE(observations1.hasSatellite(sat));
-      ASSERT_TRUE(observations2.hasSatellite(sat));
+    ASSERT_TRUE(observations1.hasSatellite(sat));
+    ASSERT_TRUE(observations2.hasSatellite(sat));
   }
 }
 
 TEST(ObservationsTest, FindCommonObservationsChangeTime)
 {
   Observations observations1;
-  observations1.loadFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+
+  std::vector<Observations> obs_rinex = loadObservationsFromRinex(rnx_file.c_str(), t_start, t_end, dt, opt);
+  observations1                       = obs_rinex[0];
+
   Observations observations2(observations1);
 
   // Change time
@@ -184,8 +202,8 @@ TEST(ObservationsTest, FindCommonObservationsChangeTime)
   ASSERT_TRUE(common_sats.size() == observations2.size());
   for (auto sat : common_sats)
   {
-      ASSERT_TRUE(observations1.hasSatellite(sat));
-      ASSERT_TRUE(observations2.hasSatellite(sat));
+    ASSERT_TRUE(observations1.hasSatellite(sat));
+    ASSERT_TRUE(observations2.hasSatellite(sat));
   }
 }