diff --git a/Code/Radiometry/otbAtmosphericCorrectionParameters2.cxx b/Code/Radiometry/otbAtmosphericCorrectionParameters2.cxx new file mode 100644 index 0000000000000000000000000000000000000000..a0d49ede6d7a428940769316679403a96b2269d3 --- /dev/null +++ b/Code/Radiometry/otbAtmosphericCorrectionParameters2.cxx @@ -0,0 +1,100 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ + +#include "otbAtmosphericCorrectionParameters2.h" + +#include <fstream> + +#include "otbAeronetFileReader.h" +#include "otbSpectralSensitivityReader.h" +#include "otbAeronetData.h" + +namespace otb +{ + + +AtmosphericCorrectionParameters2 +::AtmosphericCorrectionParameters2() +{ + /*m_SolarZenithalAngle = 361.; CHRIS + m_SolarAzimutalAngle = 361.; + m_ViewingZenithalAngle = 361.; + m_ViewingAzimutalAngle = 361.; + m_Month = 0; + m_Day = 0;*/ + m_AtmosphericPressure = 1030.; + m_WaterVaporAmount = 2.5; + m_OzoneAmount = 0.28; + m_AerosolModel = CONTINENTAL; + m_AerosolOptical = 0.2; + + /*m_WavelengthSpectralBand = InternalWavelengthSpectralBandVectorType::New(); CHRIS + m_WavelengthSpectralBand->Clear();*/ + +} + +/** Get data from aeronet file*/ +void +AtmosphericCorrectionParameters2 +::UpdateAeronetData(const std::string& file, int year, int month, int day, int hour, int minute, double epsi) +{ + if (file == "") itkExceptionMacro(<< "No Aeronet filename specified."); + + AeronetFileReader::Pointer reader = AeronetFileReader::New(); + reader->SetFileName(file); + reader->SetDay(day); + reader->SetMonth(month); + reader->SetYear(year); + reader->SetHour(hour); + reader->SetMinute(minute); + reader->SetEpsilon(epsi); + + reader->Update(); + + m_AerosolOptical = reader->GetOutput()->GetAerosolOpticalThickness(); + m_WaterVaporAmount = reader->GetOutput()->GetWater(); +} + + + +/**PrintSelf method */ +void +AtmosphericCorrectionParameters2 +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + /*os << "Solar zenithal angle : " << m_SolarZenithalAngle << std::endl; CHRIS + os << "Solar azimutal angle : " << m_SolarAzimutalAngle << std::endl; + os << "Viewing zenithal angle: " << m_ViewingZenithalAngle << std::endl; + os << "Viewing azimutal angle: " << m_ViewingAzimutalAngle << std::endl; + os << "Month : " << m_Month << std::endl; + os << "Day : " << m_Day << std::endl;*/ + os << "Atmospheric pressure : " << m_AtmosphericPressure << std::endl; + os << "Water vapor amount : " << m_WaterVaporAmount << std::endl; + os << "Ozone amount : " << m_OzoneAmount << std::endl; + os << "Aerosol model : " << m_AerosolModel << std::endl; + os << "Aerosol optical : " << m_AerosolOptical << std::endl; + + // Function values print : + /*os << "Filter function values: " << std::endl; CHRIS + for (unsigned int i = 0; i < m_WavelengthSpectralBand->Size(); ++i) + { + os << indent << "Channel " << i + 1 << " : " << std::endl; + os << indent << m_WavelengthSpectralBand->GetNthElement(i) << std::endl; + }*/ +} +} // end namespace otb diff --git a/Code/Radiometry/otbAtmosphericCorrectionParameters2.h b/Code/Radiometry/otbAtmosphericCorrectionParameters2.h index 011984139b8731eb14d6d30f920172a4c1474e14..ead49029d1f84caf8c74e3a3c96c9be1dbc6c162 100644 --- a/Code/Radiometry/otbAtmosphericCorrectionParameters2.h +++ b/Code/Radiometry/otbAtmosphericCorrectionParameters2.h @@ -101,14 +101,14 @@ public: { this->UpdateAeronetData(file, year, month, day, hour, minute, 0.4); } - void UpdateAeronetData(const std::string& file, int year, int hour, int minute, double epsi) + /*void UpdateAeronetData(const std::string& file, int year, int hour, int minute, double epsi) CHRIS { this->UpdateAeronetData(file, year, m_Month, m_Day, hour, minute, epsi); } void UpdateAeronetData(const std::string& file, int year, int hour, int minute) { this->UpdateAeronetData(file, year, m_Month, m_Day, hour, minute, 0.4); - } + }*/ /** Constructor */ AtmosphericCorrectionParameters2(); diff --git a/Code/Radiometry/otbImageMetadataCorrectionParameters.h b/Code/Radiometry/otbImageMetadataCorrectionParameters.h index 9041120f50a5e390bd056a92cb1e2ebcb375b72d..2cc1ec8c1971dc18e193b3e45013ec697adfcde5 100644 --- a/Code/Radiometry/otbImageMetadataCorrectionParameters.h +++ b/Code/Radiometry/otbImageMetadataCorrectionParameters.h @@ -38,7 +38,7 @@ namespace otb * */ -class ITK_EXPORT ImageMetadataCorrectionParameters : public itk::Object +class ITK_EXPORT ImageMetadataCorrectionParameters : public itk::DataObject { public: /** Standard typedefs */ diff --git a/Code/Radiometry/otbRadiometryCorrectionParametersToAtmosphericRadiativeTerms.h b/Code/Radiometry/otbRadiometryCorrectionParametersToAtmosphericRadiativeTerms.h index 3ba45922f405a6cb1953184f2b99fd545ee4d8fa..d71626d946f146833539f4844deb0d35f7addc02 100644 --- a/Code/Radiometry/otbRadiometryCorrectionParametersToAtmosphericRadiativeTerms.h +++ b/Code/Radiometry/otbRadiometryCorrectionParametersToAtmosphericRadiativeTerms.h @@ -19,8 +19,9 @@ #define __otbRadiometryCorrectionParametersToAtmosphericRadiativeTerms_h #include "otbAtmosphericRadiativeTerms2.h" -#include "otbAtmosphericCorrectionParameters.h" -#include "otbSIXSTraits.h" +#include "otbAtmosphericCorrectionParameters2.h" +#include "otbImageMetadataCorrectionParameters.h" +#include "otbSIXSTraits2.h" namespace otb { @@ -37,12 +38,12 @@ namespace otb public: /** Call the varSol function*/ - static AtmosphericRadiativeTerms2::Pointer Compute(AtmosphericCorrectionParameters* paramIn) + static AtmosphericRadiativeTerms2::Pointer Compute(AtmosphericCorrectionParameters2* paramAtmo, ImageMetadataCorrectionParameters* paramAcqui) { AtmosphericRadiativeTerms2::Pointer radTermsOut = AtmosphericRadiativeTerms2::New(); - typedef AtmosphericCorrectionParameters::WavelengthSpectralBandVectorType WavelengthSpectralBandVectorType; - WavelengthSpectralBandVectorType WavelengthSpectralBandVector = paramIn->GetWavelengthSpectralBand(); + typedef AtmosphericCorrectionParameters2::WavelengthSpectralBandVectorType WavelengthSpectralBandVectorType; + WavelengthSpectralBandVectorType WavelengthSpectralBandVector = paramAcqui->GetWavelengthSpectralBand(); unsigned int NbBand = WavelengthSpectralBandVector->Size(); radTermsOut->ValuesInitialization(NbBand); @@ -68,19 +69,19 @@ namespace otb upwardDirectTransmittance = 0.; upwardDiffuseTransmittanceForRayleigh = 0.; upwardDiffuseTransmittanceForAerosol = 0.; - SIXSTraits::ComputeAtmosphericParameters( - paramIn->GetSolarZenithalAngle(), /** The Solar zenithal angle */ - paramIn->GetSolarAzimutalAngle(), /** The Solar azimutal angle */ - paramIn->GetViewingZenithalAngle(), /** The Viewing zenithal angle */ - paramIn->GetViewingAzimutalAngle(), /** The Viewing azimutal angle */ - paramIn->GetMonth(), /** The Month */ - paramIn->GetDay(), /** The Day (in the month) */ - paramIn->GetAtmosphericPressure(), /** The Atmospheric pressure */ - paramIn->GetWaterVaporAmount(), /** The Water vapor amount (Total water vapor content over vertical atmospheric column) */ - paramIn->GetOzoneAmount(), /** The Ozone amount (Stratospheric ozone layer content) */ - paramIn->GetAerosolModel(), /** The Aerosol model */ - paramIn->GetAerosolOptical(), /** The Aerosol optical (radiative impact of aerosol for the reference wavelength 550-nm) */ - paramIn->GetWavelengthSpectralBand()->GetNthElement(i), /** Wavelength for the spectral band definition */ + SIXSTraits2::ComputeAtmosphericParameters( + paramAcqui->GetSolarZenithalAngle(), /** The Solar zenithal angle */ + paramAcqui->GetSolarAzimutalAngle(), /** The Solar azimutal angle */ + paramAcqui->GetViewingZenithalAngle(), /** The Viewing zenithal angle */ + paramAcqui->GetViewingAzimutalAngle(), /** The Viewing azimutal angle */ + paramAcqui->GetMonth(), /** The Month */ + paramAcqui->GetDay(), /** The Day (in the month) */ + paramAtmo->GetAtmosphericPressure(), /** The Atmospheric pressure */ + paramAtmo->GetWaterVaporAmount(), /** The Water vapor amount (Total water vapor content over vertical atmospheric column) */ + paramAtmo->GetOzoneAmount(), /** The Ozone amount (Stratospheric ozone layer content) */ + paramAtmo->GetAerosolModel(), /** The Aerosol model */ + paramAtmo->GetAerosolOptical(), /** The Aerosol optical (radiative impact of aerosol for the reference wavelength 550-nm) */ + paramAcqui->GetWavelengthSpectralBand()->GetNthElement(i), /** Wavelength for the spectral band definition */ /** Note : The Max wavelength spectral band value must be updated ! */ atmosphericReflectance, /** Atmospheric reflectance */ atmosphericSphericalAlbedo, /** atmospheric spherical albedo */ @@ -102,7 +103,7 @@ namespace otb radTermsOut->SetUpwardDirectTransmittance(i, upwardDirectTransmittance); radTermsOut->SetUpwardDiffuseTransmittanceForRayleigh(i, upwardDiffuseTransmittanceForRayleigh); radTermsOut->SetUpwardDiffuseTransmittanceForAerosol(i, upwardDiffuseTransmittanceForAerosol); - radTermsOut->SetWavelengthSpectralBand(i, paramIn->GetWavelengthSpectralBand()->GetNthElement(i)->GetCenterSpectralValue()); + radTermsOut->SetWavelengthSpectralBand(i, paramAcqui->GetWavelengthSpectralBand()->GetNthElement(i)->GetCenterSpectralValue()); } return radTermsOut; diff --git a/Code/Radiometry/otbSIXSTraits2.cxx b/Code/Radiometry/otbSIXSTraits2.cxx new file mode 100644 index 0000000000000000000000000000000000000000..d8e052615c298bfb8acecf79c891825110d46bb6 --- /dev/null +++ b/Code/Radiometry/otbSIXSTraits2.cxx @@ -0,0 +1,257 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: 2008-05-30 + Version: 2.2.0 + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#include "otbSIXSTraits2.h" + +#include "otb_6S.h" +#include "main_6s.h" +#include "otbMacro.h" + +#include <iomanip> + +namespace otb +{ + +void +SIXSTraits2::ComputeAtmosphericParameters( + const double SolarZenithalAngle, /** The Solar zenithal angle */ + const double SolarAzimutalAngle, /** The Solar azimutal angle */ + const double ViewingZenithalAngle, /** The Viewing zenithal angle */ + const double ViewingAzimutalAngle, /** The Viewing azimutal angle */ + const unsigned int Month, /** The Month */ + const unsigned int Day, /** The Day (in the month) */ + const double AtmosphericPressure, /** The Atmospheric pressure */ + const double WaterVaporAmount, /** The Water vapor amount (Total water vapor content over vertical atmospheric column) */ + const double OzoneAmount, /** The Ozone amount (Stratospheric ozone layer content) */ + const AerosolModelType& AerosolModel, /** The Aerosol model */ + const double AerosolOptical, /** The Aerosol optical (radiative impact of aerosol for the reference wavelength 550-nm) */ + WavelengthSpectralType* WavelengthSpectralBand, /** Wavelength for the spectral band definition */ + /** Note : The Max wavelength spectral band value must be updated ! */ + double& AtmosphericReflectance, /** Atmospheric reflectance */ + double& AtmosphericSphericalAlbedo, /** atmospheric spherical albedo */ + double& TotalGaseousTransmission, /** Total gaseous transmission */ + double& DownwardTransmittance, /** downward transmittance */ + double& UpwardTransmittance, /** upward transmittance */ + double& UpwardDiffuseTransmittance, /** upward diffuse transmittance */ + double& UpwardDirectTransmittance, /** Upward direct transmittance */ + double& UpwardDiffuseTransmittanceForRayleigh, /** upward diffuse transmittance for rayleigh */ + double& UpwardDiffuseTransmittanceForAerosol /** supward diffuse transmittance for aerosols */ + ) +{ +// geometrical conditions + otb_6s_doublereal asol(static_cast<otb_6s_doublereal>(SolarZenithalAngle)); + otb_6s_doublereal phi0(static_cast<otb_6s_doublereal>(SolarAzimutalAngle)); + otb_6s_doublereal avis(static_cast<otb_6s_doublereal>(ViewingZenithalAngle)); + otb_6s_doublereal phiv(static_cast<otb_6s_doublereal>(ViewingAzimutalAngle)); + otb_6s_integer month(static_cast<otb_6s_integer>(Month)); + otb_6s_integer jday(static_cast<otb_6s_integer>(Day)); + otb_6s_doublereal pressure(static_cast<otb_6s_doublereal>(AtmosphericPressure)); + otb_6s_doublereal uw(static_cast<otb_6s_doublereal>(WaterVaporAmount)); + otb_6s_doublereal uo3(static_cast<otb_6s_doublereal>(OzoneAmount)); +// atmospheric model + otb_6s_integer iaer(static_cast<otb_6s_integer>(AerosolModel)); + otb_6s_doublereal taer55(static_cast<otb_6s_doublereal>(AerosolOptical)); + + // Init output parameters + AtmosphericReflectance = 0.; + AtmosphericSphericalAlbedo = 0.; + TotalGaseousTransmission = 0.; + DownwardTransmittance = 0.; + UpwardTransmittance = 0.; + UpwardDiffuseTransmittance = 0.; + UpwardDirectTransmittance = 0.; + UpwardDiffuseTransmittanceForRayleigh = 0.; + UpwardDiffuseTransmittanceForAerosol = 0.; + + otb_6s_doublereal wlinf(0.), wlsup(0.); + otb_6s_doublereal otb_ratm__(0.), sast(0.), tgasm(0.), sdtott(0.), sutott(0.); + otb_6s_doublereal tdif_up(0.), tdir_up(0.), tdif_up_ray(0.), tdif_up_aer(0.); + + // 6S official Wavelength Spectral Band step value + const float SIXSStepOfWavelengthSpectralBandValues = .0025; + // Generate 6s Wavelength Spectral Band with the offcicial step value + ComputeWavelengthSpectralBandValuesFor6S(SIXSStepOfWavelengthSpectralBandValues, + WavelengthSpectralBand // Update + ); + try + { + + // 6S official tab size Wavelength Spectral + const otb_6s_integer S_6S_SIZE = 1501; + // Generate WavelengthSpectralBand in 6S compatible buffer s[1501] + wlinf = static_cast<otb_6s_doublereal>(WavelengthSpectralBand->GetMinSpectralValue()); + wlsup = static_cast<otb_6s_doublereal>(WavelengthSpectralBand->GetMaxSpectralValue()); + otb_6s_integer iinf = + static_cast<otb_6s_integer>((wlinf - (float) .25) / SIXSStepOfWavelengthSpectralBandValues + (float) 1.5); + otb_6s_integer cpt = iinf - 1; + otb_6s_doublereal * s(NULL); + s = new otb_6s_doublereal[S_6S_SIZE]; + memset(s, 0, S_6S_SIZE * sizeof(otb_6s_doublereal)); + const ValuesVectorType& FilterFunctionValues6S = WavelengthSpectralBand->GetFilterFunctionValues6S(); + // Set the values of FilterFunctionValues6S in s between [iinf-1; isup] + for (unsigned int i = 0; i < FilterFunctionValues6S.size() && cpt < S_6S_SIZE; ++i) + { + s[cpt] = FilterFunctionValues6S[i]; + ++cpt; + } + + // Call 6s main function + otbMsgDevMacro(<< "Start call 6S main function ..."); + otb_6s_ssssss_otb_main_function(&asol, &phi0, &avis, &phiv, &month, &jday, + &pressure, &uw, &uo3, + &iaer, + &taer55, + &wlinf, &wlsup, + s, + &otb_ratm__, + &sast, + &tgasm, + &sdtott, + &sutott, + &tdif_up, + &tdir_up, + &tdif_up_ray, + &tdif_up_aer); + otbMsgDevMacro(<< "Done call 6S main function!"); + delete[] s; + s = NULL; + } + catch (std::bad_alloc& err) + { + itkGenericExceptionMacro(<< "Exception bad_alloc in SIXSTraits2 class: " << (char*) err.what()); + } + catch (...) + { + itkGenericExceptionMacro(<< "Unknown exception in SIXSTraits2 class (catch(...)"); + } + + // Set outputs parameters + AtmosphericReflectance = static_cast<double>(otb_ratm__); + AtmosphericSphericalAlbedo = static_cast<double>(sast); + TotalGaseousTransmission = static_cast<double>(tgasm); + DownwardTransmittance = static_cast<double>(sdtott); + UpwardTransmittance = static_cast<double>(sutott); + UpwardDiffuseTransmittance = static_cast<double>(tdif_up); + UpwardDirectTransmittance = static_cast<double>(tdir_up); + UpwardDiffuseTransmittanceForRayleigh = static_cast<double>(tdif_up_ray); + UpwardDiffuseTransmittanceForAerosol = static_cast<double>(tdif_up_aer); +} + +void +SIXSTraits2::ComputeWavelengthSpectralBandValuesFor6S( + const double SIXSStepOfWavelengthSpectralBandValues, + WavelengthSpectralType* WavelengthSpectralBand + ) +{ + const double epsilon(.000001); + const double L_min = static_cast<double>(WavelengthSpectralBand->GetMinSpectralValue()); + const double L_max = static_cast<double>(WavelengthSpectralBand->GetMaxSpectralValue()); + const double L_userStep = static_cast<double>(WavelengthSpectralBand->GetUserStep()); + const ValuesVectorType& FilterFunctionValues = WavelengthSpectralBand->GetFilterFunctionValues(); + unsigned int i = 1; + unsigned int j = 1; + const double invStep = static_cast<double>(1. / L_userStep); + double value(0.); + + if (FilterFunctionValues.size() <= 1) + { + itkGenericExceptionMacro(<< "The FilterFunctionValues vector must have more than 1 values !"); + } + if (! (L_min + static_cast<double>(FilterFunctionValues.size() - 1) * L_userStep < (L_max + epsilon) ) ) + { + itkGenericExceptionMacro( + << "The following condition: " << L_min << "+(" << FilterFunctionValues.size() << "-1)*" << L_userStep << + " < (" << L_max << "+" << epsilon << ") is not respected !" << "val1 " << L_min + static_cast<double>(FilterFunctionValues.size() - 1) * L_userStep << " val2 " << L_max - epsilon); + } + + // Generate WavelengthSpectralBand if the step is not the offical 6S step value + if (vcl_abs(L_userStep - SIXSStepOfWavelengthSpectralBandValues) > epsilon) + { + ValuesVectorType values(1, FilterFunctionValues[0]); //vector size 1 with the value vect[0] + + // Stop the interpolation at the max spectral value. + value = SIXSStepOfWavelengthSpectralBandValues; + while (L_min + value <= L_max) + { + // Search the User interval that surround the StepOfWavelengthSpectralBandValues current value. + + // removed the <= here, might be wrong + while (j * L_userStep < value) + { + ++j; + } + + // Check if we are not out of bound + if (j >= FilterFunctionValues.size()) + { + itkGenericExceptionMacro( + << "Index " << j << " out of bound for FilterFunctionValues vector (size: " << FilterFunctionValues.size() << + ")." << " and value is " << value << " and SIXSStepOfWavelengthSpectralBandValues is " << SIXSStepOfWavelengthSpectralBandValues<< " and i is " << i); + } + + double valueTemp; + valueTemp = static_cast<double>(FilterFunctionValues[j - 1]) + + ((static_cast<double>(FilterFunctionValues[j]) - + static_cast<double>(FilterFunctionValues[j - 1])) * invStep) + * (value - L_userStep * (j - 1)); + values.push_back(static_cast<WavelengthSpectralBandType>(valueTemp)); + + ++i; + value += SIXSStepOfWavelengthSpectralBandValues; + } + + if (L_min + (i - 1) * SIXSStepOfWavelengthSpectralBandValues != L_max) + { + values.push_back(0); + } + // Store this values + WavelengthSpectralBand->SetFilterFunctionValues6S(values); + // Store the new Max MaxSpectralValue + WavelengthSpectralBand->SetMaxSpectralValue(static_cast<WavelengthSpectralBandType>(L_min + i * + SIXSStepOfWavelengthSpectralBandValues)); + + } + else + { + // Init with copy of FilterFunctionValues input vector values + WavelengthSpectralBand->SetFilterFunctionValues6S(FilterFunctionValues); + } +} + +void +SIXSTraits2::ComputeEnvironmentalContribution(const double diffuseTransmittanceForRayleighScattering, + const double diffuseTransmittanceForAerosolScattering, + const double radiusInKilometers, + const double altitude, + const double cosineOfViewingAngle, + double& rayleighEstimation, + double& aerosolEstimation, + double& globalEstimation) +{ + otb_6s_doublereal difr(static_cast<otb_6s_doublereal>(diffuseTransmittanceForRayleighScattering)); + otb_6s_doublereal difa(static_cast<otb_6s_doublereal>(diffuseTransmittanceForAerosolScattering)); + otb_6s_doublereal rad(static_cast<otb_6s_doublereal>(radiusInKilometers)); + otb_6s_doublereal palt(static_cast<otb_6s_doublereal>(altitude)); + otb_6s_doublereal xmuv(static_cast<otb_6s_doublereal>(cosineOfViewingAngle)); + otb_6s_doublereal fra(0.), fae(0.), fr(0.); + otb_6s_enviro_(&difr, &difa, &rad, &palt, &xmuv, &fra, &fae, &fr); + rayleighEstimation = static_cast<double>(fra); + aerosolEstimation = static_cast<double>(fae); + globalEstimation = static_cast<double>(fr); +} + +} // namespace otb diff --git a/Code/Radiometry/otbSIXSTraits2.h b/Code/Radiometry/otbSIXSTraits2.h new file mode 100644 index 0000000000000000000000000000000000000000..65c1382f5527ab98dcea64bc3c535af2daf67fd8 --- /dev/null +++ b/Code/Radiometry/otbSIXSTraits2.h @@ -0,0 +1,98 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __otbSIXSTraits2_h +#define __otbSIXSTraits2_h + +#include "otbAtmosphericCorrectionParameters2.h" + +namespace otb +{ + +/** \class SIXSTraits2 + * \brief SIXSTraits2 operations. + * + * Call 6S main function. The main method call 6S to calculate atmospheric correction parameters. + * It use by the OTB Atmospheric correction framework. + * + * \ingroup Radiometry + * + */ +class ITK_EXPORT SIXSTraits2 +{ +public: + + /** Standard class typedefs. */ + typedef SIXSTraits2 Self; + + typedef FilterFunctionValues WavelengthSpectralType; + typedef AtmosphericCorrectionParameters2::AerosolModelType AerosolModelType; + typedef WavelengthSpectralType::WavelengthSpectralBandType WavelengthSpectralBandType; + typedef WavelengthSpectralType::ValuesVectorType ValuesVectorType; + + /** Call 6S main function */ + static void ComputeAtmosphericParameters( + const double SolarZenithalAngle, /** The Solar zenithal angle */ + const double SolarAzimutalAngle, /** The Solar azimutal angle */ + const double ViewingZenithalAngle, /** The Viewing zenithal angle */ + const double ViewingAzimutalAngle, /** The Viewing azimutal angle */ + const unsigned int Month, /** The Month */ + const unsigned int Day, /** The Day (in the month) */ + const double AtmosphericPressure, /** The Atmospheric pressure */ + const double WaterVaporAmount, /** The Water vapor amount (Total water vapor content over vertical atmospheric column) */ + const double OzoneAmount, /** The Ozone amount (Stratospheric ozone layer content) */ + const AerosolModelType& AerosolModel, /** The Aerosol model */ + const double AerosolOptical, /** The Aerosol optical (radiative impact of aerosol for the reference wavelength 550-nm) */ + WavelengthSpectralType* WavelengthSpectralBand, /** Wavelength for the spectral band definition */ + /** Note : The Max wavelength spectral band value must be updated ! */ + double& AtmosphericReflectance, /** Atmospheric reflectance */ + double& AtmosphericSphericalAlbedo, /** atmospheric spherical albedo */ + double& TotalGaseousTransmission, /** Total gaseous transmission */ + double& DownwardTransmittance, /** downward transmittance */ + double& UpwardTransmittance, /** upward transmittance */ + double& UpwardDiffuseTransmittance, /** upward diffuse transmittance */ + double& UpwardDirectTransmittance, /** Upward direct transmittance */ + double& UpwardDiffuseTransmittanceForRayleigh, /** upward diffuse transmittance for rayleigh */ + double& UpwardDiffuseTransmittanceForAerosol /** upward diffuse transmittance for aerosols */ + ); + + /** + * Check the correpondance between the vector value size and + * the interval number between min and max. + * If the vector step is not at 0.0025, the new values are computed. + * The output vector values is store in the m_FilterFunctionValues6S + * of WavelengthSpectralBand + * + */ + static void ComputeWavelengthSpectralBandValuesFor6S( + const double SIXSStepOfWavelengthSpectralBandValues, + WavelengthSpectralType* WavelengthSpectralBand + ); + + static void ComputeEnvironmentalContribution(const double diffuseTransmittanceForRayleighScattering, + const double diffuseTransmittanceForAerosolScattering, + const double radiusInKilometers, + const double altitude, + const double cosineOfViewingAngle, + double& rayleighEstimation, + double& aerosolEstimation, + double& globalEstimation); +}; + +} // namespace otb + +#endif diff --git a/Testing/Code/Radiometry/otbAtmosphericCorrectionParametersToAtmosphericRadiativeTerms.cxx b/Testing/Code/Radiometry/otbAtmosphericCorrectionParametersToAtmosphericRadiativeTerms.cxx index 8c854a55820895bbbcf2ad4c5218c705b623de5e..94b06192a1a8e69ef6d84e185503904efa430258 100644 --- a/Testing/Code/Radiometry/otbAtmosphericCorrectionParametersToAtmosphericRadiativeTerms.cxx +++ b/Testing/Code/Radiometry/otbAtmosphericCorrectionParametersToAtmosphericRadiativeTerms.cxx @@ -17,7 +17,8 @@ =========================================================================*/ #include "otbRadiometryCorrectionParametersToAtmosphericRadiativeTerms.h" -#include "otbAtmosphericCorrectionParameters.h" +#include "otbAtmosphericCorrectionParameters2.h" +#include "otbImageMetadataCorrectionParameters.h" #include "otbAtmosphericRadiativeTerms2.h" #include <vector> #include <fstream> @@ -31,15 +32,17 @@ int otbAtmosphericCorrectionParametersToAtmosphericRadiativeTerms(int argc, char const char * outputFile = argv[2]; typedef otb::RadiometryCorrectionParametersToAtmosphericRadiativeTerms CorrectionParametersToRadiativeTermsType; - typedef otb::AtmosphericCorrectionParameters CorrectionParametersType; + typedef otb::AtmosphericCorrectionParameters2 AtmoCorrectionParametersType; + typedef otb::ImageMetadataCorrectionParameters AcquiCorrectionParametersType; typedef otb::AtmosphericRadiativeTerms2 RadiativeTermsType; - typedef CorrectionParametersType::AerosolModelType AerosolModelType; + typedef AtmoCorrectionParametersType::AerosolModelType AerosolModelType; typedef otb::FilterFunctionValues FilterFunctionValuesType; typedef FilterFunctionValuesType::WavelengthSpectralBandType ValueType; typedef FilterFunctionValuesType::ValuesVectorType ValuesVectorType; // Instantiating object - CorrectionParametersType::Pointer param = CorrectionParametersType::New(); + AtmoCorrectionParametersType::Pointer paramAtmo = AtmoCorrectionParametersType::New(); + AcquiCorrectionParametersType::Pointer paramAcqui = AcquiCorrectionParametersType::New(); AerosolModelType aerosolModel; FilterFunctionValuesType::Pointer functionValues = FilterFunctionValuesType::New(); ValuesVectorType vect; @@ -97,22 +100,22 @@ int otbAtmosphericCorrectionParametersToAtmosphericRadiativeTerms(int argc, char functionValues->SetMinSpectralValue(minSpectralValue); functionValues->SetMaxSpectralValue(maxSpectralValue); functionValues->SetUserStep(val); - param->SetWavelengthSpectralBandWithIndex(0, functionValues); + paramAcqui->SetWavelengthSpectralBandWithIndex(0, functionValues); //CHRIS // Set parameters - param->SetSolarZenithalAngle(static_cast<double>(solarZenithalAngle)); - param->SetSolarAzimutalAngle(static_cast<double>(solarAzimutalAngle)); - param->SetViewingZenithalAngle(static_cast<double>(viewingZenithalAngle)); - param->SetViewingAzimutalAngle(static_cast<double>(viewingAzimutalAngle)); - param->SetMonth(month); - param->SetDay(day); - param->SetAtmosphericPressure(static_cast<double>(atmosphericPressure)); - param->SetWaterVaporAmount(static_cast<double>(waterVaporAmount)); - param->SetOzoneAmount(static_cast<double>(ozoneAmount)); - param->SetAerosolModel(aerosolModel); - param->SetAerosolOptical(static_cast<double>(aerosolOptical)); - - RadiativeTermsType::Pointer radiative = CorrectionParametersToRadiativeTermsType::Compute(param); + paramAcqui->SetSolarZenithalAngle(static_cast<double>(solarZenithalAngle)); //CHRIS + paramAcqui->SetSolarAzimutalAngle(static_cast<double>(solarAzimutalAngle)); + paramAcqui->SetViewingZenithalAngle(static_cast<double>(viewingZenithalAngle)); + paramAcqui->SetViewingAzimutalAngle(static_cast<double>(viewingAzimutalAngle)); + paramAcqui->SetMonth(month); + paramAcqui->SetDay(day); + paramAtmo->SetAtmosphericPressure(static_cast<double>(atmosphericPressure)); + paramAtmo->SetWaterVaporAmount(static_cast<double>(waterVaporAmount)); + paramAtmo->SetOzoneAmount(static_cast<double>(ozoneAmount)); + paramAtmo->SetAerosolModel(aerosolModel); + paramAtmo->SetAerosolOptical(static_cast<double>(aerosolOptical)); + + RadiativeTermsType::Pointer radiative = CorrectionParametersToRadiativeTermsType::Compute(paramAtmo,paramAcqui); fout.open(outputFile);