RFT ensemble refactoring

* Compute average MD for intersections with a cell
* Create extractor for simulation well
* Remove rftReader from RifDataSourceForRftPlt
* Add function compute measured depth for RFT cells based on well path geometry
* Move statistics reader to well ensemble curve set
* Make sure both TVD and MD are cached if possible
* Add selection of grid case to use for estimation of measured depth (MD)
Add "Grid Model For MD" where the user can select a grid model. This grid model is propagated to a hidden field in EnsembleCurveSet. The grid model is then applied to RifReaderEnsembleStatisticsRft owned by EnsembleCurveSet

ApplicationLibCode/FileInterface/RifDataSourceForRftPlt.cpp

@@ -18,7 +18,6 @@
 
 #include "RifDataSourceForRftPlt.h"
 
-#include "RifReaderEclipseRft.h"
 #include "RigEclipseCaseData.h"
 
 #include "RimEclipseCase.h"
@@ -230,37 +229,6 @@ auto RifDataSourceForRftPlt::operator<=>( const RifDataSourceForRftPlt& addr2 )
     }
 }
 
-//--------------------------------------------------------------------------------------------------
-///
-//--------------------------------------------------------------------------------------------------
-RifReaderRftInterface* RifDataSourceForRftPlt::rftReader() const
-{
-    if ( m_sourceType == SourceType::GRID_MODEL_CELL_DATA || m_sourceType == SourceType::RFT_SIM_WELL_DATA )
-    {
-        // TODO: Consider changing to RimEclipseResultCase to avoid casting
-        auto eclResCase = dynamic_cast<RimEclipseResultCase*>( eclCase() );
-
-        if ( eclResCase ) return eclResCase->rftReader();
-    }
-    else if ( m_sourceType == SourceType::SUMMARY_RFT )
-    {
-        if ( m_summaryCase ) return m_summaryCase->rftReader();
-    }
-    else if ( m_sourceType == SourceType::ENSEMBLE_RFT )
-    {
-        if ( m_ensemble ) return m_ensemble->rftStatisticsReader();
-    }
-    else if ( m_sourceType == SourceType::OBSERVED_FMU_RFT )
-    {
-        if ( m_observedFmuRftData ) return m_observedFmuRftData->rftReader();
-    }
-    else if ( m_sourceType == SourceType::OBSERVED_PRESSURE_DEPTH )
-    {
-        if ( m_pressureDepthData ) return m_pressureDepthData->rftReader();
-    }
-    return nullptr;
-}
-
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationLibCode/FileInterface/RifDataSourceForRftPlt.h

@@ -31,7 +31,6 @@
 
 class RimWellLogFile;
 class RimEclipseCase;
-class RifReaderRftInterface;
 class RimSummaryCase;
 class RimSummaryCaseCollection;
 class RimObservedFmuRftData;
@@ -64,8 +63,7 @@ public:
     RifDataSourceForRftPlt( RimObservedFmuRftData* observedFmuRftData );
     RifDataSourceForRftPlt( RimPressureDepthData* pressureDepthData );
 
-    SourceType             sourceType() const;
-    RifReaderRftInterface* rftReader() const;
+    SourceType sourceType() const;
 
     RimEclipseCase*           eclCase() const;
     RimSummaryCase*           summaryCase() const;

ApplicationLibCode/FileInterface/RifReaderEclipseRft.cpp

