GOESDB.cc

Go to the documentation of this file.

#include <det/Detector.h>

#include <utl/ErrorLogger.h>
#include <utl/Point.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/AugerException.h>

#include <atm/Atmosphere.h>
#include <atm/GOESDB.h>

using namespace det;
using namespace utl;
using namespace std;


namespace atm {

  // some hardcoded parameters related to the GOES pixelization of the Auger array

  namespace {
    const int gNumEasting = 30; // # of pixels in easting
    const int gNumNorthing = 12; // # of pixels in northing
    const int gNumCloudProbabilityLevels = 5; // DB returns from [0, 4]
    const double gStartEasting = (440000 + 553.2)*m; // start of easting coverage
    const double gStartNorthing = (6060000 + 9736.34)*m; // start of northing coverage
    const double gDeltaEasting = 2394.9*m; // size of pixel in easting
    const double gDeltaNorthing = 5493.03*m; // size of pixel in northing
  }


  double
  GOESDB::GetPixelWidthEasting()
    const
  {
    return gDeltaEasting;
  }


  double
  GOESDB::GetPixelWidthNorthing()
    const
  {
    return gDeltaNorthing;
  }


  unsigned int
  GOESDB::GetNumberOfPixels()
    const
  {
    return gNumEasting * gNumNorthing;
  }


  int
  GOESDB::GetPixelId(const double easting, const double northing)
    const
  {
    const int northPixel = GetNorthPixel(northing);
    const int eastPixel = GetEastPixel(easting, northPixel);
    return gNumEasting * northPixel + eastPixel;
  }


  int
  GOESDB::GetNorthPixel(const double northing)
    const
  {
    const double deltaNorthing = GetPixelWidthNorthing();
    return (northing - gStartNorthing) / deltaNorthing;
  }


  int
  GOESDB::GetEastPixel(const double easting, const int northPixel)
    const
  {
    const double deltaEasting = GetPixelWidthEasting();
    return (easting - (gStartEasting - (northPixel%6) *
                       deltaEasting / 6)) / deltaEasting;
  }


  const utl::UTMPoint&
  GOESDB::GetPixelCenter(const unsigned int pixelId)
    const
  {
    if (pixelId >= GetNumberOfPixels()) {
      ostringstream errMsg;
      errMsg << " pixel id " << pixelId << " out of bound ";
      throw OutOfBoundException(errMsg.str());
    }

    if (fPixelCenters.empty()) {

      const double deltaEasting = GetPixelWidthEasting();
      const double deltaNorthing = GetPixelWidthNorthing();
      const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);

      for (unsigned int i = 0; i < GetNumberOfPixels(); ++i) {
        const unsigned int eastPixel = i % gNumEasting;
        const unsigned int northPixel = i / gNumEasting;
        const double easting = gStartEasting + eastPixel * deltaEasting
          - (northPixel%6) * deltaEasting / 6 + deltaEasting / 2;
        const double northing = gStartNorthing + northPixel * deltaNorthing
          + deltaNorthing / 2;
        const UTMPoint thisPoint(northing, easting, 1400*m, 19, 'H', wgs84);
        fPixelCenters.push_back(thisPoint);
      }

    }

