Commit 023b7f64 authored by Guillaume Pasero's avatar Guillaume Pasero

Merge branch '1937-keep-gcp' into 'develop'

Resolve "ExtractROI: keep GCPs"

Closes #1937

See merge request orfeotoolbox/otb!642
parents 43f4def1 5774cb13
Pipeline #3147 passed with stage
in 10 minutes and 53 seconds
......@@ -203,6 +203,20 @@ otb_test_application(NAME apTvUtQuicklookROI1Channel
${TEMP}/apTvUtQuicklookROI1Channel.tif
)
otb_test_application(NAME apTvUtQuicklookWithGCP
APP Quicklook
OPTIONS -in ${INPUTDATA}/spot5SubWithGcps.tif
-out ${TEMP}/apTvUtQuicklookWithGCP.tif uint8
-sr 4
-rox 100
-roy 100
-rsx 200
-rsy 200
VALID --compare-metadata ${NOTOL}
${BASELINE}/apTvUtQuicklookWithGCP.tif
${TEMP}/apTvUtQuicklookWithGCP.tif
)
otb_test_application(NAME apTvUtQuicklookSpot5
APP Quicklook
OPTIONS -in LARGEINPUT{SPOT5/TEHERAN/IMAGERY.TIF}
......
......@@ -30,8 +30,6 @@
#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
#include "otb_boost_string_header.h"
int otbExtractROITestMetaData(int itkNotUsed(argc), char* argv[])
{
typedef float PixelType;
......@@ -100,17 +98,20 @@ int otbExtractROITestMetaData(int itkNotUsed(argc), char* argv[])
reader00->SetFileName(argv[2]);
reader00->GenerateOutputInformation();
if (reader00->GetOutput()->GetProjectionRef() != "" || boost::algorithm::istarts_with(reader00->GetOutput()->GetProjectionRef(), "LOCAL_CS"))
// The input image should have a sensor model and GCP, but the output images
// should only have the sensor model (priority over GCP). This behaviour
// must be consistent regardless of the ROI.
if (reader00->GetOutput()->GetProjectionRef().size())
{
std::cout << "The read generated extract from index (0, 0) must NOT contain a ProjectionReference." << std::endl;
std::cout << "Found ProjectionReference: " << reader00->GetOutput()->GetProjectionRef() << std::endl;
return EXIT_FAILURE;
}
if (reader00->GetOutput()->GetGCPCount() == 0)
if (reader00->GetOutput()->GetGCPCount())
{
std::cout << "The read generated extract from index (0, 0) must contain a list a GCPs.." << std::endl;
std::cout << "The read generated extract from index (0, 0) must NOT contain a list a GCPs.." << std::endl;
return EXIT_FAILURE;
}
......@@ -118,14 +119,14 @@ int otbExtractROITestMetaData(int itkNotUsed(argc), char* argv[])
reader57->SetFileName(argv[3]);
reader57->GenerateOutputInformation();
if (reader57->GetOutput()->GetProjectionRef() != "" || boost::algorithm::istarts_with(reader57->GetOutput()->GetProjectionRef(), "LOCAL_CS"))
if (reader57->GetOutput()->GetProjectionRef().size())
{
std::cout << "The read generated extract from index (x, y) must NOT contain a ProjectionReference." << std::endl;
std::cout << "Found ProjectionReference: " << reader57->GetOutput()->GetProjectionRef() << std::endl;
return EXIT_FAILURE;
}
if (reader57->GetOutput()->GetGCPCount() != 0)
if (reader57->GetOutput()->GetGCPCount())
{
std::cout << "The read generated extract from index (x, y) must NOT contain a list a GCPs.." << std::endl;
return EXIT_FAILURE;
......
......@@ -743,14 +743,18 @@ void GDALImageIO::InternalReadImageInformation()
/* -------------------------------------------------------------------- */
/* Get the projection coordinate system of the image : ProjectionRef */
/* -------------------------------------------------------------------- */
if (dataset->GetProjectionRef() != nullptr && !std::string(dataset->GetProjectionRef()).empty())
const char* pszProjection = dataset->GetProjectionRef();
if (pszProjection != nullptr && !std::string(pszProjection).empty())
{
OGRSpatialReferenceH pSR = OSRNewSpatialReference(nullptr);
const char* pszProjection = nullptr;
pszProjection = dataset->GetProjectionRef();
if (OSRImportFromWkt(pSR, (char**)(&pszProjection)) == OGRERR_NONE)
if (strncmp(pszProjection, "LOCAL_CS",8) == 0)
{
// skip local coordinate system as they will cause crashed later
// In GDAL 3, they begin to do special processing for Transmercator local
// coordinate system
otbLogMacro(Debug, << "Skipping LOCAL_CS projection")
}
else if (OSRImportFromWkt(pSR, (char**)(&pszProjection)) == OGRERR_NONE)
{
char* pszPrettyWkt = nullptr;
OSRExportToPrettyWkt(pSR, &pszPrettyWkt, FALSE);
......@@ -761,7 +765,7 @@ void GDALImageIO::InternalReadImageInformation()
}
else
{
itk::EncapsulateMetaData<std::string>(dict, MetaDataKey::ProjectionRefKey, static_cast<std::string>(dataset->GetProjectionRef()));
itk::EncapsulateMetaData<std::string>(dict, MetaDataKey::ProjectionRefKey, static_cast<std::string>(pszProjection));
}
if (pSR != nullptr)
......@@ -1399,111 +1403,120 @@ void GDALImageIO::InternalWriteImageInformation(const void* buffer)
std::string projectionRef;
itk::ExposeMetaData<std::string>(dict, MetaDataKey::ProjectionRefKey, projectionRef);
ImageKeywordlist otb_kwl;
itk::ExposeMetaData<ImageKeywordlist>(dict, MetaDataKey::OSSIMKeywordlistKey, otb_kwl);
unsigned int gcpCount = 0;
itk::ExposeMetaData<unsigned int>(dict, MetaDataKey::GCPCountKey, gcpCount);
// In OTB we can have simultaneously projection ref, sensor keywordlist, GCPs
// but GDAL can't handle both geotransform and GCP (see issue #303). Here is the priority
// order we will be using in OTB:
// [ProjRef+geotransform] > [SensorKeywordlist+geotransform] > [GCPs]
/* -------------------------------------------------------------------- */
/* Set the GCPs */
/* Pre-compute geotransform */
/* -------------------------------------------------------------------- */
const double Epsilon = 1E-10;
if (projectionRef.empty() && (std::abs(m_Origin[0] - 0.5) > Epsilon || std::abs(m_Origin[1] - 0.5) > Epsilon ||
std::abs(m_Spacing[0] * m_Direction[0][0] - 1.0) > Epsilon || std::abs(m_Spacing[1] * m_Direction[1][1] - 1.0) > Epsilon))
{
// See issue #303 :
// If there is no ProjectionRef, and the GeoTransform is not the identity,
// then saving also GCPs is undefined behavior for GDAL, and a WGS84 projection crs
// is assigned arbitrarily
otbLogMacro(Warning, << "Skipping GCPs saving to prevent GDAL from assigning a WGS84 projref to file (" << m_FileName << ")")
}
else
{
unsigned int gcpCount = 0;
itk::ExposeMetaData<unsigned int>(dict, MetaDataKey::GCPCountKey, gcpCount);
if (gcpCount > 0)
{
GDAL_GCP* gdalGcps = new GDAL_GCP[gcpCount];
for (unsigned int gcpIndex = 0; gcpIndex < gcpCount; ++gcpIndex)
{
// Build the GCP string in the form of GCP_n
std::ostringstream lStream;
lStream << MetaDataKey::GCPParametersKey << gcpIndex;
std::string key = lStream.str();
OTB_GCP gcp;
itk::ExposeMetaData<OTB_GCP>(dict, key, gcp);
gdalGcps[gcpIndex].pszId = const_cast<char*>(gcp.m_Id.c_str());
gdalGcps[gcpIndex].pszInfo = const_cast<char*>(gcp.m_Info.c_str());
gdalGcps[gcpIndex].dfGCPPixel = gcp.m_GCPCol;
gdalGcps[gcpIndex].dfGCPLine = gcp.m_GCPRow;
gdalGcps[gcpIndex].dfGCPX = gcp.m_GCPX;
gdalGcps[gcpIndex].dfGCPY = gcp.m_GCPY;
gdalGcps[gcpIndex].dfGCPZ = gcp.m_GCPZ;
}
std::string gcpProjectionRef;
itk::ExposeMetaData<std::string>(dict, MetaDataKey::GCPProjectionKey, gcpProjectionRef);
dataset->SetGCPs(gcpCount, gdalGcps, gcpProjectionRef.c_str());
delete[] gdalGcps;
}
}
double geoTransform[6];
geoTransform[0] = m_Origin[0] - 0.5 * m_Spacing[0] * m_Direction[0][0];
geoTransform[3] = m_Origin[1] - 0.5 * m_Spacing[1] * m_Direction[1][1];
geoTransform[1] = m_Spacing[0] * m_Direction[0][0];
geoTransform[5] = m_Spacing[1] * m_Direction[1][1];
// FIXME: Here component 1 and 4 should be replaced by the orientation parameters
geoTransform[2] = 0.;
geoTransform[4] = 0.;
// only write geotransform if it has non-default values
bool writeGeotransform =
std::abs(geoTransform[0]) > Epsilon ||
std::abs(geoTransform[1] - 1.0) > Epsilon ||
std::abs(geoTransform[2]) > Epsilon ||
std::abs(geoTransform[3]) > Epsilon ||
std::abs(geoTransform[4]) > Epsilon ||
std::abs(geoTransform[5] - 1.0) > Epsilon;
/* -------------------------------------------------------------------- */
/* Set the projection coordinate system of the image : ProjectionRef */
/* Case 1: Set the projection coordinate system of the image */
/* -------------------------------------------------------------------- */
if (!projectionRef.empty())
{
{
dataset->SetProjection(projectionRef.c_str());
}
else
{
}
/* -------------------------------------------------------------------- */
/* Case 2: Sensor keywordlist */
/* -------------------------------------------------------------------- */
else if (otb_kwl.GetSize())
{
/* -------------------------------------------------------------------- */
/* Set the RPC coeffs if no projection available (since GDAL 1.10.0) */
/* Set the RPC coeffs (since GDAL 1.10.0) */
/* -------------------------------------------------------------------- */
ImageKeywordlist otb_kwl;
itk::ExposeMetaData<ImageKeywordlist>(dict, MetaDataKey::OSSIMKeywordlistKey, otb_kwl);
if (m_WriteRPCTags && otb_kwl.GetSize() != 0)
{
if (m_WriteRPCTags)
{
GDALRPCInfo gdalRpcStruct;
if (otb_kwl.convertToGDALRPC(gdalRpcStruct))
{
{
otbLogMacro(Debug, << "Saving RPC to file (" << m_FileName << ")")
char** rpcMetadata = RPCInfoToMD(&gdalRpcStruct);
dataset->SetMetadata(rpcMetadata, "RPC");
CSLDestroy(rpcMetadata);
}
}
}
}
/* -------------------------------------------------------------------- */
/* Case 3: Set the GCPs */
/* -------------------------------------------------------------------- */
else if (gcpCount)
{
otbLogMacro(Debug, << "Saving GCPs to file (" << m_FileName << ")")
std::vector<GDAL_GCP> gdalGcps(gcpCount);
std::vector<OTB_GCP> otbGcps(gcpCount);
for (unsigned int gcpIndex = 0; gcpIndex < gcpCount; ++gcpIndex)
{
// Build the GCP string in the form of GCP_n
std::ostringstream lStream;
lStream << MetaDataKey::GCPParametersKey << gcpIndex;
std::string key = lStream.str();
OTB_GCP &gcp = otbGcps[gcpIndex];
itk::ExposeMetaData<OTB_GCP>(dict, key, gcp);
gdalGcps[gcpIndex].pszId = const_cast<char*>(gcp.m_Id.c_str());
gdalGcps[gcpIndex].pszInfo = const_cast<char*>(gcp.m_Info.c_str());
gdalGcps[gcpIndex].dfGCPPixel = gcp.m_GCPCol;
gdalGcps[gcpIndex].dfGCPLine = gcp.m_GCPRow;
gdalGcps[gcpIndex].dfGCPX = gcp.m_GCPX;
gdalGcps[gcpIndex].dfGCPY = gcp.m_GCPY;
gdalGcps[gcpIndex].dfGCPZ = gcp.m_GCPZ;
if (writeGeotransform)
{
// we need to transform GCP col and row accordingly
// WARNING: assume no rotation is there
gdalGcps[gcpIndex].dfGCPPixel = (gcp.m_GCPCol - geoTransform[0]) / geoTransform[1];
gdalGcps[gcpIndex].dfGCPLine = (gcp.m_GCPRow - geoTransform[3]) / geoTransform[5];
}
}
std::string gcpProjectionRef;
itk::ExposeMetaData<std::string>(dict, MetaDataKey::GCPProjectionKey, gcpProjectionRef);
dataset->SetGCPs(gcpCount, gdalGcps.data(), gcpProjectionRef.c_str());
// disable geotransform with GCP
writeGeotransform = false;
}
/* -------------------------------------------------------------------- */
/* Set the six coefficients of affine geoTransform */
/* Save geotransform if needed. */
/* -------------------------------------------------------------------- */
if (std::abs(m_Origin[0] - 0.5) > Epsilon || std::abs(m_Origin[1] - 0.5) > Epsilon || std::abs(m_Spacing[0] * m_Direction[0][0] - 1.0) > Epsilon ||
std::abs(m_Spacing[1] * m_Direction[1][1] - 1.0) > Epsilon)
{
// Only set the geotransform if it is not identity (it may erase GCP)
itk::VariableLengthVector<double> geoTransform(6);
/// Reporting origin and spacing
// Beware : GDAL origin is at the corner of the top-left pixel
// whereas OTB/ITK origin is at the centre of the top-left pixel
geoTransform[0] = m_Origin[0] - 0.5 * m_Spacing[0] * m_Direction[0][0];
geoTransform[3] = m_Origin[1] - 0.5 * m_Spacing[1] * m_Direction[1][1];
geoTransform[1] = m_Spacing[0] * m_Direction[0][0];
geoTransform[5] = m_Spacing[1] * m_Direction[1][1];
// FIXME: Here component 1 and 4 should be replaced by the orientation parameters
geoTransform[2] = 0.;
geoTransform[4] = 0.;
dataset->SetGeoTransform(const_cast<double*>(geoTransform.GetDataPointer()));
// Error if writing to a ENVI file with a positive Y spacing
if (driverShortName == "ENVI" && geoTransform[5] > 0.)
if (writeGeotransform)
{
if ( driverShortName == "ENVI" && geoTransform[5] > 0.)
{
// Error if writing to a ENVI file with a positive Y spacing
itkExceptionMacro(<< "Can not write to ENVI file format with a positive Y spacing (" << m_FileName << ")");
}
dataset->SetGeoTransform(geoTransform);
}
}
/* -------------------------------------------------------------------- */
/* Report metadata. */
......
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