Unverified Commit a3b0598b authored by Oleg Alexandrov's avatar Oleg Alexandrov Committed by GitHub
Browse files

Update SAR to the newest format (#389) 



* Tests passing

* Add new test

* Update SAR to the newest format

* Rm stray debug statement

* sarTests: Minor adjustment

* Let a test print its output on failure, as cannot figure out the problem

* Fix a test, not sure about the reason for the divergence

* Use gmttime as suggested by the Z timestamp, and not local time

* Try setting the daylight savings time bit

* Use daylight savings time

* Try avoiding local time altogether

* Modularize code

* Include ctime for Windows to build

* Must be std::ctime

Co-authored-by: default avatarJesse Mapel <jmapel@usgs.gov>
parent 6e8db600

include/usgscsm/Utilities.h

+4 −0
Original line number Diff line number Diff line
@@ -176,4 +176,8 @@ void applyRotationToQuatVec(ale::Rotation const& r, std::vector<double> & quater 
void applyRotationTranslationToXyzVec(ale::Rotation const& r, ale::Vec3d const& t,

                                      std::vector<double> & xyz);



// Convert ephemeris time, in seconds since January 1, 2000 12:00:00 AM GMT,

// to a calendar time string, such as 2000-01-01T00:16:40Z.

std::string ephemTimeToCalendarTime(double ephemTime);



#endif  // INCLUDE_USGSCSM_UTILITIES_H_

src/UsgsAstroFrameSensorModel.cpp

+2 −10
Original line number Diff line number Diff line
@@ -693,16 +693,8 @@ std::string UsgsAstroFrameSensorModel::getSensorMode() const { 
std::string UsgsAstroFrameSensorModel::getReferenceDateAndTime() const {

  MESSAGE_LOG("Accessing reference data and time");

  time_t ephemTime = m_ephemerisTime;

  struct tm t = {0};  // Initalize to all 0's

  t.tm_year = 100;    // This is year-1900, so 100 = 2000

  t.tm_mday = 1;

  time_t timeSinceEpoch = mktime(&t);

  time_t finalTime = ephemTime + timeSinceEpoch;

  char buffer[22];

  strftime(buffer, 22, "%Y-%m-%dT%H:%M:%SZ", localtime(&finalTime));

  buffer[21] = '\0';



  return buffer;



  return ephemTimeToCalendarTime(ephemTime);

}



std::string UsgsAstroFrameSensorModel::getModelState() const {

src/UsgsAstroLsSensorModel.cpp

+6 −16
Original line number Diff line number Diff line
@@ -1212,16 +1212,8 @@ std::string UsgsAstroLsSensorModel::getReferenceDateAndTime() const { 
  double relativeTime =

    UsgsAstroLsSensorModel::getImageTime(referencePointImage);

  time_t ephemTime = m_centerEphemerisTime + relativeTime;

  struct tm t = {0};  // Initalize to all 0's

  t.tm_year = 100;    // This is year-1900, so 100 = 2000

  t.tm_mday = 1;

  time_t timeSinceEpoch = mktime(&t);

  time_t finalTime = ephemTime + timeSinceEpoch;

  char buffer[22];

  strftime(buffer, 22, "%Y-%m-%dT%H:%M:%SZ", localtime(&finalTime));

  buffer[21] = '\0';



  return buffer;

  return ephemTimeToCalendarTime(ephemTime);

}



//***************************************************************************
@@ -2380,8 +2372,6 @@ std::string UsgsAstroLsSensorModel::constructStateFromIsd( 
      state["m_detectorSampleOrigin"].dump(),

      state["m_detectorLineOrigin"].dump())







  ale::Orientations j2000_to_sensor = ale::getInstrumentPointing(jsonIsd);

  ale::State interpInstState, rotatedInstState;

  std::vector<double> positions = {};
@@ -2481,8 +2471,8 @@ std::string UsgsAstroLsSensorModel::constructStateFromIsd( 
  MESSAGE_LOG("m_currentParameterValue: {}",

              state["m_currentParameterValue"].dump())



  // get radii

  // ALE operates in Km and we want m

  // Get radii

  // ALE operates in km and we want m

  state["m_minorAxis"] = ale::getSemiMinorRadius(jsonIsd) * 1000;

  state["m_majorAxis"] = ale::getSemiMajorRadius(jsonIsd) * 1000;

  MESSAGE_LOG(

src/UsgsAstroPushFrameSensorModel.cpp

+2 −10
Original line number Diff line number Diff line
@@ -1191,16 +1191,8 @@ std::string UsgsAstroPushFrameSensorModel::getReferenceDateAndTime() const { 
  double relativeTime =

      UsgsAstroPushFrameSensorModel::getImageTime(referencePointImage);

  time_t ephemTime = m_centerEphemerisTime + relativeTime;

  struct tm t = {0};  // Initalize to all 0's

  t.tm_year = 100;    // This is year-1900, so 100 = 2000

  t.tm_mday = 1;

  time_t timeSinceEpoch = mktime(&t);

  time_t finalTime = ephemTime + timeSinceEpoch;

  char buffer[22];

  strftime(buffer, 22, "%Y-%m-%dT%H:%M:%SZ", localtime(&finalTime));

  buffer[21] = '\0';



  return buffer;

  

  return ephemTimeToCalendarTime(ephemTime);

}



//***************************************************************************

src/UsgsAstroSarSensorModel.cpp

+138 −63
Original line number Diff line number Diff line
@@ -8,6 +8,8 @@ 


#include <nlohmann/json.hpp>



#include "ale/Util.h"



using json = nlohmann::json;

using namespace std;
@@ -44,64 +46,85 @@ string UsgsAstroSarSensorModel::constructStateFromIsd( 


  std::shared_ptr<csm::WarningList> parsingWarnings(new csm::WarningList);



  state["m_modelName"] = getSensorModelName(isd, parsingWarnings.get());

  state["m_imageIdentifier"] = getImageId(isd, parsingWarnings.get());

  state["m_sensorName"] = getSensorName(isd, parsingWarnings.get());

  state["m_platformName"] = getPlatformName(isd, parsingWarnings.get());



  state["m_nLines"] = getTotalLines(isd, parsingWarnings.get());

  state["m_nSamples"] = getTotalSamples(isd, parsingWarnings.get());



  // Zero computed state values

  state["m_referencePointXyz"] = vector<double>(3, 0.0);



  // sun_position and velocity are required for getIlluminationDirection

  state["m_sunPosition"] = getSunPositions(isd, parsingWarnings.get());

  state["m_sunVelocity"] = getSunVelocities(isd, parsingWarnings.get());

  state["m_modelName"] = ale::getSensorModelName(isd);

  state["m_imageIdentifier"] = ale::getImageId(isd);

  state["m_sensorName"] = ale::getSensorName(isd);

  state["m_platformName"] = ale::getPlatformName(isd);



  state["m_centerEphemerisTime"] = getCenterTime(isd, parsingWarnings.get());

  state["m_startingEphemerisTime"] = getStartingTime(isd, parsingWarnings.get());

  MESSAGE_LOG(

      "m_modelName: {} "

      "m_imageIdentifier: {} "

      "m_sensorName: {} "

      "m_platformName: {} ",

      state["m_modelName"].dump(), state["m_imageIdentifier"].dump(),

      state["m_sensorName"].dump(), state["m_platformName"].dump());

  

  state["m_nLines"] = ale::getTotalLines(isd);

  state["m_nSamples"] = ale::getTotalSamples(isd);

  MESSAGE_LOG(

      "m_nLines: {} "

      "m_nSamples: {} ",

      state["m_nLines"].dump(), state["m_nSamples"].dump());



  ale::States inst_state = ale::getInstrumentPosition(isd);

  std::vector<double> ephemTime = inst_state.getTimes();

  if (ephemTime.size() < 2) {

    parsingWarnings->push_back(csm::Warning(

        csm::Warning::DATA_NOT_AVAILABLE, "Instrument position has less than 2 states.",

        "UsgsAstroSarSensorModel::constructStateFromIsd()"));

    MESSAGE_LOG("m_dtEphem not in ISD")

  }

  else {

    // Re-compute times so that they are linearly spaced

    double startTime = ephemTime.front();

    double timeStep = (ephemTime.back() - ephemTime.front()) / (ephemTime.size() - 1);

    state["m_t0Ephem"] = startTime;

    state["m_dtEphem"] = timeStep;

    for (int i = 0; i < ephemTime.size(); i++) {

      ephemTime[i] = startTime + i * timeStep;

    }

  }



  // Make the Sun positions and velocities be in target body coordinates

  ale::States sunState = ale::getSunPosition(isd);

  ale::Orientations j2000_to_target = ale::getBodyRotation(isd);

  

  std::vector<double> sunPositions = {};

  std::vector<double> sunVelocities = {};

  for (int i = 0; i < ephemTime.size(); i++) {



    ale::State j2000SunState = sunState.getState(ephemTime[i], ale::SPLINE);

    ale::State rotatedSunState =

        j2000_to_target.rotateStateAt(ephemTime[i], j2000SunState, ale::SLERP);

    // ALE operates in Km and we want m

    sunPositions.push_back(rotatedSunState.position.x * 1000);

    sunPositions.push_back(rotatedSunState.position.y * 1000);

    sunPositions.push_back(rotatedSunState.position.z * 1000);

    sunVelocities.push_back(rotatedSunState.velocity.x * 1000);

    sunVelocities.push_back(rotatedSunState.velocity.y * 1000);

    sunVelocities.push_back(rotatedSunState.velocity.z * 1000);

  }

  state["m_sunPosition"] = sunPositions;

  state["m_sunVelocity"] = sunVelocities;



  state["m_startingEphemerisTime"] = ale::getStartingTime(isd);

  state["m_centerEphemerisTime"] = ale::getCenterTime(isd);

  state["m_endingEphemerisTime"] = getEndingTime(isd, parsingWarnings.get());

  MESSAGE_LOG("m_startingEphemerisTime: {} "

              "m_centerEphemerisTime: {} "

              "m_endingEphemerisTime: {} ",

              state["m_startingEphemerisTime"].dump(),

              state["m_centerEphemerisTime"].dump(),

              state["m_endingEphemerisTime"].dump()

              )



  state["m_exposureDuration"] = getExposureDuration(isd, parsingWarnings.get());



  try {

    state["m_dtEphem"] = isd.at("dt_ephemeris");

  } catch (...) {

    std::string msg = "dt_ephemeris not in ISD";

    MESSAGE_LOG(msg);

    parsingWarnings->push_back(

        csm::Warning(csm::Warning::DATA_NOT_AVAILABLE, msg,

                     "UsgsAstroSarSensorModel::constructStateFromIsd()"));

  }



  try {

    state["m_t0Ephem"] = isd.at("t0_ephemeris");

  } catch (...) {

    std::string msg = "t0_ephemeris not in ISD";

    MESSAGE_LOG(msg);

    parsingWarnings->push_back(

        csm::Warning(csm::Warning::DATA_NOT_AVAILABLE, msg,

                     "UsgsAstroSarSensorModel::constructStateFromIsd()"));

  }



  state["m_positions"] = getSensorPositions(isd, parsingWarnings.get());

  state["m_velocities"] = getSensorVelocities(isd, parsingWarnings.get());



  state["m_currentParameterValue"] = vector<double>(NUM_PARAMETERS, 0.0);



  // get radii

  state["m_minorAxis"] = getSemiMinorRadius(isd, parsingWarnings.get());

  state["m_majorAxis"] = getSemiMajorRadius(isd, parsingWarnings.get());



  // set identifiers

  state["m_platformIdentifier"] = getPlatformName(isd, parsingWarnings.get());

  state["m_sensorIdentifier"] = getSensorName(isd, parsingWarnings.get());



  // get reference_height

  state["m_minElevation"] = -1000;

  state["m_maxElevation"] = 1000;



  // SAR specific values

  state["m_scaledPixelWidth"] = getScaledPixelWidth(isd, parsingWarnings.get());

  state["m_scaleConversionCoefficients"] =
@@ -111,6 +134,52 @@ string UsgsAstroSarSensorModel::constructStateFromIsd( 
  state["m_wavelength"] = getWavelength(isd, parsingWarnings.get());

  state["m_lookDirection"] = getLookDirection(isd, parsingWarnings.get());



  // Process positions and velocities

  ale::State interpInstState, rotatedInstState;

  std::vector<double> instPositions = {};

  std::vector<double> instVelocities = {};

  for (int i = 0; i < ephemTime.size(); i++) {

    

    interpInstState = inst_state.getState(ephemTime[i], ale::SPLINE);

    rotatedInstState =

      j2000_to_target.rotateStateAt(ephemTime[i], interpInstState, ale::SLERP);

    

    // ALE operates in Km and we want m

    instPositions.push_back(rotatedInstState.position.x * 1000);

    instPositions.push_back(rotatedInstState.position.y * 1000);

    instPositions.push_back(rotatedInstState.position.z * 1000);

    instVelocities.push_back(rotatedInstState.velocity.x * 1000);

    instVelocities.push_back(rotatedInstState.velocity.y * 1000);

    instVelocities.push_back(rotatedInstState.velocity.z * 1000);



  }

  state["m_positions"] = instPositions;

  state["m_velocities"] = instVelocities;



  // Get radii

  // ALE operates in km and we want m

  state["m_minorAxis"] = ale::getSemiMinorRadius(isd) * 1000;

  state["m_majorAxis"] = ale::getSemiMajorRadius(isd) * 1000;

  MESSAGE_LOG("m_minorAxis: {}"

              "m_majorAxis: {}",

              state["m_minorAxis"].dump(), state["m_majorAxis"].dump());

  

  // set identifiers

  state["m_platformIdentifier"] = ale::getPlatformName(isd);

  state["m_sensorIdentifier"] = ale::getSensorName(isd);

  MESSAGE_LOG(

      "m_platformIdentifier: {}"

      "m_sensorIdentifier: {}",

      state["m_platformIdentifier"].dump(), state["m_sensorIdentifier"].dump());



  // Get reference_height

  state["m_minElevation"] = ale::getMinHeight(isd);

  state["m_maxElevation"] = ale::getMaxHeight(isd);

  MESSAGE_LOG(

      "m_minElevation: {}"

      "m_maxElevation: {}",

      state["m_minElevation"].dump(), state["m_maxElevation"].dump());



  // Default to identity covariance

  state["m_covariance"] = vector<double>(NUM_PARAMETERS * NUM_PARAMETERS, 0.0);

  for (int i = 0; i < NUM_PARAMETERS; i++) {
@@ -501,8 +570,11 @@ csm::EcefCoord UsgsAstroSarSensorModel::imageToGround( 
  double time =

    m_startingEphemerisTime + (imagePt.line - 0.5) * m_exposureDuration;

  double groundRange = (imagePt.samp - 0.5) * m_scaledPixelWidth;

  MESSAGE_LOG("Ground range: {}", groundRange);



  std::vector<double> coeffs = getRangeCoefficients(time);

  double slantRange = groundRangeToSlantRange(groundRange, coeffs);

  MESSAGE_LOG("Slant range: {}", slantRange);



  // Compute the in-track, cross-track, nadir, coordinate system to solve in

  csm::EcefVector spacecraftPosition = getSpacecraftPosition(time);
@@ -514,10 +586,14 @@ csm::EcefCoord UsgsAstroSarSensorModel::imageToGround( 
  csm::EcefVector tHat = normalized(rejection(spacecraftPosition, vHat));

  // Cross-track unit vector

  csm::EcefVector uHat = cross(vHat, tHat);

  MESSAGE_LOG("In track vector: {}, {}, {}", vHat.x, vHat.y, vHat.z);

  MESSAGE_LOG("Nadir vector: {}, {}, {}", tHat.x, tHat.y, tHat.z);

  MESSAGE_LOG("In track vector: {}, {}, {}", uHat.x, uHat.y, uHat.z);



  // Compute the spacecraft position in the new coordinate system

  //   The cross-track unit vector is orthogonal to the position so we ignore it

  double nadirComp = dot(spacecraftPosition, tHat);

  MESSAGE_LOG("Nadir position component: {}", nadirComp);

  

  // Iterate to find proper radius value. Limit the number of iterations

  // because the desired precision may not drop beyond 1e-10 no matter
@@ -533,10 +609,12 @@ csm::EcefCoord UsgsAstroSarSensorModel::imageToGround( 
    double alpha =

        (radiusSqr - slantRange * slantRange - positionMag * positionMag) /

        (2 * nadirComp);

    MESSAGE_LOG("Iteration alpha: {}", alpha);

    double beta = sqrt(slantRange * slantRange - alpha * alpha);

    if (m_lookDirection == LEFT) {

      beta *= -1;

    }

    MESSAGE_LOG("Iteration beta: {}", beta);

    groundVec = alpha * tHat + beta * uHat + spacecraftPosition;

    pointHeight = computeEllipsoidElevation(

        groundVec.x, groundVec.y, groundVec.z, m_majorAxis, m_minorAxis);
@@ -546,9 +624,10 @@ csm::EcefCoord UsgsAstroSarSensorModel::imageToGround( 
      break;

  }



  if (achievedPrecision) {

  if (achievedPrecision)

    *achievedPrecision = fabs(pointHeight - height);

  }



  MESSAGE_LOG("Iteration ground point: {}, {}, {}", groundVec.x, groundVec.y, groundVec.z);



  csm::EcefCoord groundPt(groundVec.x, groundVec.y, groundVec.z);
@@ -665,9 +744,8 @@ csm::EcefLocus UsgsAstroSarSensorModel::imageToProximateImagingLocus( 
  }

  tangent = normalized(tangent);



  if (achievedPrecision) {

  if (achievedPrecision)

    *achievedPrecision = 0.0;

  }



  return csm::EcefLocus(closestVec.x, closestVec.y, closestVec.z, tangent.x,

                        tangent.y, tangent.z);
@@ -779,6 +857,8 @@ csm::EcefCoord UsgsAstroSarSensorModel::getAdjustedSensorPosition( 
    spacecraftPosition.y = m_positions[1] + getValue(1, adj);

    spacecraftPosition.z = m_positions[2] + getValue(2, adj);

  }

  MESSAGE_LOG("Adjusted sensor positions is {}, {}, {}.",

      spacecraftPosition.x, spacecraftPosition.y, spacecraftPosition.z);

  return csm::EcefCoord(spacecraftPosition.x, spacecraftPosition.y,

                        spacecraftPosition.z);

}
@@ -817,6 +897,8 @@ csm::EcefVector UsgsAstroSarSensorModel::getAdjustedSensorVelocity( 
    spacecraftVelocity.y = m_velocities[1] + getValue(4, adj);

    spacecraftVelocity.z = m_velocities[2] + getValue(5, adj);

  }

  MESSAGE_LOG("Adjusted sensor velocity is {}, {}, {}",

      spacecraftVelocity.x, spacecraftVelocity.y, spacecraftVelocity.z);

  return spacecraftVelocity;

}
@@ -1042,18 +1124,11 @@ string UsgsAstroSarSensorModel::getSensorMode() const { return "STRIP"; } 


string UsgsAstroSarSensorModel::getReferenceDateAndTime() const {

  csm::ImageCoord referencePointImage = groundToImage(m_referencePointXyz);

  double relativeTime = getImageTime(referencePointImage);

  double relativeTime =

    UsgsAstroSarSensorModel::getImageTime(referencePointImage);

  time_t ephemTime = m_centerEphemerisTime + relativeTime;

  struct tm t = {0};  // Initalize to all 0's

  t.tm_year = 100;    // This is year-1900, so 100 = 2000

  t.tm_mday = 1;

  time_t timeSinceEpoch = mktime(&t);

  time_t finalTime = ephemTime + timeSinceEpoch;

  char buffer[22];

  strftime(buffer, 22, "%Y-%m-%dT%H:%M:%SZ", localtime(&finalTime));

  buffer[21] = '\0';



  return buffer;



  return ephemTimeToCalendarTime(ephemTime);

}



csm::Ellipsoid UsgsAstroSarSensorModel::getEllipsoid() const {