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