mirror of
https://github.com/OPM/ResInsight.git
synced 2025-02-25 18:55:39 -06:00
Create depth adjusted LAS files
* Add support for creating depth adjustment LAS files * Add RiaDefine for LAS file depth property names * Remove incorrect check in K-index Calculator Co-authored-by: jorgenherje <jorgenherje@users.noreply.github.com> Co-authored-by: Magne Sjaastad <magne.sjaastad@ceetronsolutions.com>
This commit is contained in:
Jørgen Herje
committed by
Magne Sjaastad
Magne Sjaastad
parent
264fbbfadc
commit
d634fcd4b4
@@ -0,0 +1,411 @@
|
||||
#include "RicCreateDepthAdjustedLasFilesFeature.h"
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2023- Equinor ASA
|
||||
//
|
||||
// ResInsight is free software: you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation, either version 3 of the License, or
|
||||
// (at your option) any later version.
|
||||
//
|
||||
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
// FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
|
||||
// for more details.
|
||||
//
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
#include "RicCreateDepthAdjustedLasFilesImpl.h"
|
||||
|
||||
#include "RiaLogging.h"
|
||||
|
||||
#include "RicCreateDepthAdjustedLasFilesUi.h"
|
||||
|
||||
#include "RigCaseCellResultsData.h"
|
||||
#include "RigEclipseCaseData.h"
|
||||
#include "RigEclipseResultAddress.h"
|
||||
#include "RigEclipseWellLogExtractor.h"
|
||||
#include "RigGeoMechWellLogExtractor.h"
|
||||
#include "RigResultAccessorFactory.h"
|
||||
|
||||
#include "RimEclipseCase.h"
|
||||
#include "RimGeoMechCase.h"
|
||||
#include "RimMainPlotCollection.h"
|
||||
#include "RimWellLogFile.h"
|
||||
#include "RimWellLogPlotCollection.h"
|
||||
#include "RimWellPath.h"
|
||||
|
||||
#include "NRLib/nrlib/well/laswell.hpp"
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void LasDepthValueAndIndexPerKLayer::insertIndexAndValue( int kLayer, size_t index, double value )
|
||||
{
|
||||
m_kLayerIndexAndValuePairsMap[kLayer][index] = value;
|
||||
}
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
bool LasDepthValueAndIndexPerKLayer::hasKLayer( int kLayer ) const
|
||||
{
|
||||
return m_kLayerIndexAndValuePairsMap.find( kLayer ) != m_kLayerIndexAndValuePairsMap.end();
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::map<size_t, double> LasDepthValueAndIndexPerKLayer::indexAndValuePairs( int kLayer ) const
|
||||
{
|
||||
if ( !hasKLayer( kLayer ) ) return std::map<size_t, double>();
|
||||
|
||||
return m_kLayerIndexAndValuePairsMap.at( kLayer );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
cvf::ref<RigResultAccessor> RicCreateDepthAdjustedLasFilesImpl::createIndexKResultAccessor( RimEclipseCase* eclipseCase )
|
||||
{
|
||||
const int firstTimeStep = 0;
|
||||
const int gridIndex = 0;
|
||||
RigEclipseResultAddress indexKResAdr( RiaDefines::ResultCatType::STATIC_NATIVE, RiaResultNames::indexKResultName() );
|
||||
eclipseCase->eclipseCaseData()->results( RiaDefines::PorosityModelType::MATRIX_MODEL )->ensureKnownResultLoaded( indexKResAdr );
|
||||
|
||||
return RigResultAccessorFactory::createFromResultAddress( eclipseCase->eclipseCaseData(),
|
||||
gridIndex,
|
||||
RiaDefines::PorosityModelType::MATRIX_MODEL,
|
||||
firstTimeStep,
|
||||
indexKResAdr );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
LasDepthValueAndIndexPerKLayer RicCreateDepthAdjustedLasFilesImpl::createLasDepthIndexAndPercValuePerKLayerFromMap(
|
||||
const std::vector<double>& lasWellDepths,
|
||||
const std::map<int, IndexKDepthData>& indexKDepthDataMap )
|
||||
{
|
||||
// Create container of depth value (in percent) and its original index in a LAS file vector
|
||||
// categorized by K-layer. Depth value as percentage value between MD top and MD bottom for K-layer.
|
||||
auto lasWellDepthValueAndIndexPerKLayer = LasDepthValueAndIndexPerKLayer();
|
||||
for ( size_t i = 0; i < lasWellDepths.size(); ++i )
|
||||
{
|
||||
const double depth = lasWellDepths[i];
|
||||
for ( const auto& [indexK, depthData] : indexKDepthDataMap )
|
||||
{
|
||||
if ( depthData.mdTop <= depth && depth <= depthData.mdBottom )
|
||||
{
|
||||
const double percentage = ( depth - depthData.mdTop ) / ( depthData.mdBottom - depthData.mdTop );
|
||||
lasWellDepthValueAndIndexPerKLayer.insertIndexAndValue( indexK, i, percentage );
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
return lasWellDepthValueAndIndexPerKLayer;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// NOTE: map createIndexKDepthDataMapFromCase is created using well extractor, while sourceWellLogData depth
|
||||
/// values are from LAS file. Floating point rounding in LAS file can occur, thus depth values might be placed
|
||||
/// outside of K-layer close to top/bottom due to inaccuracy.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RicCreateDepthAdjustedLasFilesImpl::createDestinationWellsLasFiles( RimCase* selectedCase,
|
||||
RimWellPath* sourceWell,
|
||||
const std::vector<RimWellPath*> destinationWells,
|
||||
const std::vector<QString>& selectedResultProperties,
|
||||
const QString& exportFolder,
|
||||
double rkbDiff )
|
||||
{
|
||||
auto* sourceWellLogData = sourceWell->wellLogFiles()[0]->wellLogFileData();
|
||||
const auto defaultPropertyMap = createDefaultPropertyMap( selectedResultProperties, sourceWellLogData );
|
||||
|
||||
// NOTE: map createIndexKDepthDataMapFromCase is created using well extractor, while sourceWellLogData depth
|
||||
// values are from LAS file. Floating point rounding in LAS file can occur, thus depth values might be placed
|
||||
// outside of K-layer close to top/bottom due to inaccuracy.
|
||||
const auto sourceWellDepthIndexAndPercValuePerKLayer =
|
||||
createLasDepthIndexAndPercValuePerKLayerFromMap( sourceWellLogData->depthValues(),
|
||||
createIndexKDepthDataMapFromCase( selectedCase, sourceWell ) );
|
||||
for ( RimWellPath* well : destinationWells )
|
||||
{
|
||||
const std::map<int, IndexKDepthData> destinationWellIndexKDepthsMap =
|
||||
createIndexKDepthDataMapFromCase( selectedCase, well );
|
||||
|
||||
if ( destinationWellIndexKDepthsMap.empty() ) continue;
|
||||
|
||||
std::vector<double> mdValues;
|
||||
std::vector<double> tvdMslValues;
|
||||
std::vector<double> tvdRkbValues;
|
||||
std::map<QString, std::vector<double>> propertyMap = defaultPropertyMap;
|
||||
for ( const auto& [indexK, depthData] : destinationWellIndexKDepthsMap )
|
||||
{
|
||||
if ( !sourceWellDepthIndexAndPercValuePerKLayer.hasKLayer( indexK ) ) continue;
|
||||
|
||||
for ( const auto& [index, depthPerc] : sourceWellDepthIndexAndPercValuePerKLayer.indexAndValuePairs( indexK ) )
|
||||
{
|
||||
if ( sourceWellLogData->hasTvdMslChannel() )
|
||||
{
|
||||
const double tvdMslValue = depthPerc * ( depthData.tvdBottom - depthData.tvdTop ) + depthData.tvdTop;
|
||||
tvdMslValues.push_back( tvdMslValue );
|
||||
}
|
||||
if ( sourceWellLogData->hasTvdRkbChannel() )
|
||||
{
|
||||
const double tvdRkbValue = depthPerc * ( depthData.tvdBottom - depthData.tvdTop ) +
|
||||
depthData.tvdTop + rkbDiff;
|
||||
tvdRkbValues.push_back( tvdRkbValue );
|
||||
}
|
||||
|
||||
const double mdValue = depthPerc * ( depthData.mdBottom - depthData.mdTop ) + depthData.mdTop;
|
||||
mdValues.push_back( mdValue );
|
||||
for ( auto& [propertyName, values] : propertyMap )
|
||||
{
|
||||
double value = sourceWellLogData->values( propertyName )[index];
|
||||
value = value == HUGE_VAL ? sourceWellLogData->getMissingValue() : value;
|
||||
values.push_back( value );
|
||||
}
|
||||
}
|
||||
}
|
||||
createDestinationWellLasFile( well->name(),
|
||||
selectedCase->caseUserDescription(),
|
||||
mdValues,
|
||||
tvdMslValues,
|
||||
tvdRkbValues,
|
||||
propertyMap,
|
||||
sourceWellLogData,
|
||||
exportFolder );
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RicCreateDepthAdjustedLasFilesImpl::createDestinationWellLasFile( const QString& wellName,
|
||||
const QString& caseDescription,
|
||||
const std::vector<double>& mdValues,
|
||||
const std::vector<double>& tvdMslValues,
|
||||
const std::vector<double>& tvdRkbValues,
|
||||
const std::map<QString, std::vector<double>>& propertyMap,
|
||||
const RigWellLogFile* sourceWellLogData,
|
||||
const QString& exportFolder )
|
||||
{
|
||||
const auto depthUnitText = createDepthUnitText( sourceWellLogData->depthUnit() );
|
||||
const auto depthUnitComment = createDepthUnitComment( sourceWellLogData->depthUnit() );
|
||||
|
||||
const auto deptUnit = sourceWellLogData->depthUnit();
|
||||
|
||||
// Build LAS file
|
||||
NRLib::LasWell lasFile;
|
||||
|
||||
lasFile.setVersionInfo( "2.0" );
|
||||
lasFile.setDepthUnit( depthUnitText );
|
||||
lasFile.SetMissing( sourceWellLogData->getMissingValue() );
|
||||
lasFile.setStartDepth( *std::min_element( mdValues.begin(), mdValues.end() ) );
|
||||
lasFile.setStopDepth( *std::max_element( mdValues.begin(), mdValues.end() ) );
|
||||
|
||||
lasFile.addWellInfo( "WELL", wellName.toStdString() );
|
||||
lasFile.addWellInfo( "DATE", sourceWellLogData->date().toStdString() );
|
||||
|
||||
// Add Measured depth
|
||||
lasFile.AddLog( RiaDefines::propertyNameMeasuredDepth().toStdString(), depthUnitText, "Depth " + depthUnitComment, mdValues );
|
||||
|
||||
// Add tvd msl values if existing
|
||||
if ( !tvdMslValues.empty() )
|
||||
{
|
||||
const auto unitText =
|
||||
sourceWellLogData->wellLogChannelUnitString( RiaDefines::propertyNameTvdMslDepth(), deptUnit ).toStdString();
|
||||
lasFile.AddLog( RiaDefines::propertyNameTvdMslDepth().toStdString(),
|
||||
unitText,
|
||||
"True vertical depth " + depthUnitComment,
|
||||
tvdMslValues );
|
||||
}
|
||||
|
||||
// Add tvd rkb values if existing
|
||||
if ( !tvdRkbValues.empty() )
|
||||
{
|
||||
const auto unitText =
|
||||
sourceWellLogData->wellLogChannelUnitString( RiaDefines::propertyNameTvdRkbDepth(), deptUnit ).toStdString();
|
||||
lasFile.AddLog( RiaDefines::propertyNameTvdRkbDepth().toStdString(),
|
||||
unitText,
|
||||
"True vertical depth (Rotary Kelly Bushing)",
|
||||
tvdRkbValues );
|
||||
}
|
||||
|
||||
// Add property values
|
||||
for ( auto& [name, values] : propertyMap )
|
||||
{
|
||||
std::string unitText = sourceWellLogData->wellLogChannelUnitString( name ).toStdString();
|
||||
lasFile.AddLog( name.toUpper().toStdString(), unitText, "", values );
|
||||
}
|
||||
|
||||
std::vector<std::string> commentHeader;
|
||||
QString fullPathName = exportFolder + "/" + wellName + "_" + caseDescription + "_" + sourceWellLogData->date() + ".las";
|
||||
lasFile.WriteToFile( fullPathName.toStdString(), commentHeader );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::string RicCreateDepthAdjustedLasFilesImpl::createDepthUnitText( RiaDefines::DepthUnitType depthUnitType )
|
||||
{
|
||||
return depthUnitType == RiaDefines::DepthUnitType::UNIT_METER
|
||||
? "M"
|
||||
: depthUnitType == RiaDefines::DepthUnitType::UNIT_FEET ? "FT" : "";
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::string RicCreateDepthAdjustedLasFilesImpl::createDepthUnitComment( RiaDefines::DepthUnitType depthUnitType )
|
||||
{
|
||||
return depthUnitType == RiaDefines::DepthUnitType::UNIT_METER
|
||||
? "in meters"
|
||||
: depthUnitType == RiaDefines::DepthUnitType::UNIT_FEET ? "in feet" : "in Connection number";
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::map<int, RicCreateDepthAdjustedLasFilesImpl::IndexKDepthData>
|
||||
RicCreateDepthAdjustedLasFilesImpl::createIndexKDepthDataMapFromCase( RimCase* selectedCase, RimWellPath* wellPath )
|
||||
{
|
||||
RimEclipseCase* eclipseCase = dynamic_cast<RimEclipseCase*>( selectedCase );
|
||||
RimGeoMechCase* geomCase = dynamic_cast<RimGeoMechCase*>( selectedCase );
|
||||
|
||||
RimWellLogPlotCollection* wellLogCollection = RimMainPlotCollection::current()->wellLogPlotCollection();
|
||||
if ( eclipseCase != nullptr )
|
||||
{
|
||||
cvf::ref<RigEclipseWellLogExtractor> wellExtractor =
|
||||
wellLogCollection->findOrCreateExtractor( wellPath, eclipseCase );
|
||||
if ( wellExtractor.isNull() )
|
||||
{
|
||||
RiaLogging::info( QString( "Could not create RigEclipseWellLogExtractor for %1" ).arg( wellPath->name() ) );
|
||||
}
|
||||
const auto result = createIndexKDepthDataMap( wellExtractor, createIndexKResultAccessor( eclipseCase ) );
|
||||
if ( result.empty() )
|
||||
{
|
||||
RiaLogging::info( QString( "Not able to create Index-K depth map for %1" ).arg( wellPath->name() ) );
|
||||
}
|
||||
return result;
|
||||
}
|
||||
else if ( geomCase != nullptr )
|
||||
{
|
||||
cvf::ref<RigGeoMechWellLogExtractor> wellExtractor = wellLogCollection->findOrCreateExtractor( wellPath, geomCase );
|
||||
if ( wellExtractor.isNull() )
|
||||
{
|
||||
RiaLogging::info( QString( "Could not create RigGeoMechWellLogExtractor for %1" ).arg( wellPath->name() ) );
|
||||
}
|
||||
const auto result = createIndexKDepthDataMap( wellExtractor );
|
||||
if ( result.empty() )
|
||||
{
|
||||
RiaLogging::info( QString( "Not able to create Index-K depth map for %1" ).arg( wellPath->name() ) );
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
RiaLogging::info( QString( "Invalid case when creating Index-K depth map for %1" ).arg( wellPath->name() ) );
|
||||
return std::map<int, RicCreateDepthAdjustedLasFilesImpl::IndexKDepthData>();
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::map<int, RicCreateDepthAdjustedLasFilesImpl::IndexKDepthData>
|
||||
RicCreateDepthAdjustedLasFilesImpl::createIndexKDepthDataMap( cvf::ref<RigEclipseWellLogExtractor> wellExtractor,
|
||||
cvf::ref<RigResultAccessor> indexKResAcc )
|
||||
{
|
||||
std::vector<double> wellIndexKValues;
|
||||
wellExtractor->curveData( indexKResAcc.p(), &wellIndexKValues );
|
||||
return createIndexKDepthDataMapFromVectors( wellExtractor->cellIntersectionMDs(),
|
||||
wellExtractor->cellIntersectionTVDs(),
|
||||
wellIndexKValues );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::map<int, RicCreateDepthAdjustedLasFilesImpl::IndexKDepthData>
|
||||
RicCreateDepthAdjustedLasFilesImpl::createIndexKDepthDataMap( cvf::ref<RigGeoMechWellLogExtractor> wellExtractor )
|
||||
{
|
||||
const int frameIdx = 0;
|
||||
const int timeStepIdx = 0;
|
||||
RigFemResultAddress indexKResAdr( RigFemResultPosEnum::RIG_ELEMENT_NODAL, "INDEX", "INDEX_K" );
|
||||
std::vector<double> wellIndexKValues;
|
||||
wellExtractor->curveData( indexKResAdr, timeStepIdx, frameIdx, &wellIndexKValues );
|
||||
return createIndexKDepthDataMapFromVectors( wellExtractor->cellIntersectionMDs(),
|
||||
wellExtractor->cellIntersectionTVDs(),
|
||||
wellIndexKValues );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::map<int, RicCreateDepthAdjustedLasFilesImpl::IndexKDepthData>
|
||||
RicCreateDepthAdjustedLasFilesImpl::createIndexKDepthDataMapFromVectors( const std::vector<double>& wellMdValues,
|
||||
const std::vector<double>& wellTvdValues,
|
||||
const std::vector<double>& wellIndexKValues )
|
||||
{
|
||||
std::map<int, IndexKDepthData> indexKDepthsMap;
|
||||
|
||||
// Must have non-empty equal length vectors!
|
||||
if ( wellIndexKValues.empty() )
|
||||
{
|
||||
RiaLogging::info( QString( "Empty vector of index-K values" ) );
|
||||
return indexKDepthsMap;
|
||||
}
|
||||
if ( wellMdValues.size() != wellTvdValues.size() || wellMdValues.size() != wellIndexKValues.size() )
|
||||
{
|
||||
return indexKDepthsMap;
|
||||
}
|
||||
|
||||
int prevKLayer = -1;
|
||||
for ( size_t i = 0; i < wellIndexKValues.size(); ++i )
|
||||
{
|
||||
// Asymptotically increasing k-indexes!
|
||||
const auto kLayer = static_cast<int>( wellIndexKValues[i] );
|
||||
if ( kLayer < prevKLayer ) break;
|
||||
|
||||
if ( indexKDepthsMap.find( kLayer ) == indexKDepthsMap.end() )
|
||||
{
|
||||
indexKDepthsMap[kLayer] = IndexKDepthData();
|
||||
indexKDepthsMap[kLayer].mdTop = wellMdValues[i];
|
||||
indexKDepthsMap[kLayer].mdBottom = wellMdValues[i];
|
||||
indexKDepthsMap[kLayer].tvdTop = wellTvdValues[i];
|
||||
indexKDepthsMap[kLayer].tvdBottom = wellTvdValues[i];
|
||||
}
|
||||
else
|
||||
{
|
||||
indexKDepthsMap[kLayer].mdBottom = wellMdValues[i];
|
||||
indexKDepthsMap[kLayer].tvdBottom = wellTvdValues[i];
|
||||
}
|
||||
}
|
||||
|
||||
return indexKDepthsMap;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::map<QString, std::vector<double>>
|
||||
RicCreateDepthAdjustedLasFilesImpl::createDefaultPropertyMap( const std::vector<QString>& selectedProperties,
|
||||
const RigWellLogFile* wellLogFile )
|
||||
{
|
||||
const QStringList lasDepthNames = QStringList( { RiaDefines::propertyNameMeasuredDepth(),
|
||||
RiaDefines::propertyNameTvdMslDepth(),
|
||||
RiaDefines::propertyNameTvdRkbDepth() } );
|
||||
std::vector<QString> validPropertyNames;
|
||||
for ( const auto& propertyName : selectedProperties )
|
||||
{
|
||||
if ( !lasDepthNames.contains( propertyName ) && wellLogFile->wellLogChannelNames().contains( propertyName ) )
|
||||
{
|
||||
validPropertyNames.push_back( propertyName );
|
||||
}
|
||||
}
|
||||
std::map<QString, std::vector<double>> defaultPropertyMap;
|
||||
for ( const auto& name : validPropertyNames )
|
||||
{
|
||||
defaultPropertyMap[name];
|
||||
}
|
||||
return defaultPropertyMap;
|
||||
}
|
||||
Reference in New Issue
Block a user