OTB  6.7.0
Orfeo Toolbox
otbGeographicalDistance.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbGeographicalDistance_hxx
22 #define otbGeographicalDistance_hxx
23 
25 #include "otbMath.h"
26 
27 namespace otb
28 {
29 
30 template <class TVector>
32 ::GeographicalDistance() : m_EarthRadius(6371000)
33 {}
34 
35 
36 template <class TVector>
37 double
39 ::Evaluate(const VectorType & x) const
40 {
41  // First check if vector length is sufficient
42  if(x.Size()<2)
43  itkExceptionMacro(<<"Vector length must be at least 2 to compute geographical distance.");
44 
45  // Call evaluate implementation with the first point being the
46  // origin
47  VectorType origin(x);
48  origin[0]=this->GetOrigin()[0];
49  origin[1]=this->GetOrigin()[1];
50 
51  return this->Evaluate(origin, x);
52 }
53 
54 template <class TVector>
55 double
57 ::Evaluate(const VectorType & x, const VectorType & y) const
58 {
59  // First check if vector length is sufficient
60  if(x.Size()<2 || y.Size()<2)
61  itkExceptionMacro(<<"Vector length must be at least 2 to compute geographical distance.");
62 
63  // Build some const variables
64  const double One = itk::NumericTraits<double>::One;
65  const double Two = One + One;
66  const double Deg2Rad = CONST_PI/180.;
67 
68  // Compute latitude and longitude differences
69  double dLat = (std::fabs(x[1] - y[1])) * Deg2Rad;
70  double dLon = (std::fabs(x[0] - y[0])) * Deg2Rad;
71 
72  // Compute dx in meters
73  double a = std::sin(dLat / Two) * std::sin(dLat / Two) + std::cos(y[1] * Deg2Rad) * std::cos(
74  x[1] * Deg2Rad) * std::sin(dLon / Two) * std::sin(dLon / Two);
75  double c = Two * std::atan2(std::sqrt(a), std::sqrt(One - a));
76  double d = m_EarthRadius * c;
77 
78  // Return result
79  return d;
80 }
81 
82 template <class TVector>
83 void
85 ::PrintSelf(std::ostream & os, itk::Indent indent) const
86 {
87  // Call superclass implementation
88  Superclass::PrintSelf(os, indent);
89 
90  // Earth radius
91  os<<indent<<"Earth radius: "<<m_EarthRadius<<std::endl;
92 }
93 } // End namespace otb
94 
95 #endif
constexpr double CONST_PI
Definition: otbMath.h:48
double Evaluate(const VectorType &x) const override
void PrintSelf(std::ostream &os, itk::Indent indent) const override