Commit 7b09c7e4 authored by Cédric Traizet's avatar Cédric Traizet

Merge branch 'compatibility_gdal3' into 'develop'

Compatibility with gdal 3

See merge request !511
parents f1f29b3a a9ad9c59
Pipeline #1747 passed with stages
in 7 minutes and 51 seconds
PROJCS["UTM Zone 31, Northern Hemisphere",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]
Forward projection:
[1.44, 43.605] -> [374100.828760369, 4829184.80636283]
PROJCS["UTM Zone 31, Northern Hemisphere",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]
Inverse projection:
[374100.8, 4829184.8] -> [1.43999964524438, 43.6049999378673]
PROJCS["WGS 84 / UTM zone 31N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32631"]]
Forward projection:
[1.44, 43.605] -> [374100.828760369, 4829184.80636283]
PROJCS["WGS 84 / UTM zone 31N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32631"]]
Inverse projection:
[374100.8, 4829184.8] -> [1.43999964524438, 43.6049999378673]
PROJCS["UTM Zone 31, Northern Hemisphere",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]
Forward projection:
[1.44, 43.605] -> [374100.8287608, 4829184.80641377]
PROJCS["UTM Zone 31, Northern Hemisphere",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]
Inverse projection:
[374100.8, 4829184.8] -> [1.439999645251, 43.6049999374087]
ProjRef: GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AXIS["Latitude",NORTH],
AXIS["Longitude",EAST],
AUTHORITY["EPSG","4326"]]
Origin: [1943.5, 2364.5]
Spacing: [1, 1]
......@@ -21,6 +21,7 @@
#define otbSpatialReference_h
#include "OTBGdalAdaptersExport.h"
#include "ogr_spatialref.h"
#include <memory>
#include <string>
......@@ -150,6 +151,14 @@ public:
*/
unsigned int ToEPSG() const;
#if GDAL_VERSION_NUM >= 3000000
/** Set the Axis mapping strategy
* proxy to the eponym ogr spatial reference method
* \param strategy Axis mapping stategy
*/
void SetAxisMappingStrategy(OSRAxisMappingStrategy strategy);
#endif
/**
* Find which UTM zone a given (lat,lon) point falls in.
* \pre -180<=lon<=180
......
......@@ -291,4 +291,12 @@ void SpatialReference::UTMFromGeoPoint(double lon, double lat, unsigned int & zo
// post conditions
assert(zone<=60);
}
#if GDAL_VERSION_NUM >= 3000000
void SpatialReference::SetAxisMappingStrategy(OSRAxisMappingStrategy strategy)
{
m_SR.get()->SetAxisMappingStrategy(strategy);
};
#endif
}
......@@ -70,6 +70,11 @@ GenericMapProjection<TDirectionOfMapping, TScalarType, NInputDimensions, NOutput
SpatialReference wgs84 = SpatialReference::FromWGS84();
SpatialReference wktSpatialReference = SpatialReference::FromDescription(projectionRefWkt);
#if GDAL_VERSION_NUM >= 3000000
wgs84.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
wktSpatialReference.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
#endif
if(DirectionOfMapping == TransformDirection::INVERSE)
{
std::unique_ptr<CoordinateTransformation> newMapProjection(new CoordinateTransformation(wktSpatialReference,wgs84));
......
......@@ -23,7 +23,7 @@
#include <iomanip>
#include "otbGenericMapProjection.h"
#include "ogr_spatialref.h"
int otbGenericMapProjection(int itkNotUsed(argc), char* argv[])
{
const char * outFileName = argv[1];
......@@ -35,11 +35,19 @@ int otbGenericMapProjection(int itkNotUsed(argc), char* argv[])
/** Test the ability to instantiate a projection from a string*/
std::string projectionRefWkt =
"PROJCS[\"UTM Zone 31, Northern Hemisphere\", GEOGCS[\"WGS 84\", DATUM[\"WGS_1984\", SPHEROID[\"WGS 84\", 6378137, 298.257223563, AUTHORITY[\"EPSG\",\"7030\"]], TOWGS84[0, 0, 0, 0, 0, 0, 0], AUTHORITY[\"EPSG\",\"6326\"]], PRIMEM[\"Greenwich\", 0, AUTHORITY[\"EPSG\",\"8901\"]], UNIT[\"degree\", 0.0174532925199433, AUTHORITY[\"EPSG\",\"9108\"]], AXIS[\"Lat\", NORTH], AXIS[\"Long\", EAST], AUTHORITY[\"EPSG\",\"4326\"]], PROJECTION[\"Transverse_Mercator\"], PARAMETER[\"latitude_of_origin\", 0], PARAMETER[\"central_meridian\", 3], PARAMETER[\"scale_factor\", 0.9996], PARAMETER[\"false_easting\", 500000], PARAMETER[\"false_northing\", 0], UNIT[\"Meter\", 1]]";
auto baselineSpatialReference = otb::SpatialReference::FromDescription(projectionRefWkt);
typedef otb::GenericMapProjection<otb::TransformDirection::FORWARD> GenericMapProjection;
GenericMapProjection::Pointer genericMapProjection = GenericMapProjection::New();
genericMapProjection->SetWkt(projectionRefWkt);
file << genericMapProjection->GetWkt() << std::endl << std::endl;
auto testSpatialReferenceForward = otb::SpatialReference::FromDescription(genericMapProjection->GetWkt());
if (testSpatialReferenceForward!=baselineSpatialReference)
{
std::cerr << "The spatial reference used in the forward test is different from the input spatial reference" << std::endl;
return EXIT_FAILURE;
}
itk::Point<double, 2> point;
point[0] = 1.44;
......@@ -52,7 +60,14 @@ int otbGenericMapProjection(int itkNotUsed(argc), char* argv[])
typedef otb::GenericMapProjection<otb::TransformDirection::INVERSE> GenericMapProjectionInverse;
GenericMapProjectionInverse::Pointer genericMapProjectionInverse = GenericMapProjectionInverse::New();
genericMapProjectionInverse->SetWkt(projectionRefWkt);
file << genericMapProjectionInverse->GetWkt() << std::endl << std::endl;
auto testSpatialReferenceReverse = otb::SpatialReference::FromDescription(genericMapProjectionInverse->GetWkt());
if (testSpatialReferenceForward!=baselineSpatialReference)
{
std::cerr << "The spatial reference used in the reverse test is different from the input spatial reference" << std::endl;
return EXIT_FAILURE;
}
point[0] = 374100.8;
point[1] = 4829184.8;
......
......@@ -289,6 +289,9 @@ void OGRVectorDataIO::Write(const itk::DataObject* datag, char ** /** unused */)
if (projectionInformationAvailable)
{
oSRS = static_cast<OGRSpatialReference *>(OSRNewSpatialReference(projectionRefWkt.c_str()));
#if GDAL_VERSION_NUM >= 3000000
oSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
#endif
}
// Retrieving root node
......
......@@ -2613,26 +2613,37 @@ void TestHelper::ogrReportOnLayer(OGRLayer * ref_poLayer,
CheckValueTolerance("Extent: MaxX", ref_oExt.MaxX, test_oExt.MaxX, nbdiff, m_ReportErrors, epsilon);
CheckValueTolerance("Extent: MaxY", ref_oExt.MaxY, test_oExt.MaxY, nbdiff, m_ReportErrors, epsilon);
}
/// Spatial reference comparison :
/// We use the IsSame GDAL method to determine if the spatial references are equivalent,
/// if they're not, a wkt string comparison is done for error reporting.
/// Note that two spatial references might be equivalent without having the same
/// wkt representation.
if ((ref_poLayer->GetSpatialRef() == nullptr) ||
(test_poLayer->GetSpatialRef() == nullptr) ||
( ! ref_poLayer->GetSpatialRef()->IsSame( test_poLayer->GetSpatialRef() ) ))
{
char *ref_pszWKT;
char *test_pszWKT;
char *ref_pszWKT;
char *test_pszWKT;
if (ref_poLayer->GetSpatialRef() == nullptr) ref_pszWKT = CPLStrdup("(unknown)");
else
{
ref_poLayer->GetSpatialRef()->exportToPrettyWkt(&ref_pszWKT);
}
if (test_poLayer->GetSpatialRef() == nullptr) test_pszWKT = CPLStrdup("(unknown)");
else
{
test_poLayer->GetSpatialRef()->exportToPrettyWkt(&test_pszWKT);
}
otbCheckStringValue("Layer SRS WKT", ref_pszWKT, test_pszWKT, nbdiff, m_ReportErrors);
CPLFree(ref_pszWKT);
CPLFree(test_pszWKT);
if (ref_poLayer->GetSpatialRef() == nullptr) ref_pszWKT = CPLStrdup("(unknown)");
else
{
ref_poLayer->GetSpatialRef()->exportToPrettyWkt(&ref_pszWKT);
}
if (test_poLayer->GetSpatialRef() == nullptr) test_pszWKT = CPLStrdup("(unknown)");
else
{
test_poLayer->GetSpatialRef()->exportToPrettyWkt(&test_pszWKT);
}
otbCheckStringValue("Layer SRS WKT", ref_pszWKT, test_pszWKT, nbdiff, m_ReportErrors);
CPLFree(ref_pszWKT);
CPLFree(test_pszWKT);
}
otbCheckStringValue("FID Column", ref_poLayer->GetFIDColumn(), test_poLayer->GetFIDColumn(), nbdiff, m_ReportErrors);
otbCheckStringValue("Geometry Column", ref_poLayer->GetGeometryColumn(),
test_poLayer->GetGeometryColumn(), nbdiff, m_ReportErrors);
......
......@@ -128,9 +128,16 @@ PersistentImageSampleExtractorFilter<TInputImage>
if(projectionInformationAvailable)
{
OGRSpatialReference imgSRS;
#if GDAL_VERSION_NUM >= 3000000 // importFromWkt is const-correct in GDAL 3
const char *projWktCstr = projectionRefWkt.c_str();
OGRErr err = imgSRS.importFromWkt( &projWktCstr );
#else
const char *projWktCstr = projectionRefWkt.c_str();
char **projWktPointer = const_cast<char**>(&projWktCstr);
OGRErr err = imgSRS.importFromWkt( projWktPointer );
#endif
if (err == OGRERR_NONE)
{
// get input layer
......
......@@ -335,7 +335,11 @@ void Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::GenerateOutpu
if (!m_ProjectionRef.empty())
{
OGRSpatialReference oSRS;
#if GDAL_VERSION_NUM >= 3000000 // importFromWkt is const-correct in GDAL 3
const char *wkt = m_ProjectionRef.c_str();
#else
char *wkt = const_cast<char *> (m_ProjectionRef.c_str());
#endif
oSRS.importFromWkt(&wkt);
m_IsGeographic = oSRS.IsGeographic(); // TODO check if this test is valid for all projection systems
}
......
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