@@ -246,19 +246,17 @@ void RifReaderEclipseRft::values( const RifEclipseRftAddress& rftAddress, std::v
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-void RifReaderEclipseRft::cellIndices( const RifEclipseRftAddress& rftAddress, std::vector<caf::VecIjk>* indices )
+std::vector<caf::VecIjk> RifReaderEclipseRft::cellIndices( const QString& wellName, const QDateTime& timeStep )
 {
-    CVF_ASSERT( indices );
-
     if ( !m_ecl_rft_file )
     {
         open();
     }
 
-    indices->clear();
+    int index = indexFromAddress( wellName, timeStep );
+    if ( index < 0 ) return {};
 
-    int index = indexFromAddress( rftAddress );
-    if ( index < 0 ) return;
+    std::vector<caf::VecIjk> indices;
 
     ecl_rft_node_type* node = ecl_rft_file_iget_node( m_ecl_rft_file, index );
 
@@ -268,8 +266,10 @@ void RifReaderEclipseRft::cellIndices( const RifEclipseRftAddress& rftAddress, s
         ecl_rft_node_iget_ijk( node, cellIdx, &i, &j, &k );
 
         caf::VecIjk ijk( (size_t)i, (size_t)j, (size_t)k );
-        indices->push_back( ijk );
+        indices.push_back( ijk );
     }
+
+    return indices;
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -453,3 +453,19 @@ int RifReaderEclipseRft::indexFromAddress( const RifEclipseRftAddress& rftAddres
 
     return -1;
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+int RifReaderEclipseRft::indexFromAddress( const QString& wellName, const QDateTime& timeStep ) const
+{
+    for ( const auto& [rftAddr, nodeIdx] : m_rftAddressToLibeclNodeIdx )
+    {
+        if ( rftAddr.wellName() == wellName && rftAddr.timeStep() == timeStep )
+        {
+            return nodeIdx;
+        }
+    }
+
+    return -1;
+}

ApplicationLibCode/FileInterface/RifReaderEclipseRft.h

@@ -45,7 +45,7 @@ public:
 
     std::set<RifEclipseRftAddress> eclipseRftAddresses() override;
     void                           values( const RifEclipseRftAddress& rftAddress, std::vector<double>* values ) override;
-    void                           cellIndices( const RifEclipseRftAddress& rftAddress, std::vector<caf::VecIjk>* indices ) override;
+    std::vector<caf::VecIjk>       cellIndices( const QString& wellName, const QDateTime& timeStep ) override;
 
     std::set<QDateTime> availableTimeSteps( const QString& wellName ) override;
     std::set<QDateTime> availableTimeSteps( const QString&                                               wellName,
@@ -60,6 +60,7 @@ public:
 private:
     void open();
     int  indexFromAddress( const RifEclipseRftAddress& rftAddress ) const;
+    int  indexFromAddress( const QString& wellName, const QDateTime& timeStep ) const;
 
 private:
     // NOLINTBEGIN(modernize-use-using)

ApplicationLibCode/FileInterface/RifReaderEnsembleStatisticsRft.cpp

@@ -18,20 +18,23 @@
 
 #include "RifReaderEnsembleStatisticsRft.h"
 
-#include "RiaCurveMerger.h"
-#include "RiaWeightedMeanCalculator.h"
+#include "RiaExtractionTools.h"
+
 #include "RigStatisticsMath.h"
 
 #include "RimSummaryCase.h"
 #include "RimSummaryCaseCollection.h"
+#include "RimTools.h"
 
 #include "cafAssert.h"
 
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-RifReaderEnsembleStatisticsRft::RifReaderEnsembleStatisticsRft( const RimSummaryCaseCollection* summaryCaseCollection )
+RifReaderEnsembleStatisticsRft::RifReaderEnsembleStatisticsRft( const RimSummaryCaseCollection* summaryCaseCollection,
+                                                                RimEclipseCase*                 eclipseCase )
     : m_summaryCaseCollection( summaryCaseCollection )
+    , m_eclipseCase( eclipseCase )
 {
 }
 
@@ -40,6 +43,8 @@ RifReaderEnsembleStatisticsRft::RifReaderEnsembleStatisticsRft( const RimSummary
 //--------------------------------------------------------------------------------------------------
 std::set<RifEclipseRftAddress> RifReaderEnsembleStatisticsRft::eclipseRftAddresses()
 {
+    if ( !m_summaryCaseCollection ) return {};
+
     std::set<RifEclipseRftAddress> allAddresses;
     for ( auto summaryCase : m_summaryCaseCollection->allSummaryCases() )
     {
@@ -78,7 +83,10 @@ std::set<RifEclipseRftAddress> RifReaderEnsembleStatisticsRft::eclipseRftAddress
 //--------------------------------------------------------------------------------------------------
 void RifReaderEnsembleStatisticsRft::values( const RifEclipseRftAddress& rftAddress, std::vector<double>* values )
 {
-    CAF_ASSERT( rftAddress.wellLogChannel() == RifEclipseRftAddress::RftWellLogChannelType::TVD ||
+    if ( !m_summaryCaseCollection ) return;
+
+    CAF_ASSERT( rftAddress.wellLogChannel() == RifEclipseRftAddress::RftWellLogChannelType::MD ||
+                rftAddress.wellLogChannel() == RifEclipseRftAddress::RftWellLogChannelType::TVD ||
                 rftAddress.wellLogChannel() == RifEclipseRftAddress::RftWellLogChannelType::PRESSURE_MEAN ||
                 rftAddress.wellLogChannel() == RifEclipseRftAddress::RftWellLogChannelType::PRESSURE_P10 ||
                 rftAddress.wellLogChannel() == RifEclipseRftAddress::RftWellLogChannelType::PRESSURE_P50 ||
@@ -88,7 +96,7 @@ void RifReaderEnsembleStatisticsRft::values( const RifEclipseRftAddress& rftAddr
     auto it = m_cachedValues.find( rftAddress );
     if ( it == m_cachedValues.end() )
     {
-        calculateStatistics( rftAddress );
+        calculateStatistics( rftAddress.wellName(), rftAddress.timeStep() );
     }
     *values = m_cachedValues[rftAddress];
 }
@@ -98,6 +106,8 @@ void RifReaderEnsembleStatisticsRft::values( const RifEclipseRftAddress& rftAddr
 //--------------------------------------------------------------------------------------------------
 std::set<QDateTime> RifReaderEnsembleStatisticsRft::availableTimeSteps( const QString& wellName )
 {
+    if ( !m_summaryCaseCollection ) return {};
+
     std::set<QDateTime> allTimeSteps;
     for ( auto summaryCase : m_summaryCaseCollection->allSummaryCases() )
     {
@@ -116,6 +126,8 @@ std::set<QDateTime> RifReaderEnsembleStatisticsRft::availableTimeSteps( const QS
 std::set<QDateTime> RifReaderEnsembleStatisticsRft::availableTimeSteps( const QString& wellName,
                                                                         const RifEclipseRftAddress::RftWellLogChannelType& wellLogChannelName )
 {
+    if ( !m_summaryCaseCollection ) return {};
+
     std::set<QDateTime> allTimeSteps;
     for ( auto summaryCase : m_summaryCaseCollection->allSummaryCases() )
     {
@@ -135,6 +147,8 @@ std::set<QDateTime>
     RifReaderEnsembleStatisticsRft::availableTimeSteps( const QString&                                               wellName,
                                                         const std::set<RifEclipseRftAddress::RftWellLogChannelType>& relevantChannels )
 {
+    if ( !m_summaryCaseCollection ) return {};
+
     std::set<QDateTime> allTimeSteps;
     for ( auto summaryCase : m_summaryCaseCollection->allSummaryCases() )
     {
@@ -152,6 +166,8 @@ std::set<QDateTime>
 //--------------------------------------------------------------------------------------------------
 std::set<RifEclipseRftAddress::RftWellLogChannelType> RifReaderEnsembleStatisticsRft::availableWellLogChannels( const QString& wellName )
 {
+    if ( !m_summaryCaseCollection ) return {};
+
     std::set<RifEclipseRftAddress::RftWellLogChannelType> allWellLogChannels;
     for ( auto summaryCase : m_summaryCaseCollection->allSummaryCases() )
     {
@@ -170,6 +186,8 @@ std::set<RifEclipseRftAddress::RftWellLogChannelType> RifReaderEnsembleStatistic
 //--------------------------------------------------------------------------------------------------
 std::set<QString> RifReaderEnsembleStatisticsRft::wellNames()
 {
+    if ( !m_summaryCaseCollection ) return {};
+
     std::set<QString> allWellNames;
     for ( auto summaryCase : m_summaryCaseCollection->allSummaryCases() )
     {
@@ -185,47 +203,110 @@ std::set<QString> RifReaderEnsembleStatisticsRft::wellNames()
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-void RifReaderEnsembleStatisticsRft::calculateStatistics( const RifEclipseRftAddress& rftAddress )
+void RifReaderEnsembleStatisticsRft::calculateStatistics( const QString& wellName, const QDateTime& timeStep )
 {
-    const QString&       wellName = rftAddress.wellName();
-    const QDateTime&     timeStep = rftAddress.timeStep();
-    RifEclipseRftAddress depthAddress =
-        RifEclipseRftAddress::createAddress( wellName, timeStep, RifEclipseRftAddress::RftWellLogChannelType::TVD );
-    RifEclipseRftAddress pressAddress =
-        RifEclipseRftAddress::createAddress( wellName, timeStep, RifEclipseRftAddress::RftWellLogChannelType::PRESSURE );
+    if ( !m_summaryCaseCollection ) return;
 
-    RifEclipseRftAddress p10Address =
-        RifEclipseRftAddress::createAddress( wellName, timeStep, RifEclipseRftAddress::RftWellLogChannelType::PRESSURE_P10 );
-    RifEclipseRftAddress p50Address =
-        RifEclipseRftAddress::createAddress( wellName, timeStep, RifEclipseRftAddress::RftWellLogChannelType::PRESSURE_P50 );
-    RifEclipseRftAddress p90Address =
-        RifEclipseRftAddress::createAddress( wellName, timeStep, RifEclipseRftAddress::RftWellLogChannelType::PRESSURE_P90 );
-    RifEclipseRftAddress meanAddress =
-        RifEclipseRftAddress::createAddress( wellName, timeStep, RifEclipseRftAddress::RftWellLogChannelType::PRESSURE_MEAN );
+    using ChannelType                 = RifEclipseRftAddress::RftWellLogChannelType;
+    RifEclipseRftAddress pressAddress = RifEclipseRftAddress::createAddress( wellName, timeStep, ChannelType::PRESSURE );
+    RifEclipseRftAddress tvdAddress   = RifEclipseRftAddress::createAddress( wellName, timeStep, ChannelType::TVD );
 
-    RiaCurveMerger<double> curveMerger;
-
-    RiaWeightedMeanCalculator<size_t> dataSetSizeCalc;
-
-    for ( RimSummaryCase* summaryCase : m_summaryCaseCollection->allSummaryCases() )
+    RigEclipseWellLogExtractor* extractor = RiaExtractionTools::findOrCreateWellLogExtractor( wellName, m_eclipseCase );
+    if ( extractor )
     {
-        RifReaderRftInterface* reader = summaryCase->rftReader();
-        if ( reader )
+        // Create a well log extractor if a well path and an Eclipse case is defined
+        // Use the extractor to compute measured depth for RFT cells
+        // The TVD values is extracted from the first summary case
+
+        RifEclipseRftAddress              mdAddress = RifEclipseRftAddress::createAddress( wellName, timeStep, ChannelType::MD );
+        RiaCurveMerger<double>            curveMerger;
+        RiaWeightedMeanCalculator<size_t> dataSetSizeCalc;
+        std::vector<double>               tvdDepthsForFirstRftCase;
+
+        for ( RimSummaryCase* summaryCase : m_summaryCaseCollection->allSummaryCases() )
         {
-            std::vector<double> depths;
+            auto reader = summaryCase->rftReader();
+            if ( !reader ) continue;
+
             std::vector<double> pressures;
-            reader->values( depthAddress, &depths );
             reader->values( pressAddress, &pressures );
-            if ( !depths.empty() && !pressures.empty() )
+
+            if ( tvdDepthsForFirstRftCase.empty() )
             {
-                dataSetSizeCalc.addValueAndWeight( depths.size(), 1.0 );
-                curveMerger.addCurveData( depths, pressures );
+                reader->values( tvdAddress, &tvdDepthsForFirstRftCase );
+            }
+
+            std::vector<double> measuredDepths = reader->computeMeasuredDepth( wellName, timeStep, extractor );
+
+            if ( !measuredDepths.empty() && !pressures.empty() )
+            {
+                dataSetSizeCalc.addValueAndWeight( measuredDepths.size(), 1.0 );
+                curveMerger.addCurveData( measuredDepths, pressures );
             }
         }
-    }
-    curveMerger.computeInterpolatedValues( false );
 
-    clearData( wellName, timeStep );
+        extractStatisticsFromCurveMerger( wellName, timeStep, mdAddress, curveMerger, dataSetSizeCalc );
+
+        if ( m_cachedValues[mdAddress].size() == tvdDepthsForFirstRftCase.size() )
+        {
+            // The number of RFT cells can vary between realizations in some cases. TVD depth is only given if the number of RFT cells is
+            // identical between realizations.
+
+            m_cachedValues[tvdAddress] = tvdDepthsForFirstRftCase;
+        }
+    }
+    else
+    {
+        // Compute statistics based on TVD depths. No measured depth can be estimated.
+        // This concept works well for vertical wells, but does not work for horizontal wells.
+
+        RiaCurveMerger<double>            curveMerger;
+        RiaWeightedMeanCalculator<size_t> dataSetSizeCalc;
+        RifEclipseRftAddress              tvdAddress = RifEclipseRftAddress::createAddress( wellName, timeStep, ChannelType::TVD );
+
+        for ( RimSummaryCase* summaryCase : m_summaryCaseCollection->allSummaryCases() )
+        {
+            auto reader = summaryCase->rftReader();
+            if ( !reader ) continue;
+
+            std::vector<double> pressures;
+            reader->values( pressAddress, &pressures );
+
+            std::vector<double> tvdDepths;
+            reader->values( tvdAddress, &tvdDepths );
+
+            if ( !tvdDepths.empty() && !pressures.empty() )
+            {
+                dataSetSizeCalc.addValueAndWeight( tvdDepths.size(), 1.0 );
+                curveMerger.addCurveData( tvdDepths, pressures );
+            }
+        }
+
+        extractStatisticsFromCurveMerger( wellName, timeStep, tvdAddress, curveMerger, dataSetSizeCalc );
+    }
+}
+
+//--------------------------------------------------------------------------------------------------
+/// Compute statistics for values, either based on measured depth or TVD
+//--------------------------------------------------------------------------------------------------
+void RifReaderEnsembleStatisticsRft::extractStatisticsFromCurveMerger( const QString&                     wellName,
+                                                                       const QDateTime&                   timeStep,
+                                                                       RifEclipseRftAddress               depthAddress,
+                                                                       RiaCurveMerger<double>&            curveMerger,
+                                                                       RiaWeightedMeanCalculator<size_t>& dataSetSizeCalc )
+{
+    using ChannelType = RifEclipseRftAddress::RftWellLogChannelType;
+
+    CAF_ASSERT( depthAddress.wellLogChannel() == ChannelType::MD || depthAddress.wellLogChannel() == ChannelType::TVD );
+
+    auto p10Address  = RifEclipseRftAddress::createAddress( wellName, timeStep, ChannelType::PRESSURE_P10 );
+    auto p50Address  = RifEclipseRftAddress::createAddress( wellName, timeStep, ChannelType::PRESSURE_P50 );
+    auto p90Address  = RifEclipseRftAddress::createAddress( wellName, timeStep, ChannelType::PRESSURE_P90 );
+    auto meanAddress = RifEclipseRftAddress::createAddress( wellName, timeStep, ChannelType::PRESSURE_MEAN );
+
+    clearCache( wellName, timeStep );
+
+    curveMerger.computeInterpolatedValues( false );
 
     const std::vector<double>& allDepths = curveMerger.allXValues();
     if ( !allDepths.empty() )
@@ -259,7 +340,7 @@ void RifReaderEnsembleStatisticsRft::calculateStatistics( const RifEclipseRftAdd
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-void RifReaderEnsembleStatisticsRft::clearData( const QString& wellName, const QDateTime& timeStep )
+void RifReaderEnsembleStatisticsRft::clearCache( const QString& wellName, const QDateTime& timeStep )
 {
     for ( auto it = m_cachedValues.begin(); it != m_cachedValues.end(); )
     {

ApplicationLibCode/FileInterface/RifReaderEnsembleStatisticsRft.h

@@ -20,14 +20,18 @@
 
 #include "RifReaderEclipseRft.h"
 
+#include "RiaCurveMerger.h"
+#include "RiaWeightedMeanCalculator.h"
+
 #include "cvfObject.h"
 
 class RimSummaryCaseCollection;
+class RimEclipseCase;
 
 class RifReaderEnsembleStatisticsRft : public RifReaderRftInterface, public cvf::Object
 {
 public:
-    RifReaderEnsembleStatisticsRft( const RimSummaryCaseCollection* summaryCaseCollection );
+    RifReaderEnsembleStatisticsRft( const RimSummaryCaseCollection* summaryCaseCollection, RimEclipseCase* eclipseCase );
 
     std::set<RifEclipseRftAddress> eclipseRftAddresses() override;
     void                           values( const RifEclipseRftAddress& rftAddress, std::vector<double>* values ) override;
@@ -41,11 +45,19 @@ public:
     std::set<QString>                                     wellNames() override;
 
 private:
-    void calculateStatistics( const RifEclipseRftAddress& rftAddress );
-    void clearData( const QString& wellName, const QDateTime& timeStep );
+    void calculateStatistics( const QString& wellName, const QDateTime& timeStep );
+
+    void extractStatisticsFromCurveMerger( const QString&                     wellName,
+                                           const QDateTime&                   timeStep,
+                                           RifEclipseRftAddress               depthAddress,
+                                           RiaCurveMerger<double>&            curveMerger,
+                                           RiaWeightedMeanCalculator<size_t>& dataSetSizeCalc );
+
+    void clearCache( const QString& wellName, const QDateTime& timeStep );
 
 private:
     const RimSummaryCaseCollection* m_summaryCaseCollection;
+    RimEclipseCase*                 m_eclipseCase;
 
     std::map<RifEclipseRftAddress, std::vector<double>> m_cachedValues;
 };

ApplicationLibCode/FileInterface/RifReaderOpmRft.cpp

@@ -272,40 +272,44 @@ std::set<QString> RifReaderOpmRft::wellNames()
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-void RifReaderOpmRft::cellIndices( const RifEclipseRftAddress& rftAddress, std::vector<caf::VecIjk>* indices )
+std::vector<caf::VecIjk> RifReaderOpmRft::cellIndices( const QString& wellName, const QDateTime& timeStep )
 {
-    if ( !openFiles() ) return;
+    if ( !openFiles() ) return {};
 
-    auto wellName = rftAddress.wellName().toStdString();
+    auto stdWellName = wellName.toStdString();
 
-    auto date = rftAddress.timeStep().date();
+    auto date = timeStep.date();
     int  y    = date.year();
     int  m    = date.month();
     int  d    = date.day();
 
+    std::vector<caf::VecIjk> indices;
+
     try
     {
         auto resultNameI = "CONIPOS";
-        auto dataI       = m_opm_rft->getRft<int>( resultNameI, wellName, y, m, d );
+        auto dataI       = m_opm_rft->getRft<int>( resultNameI, stdWellName, y, m, d );
 
         auto resultNameJ = "CONJPOS";
-        auto dataJ       = m_opm_rft->getRft<int>( resultNameJ, wellName, y, m, d );
+        auto dataJ       = m_opm_rft->getRft<int>( resultNameJ, stdWellName, y, m, d );
 
         auto resultNameK = "CONKPOS";
-        auto dataK       = m_opm_rft->getRft<int>( resultNameK, wellName, y, m, d );
+        auto dataK       = m_opm_rft->getRft<int>( resultNameK, stdWellName, y, m, d );
 
         if ( !dataI.empty() && ( dataI.size() == dataJ.size() ) && ( dataI.size() == dataK.size() ) )
         {
             for ( size_t n = 0; n < dataI.size(); n++ )
             {
                 // NB: Transform to zero-based cell indices
-                indices->push_back( caf::VecIjk( dataI[n] - 1, dataJ[n] - 1, dataK[n] - 1 ) );
+                indices.push_back( caf::VecIjk( dataI[n] - 1, dataJ[n] - 1, dataK[n] - 1 ) );
             }
         }
     }
     catch ( ... )
     {
     }
+
+    return indices;
 }
 
 //--------------------------------------------------------------------------------------------------

ApplicationLibCode/FileInterface/RifReaderOpmRft.h

@@ -50,7 +50,7 @@ public:
     std::set<RifEclipseRftAddress::RftWellLogChannelType> availableWellLogChannels( const QString& wellName ) override;
     std::set<QString>                                     wellNames() override;
 
-    void cellIndices( const RifEclipseRftAddress& rftAddress, std::vector<caf::VecIjk>* indices ) override;
+    std::vector<caf::VecIjk> cellIndices( const QString& wellName, const QDateTime& timeStep ) override;
 
     std::map<int, int> branchIdsAndOneBasedIndices( const QString& wellName, const QDateTime& timeStep, RiaDefines::RftBranchType branchType );
 

ApplicationLibCode/FileInterface/RifReaderRftInterface.cpp

@@ -15,8 +15,17 @@
 //  for more details.
 //
 /////////////////////////////////////////////////////////////////////////////////
+
 #include "RifReaderRftInterface.h"
 
+#include "RigEclipseCaseData.h"
+#include "RigEclipseWellLogExtractor.h"
+#include "RigMainGrid.h"
+#include "RigWellPath.h"
+#include "RigWellPathGeometryTools.h"
+
+#include "cafVecIjk.h"
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
@@ -34,9 +43,104 @@ std::set<RifEclipseRftAddress> RifReaderRftInterface::eclipseRftAddresses( const
     return matchingAddresses;
 }
 
+//--------------------------------------------------------------------------------------------------
+// Compute average measured depth for cell based on grid intersections for cells. If the well path geometry do not contain measured depth
+// for a grid cell, the measured depth is estimated based on existing geometry for the well path.
+//--------------------------------------------------------------------------------------------------
+std::vector<double>
+    RifReaderRftInterface::computeMeasuredDepth( const QString& wellName, const QDateTime& timeStep, RigEclipseWellLogExtractor* eclExtractor )
+{
+    if ( !eclExtractor ) return {};
+    if ( !eclExtractor->caseData() ) return {};
+
+    if ( eclExtractor->cellIntersectionMDs().empty() ) return {};
+
+    auto mainGrid = eclExtractor->caseData()->mainGrid();
+    if ( !mainGrid ) return {};
+
+    std::vector<double> avgMeasuredDepthForCells;
+    std::vector<double> tvdValuesToEstimate;
+
+    auto cellIjk = cellIndices( wellName, timeStep );
+    for ( const caf::VecIjk& ijk : cellIjk )
+    {
+        auto globalCellIndex = mainGrid->cellIndexFromIJK( ijk.i(), ijk.j(), ijk.k() );
+
+        auto avgMd = eclExtractor->averageMdForCell( globalCellIndex );
+        if ( avgMd.has_value() )
+        {
+            avgMeasuredDepthForCells.push_back( avgMd.value() );
+        }
+        else
+        {
+            // The RFT cell is not part of cells intersected by well path
+            // Use the TVD of cell center to estimate measured depth
+
+            avgMeasuredDepthForCells.push_back( std::numeric_limits<double>::infinity() );
+            auto center = mainGrid->cell( globalCellIndex ).center();
+            tvdValuesToEstimate.push_back( -center.z() );
+        }
+    }
+
+    if ( !tvdValuesToEstimate.empty() )
+    {
+        auto wellPathMd  = eclExtractor->wellPathGeometry()->measuredDepths();
+        auto wellPathTvd = eclExtractor->wellPathGeometry()->trueVerticalDepths();
+
+        // Estimate measured depth for cells that do not have  measured depth
+        auto estimatedMeasuredDepth = RigWellPathGeometryTools::interpolateMdFromTvd( wellPathMd, wellPathTvd, tvdValuesToEstimate );
+
+        double previousMd = std::numeric_limits<double>::infinity();
+
+        // Replace infinity MD values with estimated MD values based on well path geometry
+        // previousMd contains the last known measured depth value as we move along the well path
+
+        size_t estimatedIndex = 0;
+        for ( auto& measuredDepth : avgMeasuredDepthForCells )
+        {
+            if ( measuredDepth == std::numeric_limits<double>::infinity() )
+            {
+                // No measured depth for cell is found, try to estimate MD based on MD for previous cell
+                if ( estimatedIndex < estimatedMeasuredDepth.size() )
+                {
+                    double estimatedMd = estimatedMeasuredDepth[estimatedIndex++];
+                    if ( previousMd != std::numeric_limits<double>::infinity() )
+                    {
+                        if ( previousMd < estimatedMd )
+                        {
+                            // The estimated MD is larger than previous MD, use the estimated MD
+                            measuredDepth = estimatedMd;
+                        }
+                        else
+                        {
+                            // The estimated MD is smaller than previous MD, use the previous MD + 1.0
+                            measuredDepth = previousMd + 1.0;
+                        }
+                    }
+                    else
+                    {
+                        // We do not have a valid previous MD, use the estimated MD
+                        measuredDepth = estimatedMd;
+                    }
+                }
+                else
+                {
+                    measuredDepth = previousMd + 1.0;
+                }
+            }
+
+            // Assign the current MD as previous MD to be used for next estimated MD
+            previousMd = measuredDepth;
+        }
+    }
+
+    return avgMeasuredDepthForCells;
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-void RifReaderRftInterface::cellIndices( const RifEclipseRftAddress& rftAddress, std::vector<caf::VecIjk>* indices )
+std::vector<caf::VecIjk> RifReaderRftInterface::cellIndices( const QString& wellName, const QDateTime& timeStep )
 {
+    return {};
 }

ApplicationLibCode/FileInterface/RifReaderRftInterface.h

@@ -31,6 +31,8 @@ namespace caf
 class VecIjk;
 };
 
+class RigEclipseWellLogExtractor;
+
 class RifReaderRftInterface
 {
 public:
@@ -46,5 +48,9 @@ public:
     virtual std::set<RifEclipseRftAddress::RftWellLogChannelType> availableWellLogChannels( const QString& wellName ) = 0;
     virtual std::set<QString>                                     wellNames()                                         = 0;
 
-    virtual void cellIndices( const RifEclipseRftAddress& rftAddress, std::vector<caf::VecIjk>* indices );
+    // To be moved into Rig structures
+    std::vector<double> computeMeasuredDepth( const QString& wellName, const QDateTime& timeStep, RigEclipseWellLogExtractor* extractor );
+
+    // TODO: Move to protected or private
+    virtual std::vector<caf::VecIjk> cellIndices( const QString& wellName, const QDateTime& timeStep );
 };