diff --git a/Code/Radiometry/otbAtmosphericCorrectionParameters.cxx b/Code/Radiometry/otbAtmosphericCorrectionParameters.cxx index 1c37a3efc2194cf1fa4926073ba6b6979dab59d1..32c05266049c20bc6ee50084ae3ecb50d830124f 100644 --- a/Code/Radiometry/otbAtmosphericCorrectionParameters.cxx +++ b/Code/Radiometry/otbAtmosphericCorrectionParameters.cxx @@ -37,6 +37,9 @@ AtmosphericCorrectionParameters m_OzoneAmount = 0.28; m_AerosolModel = CONTINENTAL; m_AerosolOptical = 0.2; + m_AeronetFileName = ""; + m_Day = 1; + m_Month = 1; } @@ -69,24 +72,12 @@ void AtmosphericCorrectionParameters ::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/otbAtmosphericCorrectionParameters.h b/Code/Radiometry/otbAtmosphericCorrectionParameters.h index 3d7d0409febdc2eb19978943a8a25d5c7505a26e..792b681e094372511acd005dfd2cdd7c4b2f999c 100644 --- a/Code/Radiometry/otbAtmosphericCorrectionParameters.h +++ b/Code/Radiometry/otbAtmosphericCorrectionParameters.h @@ -90,17 +90,37 @@ public: itkGetMacro(AerosolOptical, double); + /** Get/Set Aeronet file name. */ + itkSetMacro(AeronetFileName, std::string); + itkGetMacro(AeronetFileName, std::string); + /** Read the aeronet data and extract aerosol optical and water vapor amount. */ void ReadAeronetData(const std::string& file, int year, int month, int day, int hour, int minute, double epsi); - void UpdateAeronetData(const std::string& file, int year, int month, int day, int hour, int minute, double epsi); void UpdateAeronetData(const std::string& file, int year, int month, int day, int hour, int minute) { this->UpdateAeronetData(file, year, month, day, hour, minute, 0.4); } + void UpdateAeronetData(int year, int month, int day, int hour, int minute, double epsi) + { + this->UpdateAeronetData(m_AeronetFileName, year, month, day, hour, minute, epsi); + } + void UpdateAeronetData(int year, int month, int day, int hour, int minute) + { + this->UpdateAeronetData(m_AeronetFileName, year, month, day, hour, minute, 0.4); + } + void UpdateAeronetData(const std::string& file, int year, int hour, int minute) + { + this->UpdateAeronetData(file, year, m_Month, m_Day, hour, minute, 0.4); + } + void UpdateAeronetData(int year, int hour, int minute) + { + this->UpdateAeronetData(m_AeronetFileName, year, m_Month, m_Day, hour, minute, 0.4); + } + /*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); @@ -123,7 +143,12 @@ protected: private: AtmosphericCorrectionParameters(const Self &); //purposely not implemented void operator =(const Self&); //purposely not implemented - + /** Path to an Aeronet data file, allows to compute aerosol optical and water vapor amounts. */ + std::string m_AeronetFileName; + /** Day */ + int m_Day; + /** Month */ + int m_Month; /** The Atmospheric pressure */ double m_AtmosphericPressure; /** The Water vapor amount (Total water vapor content over vertical atmospheric column) */ diff --git a/Code/Radiometry/otbImageMetadataCorrectionParameters.cxx b/Code/Radiometry/otbImageMetadataCorrectionParameters.cxx index 1e7ad8cedd81dee13aa31a9beb753b2e85a3484e..30a79df1855811ef0de7c67a9b6e6ae985062584 100644 --- a/Code/Radiometry/otbImageMetadataCorrectionParameters.cxx +++ b/Code/Radiometry/otbImageMetadataCorrectionParameters.cxx @@ -38,6 +38,7 @@ ImageMetadataCorrectionParameters m_Month = 0; m_Day = 0; m_Year = 0; + m_FilterFunctionValuesFileName = ""; m_WavelengthSpectralBand = InternalWavelengthSpectralBandVectorType::New(); m_WavelengthSpectralBand->Clear(); diff --git a/Code/Radiometry/otbImageMetadataCorrectionParameters.h b/Code/Radiometry/otbImageMetadataCorrectionParameters.h index 2cc1ec8c1971dc18e193b3e45013ec697adfcde5..0bde6a032577109da60efea781abe2beb4896419 100644 --- a/Code/Radiometry/otbImageMetadataCorrectionParameters.h +++ b/Code/Radiometry/otbImageMetadataCorrectionParameters.h @@ -95,6 +95,10 @@ public: */ itkSetMacro(Year, unsigned int); itkGetMacro(Year, unsigned int); + + /** Get/Set FilterFunction file name. */ + itkSetMacro(FilterFunctionValuesFileName, std::string); + itkGetMacro(FilterFunctionValuesFileName, std::string); /** * Set/Get the wavelength spectral band. @@ -128,6 +132,10 @@ public: * Read a file that contains filter function values on the 6S format. */ void LoadFilterFunctionValue(const std::string& filename); + void LoadFilterFunctionValue() + { + this->LoadFilterFunctionValue(m_FilterFunctionValuesFileName); + } /** Constructor */ ImageMetadataCorrectionParameters(); @@ -157,6 +165,7 @@ private: unsigned int m_Day; /** The Year */ unsigned int m_Year; + std::string m_FilterFunctionValuesFileName; /** Wavelength for the each spectral band definition */ WavelengthSpectralBandVectorType m_WavelengthSpectralBand; diff --git a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h index 0362bc1a70fc9330c97b674092e3518608cc63d4..b880a23950323a908117cc094a5b9e87cf2cf549 100644 --- a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h +++ b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h @@ -170,13 +170,13 @@ public: typedef otb::RadiometryCorrectionParametersToAtmosphericRadiativeTerms CorrectionParametersToRadiativeTermsType; - typedef otb::AtmosphericCorrectionParameters AtmoCorrectionParametersType; + typedef otb::AtmosphericCorrectionParameters AtmoCorrectionParametersType; typedef typename AtmoCorrectionParametersType::Pointer AtmoCorrectionParametersPointerType; typedef otb::ImageMetadataCorrectionParameters AcquiCorrectionParametersType; typedef typename AcquiCorrectionParametersType::Pointer AcquiCorrectionParametersPointerType; - typedef otb::AtmosphericRadiativeTerms AtmosphericRadiativeTermsType; + typedef otb::AtmosphericRadiativeTerms AtmosphericRadiativeTermsType; typedef typename AtmosphericRadiativeTermsType::Pointer AtmosphericRadiativeTermsPointerType; diff --git a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx index 900e5901ef66d7e0267cb0275fba10653b906f26..ff777e6fc7ccb9c6abc10fef0ef55f2006657a68 100644 --- a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx +++ b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx @@ -34,12 +34,11 @@ ReflectanceToSurfaceReflectanceImageFilter<TInputImage, TOutputImage> m_IsSetAtmosphericRadiativeTerms(false), m_IsSetAtmoCorrectionParameters(false), m_IsSetAcquiCorrectionParameters(false), - m_UseGenerateParameters(true), - m_AtmosphericRadiativeTerms(NULL), - m_AtmoCorrectionParameters(NULL), - m_AcquiCorrectionParameters(NULL) + m_UseGenerateParameters(true) { - + m_AtmosphericRadiativeTerms = AtmosphericRadiativeTermsType::New(); + m_AtmoCorrectionParameters = AtmoCorrectionParametersType::New(); + m_AcquiCorrectionParameters = AcquiCorrectionParametersType::New(); } template <class TInputImage, class TOutputImage> @@ -97,14 +96,30 @@ ReflectanceToSurfaceReflectanceImageFilter<TInputImage, TOutputImage> { itkExceptionMacro(<< "Atmospheric correction parameters must be provided before updating the atmospheric radiative terms"); } + + + // load filter function values + bool SetFilterFunctionValuesFileName=false; + if (m_AcquiCorrectionParameters->GetFilterFunctionValuesFileName() != "") + { + m_AcquiCorrectionParameters->LoadFilterFunctionValue(); + SetFilterFunctionValuesFileName=true; + } + + + MetaDataDictionaryType dict = this->GetInput()->GetMetaDataDictionary(); + OpticalImageMetadataInterface::Pointer imageMetadataInterface = OpticalImageMetadataInterfaceFactory::CreateIMI(dict); + + + if (m_AtmoCorrectionParameters->GetAeronetFileName() != "") + m_AtmoCorrectionParameters->UpdateAeronetData(imageMetadataInterface->GetYear(), + imageMetadataInterface->GetHour(), + imageMetadataInterface->GetMinute()); // Acquisition parameters if (!m_IsSetAcquiCorrectionParameters) // Get info from image metadata interface { - MetaDataDictionaryType dict = this->GetInput()->GetMetaDataDictionary(); - OpticalImageMetadataInterface::Pointer imageMetadataInterface = OpticalImageMetadataInterfaceFactory::CreateIMI(dict); - m_AcquiCorrectionParameters = AcquiCorrectionParametersType::New(); m_AcquiCorrectionParameters->SetSolarZenithalAngle(90. - imageMetadataInterface->GetSunElevation()); @@ -115,25 +130,30 @@ ReflectanceToSurfaceReflectanceImageFilter<TInputImage, TOutputImage> m_AcquiCorrectionParameters->SetDay(imageMetadataInterface->GetDay()); m_AcquiCorrectionParameters->SetMonth(imageMetadataInterface->GetMonth()); - if (imageMetadataInterface->GetSpectralSensitivity()->Capacity() > 0) - { - m_AcquiCorrectionParameters->SetWavelengthSpectralBand(imageMetadataInterface->GetSpectralSensitivity()); - - } - else - { - otbMsgDevMacro(<< "use dummy filter"); - WavelengthSpectralBandVectorType spectralDummy; - spectralDummy->Clear(); - for (unsigned int i = 0; i < this->GetInput()->GetNumberOfComponentsPerPixel(); ++i) - { - spectralDummy->PushBack(FilterFunctionValuesType::New()); - } - m_AcquiCorrectionParameters->SetWavelengthSpectralBand(spectralDummy); - } + + if (!SetFilterFunctionValuesFileName) + { + if (imageMetadataInterface->GetSpectralSensitivity()->Capacity() > 0) + { + m_AcquiCorrectionParameters->SetWavelengthSpectralBand(imageMetadataInterface->GetSpectralSensitivity()); + } + else + { + otbMsgDevMacro(<< "use dummy filter"); + WavelengthSpectralBandVectorType spectralDummy; + spectralDummy->Clear(); + for (unsigned int i = 0; i < this->GetInput()->GetNumberOfComponentsPerPixel(); ++i) + { + spectralDummy->PushBack(FilterFunctionValuesType::New()); + } + m_AcquiCorrectionParameters->SetWavelengthSpectralBand(spectralDummy); + } + } } + + m_AtmosphericRadiativeTerms = CorrectionParametersToRadiativeTermsType::Compute(m_AtmoCorrectionParameters,m_AcquiCorrectionParameters); }