PlotGOES.cc

Go to the documentation of this file.

// Keep all further includes in the header
#include "PlotGOES.h"

using namespace PlotGOESNS;

PlotGOES::PlotGOES()
{}

PlotGOES::~PlotGOES()
{}

fwk::VModule::ResultFlag
PlotGOES::Init()
{
  const std::string bTopStr("PlotGOES");
  const utl::Branch bTop(fwk::CentralConfig::GetInstance()->GetTopBranch(bTopStr));

  if (!bTop) {
    FATAL("Could not find branch " + bTopStr);
    return eFailure;
  }

  bTop.GetChild("TimeStart").GetData(fTimeStart);
  bTop.GetChild("TimeFinish").GetData(fTimeFinish);
  fTimeStep = bTop.GetChild("TimeStep").Get<double>();

  bTop.GetChild("PlotFormats").GetData(fPlotFormats);

  {
    std::string str;
    bTop.GetChild("PlotFDFOV").GetData(str);
    if (str == "on")
      fPlotFDFOV = true;
    else if (str == "off")
      fPlotFDFOV = false;
    else {
      FATAL("Unrecognised setting \"" + str + "\" to <PlotFDFOV>");
      return eFailure;
    }
  }

  return eSuccess;
}

fwk::VModule::ResultFlag
PlotGOES::Run(evt::Event& /*event*/)
{
  // Prepare ROOT plotting
  //

  // For paper
  const int textFont(132); // Times new roman
  const double textSize(0.06);
  const bool showTitle(true);

  // For poster
  // const int textFont(42); // Helvetica
  // const double textSize(0.055);
  // const bool showTitle(false);

  gStyle->SetTitleFont(textFont, "");    // Hist title
  gStyle->SetTitleFont(textFont, "XYZ"); // Axis titles
  gStyle->SetLabelFont(textFont, "XYZ");

  const int canvasW(800);
  const int canvasH(600);

  // Set margins
  const double marginH(0.12); // 0.1 is default
  const double marginW((1 - (1 - 2*marginH)*canvasH/canvasW)/2); // Square plot
  gStyle->SetPadLeftMargin(marginW);
  gStyle->SetPadRightMargin(marginW);
  gStyle->SetPadTopMargin(marginH);
  gStyle->SetPadBottomMargin(marginH);

  // Palette
  {
    const int nStops(6); // Points between which to linearly interpolate
    const int nColours(255); // True depth of palette

    double stops[nStops] = {0.0, 0.2, 0.4, 0.6, 0.8, 1.0};
    double r[nStops]     = {1.0, 1.0, 0.9, 0.8, 0.7, 0.6};
    double g[nStops]     = {1.0, 1.0, 0.9, 0.8, 0.7, 0.6};
    double b[nStops]     = {1.0, 1.0, 0.9, 0.8, 0.7, 0.6};

    TColor::CreateGradientColorTable(nStops, stops, r, g, b, nColours);
  }

  // Constants
  //

  // The standard point on tile 19H which we tend to use as an origin when drawing the array
  const double augerEasting(440000.0);
  const double augerNorthing(6060000.0);

  // Use https://maps.auger.org.ar to select station IDs around the edge
  // This list starts at the top left (near Coihueco) and goes clockwise
  int nStations(35);
  unsigned int edgeIDs[nStations] = {
    1047, 1041, 1051,  979,  985, 1549, 1506, 1483, 1501, 1131,
    1133, 1145, 1148, 1666, 1668, 1676,  790, 1175,  768,  496,
     572,  531,  462,  476,  321,  379,  388,  377, 1260, 1247,
     945,  618, 1314,  615,  787
  };

  det::Detector& det( det::Detector::GetInstance() );
  const sdet::SDetector& sdet( det.GetSDetector() );
  const fdet::FDetector& fdet( det.GetFDetector() );
  static const utl::ReferenceEllipsoid& wgs84( utl::ReferenceEllipsoid::GetWGS84() );

  // Choose a time when all the stations in the above list exist
  det.Update( utl::UTCDateTime(2020, 1, 1).GetTimeStamp() );

  // Fill a graph with the locations of the stations, wrapped back around
  TGraph sdArray;
  double pair[2];
  for (int iSDID(0); iSDID < nStations; ++iSDID) {
    utl::Point p( sdet.GetStation(edgeIDs[iSDID]).GetPosition() );
    utl::UTMPoint utmP(p, utl::ReferenceEllipsoid::GetWGS84());

    const double utmX((utmP.GetEasting() - augerEasting)/utl::km);
    const double utmY((utmP.GetNorthing() - augerNorthing)/utl::km);

    if (iSDID == 0) {
      pair[0] = utmX;
      pair[1] = utmY;
    }

    sdArray.SetPoint(sdArray.GetN(), utmX, utmY);
  }
  sdArray.SetPoint(sdArray.GetN(), pair[0], pair[1]);

  // Fill a vector with polygons showing the outline of the FD telescope fields of view
  std::vector<TPolyLine> telFOVs;
  for (auto itEye( fdet.EyesBegin() ); itEye != fdet.EyesEnd(); ++itEye) {
    // Skip HEAT, virtual eyes, etc
    if (itEye->GetId() > 4)
      continue;

    const utl::CoordinateSystemPtr eyeCS( itEye->GetEyeCoordinateSystem() );

    /* Because each telescope has an independent pointing which isn't guaranteed to be
      exactly 30 degrees apart from the previous, we need to smooth things over a little bit by
      forcing each telescope FOV after the first one to overlap with the previous one. The
      consequence of not doing this is that some borders between telescopes clearly have two
      lines that are very close together, which clutters the plot */
    bool usePrevPoint(false);
    utl::UTMPoint prevC(0.0, 0.0, 0.0, wgs84); // unset

    for (auto itTel( itEye->TelescopesBegin() ); itTel != itEye->TelescopesEnd(); ++itTel) {
      const utl::Vector& telAxis( itTel->GetAxis() );
      const utl::Point& telBase( itTel->GetPosition() );

      const utl::UTMPoint ptB(telBase, wgs84);

      // Positive and negative rotations about a vertical Z axis (think right-hand rule)
      utl::TransformationPolicy rotP(utl::TransformationMatrix::RotationZ(15*utl::degree), eyeCS);
      utl::TransformationPolicy rotN(utl::TransformationMatrix::RotationZ(-15*utl::degree), eyeCS);

      const utl::Vector ptBtoPtA(rotN*telAxis);
      utl::UTMPoint ptA(telBase + 20*utl::km*ptBtoPtA, wgs84);

      const utl::Vector ptBtoPtC(rotP*telAxis);
      const utl::UTMPoint ptC(telBase + 20*utl::km*ptBtoPtC, wgs84);

      // Swap out point A for the previous point C, if this isn't the first telescope for an eye
      if (usePrevPoint)
        ptA = prevC;
      else
        usePrevPoint = true;
      prevC = ptC;

      const int v(3); // 3 vertices for an open triangle

      double xPts[v] = {
        ptA.GetEasting() ,
        ptB.GetEasting() ,
        ptC.GetEasting()
      };

      double yPts[v] = {
        ptA.GetNorthing() ,
        ptB.GetNorthing() ,
        ptC.GetNorthing()
      };

      // Reference to relative origin; convert to km
      for (int i(0); i < v; ++i) {
        xPts[i] = (xPts[i] - augerEasting)/utl::km;
        yPts[i] = (yPts[i] - augerNorthing)/utl::km;
      } // i

      telFOVs.emplace_back(v, xPts, yPts);
    } // itTel
  } // itEye

  //

  // Loop over times
  for (utl::TimeStamp detTime(fTimeStart); detTime <= fTimeFinish; detTime += fTimeStep) {
    det.Update(detTime);
    const atm::GOESDB& goesdb( det.GetAtmosphere().GetGOESDB() );

    // Always get these again, as they may become time-dependent in future
    const double pixW( goesdb.GetPixelWidthEasting() );
    const double pixH( goesdb.GetPixelWidthNorthing() );
    const unsigned int nPix( goesdb.GetNumberOfPixels() );

    TH2Poly hist;
    bool noData(false);

    // This technique is to get all the cloud probabilities at once, because we intend to plot all of them. If your application does not need to get *all* the cloud probabilities, then you should probably loop over your desired pixels and get them individually

    atm::GOESDB::ProbabilityMap pixMap;
    try {
      pixMap = goesdb.GetAllCloudProbabilities();
    }
    catch (utl::DataNotFoundInDBException&) {
      continue;
    }

    for (unsigned int iPix(0); iPix < nPix; ++iPix) {
      double prob(0.0);
      try {
        // Cloud prob is 0, 0.2, 0.4, 0.6, or 0.8! These are lower bounds
        prob = pixMap.at(iPix);

        // Add 0.1 so the val is in the middle of the 20%-wide region it represents
        prob += 0.1;
      }
      catch (std::out_of_range&) {
        noData = true;
        break;
      }

      const utl::UTMPoint& pix( goesdb.GetPixelCenter(iPix) );
      const double pixX( pix.GetEasting() );
      const double pixY( pix.GetNorthing() );

      const int v(5); // 5 vertices for a closed rectangle

      double xPts[v] = {
        pixX - pixW / 2.0 ,
        pixX + pixW / 2.0 ,
        pixX + pixW / 2.0 ,
        pixX - pixW / 2.0 ,
        pixX - pixW / 2.0
      };

      double yPts[v] = {
        pixY + pixH / 2.0 ,
        pixY + pixH / 2.0 ,
        pixY - pixH / 2.0 ,
        pixY - pixH / 2.0 ,
        pixY + pixH / 2.0
      };

      // Reference to relative origin; convert to km
      for (int i(0); i < v; ++i) {
        xPts[i] = (xPts[i] - augerEasting)/utl::km;
        yPts[i] = (yPts[i] - augerNorthing)/utl::km;
      } // i

      // Add the bin and set its value
      hist.SetBinContent(hist.AddBin(v, xPts, yPts), prob);
    } // iPix

    if (noData) {
      std::ostringstream status;
      status << "Missing pixel data at " << detTime;
      WARNING(status);
      continue;
    }

    // Add extra points outside of range, so we can SetRangeUser() further out
    hist.AddBin(-10, -10, -10, -10);
    hist.AddBin(85, 85, 85, 85);

    // Plotting
    {
      TCanvas canv("canv", "", canvasW, canvasH);
      canv.SetCanvasSize(canvasW, canvasH);

      {
        std::ostringstream conv;
        conv << utl::UTCDateTime(detTime);
        hist.SetTitle( conv.str().c_str() );
      }

      hist.SetStats(false);
      hist.SetTitleSize(1.1*textSize);

      if (!showTitle)
        hist.SetTitle("");

      TAxis* histX( hist.GetXaxis() );
      TAxis* histY( hist.GetYaxis() );
      TAxis* histZ( hist.GetZaxis() );

      histX->SetTitle("Relative easting (km)");
      histX->SetTitleOffset(0.95);
      histX->CenterTitle();
      histX->SetTitleSize(textSize);
      histX->SetLabelSize(textSize);
      histX->SetRangeUser(-5, 75);

      histY->SetTitle("Relative northing (km)");
      histY->SetTitleOffset(0.95);
      histY->CenterTitle();
      histY->SetTitleSize(textSize);
      histY->SetLabelSize(textSize);
      histY->SetRangeUser(0, 80);

      histZ->SetTitle("Cloud probability");
      histZ->SetTitleOffset(0.95);
      histZ->CenterTitle();
      histZ->SetTitleSize(textSize);
      histZ->SetLabelSize(textSize);
      histZ->RotateTitle();
      histZ->SetRangeUser(0, 1);
      histZ->SetNdivisions(1005);

      // There are only 5 levels, so clip to those
      hist.SetContour(5);

      hist.Draw("colz");

      sdArray.SetLineColor(1);
      sdArray.SetLineWidth(2);
      sdArray.Draw("same l");

      if (fPlotFDFOV) {
        for (auto itFOV( telFOVs.begin() ); itFOV != telFOVs.end(); ++itFOV) {
          itFOV->SetLineColor(1);
          itFOV->SetLineStyle(2);
          itFOV->Draw();
        } // itFOV
      }

      for (auto itF(fPlotFormats.begin()); itF != fPlotFormats.end(); ++itF) {
        std::ostringstream conv;
        conv << detTime.GetGPSSecond() << "." << *itF;
        canv.SaveAs( conv.str().c_str() );
      } // itF
    }
  }

  return eSuccess;
}

fwk::VModule::ResultFlag
PlotGOES::Finish()
{
  return eSuccess;
}