RdStationPositionCorrection.cc

Go to the documentation of this file.

#include "RdStationPositionCorrection.h"

#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CentralConfig.h>

#include <det/VManager.h>
#include <rdet/RDetector.h>

#include <utl/ReadStream.h>
#include <utl/UTMPoint.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/PhysicalConstants.h>
#include <utl/UTCDateTime.h>

#include <evt/Event.h>
#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/StationRecData.h>

#include <vector>


using namespace evt;
using namespace revt;
using namespace std;
using namespace utl;
using namespace fwk;


namespace RdStationPositionCorrectionOG {

  struct UTM {
    double fNorthing;
    double fEasting;
    double fHeight;
  };


  inline
  istream&
  operator>>(istream& is, UTM& utm)
  {
    return is >> utm.fNorthing >> ws >> utm.fEasting >> ws >> utm.fHeight;
  }


  struct KnownPositions {
    unsigned int fStationId;
    UTM fUsedPosition;
    UTM fTruePosition;
    UTCDateTime fCorrectionValidUntil;
  };


  inline
  istream&
  operator>>(istream& is, KnownPositions& pos)
  {
    return is >> pos.fStationId >> ws >> pos.fUsedPosition >> ws >> pos.fTruePosition >> ws
              >> pos.fCorrectionValidUntil;
  }


  VModule::ResultFlag
  RdStationPositionCorrection::Init()
  {
    Branch topB = CentralConfig::GetInstance()->GetTopBranch("RdStationPositionCorrection");

    fDistanceTolerance = topB.GetChild("DistanceTolerance").Get<double>();
    const double effectiveVerticalDelay = topB.GetChild("EffectiveVerticalDelay").Get<double>();
    const string rdStationGPSPositionsFile = topB.GetChild("RdStationGPSPositionsFile").Get<string>();

    // trigger detector update to read standard station positions
    det::Detector& det = det::Detector::GetInstance();
    det.Update(GetCurrentSystemTime());
    const rdet::RDetector& rDet = det.GetRDetector();

    ifstream posFile(rdStationGPSPositionsFile.c_str());
    if (!posFile.is_open()) {
      ostringstream err;
      err << "Cannot open GPS positions file '" << rdStationGPSPositionsFile << "'!";
      ERROR(err);
      return eFailure;
    }

    vector<KnownPositions> positions;
    ReadStream(posFile).GetLines(positions);

    fOffsets.clear();
    fValidityTimes.clear();

    const int zone = 19;
    const char band = 'H';
    const ReferenceEllipsoid& ellipsoid = ReferenceEllipsoid::GetWGS84();
    const CoordinateSystemPtr siteCS = det.GetSiteCoordinateSystem();

    for (const auto& pos : positions) {
      // get position used in GPS unit
      const unsigned int sId = pos.fStationId;
      const UTM& usedpos = pos.fUsedPosition;
      const UTMPoint usedposUTM(usedpos.fNorthing, usedpos.fEasting, usedpos.fHeight, zone, band, ellipsoid);
      const Point usedPosition = usedposUTM.GetPoint(siteCS);

      try {
        // get true position which should have been used in GPS unit
        const UTM& truepos = pos.fTruePosition;
        const UTMPoint trueposUTM(truepos.fNorthing, truepos.fEasting, truepos.fHeight, zone, band, ellipsoid);
        const Point truePosition = trueposUTM.GetPoint(siteCS);
        const CoordinateSystemPtr localCS = LocalCoordinateSystem::Create(truePosition);

        const double dt = (truePosition - usedPosition).GetZ(localCS) * effectiveVerticalDelay;

        // check for double entries in GPS positions correction file
        if (fOffsets.find(sId) != fOffsets.end()) {
          ostringstream err;
          err << "Duplicate entry for station id " << sId << " found in the positions file!";
          ERROR(err);
          return eFailure;
        }

        // store offset and validity time of correction
        fOffsets.emplace(sId, dt);
        fValidityTimes.emplace(sId, pos.fCorrectionValidUntil.GetTimeStamp());

        // cross-check that position listed in GPS correction file is consistent with position in RDetector
        const rdet::Station& dStation = rDet.GetStation(sId);
        const Point detPosition = dStation.GetPosition();
        const double posDistance = (truePosition - detPosition).GetMag();

        if (posDistance > fDistanceTolerance) {
          ostringstream err;
          err << "Distance of " << posDistance / meter << " m between current station position in RDetector"
                 " and GPS corrections file for station id " << sId << " is too large!";
          ERROR(err);
          return eFailure;
        }

      } catch (NonExistentComponentException&) {
        ostringstream warn;
        warn << "Station id " << sId << " has a GPS position in the corrections file but is not found in the RDetector!";
        WARNING(warn);
      }
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdStationPositionCorrection::Run(Event& event)
  {
    if (!event.HasREvent())
      return eSuccess;

    REvent& rEvent = event.GetREvent();

    for (auto& station : rEvent.StationsRange()) {
      if (!station.HasRecData())
        continue;

      StationRecData& recData = station.GetRecData();

      if (!recData.HasParameter(eTraceStartTime)) {
        ostringstream err;
        err << "Trying to correct trace start time for station id " << station.GetId()
            << ", but no trace start time available in ParameterStorage!";
        ERROR(err);
        return eFailure;
      }

      const TimeInterval oldTime = recData.GetParameter(eTraceStartTime);
      const auto offsetIt = fOffsets.find(station.GetId());

      if (offsetIt == fOffsets.end()) // no correction value found for this station
        continue;

      // check validity time to see if correction needs to be applied, it is ensured that the entry is present in the map
      if (rEvent.GetHeader().GetTime() > fValidityTimes[station.GetId()])
        continue;

      const TimeInterval& offset = offsetIt->second;
      recData.SetParameter(eTraceStartTime, oldTime - offset, false); // set corrected start time
      recData.SetParameter(eTraceStartTimeCorrectionByGPSPosition, -offset); // store value that has been added to the trace start time

      // also change relative position of signal search window accordingly
      recData.SetParameter(eSignalSearchWindowStart, recData.GetParameter(eSignalSearchWindowStart) + offset, false);
      recData.SetParameter(eSignalSearchWindowStop, recData.GetParameter(eSignalSearchWindowStop) + offset, false);
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdStationPositionCorrection::Finish()
  {
    return eSuccess;
  }

}