diff --git a/Code/Projections/otbUtmMapProjection.h b/Code/Projections/otbUtmMapProjection.h index 03b7b985e8cff9c7536be44880f8d75e6b46990c..2a2e31258f6982bef643418be0005f3f698b76db 100644 --- a/Code/Projections/otbUtmMapProjection.h +++ b/Code/Projections/otbUtmMapProjection.h @@ -18,7 +18,6 @@ #ifndef __otbUtmMapProjection_h #define __otbUtmMapProjection_h -#include "projection/ossimMapProjection.h" #include "projection/ossimUtmProjection.h" #include "otbMapProjection.h" @@ -26,18 +25,20 @@ namespace otb { /** \class UtmMapProjection * \brief This class implements the UTM map projection. + * * It converts coordinates in longitude,latitude (WGS84) to UTM map coordinates. + * */ -template <InverseOrForwardTransformationEnum transform> -class ITK_EXPORT UtmMapProjection : public MapProjection<ossimUtmProjection, transform> +template <InverseOrForwardTransformationEnum TTransform> +class ITK_EXPORT UtmMapProjection : public MapProjection<ossimUtmProjection, TTransform> { public: /** Standard class typedefs. */ - typedef UtmMapProjection Self; - typedef MapProjection<ossimUtmProjection, transform> Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; + typedef UtmMapProjection Self; + typedef MapProjection<ossimUtmProjection, TTransform> Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; typedef typename Superclass::ScalarType ScalarType; typedef itk::Point<ScalarType, 2> InputPointType; @@ -50,20 +51,16 @@ public: itkTypeMacro(UtmMapProjection, MapProjection); virtual void SetZone(long zone); - virtual void SetZone(const InputPointType& ground); virtual void SetHemisphere(char hemisphere); - virtual long GetZone(); + virtual int GetZone() const; virtual const char GetHemisphere() const; virtual void SetZoneAndHemisphereFromGeoPoint(const InputPointType& geoPoint); -// virtual void SetZoneAndHemisphereFromCartoPoint(const OutputPointType &cartoPoint); - virtual int GetZoneFromGeoPoint(const InputPointType& geoPoint); - -// virtual void Initialize(const InputPointType& middlePoint); + virtual int GetZoneFromGeoPoint(const InputPointType& geoPoint) const; protected: - UtmMapProjection(); - virtual ~UtmMapProjection(); + UtmMapProjection() {}; + virtual ~UtmMapProjection() {}; private: UtmMapProjection(const Self &); //purposely not implemented diff --git a/Code/Projections/otbUtmMapProjection.txx b/Code/Projections/otbUtmMapProjection.txx index 6e9d3f90b42991ab5a8dc7ca323e66811462afdb..50d72ff3952f4f580832a469a0bef60a06b4fcc8 100644 --- a/Code/Projections/otbUtmMapProjection.txx +++ b/Code/Projections/otbUtmMapProjection.txx @@ -20,123 +20,72 @@ #define __otbUtmMapProjection_txx #include "otbUtmMapProjection.h" -#include "otbMapProjections.h" namespace otb { -template <InverseOrForwardTransformationEnum transform> -UtmMapProjection<transform> -::UtmMapProjection() -{ -} - -template <InverseOrForwardTransformationEnum transform> -UtmMapProjection<transform> -::~UtmMapProjection() -{ -} ///Set the zone -template <InverseOrForwardTransformationEnum transform> -void UtmMapProjection<transform> +template <InverseOrForwardTransformationEnum TTransform> +void UtmMapProjection<TTransform> ::SetZone(long zone) { this->m_MapProjection->setZone(zone); this->Modified(); } -///Set the zone -template <InverseOrForwardTransformationEnum transform> -void UtmMapProjection<transform> -::SetZone(const InputPointType& ground) -{ - ossimGpt ossimGround; - ossimGround.lon = ground[0]; - ossimGround.lat = ground[1]; - this->m_MapProjection->setZone(ossimGround); - this->Modified(); -} - ///Set the hemisphere -template <InverseOrForwardTransformationEnum transform> -void UtmMapProjection<transform> +template <InverseOrForwardTransformationEnum TTransform> +void UtmMapProjection<TTransform> ::SetHemisphere(char hemisphere) { this->m_MapProjection->setHemisphere(hemisphere); this->Modified(); } -template <InverseOrForwardTransformationEnum transform> -void UtmMapProjection<transform> +template <InverseOrForwardTransformationEnum TTransform> +void UtmMapProjection<TTransform> ::SetZoneAndHemisphereFromGeoPoint(const InputPointType& geoPoint) { - double latitude = geoPoint[1]; - char hemisphere; - int zone; + char hemisphere; - if (latitude > 0.) hemisphere = 'N'; + if (geoPoint[1] > 0.) hemisphere = 'N'; else hemisphere = 'S'; this->SetHemisphere(hemisphere); - zone = this->GetZoneFromGeoPoint(geoPoint); + int zone = this->GetZoneFromGeoPoint(geoPoint); this->SetZone(zone); } -/*template <InverseOrForwardTransformationEnum transform> -void UtmMapProjection<transform> -::SetZoneAndHemisphereFromCartoPoint(const OutputPointType &cartoPoint) -{ -InputPointType geoPoint; - - // TODO : Tester que la projection est bien inverse !!! -geoPoint = this->TransformPoint(cartoPoint); -this->SetZoneAndHemisphereFromGeoPoint(geoPoint); -} */ - ///\return the zone -template <InverseOrForwardTransformationEnum transform> -long UtmMapProjection<transform> -::GetZone() +template <InverseOrForwardTransformationEnum TTransform> +int UtmMapProjection<TTransform> +::GetZone() const { - long zone; - zone = this->m_MapProjection->getZone(); + int zone = this->m_MapProjection->getZone(); return zone; } ///\return the hemisphere -template <InverseOrForwardTransformationEnum transform> -const char UtmMapProjection<transform> +template <InverseOrForwardTransformationEnum TTransform> +const char UtmMapProjection<TTransform> ::GetHemisphere() const { - char hemisphere = 0; - hemisphere = this->m_MapProjection->getHemisphere(); + char hemisphere = this->m_MapProjection->getHemisphere(); return hemisphere; } -template <InverseOrForwardTransformationEnum transform> -int UtmMapProjection<transform> -::GetZoneFromGeoPoint(const InputPointType& geoPoint) +template <InverseOrForwardTransformationEnum TTransform> +int UtmMapProjection<TTransform> +::GetZoneFromGeoPoint(const InputPointType& geoPoint) const { - double longitude = geoPoint[0]; - //double latitude = geoPoint[1]; - int zone; - - // Each UTM zone is a narrow zone of 6 degrees in width - // Zone 31 is between 0 and 6 degrees (lon) - // There is 60 zones in each hemisphere - zone = ((static_cast<int>(floor(longitude / 6)) + 30) % 60 + 60) % 60 + 1; + //use ossim to handle the special case of UTM + ossimGpt point(geoPoint[1], geoPoint[0]); + int zone = this->m_MapProjection->computeZone(point); return zone; } -/* template <InverseOrForwardTransformationEnum transform> - void UtmMapProjection<transform> - ::Initialize(const InputPointType& middlePoint) - { - this->SetZoneAndHemisphereFromCartoPoint(middlePoint); -}*/ - } #endif diff --git a/Testing/Code/Projections/otbUtmMapProjection.cxx b/Testing/Code/Projections/otbUtmMapProjection.cxx index 589e121711327366445fe2d98e960e8453cc0d9e..a7ab9e9aa3f7fb5f5636fc1f4a80ef452db8a582 100644 --- a/Testing/Code/Projections/otbUtmMapProjection.cxx +++ b/Testing/Code/Projections/otbUtmMapProjection.cxx @@ -45,18 +45,37 @@ int otbUtmMapProjection(int argc, char* argv[]) otb::UtmForwardProjection::Pointer lUtmProjection2 = otb::UtmForwardProjection::New(); file << lUtmProjection2->GetWkt() << std::endl << std::endl; - lUtmProjection->SetHemisphere('S'); + lUtmProjection2->SetHemisphere('S'); file << "We should be in the south hemisphere now\n"; file << lUtmProjection2->GetWkt() << std::endl << std::endl; - lUtmProjection->SetZone(48); + lUtmProjection2->SetZone(48); file << "We should be in zone 48\n"; file << lUtmProjection2->GetWkt() << std::endl << std::endl; - lUtmProjection->SetHemisphere('N'); + lUtmProjection2->SetHemisphere('N'); file << "And back in the north hemisphere\n"; file << lUtmProjection2->GetWkt() << std::endl << std::endl; + lUtmProjection2->SetHemisphere('S');//Back in the south to make sure the following works + otb::UtmForwardProjection::InputPointType groundPoint; + groundPoint[0] = -4.5; + groundPoint[1] = 48.5; + lUtmProjection2->SetZoneAndHemisphereFromGeoPoint(groundPoint); + file << "We should be in zone 30 in north hemisphere\n"; + file << lUtmProjection2->GetWkt() << std::endl << std::endl; + + //And checking accesors + file << "Zone should be 30: " << lUtmProjection2->GetZone() << "\n"; + file << "Hemisphere should be N: " << lUtmProjection2->GetHemisphere() << "\n"; + + groundPoint[0] = 1.5; + groundPoint[1] = 43.5; + file << "Zone should be 31: " << lUtmProjection2->GetZoneFromGeoPoint(groundPoint) << "\n"; + + groundPoint[0] = 5; + groundPoint[1] = 60; + file << "Zone should be 32 (exception near Norway): " << lUtmProjection2->GetZoneFromGeoPoint(groundPoint) << "\n"; file.close(); return EXIT_SUCCESS;