    return fPixelCenters[pixelId];
  }


  double
  GOESDB::GetCloudProbability(const unsigned int pId)
    const
  {
    return IndexToProbability(GetCloudProbabilityIndex(pId));
  }


  double
  GOESDB::IndexToProbability(const unsigned int cpIndex)
    const
  {
    return cpIndex / double(gNumCloudProbabilityLevels);
  }


  int
  GOESDB::GetCloudProbabilityIndex(const unsigned int pId)
    const
  {
    if (!HasData(pId)) {
      VManager::IndexMap im;
      im["pixel_id"] = boost::lexical_cast<string>(pId);
      throw DataNotFoundInDBException("cloud_probability", "cloud_pixel", im,
                                      boost::lexical_cast<string>(Detector::GetInstance().GetTime()));
    }
    return fPixelMap[pId];
  }


  bool
  GOESDB::HasData(const unsigned int pId)
    const
  {
    if (fPixelMap.find(pId) != fPixelMap.end())
      return true;

    const VManager& manager = Detector::GetInstance().GetAManagerRegister().GetManager("AGOESSQLManager");

    int cloudProb;
    const VManager::IndexMap im({{"pixel_id", boost::lexical_cast<string>(pId)}});
    if (manager.GetData(cloudProb, "cloud_pixel", "cloud_probability", im) != VManager::eFound)
      return false;

    if (cloudProb < 0 || cloudProb > gNumCloudProbabilityLevels - 1) {
      ostringstream errMsg;
      errMsg << " cloud index from DB out of range! Index = " << cloudProb
             << ", expected range is [" << 0 << ", "
             << gNumCloudProbabilityLevels - 1 << "]"
             << " (maybe gNumCloudProbabilityLevels needs to be changed?)";
      throw OutOfBoundException(errMsg.str());
    }
    fPixelMap[pId] = cloudProb;
    return true;
  }


  GOESDB::ProbabilityMap
  GOESDB::GetAllCloudProbabilities()
    const
  {
    if (fPixelMap.size() < GetNumberOfPixels()) {

      const VManager& vmanager =
        Detector::GetInstance().GetAManagerRegister().GetManager("AGOESSQLManager");
      const AGOESSQLManager& manager = dynamic_cast<const AGOESSQLManager&>(vmanager);

      const VManager::IndexMap im;

      if (manager.GetData(fPixelMap, "cloud_pixel", "cloud_probability", im) == VManager::eFound) {

        auto maxCloudProb = max_element(fPixelMap.begin(), fPixelMap.end());
        if (maxCloudProb->second < 0 || maxCloudProb->second > gNumCloudProbabilityLevels - 1) {
          ostringstream errMsg;
          errMsg << " cloud index from DB out of range! Index = " << maxCloudProb->first
                 << ", expected range is [" << 0 << ", "
                 << gNumCloudProbabilityLevels - 1 << "]"
                 << " (maybe gNumCloudProbabilityLevels needs to be changed?)";
          throw OutOfBoundException(errMsg.str());
        }
      } else {
        ostringstream msg;
        msg << "request for data not found in the database for time "
            << Detector::GetInstance().GetTime() << endl;
        // ERROR(msg);
        VManager::IndexMap im;
        throw DataNotFoundInDBException("pixel_id,cloud_probability", "cloud_pixel", im,
                                        boost::lexical_cast<string>(Detector::GetInstance().GetTime()));
      }

    }

    ProbabilityMap retval;
    for (const auto& p : fPixelMap)
      retval[p.first] = IndexToProbability(p.second);

    return retval;
  }


  double
  GOESDB::GetCloudProbability(const Point& p)
    const
  {
    utl::UTMPoint utmPoint(p, ReferenceEllipsoid::GetWGS84());
    return GetCloudProbability(utmPoint);
  }


  double
  GOESDB::GetCloudProbability(const UTMPoint& p)
    const
  {
    const double x = p.GetEasting();
    const double y = p.GetNorthing();
    return GetCloudProbability(x, y);
  }


  double
  GOESDB::GetCloudProbability(const double easting, const double northing)
    const
  {
    const int pId = GetPixelId(easting, northing);
    if (pId < 0 || pId >= int(GetNumberOfPixels()))
      throw NonExistentComponentException("Request Point is not within the range covered by GOES data");

    return GetCloudProbability(pId);
  }


  // utility function for intersection of two line segments
  inline
  bool
  SegmentIntersection(double p0X, double p0Y, double p1X, double p1Y,
                      double p2X, double p2Y, double p3X, double p3Y)
  {
    const double s1X = p1X - p0X;
    const double s1Y = p1Y - p0Y;
    const double s2X = p3X - p2X;
    const double s2Y = p3Y - p2Y;

    if (s1X * s2Y == s2X * s1Y)
      return false;

    const double d = s1X * s2Y - s2X * s1Y;

    const double s = (-s1Y * (p0X - p2X) + s1X * (p0Y - p2Y)) / d;
    const double t = ( s2X * (p0Y - p2Y) - s2Y * (p0X - p2X)) / d;

    return (s >= 0 && s <= 1 && t >= 0 && t <= 1);
  }


  double
  GOESDB::GetMaximumCloudProbability(const utl::Point& pos1,
                                     const utl::Point& pos2)
    const
  {
    const double dX = GetPixelWidthEasting() / 2;
    const double dY = GetPixelWidthNorthing() / 2;

    UTMPoint p1(pos1, ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));
    UTMPoint p2(pos2, ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));
    const double x1 = p1.GetEasting();
    const double y1 = p1.GetNorthing();
    const double x2 = p2.GetEasting();
    const double y2 = p2.GetNorthing();

    // the pixels that contain point 1 and 2
    const int pixelId1 = GetPixelId(x1, y1);
    const int pixelId2 = GetPixelId(x2, y2);

    int maxCloudProbabilityIndex = -1;

    for (int i = 0; i < int(GetNumberOfPixels()); ++i) {
      bool isCloudPixel = false;
      if (i == pixelId1 || i == pixelId2)
        isCloudPixel = true;
      else {
        const utl::UTMPoint& pixelCenter = GetPixelCenter(i);
        const double x = pixelCenter.GetEasting();
        const double y = pixelCenter.GetNorthing();
        if (SegmentIntersection(x1, y1, x2, y2, x - dX, y - dY, x - dX, y + dY) ||
            SegmentIntersection(x1, y1, x2, y2, x + dX, y - dY, x + dX, y + dY) ||
            SegmentIntersection(x1, y1, x2, y2, x - dX, y - dY, x + dX, y - dY) ||
            SegmentIntersection(x1, y1, x2, y2, x - dX, y + dY, x + dX, y + dY))
          isCloudPixel = true;
      }

      if (isCloudPixel) {
        const int thisIndex = GetCloudProbabilityIndex(i);
        if (thisIndex > maxCloudProbabilityIndex)
          maxCloudProbabilityIndex = thisIndex;
      }
    }

    return IndexToProbability(maxCloudProbabilityIndex);
  }

}