Skip to content

GitLab

Commit a59b1cd7 authored by remi cresson's avatar remi cresson
Browse files

ENH: rainfall data is optional

parent 3399d36c
Hide whitespace changes
Inline Side-by-side
app/otbTimeSeriesIndexTrend.cxx
View file @ a59b1cd7
......@@ -170,9 +170,11 @@ private:
// Input rainfall time series images
AddParameter(ParameterType_InputImageList, "rainfts", "Input rainfall time series images");
MandatoryOff("rainfts");
// Input rainfall time series dates
AddParameter(ParameterType_InputFilename, "rainfdates", "Input rainfall time series dates ASCII file (Must be YYYYMMDD format)");
MandatoryOff("rainfdates");
// Parameter group for NDVI time series
AddParameter(ParameterType_Group, "ndvi", "NDVI Time series");
......@@ -374,131 +376,141 @@ private:
SetParameterOutputImage("ndvilabel", m_NDVITimeSeriesLabelFilter->GetOutput());
}
// Get the rainfall time series dates
std::vector<DateType> rainDates = GetTimeSeriesDates("rainfdates");
// Get the rainfall time series input images list
FloatVectorImageListType::Pointer inRFList = this->GetParameterImageList("rainfts");
if( inRFList->Size() < 2 )
// Check rainfall mandatory parameters
if (HasValue("rainfdates") && HasValue("rainfts"))
{
otbAppLogFATAL("At least two rainfall images are required.");
}
// Create one stack for input RF images list
m_RFConcatener = ListConcatenerFilterType::New();
m_RFImageList = ImageListType::New();
m_RFExtractorList = ExtractROIFilterListType::New();
// Split each input vector image into image
// and generate an mono channel image list
inRFList->GetNthElement(0)->UpdateOutputInformation();
FloatVectorImageType::SizeType rfSize = inRFList->GetNthElement(0)->GetLargestPossibleRegion().GetSize();
for( unsigned int i=0; i<inRFList->Size(); i++ )
{
FloatVectorImageType::Pointer vectIm = inRFList->GetNthElement(i);
vectIm->UpdateOutputInformation();
if( rfSize != vectIm->GetLargestPossibleRegion().GetSize() )
{
otbAppLogFATAL("Input NDVI image size number " << i << " mismatch");
}
// Get the rainfall time series dates
std::vector<DateType> rainDates = GetTimeSeriesDates("rainfdates");
// Get the rainfall time series input images list
FloatVectorImageListType::Pointer inRFList = this->GetParameterImageList("rainfts");
if( inRFList->Size() < 2 )
{
otbAppLogFATAL("At least two rainfall images are required.");
}
// Create one stack for input RF images list
m_RFConcatener = ListConcatenerFilterType::New();
m_RFImageList = ImageListType::New();
m_RFExtractorList = ExtractROIFilterListType::New();
// Split each input vector image into image
// and generate an mono channel image list
inRFList->GetNthElement(0)->UpdateOutputInformation();
FloatVectorImageType::SizeType rfSize = inRFList->GetNthElement(0)->GetLargestPossibleRegion().GetSize();
for( unsigned int i=0; i<inRFList->Size(); i++ )
{
FloatVectorImageType::Pointer vectIm = inRFList->GetNthElement(i);
vectIm->UpdateOutputInformation();
if( rfSize != vectIm->GetLargestPossibleRegion().GetSize() )
{
otbAppLogFATAL("Input NDVI image size number " << i << " mismatch");
}
for( unsigned int j=0; j<vectIm->GetNumberOfComponentsPerPixel(); j++)
{
ExtractROIFilterType::Pointer extractor = ExtractROIFilterType::New();
extractor->SetInput( vectIm );
extractor->SetChannel( j+1 );
extractor->UpdateOutputInformation();
m_RFExtractorList->PushBack( extractor );
m_RFImageList->PushBack( extractor->GetOutput() );
}
}
m_RFConcatener->SetInput( m_RFImageList );
m_RFConcatener->UpdateOutputInformation();
if (m_RFConcatener->GetOutput()->GetNumberOfComponentsPerPixel() != rainDates.size())
otbAppLogFATAL("Number of rainfall images and number of dates are different! There is "
<< (m_RFConcatener->GetOutput()->GetNumberOfComponentsPerPixel()) << " images "
<< "and " << rainDates.size() << " dates.");
// Check dates consistency
int firstYearForNDVI = ndviDates[0].year;
int lastYearForNDVI = ndviDates[ndviDates.size()-1].year;
int firstYearForRF = rainDates[0].year;
int lastYearForRF = rainDates[rainDates.size()-1].year;
if (firstYearForNDVI != firstYearForRF)
{
otbAppLogFATAL("NDVI time series and rainfall time series must start the same year!"
<< " NDVI start year " << firstYearForNDVI
<< "and rainfall start year " << firstYearForRF << ".");
}
if (lastYearForNDVI != lastYearForRF)
{
otbAppLogFATAL("NDVI time series and rainfall time series must end the same year!"
<< " NDVI last year " << lastYearForNDVI
<< "and rainfall last year " << lastYearForRF << ".");
}
// Reduce rainfall
m_RFCumulReduceFilter = TSCumulatedRangeReduceFilterType::New();
m_RFCumulReduceFilter->SetInput(m_RFConcatener->GetOutput());
m_RFCumulReduceFilter->GetFunctor().SetDates(ndviDates);
m_RFCumulReduceFilter->GetFunctor().SetDay1(GetParameterInt("rain.reduce.cumul.day1"));
m_RFCumulReduceFilter->GetFunctor().SetDay2(GetParameterInt("rain.reduce.cumul.day2"));
m_RFCumulReduceFilter->GetFunctor().SetMonth1(GetParameterInt("rain.reduce.cumul.month1"));
m_RFCumulReduceFilter->GetFunctor().SetMonth2(GetParameterInt("rain.reduce.cumul.month2"));
// Resample rainfall
LinearInterpolatorType::Pointer linInterpolator = LinearInterpolatorType::New();
m_ResampleFilter = ResampleImageFilterType::New();
m_ResampleFilter->SetInput(m_RFCumulReduceFilter->GetOutput());
m_ResampleFilter->SetOutputOrigin(m_NDVIConcatener->GetOutput()->GetOrigin());
m_ResampleFilter->SetOutputSpacing(m_NDVIConcatener->GetOutput()->GetSpacing());
m_ResampleFilter->SetOutputSize(m_NDVIConcatener->GetOutput()->GetLargestPossibleRegion().GetSize());
m_ResampleFilter->SetInterpolator(linInterpolator);
// Compute residues
m_ResiduesFilter = ResiduesFilterType::New();
m_ResiduesFilter->SetInput1(reducedNDVI);
m_ResiduesFilter->SetInput2(m_RFCumulReduceFilter->GetOutput());
// Output residues
if (HasValue("residues"))
{
SetParameterOutputImage("residues", m_ResiduesFilter->GetOutput());
}
// Regression (residues)
m_NDVIResiduesRegressionFilter = SlopeAndPValueFilterType::New();
m_NDVIResiduesRegressionFilter->SetInput(m_ResiduesFilter->GetOutput());
// Output residues trend
if (HasValue("restrend"))
{
SetParameterOutputImage("restrend", m_NDVIResiduesRegressionFilter->GetOutput());
}
// Output residues trend labels
if (HasValue("reslabel"))
{
m_NDVIResiduesLabelFilter = SlopeAndPValueLabelFilterType::New();
m_NDVIResiduesLabelFilter->SetInput(m_NDVIResiduesRegressionFilter->GetOutput());
SetParameterOutputImage("reslabel", m_NDVIResiduesLabelFilter->GetOutput());
}
// Correlation between NDVI and Rainfall
m_CorrelationFilter = PearsonCorrelationFilterType::New();
m_CorrelationFilter->SetInput1(reducedNDVI);
m_CorrelationFilter->SetInput2(m_RFCumulReduceFilter->GetOutput());
// Factors labeling
m_FactorsLabelFilter = FactorsLabelingFilterType::New();
m_FactorsLabelFilter->SetInput1(m_NDVITimeSeriesRegressionFilter->GetOutput());
m_FactorsLabelFilter->SetInput2(m_NDVIResiduesRegressionFilter->GetOutput());
m_FactorsLabelFilter->SetInput3(m_CorrelationFilter->GetOutput());
// Write factors image
SetParameterOutputImage("factorslabel", m_FactorsLabelFilter->GetOutput());
for( unsigned int j=0; j<vectIm->GetNumberOfComponentsPerPixel(); j++)
{
ExtractROIFilterType::Pointer extractor = ExtractROIFilterType::New();
extractor->SetInput( vectIm );
extractor->SetChannel( j+1 );
extractor->UpdateOutputInformation();
m_RFExtractorList->PushBack( extractor );
m_RFImageList->PushBack( extractor->GetOutput() );
}
}
m_RFConcatener->SetInput( m_RFImageList );
m_RFConcatener->UpdateOutputInformation();
if (m_RFConcatener->GetOutput()->GetNumberOfComponentsPerPixel() != rainDates.size())
otbAppLogFATAL("Number of rainfall images and number of dates are different! There is "
<< (m_RFConcatener->GetOutput()->GetNumberOfComponentsPerPixel()) << " images "
<< "and " << rainDates.size() << " dates.");
// Check dates consistency
int firstYearForNDVI = ndviDates[0].year;
int lastYearForNDVI = ndviDates[ndviDates.size()-1].year;
int firstYearForRF = rainDates[0].year;
int lastYearForRF = rainDates[rainDates.size()-1].year;
if (firstYearForNDVI != firstYearForRF)
{
otbAppLogFATAL("NDVI time series and rainfall time series must start the same year!"
<< " NDVI start year " << firstYearForNDVI
<< "and rainfall start year " << firstYearForRF << ".");
}
if (lastYearForNDVI != lastYearForRF)
{
otbAppLogFATAL("NDVI time series and rainfall time series must end the same year!"
<< " NDVI last year " << lastYearForNDVI
<< "and rainfall last year " << lastYearForRF << ".");
}
// Reduce rainfall
m_RFCumulReduceFilter = TSCumulatedRangeReduceFilterType::New();
m_RFCumulReduceFilter->SetInput(m_RFConcatener->GetOutput());
m_RFCumulReduceFilter->GetFunctor().SetDates(ndviDates);
m_RFCumulReduceFilter->GetFunctor().SetDay1(GetParameterInt("rain.reduce.cumul.day1"));
m_RFCumulReduceFilter->GetFunctor().SetDay2(GetParameterInt("rain.reduce.cumul.day2"));
m_RFCumulReduceFilter->GetFunctor().SetMonth1(GetParameterInt("rain.reduce.cumul.month1"));
m_RFCumulReduceFilter->GetFunctor().SetMonth2(GetParameterInt("rain.reduce.cumul.month2"));
// Resample rainfall
LinearInterpolatorType::Pointer linInterpolator = LinearInterpolatorType::New();
m_ResampleFilter = ResampleImageFilterType::New();
m_ResampleFilter->SetInput(m_RFCumulReduceFilter->GetOutput());
m_ResampleFilter->SetOutputOrigin(m_NDVIConcatener->GetOutput()->GetOrigin());
m_ResampleFilter->SetOutputSpacing(m_NDVIConcatener->GetOutput()->GetSpacing());
m_ResampleFilter->SetOutputSize(m_NDVIConcatener->GetOutput()->GetLargestPossibleRegion().GetSize());
m_ResampleFilter->SetInterpolator(linInterpolator);
// Compute residues
m_ResiduesFilter = ResiduesFilterType::New();
m_ResiduesFilter->SetInput1(reducedNDVI);
m_ResiduesFilter->SetInput2(m_RFCumulReduceFilter->GetOutput());
// Output residues
if (HasValue("residues"))
{
SetParameterOutputImage("residues", m_ResiduesFilter->GetOutput());
}
// Regression (residues)
m_NDVIResiduesRegressionFilter = SlopeAndPValueFilterType::New();
m_NDVIResiduesRegressionFilter->SetInput(m_ResiduesFilter->GetOutput());
// Output residues trend
if (HasValue("restrend"))
{
SetParameterOutputImage("restrend", m_NDVIResiduesRegressionFilter->GetOutput());
}
// Output residues trend labels
if (HasValue("reslabel"))
else
{
m_NDVIResiduesLabelFilter = SlopeAndPValueLabelFilterType::New();
m_NDVIResiduesLabelFilter->SetInput(m_NDVIResiduesRegressionFilter->GetOutput());
SetParameterOutputImage("reslabel", m_NDVIResiduesLabelFilter->GetOutput());
otbAppLogINFO("Rainfall data not set. Skipping rainfall based indices.")
}
// Correlation between NDVI and Rainfall
m_CorrelationFilter = PearsonCorrelationFilterType::New();
m_CorrelationFilter->SetInput1(reducedNDVI);
m_CorrelationFilter->SetInput2(m_RFCumulReduceFilter->GetOutput());
// Factors labeling
m_FactorsLabelFilter = FactorsLabelingFilterType::New();
m_FactorsLabelFilter->SetInput1(m_NDVITimeSeriesRegressionFilter->GetOutput());
m_FactorsLabelFilter->SetInput2(m_NDVIResiduesRegressionFilter->GetOutput());
m_FactorsLabelFilter->SetInput3(m_CorrelationFilter->GetOutput());
// Write factors image
SetParameterOutputImage("factorslabel", m_FactorsLabelFilter->GetOutput());
} // DOExecute()
void AfterExecuteAndWriteOutputs()
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment