diff --git a/Code/UtilitiesAdapters/OGRAdapters/otbOGRDataSourceWrapper.cxx b/Code/UtilitiesAdapters/OGRAdapters/otbOGRDataSourceWrapper.cxx index c6636911d999bf6c4a218a53234d9a9b7995ea0b..b6056f02eff064efccf3c5cf50a67b432447be31 100644 --- a/Code/UtilitiesAdapters/OGRAdapters/otbOGRDataSourceWrapper.cxx +++ b/Code/UtilitiesAdapters/OGRAdapters/otbOGRDataSourceWrapper.cxx @@ -402,13 +402,14 @@ int otb::ogr::DataSource::Size(bool doForceComputation) const /*=================================[ Misc ]==================================*/ /*===========================================================================*/ -void otb::ogr::DataSource::GetGlobalExtent(double & ulx, +std::string otb::ogr::DataSource::GetGlobalExtent(double & ulx, double & uly, double & lrx, double & lry, bool force) const { OGREnvelope sExtent; + const_iterator lit = this->begin(); if(lit==this->end()) @@ -416,7 +417,11 @@ void otb::ogr::DataSource::GetGlobalExtent(double & ulx, itkGenericExceptionMacro(<< "Cannot compute global extent because there are no layers in the DataSource"); } - const OGRErr res = lit->ogr().GetExtent(&sExtent,force); + std::string outwkt = lit->GetProjectionRef(); + + const OGRSpatialReference * ref_srs = lit->GetSpatialRef(); + + const OGRErr res = lit->ogr().GetExtent(&sExtent,force); if(res!= OGRERR_NONE) @@ -429,8 +434,8 @@ void otb::ogr::DataSource::GetGlobalExtent(double & ulx, for(; lit!=this->end(); ++lit) { - OGREnvelope cExtent; - + OGREnvelope cExtent; + const OGRErr cres = lit->ogr().GetExtent(&cExtent,force); if(cres!= OGRERR_NONE) @@ -438,6 +443,35 @@ void otb::ogr::DataSource::GetGlobalExtent(double & ulx, itkGenericExceptionMacro(<< "Cannot retrieve extent of layer <" <<lit->GetName()<<">: " << CPLGetLastErrorMsg()); } + + const OGRSpatialReference * current_srs = lit->GetSpatialRef(); + + // If both srs are valid and if they are different + if(ref_srs && current_srs && current_srs->IsSame(ref_srs) == 0) + { + // Reproject cExtent in ref_srs + // OGRCreateCoordinateTransformation is not const-correct + OGRCoordinateTransformation * coordTransformation = OGRCreateCoordinateTransformation( + const_cast<OGRSpatialReference *>(current_srs), + const_cast<OGRSpatialReference *>(ref_srs)); + + coordTransformation->Transform(1,&cExtent.MinX,&cExtent.MinY); + coordTransformation->Transform(1,&cExtent.MaxX,&cExtent.MaxY); + + double real_minx = std::min(cExtent.MinX,cExtent.MaxX); + double real_miny = std::min(cExtent.MinY,cExtent.MaxY); + double real_maxx = std::max(cExtent.MinX,cExtent.MaxX); + double real_maxy = std::max(cExtent.MinY,cExtent.MaxY); + + cExtent.MinX = real_minx; + cExtent.MinY = real_miny; + cExtent.MaxX = real_maxx; + cExtent.MaxY = real_maxy; + + CPLFree(coordTransformation); + } + + // If srs are invalid, we assume that extent are coherent // Merge with previous layer sExtent.Merge(cExtent); @@ -447,6 +481,8 @@ void otb::ogr::DataSource::GetGlobalExtent(double & ulx, uly = sExtent.MinY; lrx = sExtent.MaxX; lry = sExtent.MaxY; + + return outwkt; } diff --git a/Code/UtilitiesAdapters/OGRAdapters/otbOGRDataSourceWrapper.h b/Code/UtilitiesAdapters/OGRAdapters/otbOGRDataSourceWrapper.h index 1e67896b45335c6031764d9ee9667339739bd67b..4f5babf28e549ae5056ae25f227e0c3ff656fc3f 100644 --- a/Code/UtilitiesAdapters/OGRAdapters/otbOGRDataSourceWrapper.h +++ b/Code/UtilitiesAdapters/OGRAdapters/otbOGRDataSourceWrapper.h @@ -212,7 +212,9 @@ public: */ int Size(bool doForceComputation) const; - /** Allow to retrieve the union of the extents of all layers + /** Allow to retrieve the union of the extents of all + * layers. In case of multiple layers with different SRS, the + * global extent is expressed in the SRS of the first layer. * \param[out] ulx reference to upper-left x coordinate of the * extent * \param[out] uly reference to upper-left y coordinate of the @@ -224,9 +226,11 @@ public: * \param[in] force Force computation of layers extents if not * available. May force the driver to walk all geometries to * compute the extent. + * \return The Wkt of the extent projection (which is the wkt of + * the first layer SRS) * \throw itk::ExceptionObject if the layers extents can not be retrieved. */ - void GetGlobalExtent(double & ulx, double & uly, double & lrx, double & lry, bool force = false) const; + std::string GetGlobalExtent(double & ulx, double & uly, double & lrx, double & lry, bool force = false) const; /** Grafts data and information from one data source to another.