BUG: Utm projection GetZoneFromGeoPoint didn't handle special cases, const...

BUG: Utm projection GetZoneFromGeoPoint didn't handle special cases,  const correctness, coverage, fix error in test, remove useless methods

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

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

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;