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             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "ogr_geo_utils.h"
      14             : #include "cpl_port.h"
      15             : #include "cpl_error.h"
      16             : 
      17             : constexpr double DEG2RAD = M_PI / 180.0;
      18             : constexpr double RAD2DEG = 1.0 / DEG2RAD;
      19             : 
      20        1222 : static double OGR_Safe_acos(double x)
      21             : {
      22        1222 :     if (x > 1)
      23           0 :         x = 1;
      24        1222 :     else if (x < -1)
      25           0 :         x = -1;
      26        1222 :     return acos(x);
      27             : }
      28             : 
      29             : /************************************************************************/
      30             : /*                      OGR_GreatCircle_Distance()                      */
      31             : /************************************************************************/
      32             : 
      33          12 : double OGR_GreatCircle_Distance(double LatA_deg, double LonA_deg,
      34             :                                 double LatB_deg, double LonB_deg,
      35             :                                 double dfRadius)
      36             : {
      37          12 :     const double cosP = cos((LonB_deg - LonA_deg) * DEG2RAD);
      38          12 :     const double LatA_rad = LatA_deg * DEG2RAD;
      39          12 :     const double LatB_rad = LatB_deg * DEG2RAD;
      40          12 :     const double cosa = cos(LatA_rad);
      41          12 :     const double sina = sin(LatA_rad);
      42          12 :     const double cosb = cos(LatB_rad);
      43          12 :     const double sinb = sin(LatB_rad);
      44          12 :     const double cos_angle = sina * sinb + cosa * cosb * cosP;
      45          12 :     return OGR_Safe_acos(cos_angle) * dfRadius;
      46             : }
      47             : 
      48             : /************************************************************************/
      49             : /*                      OGR_GreatCircle_InitialHeading()                */
      50             : /************************************************************************/
      51             : 
      52           0 : double OGR_GreatCircle_InitialHeading(double LatA_deg, double LonA_deg,
      53             :                                       double LatB_deg, double LonB_deg)
      54             : {
      55           0 :     if (fabs(LatA_deg - 90) < 1e-10 || fabs(LatB_deg + 90) < 1e-10)
      56             :     {
      57           0 :         return 180;
      58             :     }
      59           0 :     else if (fabs(LatA_deg + 90) < 1e-10 || fabs(LatB_deg - 90) < 1e-10)
      60             :     {
      61           0 :         return 0;
      62             :     }
      63           0 :     else if (fabs(fmod(LonA_deg - LonB_deg, 360.0)) < 1e-10 &&
      64           0 :              fabs(LatA_deg - LatB_deg) < 1e-10)
      65             :     {
      66           0 :         return 0;  // Arbitrary number
      67             :     }
      68           0 :     else if (fabs(LatA_deg) < 1e-10 && fabs(LatB_deg) < 1e-10)
      69             :     {
      70           0 :         return (LonB_deg > LonA_deg) ? 90.0 : 270.0;
      71             :     }
      72           0 :     else if (fabs(fmod(LonA_deg - LonB_deg, 360.0)) < 1e-10)
      73             :     {
      74           0 :         return (LatA_deg > LatB_deg) ? 180.0 : 0.0;
      75             :     }
      76             :     else
      77             :     {
      78           0 :         const double LatA_rad = LatA_deg * DEG2RAD;
      79           0 :         const double LatB_rad = LatB_deg * DEG2RAD;
      80             : 
      81           0 :         const double cos_LatA = cos(LatA_rad);
      82           0 :         const double sin_LatA = sin(LatA_rad);
      83             : 
      84           0 :         const double diffG = (LonA_deg - LonB_deg) * DEG2RAD;
      85           0 :         const double cos_diffG = cos(diffG);
      86           0 :         const double sin_diffG = sin(diffG);
      87             : 
      88           0 :         const double denom = sin_LatA * cos_diffG - cos_LatA * tan(LatB_rad);
      89           0 :         if (denom == 0.0)
      90             :         {
      91             :             // Can be the case if Lat_A = -Lat_B and abs(LonA - LonB) = 180
      92           0 :             return 0.0;
      93             :         }
      94             : 
      95           0 :         double track = atan(sin_diffG / denom) * RAD2DEG;
      96             : 
      97           0 :         if (denom > 0.0)
      98             :         {
      99           0 :             track = 180 + track;
     100             :         }
     101           0 :         else if (track < 0)
     102             :         {
     103           0 :             track = 360 + track;
     104             :         }
     105             : 
     106           0 :         return track;
     107             :     }
     108             : }
     109             : 
     110             : /************************************************************************/
     111             : /*                     OGR_GreatCircle_ExtendPosition()                 */
     112             : /************************************************************************/
     113             : 
     114         609 : int OGR_GreatCircle_ExtendPosition(double dfLatA_deg, double dfLonA_deg,
     115             :                                    double dfDistance, double dfHeadingInA,
     116             :                                    double dfRadius, double *pdfLatB_deg,
     117             :                                    double *pdfLonB_deg)
     118             : {
     119         609 :     const double dfHeadingRad = dfHeadingInA * DEG2RAD;
     120         609 :     const double cos_Heading = cos(dfHeadingRad);
     121         609 :     const double sin_Heading = sin(dfHeadingRad);
     122             : 
     123         609 :     const double dfDistanceRad = dfDistance / dfRadius;
     124         609 :     const double cos_Distance = cos(dfDistanceRad);
     125         609 :     const double sin_Distance = sin(dfDistanceRad);
     126             : 
     127         609 :     const double dfLatA_rad = dfLatA_deg * DEG2RAD;
     128         609 :     const double cos_complement_LatA = sin(dfLatA_rad);
     129         609 :     const double sin_complement_LatA = cos(dfLatA_rad);
     130             : 
     131         609 :     if (dfDistance == 0.0)
     132             :     {
     133           0 :         *pdfLatB_deg = dfLatA_deg;
     134           0 :         *pdfLonB_deg = dfLonA_deg;
     135           0 :         return 1;
     136             :     }
     137             : 
     138         609 :     if (fabs(dfLatA_deg) >= 90.0)
     139             :     {
     140           0 :         *pdfLatB_deg = dfLatA_deg;
     141           0 :         *pdfLonB_deg = dfLonA_deg;
     142           0 :         return 0;
     143             :     }
     144             : 
     145         609 :     if (fabs(sin_Heading) < 1e-8)
     146             :     {
     147           4 :         *pdfLonB_deg = dfLonA_deg;
     148           4 :         if (fabs(fmod(dfHeadingInA + 360.0, 360.0)) < 1e-8)
     149             :         {
     150           2 :             *pdfLatB_deg = dfLatA_deg + dfDistanceRad * RAD2DEG;
     151             :         }
     152             :         else
     153             :         {
     154           2 :             *pdfLatB_deg = dfLatA_deg - dfDistanceRad * RAD2DEG;
     155             :         }
     156           4 :         return 1;
     157             :     }
     158             : 
     159         605 :     if (fabs(cos_complement_LatA) < 1e-8 && fabs(cos_Heading) < 1e-8)
     160             :     {
     161           0 :         *pdfLatB_deg = dfLatA_deg;
     162           0 :         if (fabs(dfHeadingInA - 90.0) < 1e-8)
     163             :         {
     164           0 :             *pdfLonB_deg = dfLonA_deg + dfDistanceRad * RAD2DEG;
     165             :         }
     166             :         else
     167             :         {
     168           0 :             *pdfLonB_deg = dfLonA_deg - dfDistanceRad * RAD2DEG;
     169             :         }
     170           0 :         return 1;
     171             :     }
     172             : 
     173         605 :     const double cos_complement_latB =
     174         605 :         cos_Distance * cos_complement_LatA +
     175         605 :         sin_Distance * sin_complement_LatA * cos_Heading;
     176             : 
     177         605 :     const double complement_latB = OGR_Safe_acos(cos_complement_latB);
     178             : 
     179         605 :     const double dfDenomin = sin(complement_latB) * sin_complement_LatA;
     180         605 :     if (dfDenomin == 0.0)
     181           0 :         CPLDebug("OGR", "OGR_GreatCircle_Distance: dfDenomin == 0.0");
     182         605 :     const double Cos_dG =
     183         605 :         (cos_Distance - cos_complement_latB * cos_complement_LatA) / dfDenomin;
     184         605 :     *pdfLatB_deg = 90 - complement_latB * RAD2DEG;
     185             : 
     186         605 :     const double dG_deg = OGR_Safe_acos(Cos_dG) * RAD2DEG;
     187             : 
     188         605 :     if (sin_Heading < 0)
     189         176 :         *pdfLonB_deg = dfLonA_deg - dG_deg;
     190             :     else
     191         429 :         *pdfLonB_deg = dfLonA_deg + dG_deg;
     192             : 
     193         605 :     if (*pdfLonB_deg > 180)
     194           3 :         *pdfLonB_deg -= 360;
     195         602 :     else if (*pdfLonB_deg <= -180)
     196           0 :         *pdfLonB_deg += 360;
     197             : 
     198         605 :     return 1;
     199             : }