LidarCloudDBModel.cc

Go to the documentation of this file.

#include <det/Detector.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>
#include <fdet/Pixel.h>

#include <atm/CloudResult.h>
#include <atm/LidarDB.h>
#include <atm/LidarZone.h>
#include <atm/LidarCloudDBModel.h>

#include <fwk/CoordinateSystemRegistry.h>

#include <utl/Point.h>
#include <utl/UTMPoint.h>
#include <utl/ErrorLogger.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/Vector.h>

#include <iostream>
#include <string>
#include <sstream>
#include <limits>

using namespace det;
using namespace fdet;
using namespace atm;
using namespace fwk;
using namespace utl;
using namespace std;

LidarCloudDBModel::LidarCloudDBModel() :
  VCloudModel(),
  fCloudHeightMap(0),
  fLidarDbIsEmptyNow(false)
{
}

LidarCloudDBModel::~LidarCloudDBModel()
{
  delete fCloudHeightMap;
}

void
LidarCloudDBModel::Init()
{
}

bool
LidarCloudDBModel::HasData()
  const
{
  return true;
  // return CheckForUpdates();
}

CloudResult
LidarCloudDBModel::EvaluateCloudCoverage(const unsigned int eyeId, const unsigned int telescopeId, const unsigned int pixelId, const utl::Point& x)
  const
{
  bool hasLidarData;

  // If we have lidar data, get the cloud base height from the nearest FD zone
  if ((hasLidarData = CheckForUpdates())) {

    const ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));
    const double pointHeight = (wgs84.PointToLatitudeLongitudeHeight(x)).get<2>();

    const string lidarZoneName = GetZoneName(x);
    const double cloudHeight = ((*fCloudHeightMap)[lidarZoneName]);

//#warning fix negative height in LIDAR DB
    // -999 in LIDAR DB for no clouds detected...
    if ( cloudHeight < -998 )
      return CloudResult(0, hasLidarData);

    // If the Point altitude is below the lowest cloud base altitude, then
    // treat the pixel as unobscured regardless of the cloud camera data
    if (pointHeight < cloudHeight)
      return CloudResult(0, hasLidarData);
  }

  // If we do not have lidar data, or if the cloud base is below the point
  // altitude, check the IR camera database for coverage data
  if (det::Detector::GetInstance().GetFDetector().
      GetEye(eyeId).GetTelescope(telescopeId).GetPixel(pixelId).HasCloudFraction()) {
    float cloudFraction = det::Detector::GetInstance().GetFDetector().
      GetEye(eyeId).GetTelescope(telescopeId).GetPixel(pixelId).GetCloudFraction();

    return CloudResult(cloudFraction, hasLidarData);
  }

  // No data available: return default CloudResult
  else {
    DEBUGLOG("no cloud data in pixel");
    return CloudResult();
  }
}

// Check for cached lidar data; return true if found.
bool
LidarCloudDBModel::CheckForUpdates()
  const
{
  // Minimize queries to the lidar database.
  if (fCurrentTime == det::Detector::GetInstance().GetTime() && fCurrentTime)
    return !fLidarDbIsEmptyNow;

  fCurrentTime = det::Detector::GetInstance().GetTime();
  fLidarDbIsEmptyNow = false;

  delete fCloudHeightMap;
  fCloudHeightMap = new map<string, double>;

  fZonePositions.clear();

  const LidarDB& lidarDB =
    det::Detector::GetInstance().GetAtmosphere().GetLidarDB();

  unsigned int zoneCount = 0;
  for (LidarDB::ZoneIterator zIt = lidarDB.ZonesBegin();
       zIt != lidarDB.ZonesEnd(); ++zIt, ++zoneCount)
  {
    // Fill zone name and position.
    fZonePositions[zIt->GetName()] =
      UTMPoint(zIt->GetNorthing(), zIt->GetEasting(),
               0., 19, 'H',
               ReferenceEllipsoid::GetEllipsoidIDFromString("WGS84")).GetPoint();
    // Get lidar cloud base data.
    (*fCloudHeightMap)[zIt->GetName()] = zIt->GetLowestCloudHeight();
  }

  // Mindless copy of MeasuredDBMieModel:
  // If there are no data, an exception should have been thrown by the
  // low-level interface (no?); so reaching this point would be a bug.
  if (zoneCount == 0) {
    fLidarDbIsEmptyNow = true;

    ostringstream msg;
    msg << "Cloud model requested data from lidar database, but none were found"
        << " at " << det::Detector::GetInstance().GetTime();

    DEBUGLOG(msg);
    return false;
  }
  return true;
}

// Get the LidarZone closest to the current Eye, identified by the Eye number.
// NOTE: this is just for test; it's probably not what we want.
// Instead, it's probably best to pick the LidarZone nearest the viewing point.
string
LidarCloudDBModel::GetZoneName(const unsigned int eyeId)
  const
{
  const utl::Point& eyePos =
    det::Detector::GetInstance().GetFDetector().GetEye(eyeId).GetPosition();

  string zoneName("");
  double smallestDist = numeric_limits<float>::max(); // was HUGE

  for (map<string, Point>::const_iterator zIt = fZonePositions.begin();
       zIt != fZonePositions.end(); ++zIt)
  {
    const double dist = (eyePos - zIt->second).GetMag();

    if (dist <= smallestDist) {
      smallestDist = dist;
      zoneName = zIt->first;
    }
  }

  return zoneName;
}

// Get the LidarZone nearest to the specified point.
string
LidarCloudDBModel::GetZoneName(const utl::Point& x)
  const
{
  string zoneName("");
  double smallestDist = numeric_limits<float>::max(); // was HUGE

  for (map<string, Point>::const_iterator zIt = fZonePositions.begin();
       zIt != fZonePositions.end(); ++zIt)
  {
    const double dist = (x - zIt->second).GetMag();

    if (dist <= smallestDist) {
      smallestDist = dist;
      zoneName = zIt->first;
    }
  }

  return zoneName;
}