Commit f3e847b5 authored by Cédric Traizet's avatar Cédric Traizet
Browse files

REFAC: create a in-memory vrt geoid+DEM file to manage elevation in GDALRPCtransformer

parent 09ea983e
......@@ -143,6 +143,8 @@ public:
/** Path to the in-memory vrt */
const std::string DEM_DATASET_PATH = "/vsimem/otb_dem_dataset.vrt";
const std::string DEM_WARPED_DATASET_PATH = "/vsimem/otb_dem_warped_dataset.vrt";
const std::string DEM_SHIFTED_DATASET_PATH = "/vsimem/otb_dem_shifted_dataset.vrt";
protected:
DEMHandler();
......@@ -152,7 +154,10 @@ protected:
private:
DEMHandler(const Self&) = delete;
void operator=(const Self&) = delete;
void operator=(const Self&) = delete;
void CreateShiftedDataset();
/** List of RAII capsules on all opened DEM datasets for memory management */
std::vector<otb::GDALDatasetWrapper::Pointer> m_DatasetList;
......@@ -162,6 +167,9 @@ private:
/** Pointer to the geoid dataset */
GDALDataset* m_GeoidDS;
/** Pointer to the sifted dataset */
GDALDataset* m_ShiftedDS;
/** Default height above elliposid, used when no DEM or geoid height is available. */
double m_DefaultHeightAboveEllipsoid;
......
......@@ -24,6 +24,10 @@
#include <boost/range/iterator_range.hpp>
#include "gdal_utils.h"
//Warp utility
#include "gdalwarper.h"
#include "vrtdataset.h"
//TODO C++ 17 : use std::optional instead
#include <boost/optional.hpp>
......@@ -181,7 +185,10 @@ DEMHandler::~DEMHandler()
}
ClearDEMs();
VSIUnlink(DEM_DATASET_PATH.c_str());
VSIUnlink(DEM_WARPED_DATASET_PATH.c_str());
VSIUnlink(DEM_SHIFTED_DATASET_PATH.c_str());
}
void DEMHandler::OpenDEMFile(const std::string& path)
......@@ -198,7 +205,7 @@ void DEMHandler::OpenDEMFile(const std::string& path)
}
void DEMHandler::OpenDEMDirectory(const std::string& DEMDirectory)
{
{
// TODO : RemoveOSSIM
OssimDEMHandler::Instance()->OpenDEMDirectory(DEMDirectory);
......@@ -238,10 +245,15 @@ void DEMHandler::OpenDEMDirectory(const std::string& DEMDirectory)
m_Dataset = static_cast<GDALDataset*>(GDALOpen(DEM_DATASET_PATH.c_str(), GA_ReadOnly));
m_DEMDirectories.push_back(DEMDirectory);
}
if(m_GeoidDS)
{
CreateShiftedDataset();
}
}
bool DEMHandler::OpenGeoidFile(const std::string& geoidFile)
{
{
// TODO : RemoveOSSIM
OssimDEMHandler::Instance()->OpenGeoidFile(geoidFile);
......@@ -260,9 +272,110 @@ bool DEMHandler::OpenGeoidFile(const std::string& geoidFile)
m_GeoidFilename = geoidFile;
}
if(m_Dataset)
{
CreateShiftedDataset();
}
return pbError;
}
void DEMHandler::CreateShiftedDataset()
{
// WIP : no data is not handed at the moment
double geoTransform[6];
m_Dataset->GetGeoTransform(geoTransform);
// Warp geoid dataset
auto warpOptions = GDALCreateWarpOptions();;
warpOptions->hSrcDS = m_GeoidDS;
//warpOptions->hDstDS = m_Dataset;
warpOptions->eResampleAlg = GRA_Bilinear;
warpOptions->eWorkingDataType = GDT_Float64;
warpOptions->nBandCount = 1;
warpOptions->panSrcBands =
(int *) CPLMalloc(sizeof(int) * warpOptions->nBandCount );
warpOptions->panSrcBands[0] = 1;
warpOptions->panDstBands =
(int *) CPLMalloc(sizeof(int) * warpOptions->nBandCount );
warpOptions->panDstBands[0] = 1;
// Establish reprojection transformer.
warpOptions->pTransformerArg =
GDALCreateGenImgProjTransformer( m_GeoidDS,
GDALGetProjectionRef(m_GeoidDS),
m_Dataset,
GDALGetProjectionRef(m_Dataset),
FALSE, 0.0, 1 );
warpOptions->pfnTransformer = GDALGenImgProjTransform;
auto ds = static_cast<GDALDataset *> (GDALCreateWarpedVRT(m_GeoidDS,
m_Dataset->GetRasterXSize(),
m_Dataset->GetRasterYSize(),
geoTransform,
warpOptions));
/*
auto warpedDataset = dynamic_cast<VRTDataset*>(ds);
if(warpedDataset)
{
auto xmlWarpedDs= CPLSerializeXMLTree(warpedDataset->SerializeToXML(nullptr));
std::cout << xmlWarpedDs << std::endl;
}
*/
ds->SetDescription(DEM_WARPED_DATASET_PATH.c_str());
GDALClose(ds);
GDALDriver *poDriver = (GDALDriver *) GDALGetDriverByName( "VRT" );
GDALDataset *poVRTDS;
poVRTDS = poDriver->Create( DEM_SHIFTED_DATASET_PATH.c_str(), m_Dataset->GetRasterXSize(), m_Dataset->GetRasterYSize(), 0, GDT_Float64, NULL );
poVRTDS->SetGeoTransform(geoTransform);
poVRTDS->SetSpatialRef(m_Dataset->GetSpatialRef());
char** derivedBandOptions = NULL;
derivedBandOptions = CSLAddNameValue(derivedBandOptions, "subclass", "VRTDerivedRasterBand");
derivedBandOptions = CSLAddNameValue(derivedBandOptions, "PixelFunctionType", "sum");
poVRTDS->AddBand(GDT_Float64, derivedBandOptions);
GDALRasterBand *poBand = poVRTDS->GetRasterBand( 1 );
// TODO use std string (and boost format ?) or stream
char demVrtXml[10000];
sprintf( demVrtXml,
"<SimpleSource>"
" <SourceFilename>%s</SourceFilename>"
" <SourceBand>1</SourceBand>"
"</SimpleSource>",
DEM_DATASET_PATH.c_str());
poBand->SetMetadataItem( "source_0", demVrtXml, "new_vrt_sources" );
char geoidVrtXml[10000];
sprintf( geoidVrtXml,
"<SimpleSource>"
" <SourceFilename>%s</SourceFilename>"
" <SourceBand>1</SourceBand>"
"</SimpleSource>",
DEM_WARPED_DATASET_PATH.c_str());
poBand->SetMetadataItem( "source_1", geoidVrtXml, "new_vrt_sources" );
GDALClose(poVRTDS);
}
double DEMHandler::GetHeightAboveEllipsoid(double lon, double lat) const
{
double result = 0.;
......
......@@ -54,9 +54,18 @@ GDALRPCTransformer::GDALRPCTransformer(double LineOffset, double SampleOffset, d
this->m_GDALRPCInfo.dfMAX_LAT = 90.0;
auto & demHandler = otb::DEMHandler::GetInstance();
if (demHandler.GetDEMCount() > 0)
{
this->SetOption("RPC_DEM", demHandler.DEM_DATASET_PATH);
if (demHandler.GetGeoidFile().empty())
{
this->SetOption("RPC_DEM", demHandler.DEM_DATASET_PATH);
}
else
{
this->SetOption("RPC_DEM", demHandler.DEM_SHIFTED_DATASET_PATH);
}
this->SetOption("RPC_DEM_MISSING_VALUE", std::to_string(demHandler.GetDefaultHeightAboveEllipsoid()));
}
else
......
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