Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  OGR
       4             :  * Purpose:  Geo-computations
       5             :  * Author:   Even Rouault, even dot rouault at spatialys.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2008-2010, Even Rouault <even dot rouault at spatialys.com>
       9             :  *
      10             :  * Permission is hereby granted, free of charge, to any person obtaining a
      11             :  * copy of this software and associated documentation files (the "Software"),
      12             :  * to deal in the Software without restriction, including without limitation
      13             :  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
      14             :  * and/or sell copies of the Software, and to permit persons to whom the
      15             :  * Software is furnished to do so, subject to the following conditions:
      16             :  *
      17             :  * The above copyright notice and this permission notice shall be included
      18             :  * in all copies or substantial portions of the Software.
      19             :  *
      20             :  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
      21             :  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      22             :  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
      23             :  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      24             :  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      25             :  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
      26             :  * DEALINGS IN THE SOFTWARE.
      27             :  ****************************************************************************/
      28             : 
      29             : #include "ogr_geo_utils.h"
      30             : #include "cpl_port.h"
      31             : #include "cpl_error.h"
      32             : 
      33             : constexpr double RAD2METER = (180.0 / M_PI) * 60.0 * 1852.0;
      34             : constexpr double METER2RAD = 1.0 / RAD2METER;
      35             : 
      36             : constexpr double DEG2RAD = M_PI / 180.0;
      37             : constexpr double RAD2DEG = 1.0 / DEG2RAD;
      38             : 
      39        1204 : static double OGR_Safe_acos(double x)
      40             : {
      41        1204 :     if (x > 1)
      42           0 :         x = 1;
      43        1204 :     else if (x < -1)
      44           0 :         x = -1;
      45        1204 :     return acos(x);
      46             : }
      47             : 
      48             : /************************************************************************/
      49             : /*                      OGR_GreatCircle_Distance()                      */
      50             : /************************************************************************/
      51             : 
      52          12 : double OGR_GreatCircle_Distance(double LatA_deg, double LonA_deg,
      53             :                                 double LatB_deg, double LonB_deg)
      54             : {
      55          12 :     const double cosP = cos((LonB_deg - LonA_deg) * DEG2RAD);
      56          12 :     const double LatA_rad = LatA_deg * DEG2RAD;
      57          12 :     const double LatB_rad = LatB_deg * DEG2RAD;
      58          12 :     const double cosa = cos(LatA_rad);
      59          12 :     const double sina = sin(LatA_rad);
      60          12 :     const double cosb = cos(LatB_rad);
      61          12 :     const double sinb = sin(LatB_rad);
      62          12 :     const double cos_angle = sina * sinb + cosa * cosb * cosP;
      63          12 :     return OGR_Safe_acos(cos_angle) * RAD2METER;
      64             : }
      65             : 
      66             : /************************************************************************/
      67             : /*                      OGR_GreatCircle_InitialHeading()                */
      68             : /************************************************************************/
      69             : 
      70           0 : double OGR_GreatCircle_InitialHeading(double LatA_deg, double LonA_deg,
      71             :                                       double LatB_deg, double LonB_deg)
      72             : {
      73           0 :     if (fabs(LatA_deg - 90) < 1e-10 || fabs(LatB_deg + 90) < 1e-10)
      74             :     {
      75           0 :         return 180;
      76             :     }
      77           0 :     else if (fabs(LatA_deg + 90) < 1e-10 || fabs(LatB_deg - 90) < 1e-10)
      78             :     {
      79           0 :         return 0;
      80             :     }
      81           0 :     else if (fabs(fmod(LonA_deg - LonB_deg, 360.0)) < 1e-10 &&
      82           0 :              fabs(LatA_deg - LatB_deg) < 1e-10)
      83             :     {
      84           0 :         return 0;  // Arbitrary number
      85             :     }
      86           0 :     else if (fabs(LatA_deg) < 1e-10 && fabs(LatB_deg) < 1e-10)
      87             :     {
      88           0 :         return (LonB_deg > LonA_deg) ? 90.0 : 270.0;
      89             :     }
      90           0 :     else if (fabs(fmod(LonA_deg - LonB_deg, 360.0)) < 1e-10)
      91             :     {
      92           0 :         return (LatA_deg > LatB_deg) ? 180.0 : 0.0;
      93             :     }
      94             :     else
      95             :     {
      96           0 :         const double LatA_rad = LatA_deg * DEG2RAD;
      97           0 :         const double LatB_rad = LatB_deg * DEG2RAD;
      98             : 
      99           0 :         const double cos_LatA = cos(LatA_rad);
     100           0 :         const double sin_LatA = sin(LatA_rad);
     101             : 
     102           0 :         const double diffG = (LonA_deg - LonB_deg) * DEG2RAD;
     103           0 :         const double cos_diffG = cos(diffG);
     104           0 :         const double sin_diffG = sin(diffG);
     105             : 
     106           0 :         const double denom = sin_LatA * cos_diffG - cos_LatA * tan(LatB_rad);
     107           0 :         if (denom == 0.0)
     108             :         {
     109             :             // Can be the case if Lat_A = -Lat_B and abs(LonA - LonB) = 180
     110           0 :             return 0.0;
     111             :         }
     112             : 
     113           0 :         double track = atan(sin_diffG / denom) * RAD2DEG;
     114             : 
     115           0 :         if (denom > 0.0)
     116             :         {
     117           0 :             track = 180 + track;
     118             :         }
     119           0 :         else if (track < 0)
     120             :         {
     121           0 :             track = 360 + track;
     122             :         }
     123             : 
     124           0 :         return track;
     125             :     }
     126             : }
     127             : 
     128             : /************************************************************************/
     129             : /*                     OGR_GreatCircle_ExtendPosition()                 */
     130             : /************************************************************************/
     131             : 
     132         609 : int OGR_GreatCircle_ExtendPosition(double dfLatA_deg, double dfLonA_deg,
     133             :                                    double dfDistance, double dfHeadingInA,
     134             :                                    double *pdfLatB_deg, double *pdfLonB_deg)
     135             : {
     136         609 :     const double dfHeadingRad = dfHeadingInA * DEG2RAD;
     137         609 :     const double cos_Heading = cos(dfHeadingRad);
     138         609 :     const double sin_Heading = sin(dfHeadingRad);
     139             : 
     140         609 :     const double dfDistanceRad = dfDistance * METER2RAD;
     141         609 :     const double cos_Distance = cos(dfDistanceRad);
     142         609 :     const double sin_Distance = sin(dfDistanceRad);
     143             : 
     144         609 :     const double dfLatA_rad = dfLatA_deg * DEG2RAD;
     145         609 :     const double cos_complement_LatA = sin(dfLatA_rad);
     146         609 :     const double sin_complement_LatA = cos(dfLatA_rad);
     147             : 
     148         609 :     if (dfDistance == 0.0)
     149             :     {
     150           0 :         *pdfLatB_deg = dfLatA_deg;
     151           0 :         *pdfLonB_deg = dfLonA_deg;
     152           0 :         return 1;
     153             :     }
     154             : 
     155         609 :     if (fabs(dfLatA_deg) >= 90.0)
     156             :     {
     157           0 :         *pdfLatB_deg = dfLatA_deg;
     158           0 :         *pdfLonB_deg = dfLonA_deg;
     159           0 :         return 0;
     160             :     }
     161             : 
     162         609 :     if (fabs(sin_Heading) < 1e-8)
     163             :     {
     164          13 :         *pdfLonB_deg = dfLonA_deg;
     165          13 :         if (fabs(fmod(dfHeadingInA + 360.0, 360.0)) < 1e-8)
     166             :         {
     167           8 :             *pdfLatB_deg = dfLatA_deg + dfDistanceRad * RAD2DEG;
     168             :         }
     169             :         else
     170             :         {
     171           5 :             *pdfLatB_deg = dfLatA_deg - dfDistanceRad * RAD2DEG;
     172             :         }
     173          13 :         return 1;
     174             :     }
     175             : 
     176         596 :     if (fabs(cos_complement_LatA) < 1e-8 && fabs(cos_Heading) < 1e-8)
     177             :     {
     178           0 :         *pdfLatB_deg = dfLatA_deg;
     179           0 :         if (fabs(dfHeadingInA - 90.0) < 1e-8)
     180             :         {
     181           0 :             *pdfLonB_deg = dfLonA_deg + dfDistanceRad * RAD2DEG;
     182             :         }
     183             :         else
     184             :         {
     185           0 :             *pdfLonB_deg = dfLonA_deg - dfDistanceRad * RAD2DEG;
     186             :         }
     187           0 :         return 1;
     188             :     }
     189             : 
     190         596 :     const double cos_complement_latB =
     191         596 :         cos_Distance * cos_complement_LatA +
     192         596 :         sin_Distance * sin_complement_LatA * cos_Heading;
     193             : 
     194         596 :     const double complement_latB = OGR_Safe_acos(cos_complement_latB);
     195             : 
     196         596 :     const double dfDenomin = sin(complement_latB) * sin_complement_LatA;
     197         596 :     if (dfDenomin == 0.0)
     198           0 :         CPLDebug("OGR", "OGR_GreatCircle_Distance: dfDenomin == 0.0");
     199         596 :     const double Cos_dG =
     200         596 :         (cos_Distance - cos_complement_latB * cos_complement_LatA) / dfDenomin;
     201         596 :     *pdfLatB_deg = 90 - complement_latB * RAD2DEG;
     202             : 
     203         596 :     const double dG_deg = OGR_Safe_acos(Cos_dG) * RAD2DEG;
     204             : 
     205         596 :     if (sin_Heading < 0)
     206         308 :         *pdfLonB_deg = dfLonA_deg - dG_deg;
     207             :     else
     208         288 :         *pdfLonB_deg = dfLonA_deg + dG_deg;
     209             : 
     210         596 :     if (*pdfLonB_deg > 180)
     211           0 :         *pdfLonB_deg -= 360;
     212         596 :     else if (*pdfLonB_deg <= -180)
     213           0 :         *pdfLonB_deg += 360;
     214             : 
     215         596 :     return 1;
     216             : }