1

2# -*- coding: utf-8 -*-

3

4 u'''Various utility functions.

5

6 After I{(C) Chris Veness 2011-2024} published under the same MIT Licence**, see

7 U{Latitude/Longitude<https://www.Movable-Type.co.UK/scripts/latlong.html>} and

8 U{Vector-based geodesy<https://www.Movable-Type.co.UK/scripts/latlong-vectors.html>}.

9'''

10# make sure int/int division yields float quotient, see .basics

11 from __future__ import division as _; del _ # PYCHOK semicolon

12

13 from pygeodesy.basics import _copysign, isinstanceof, isint, isstr, neg

14 from pygeodesy.constants import EPS, EPS0, INF, NAN, PI, PI2, PI_2, R_M, \

15 _M_KM, _M_NM, _M_SM, _0_0, _1__90, _0_5, _1_0, \

16 _N_1_0, _2__PI, _10_0, _90_0, _180_0, _N_180_0, \

17 _360_0, _400_0, isnan, isnear0, _copysign_0_0, \

18 _float, _isfinite, _over, _umod_360, _umod_PI2

19 from pygeodesy.errors import _ValueError, _xkwds, _xkwds_get1, _ALL_LAZY, _MODS

20 from pygeodesy.internals import _passarg, _passargs # , _MODS?

21 from pygeodesy.interns import _edge_, _radians_, _semi_circular_, _SPACE_

22# from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS # from .errors

23 from pygeodesy.units import Degrees, Degrees_, Feet, Float, Lam, Lamd, \

24 Meter, Meter2, Radians, Radians_

25

26 from math import acos, asin, atan2 as _atan2, cos, degrees, fabs, radians, \

27 sin, tan as _tan # pow

28

29__all__ = _ALL_LAZY.utily

30__version__ = '24.11.26'

31

32_G_DEG = _float(_400_0 / _360_0) # grades per degree

33_G_RAD = _float(_400_0 / PI2) # grades per radian

34_M_CHAIN = _float( 20.1168) # meter per yard2m(1) * 22

35_M_FATHOM = _float( 1.8288) # meter per yard2m(1) * 2 or _M_NM * 1e-3

36_M_FOOT = _float( 0.3048) # meter per Int'l foot, 1 / 3.2808398950131 = 10_000 / (254 * 12)

37_M_FOOT_GE = _float( 0.31608) # meter per German Fuss, 1 / 3.1637560111364

38_M_FOOT_FR = _float( 0.3248406) # meter per French Pied-du-Roi or pied, 1 / 3.0784329298739

39_M_FOOT_US = _float( 0.3048006096012192) # meter per US Survey foot, 1200 / 3937

40_M_FURLONG = _float( 201.168) # meter per furlong, 220 * yard2m(1) = 10 * m2chain(1)

41# _M_KM = _float(1000.0) # meter per kilo meter

42# _M_NM = _float(1852.0) # meter per nautical mile

43# _M_SM = _float(1609.344) # meter per statute mile

44_M_TOISE = _float( 1.9490436) # meter per French toise, 6 pieds = 6 / 3.0784329298739

45_M_YARD_UK = _float( 0.9144) # meter per yard, 254 * 12 * 3 / 10_000 = 3 * _M_FOOT

46

47

48 def _abs1nan(x):

49 '''(INTERNAL) Bracket C{x}.

50 '''

51 return _N_1_0 < x < _1_0 or isnan(x)

52

53

54 def acos1(x):

55 '''Return C{math.acos(max(-1, min(1, B{x})))}.

56 '''

57 return acos(x) if _abs1nan(x) else (PI if x < 0 else _0_0)

58

59

60 def acre2ha(acres):

61 '''Convert acres to hectare.

62

63 @arg acres: Value in acres (C{scalar}).

64

65 @return: Value in C{hectare} (C{float}).

66

67 @raise ValueError: Invalid B{C{acres}}.

68 '''

69 # 0.40468564224 == acre2m2(1) / 10_000

70 return Float(ha=Float(acres) * 0.40468564224)

71

72

73 def acre2m2(acres):

74 '''Convert acres to I{square} meter.

75

76 @arg acres: Value in acres (C{scalar}).

77

78 @return: Value in C{meter^2} (C{float}).

79

80 @raise ValueError: Invalid B{C{acres}}.

81 '''

82 # 4046.8564224 == chain2m(1) * furlong2m(1)

83 return Meter2(Float(acres) * 4046.8564224)

84

85

86 def asin1(x):

87 '''Return C{math.asin(max(-1, min(1, B{x})))}.

88 '''

89 return asin(x) if _abs1nan(x) else _copysign(PI_2, x)

90

91

92 def atan1(y, x=_1_0):

93 '''Return C{atan(B{y} / B{x})} angle in C{radians} M{[-PI/2..+PI/2]}

94 using C{atan2} for consistency and to avoid C{ZeroDivisionError}.

95 '''

96 return _atan2(-y, -x) if x < 0 else _atan2(y, x or _0_0) # -0. to 0.

97

98

99 def atan1d(y, x=_1_0):

100 '''Return C{atan(B{y} / B{x})} angle in C{degrees} M{[-90..+90]}

101 using C{atan2d} for consistency and to avoid C{ZeroDivisionError}.

102

103 @see: Function L{pygeodesy.atan2d}.

104 '''

105 return atan2d(-y, -x) if x < 0 else atan2d(y, x or _0_0) # -0. to 0.

106

107

108 def atan2(y, x):

109 '''Return C{atan2(B{y}, B{x})} in radians M{[-PI..+PI]}.

110

111 @see: I{Karney}'s C++ function U{Math.atan2d

112 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Math.html>}.

113 '''

114 return _atan2u(y, x, _passarg, PI, PI_2)

115

116

117 def atan2b(y, x):

118 '''Return C{atan2(B{y}, B{x})} in degrees M{[0..+360]}.

119

120 @see: Function L{pygeodesy.atan2d}.

121 '''

122 b = atan2d(y, x)

123 if b < 0:

124 b += _360_0

125 return b

126

127

128 def atan2d(y, x, reverse=False):

129 '''Return C{atan2(B{y}, B{x})} in degrees M{[-180..+180]},

130 optionally I{reversed} (by 180 degrees for C{azimuth}s).

131

132 @see: I{Karney}'s C++ function U{Math.atan2d

133 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Math.html>}.

134 '''

135 d = _atan2u(y, x, degrees, _180_0, _90_0)

136 return _azireversed(d) if reverse else d

137

138

139 def _atan2u(y, x, _2u, H, Q): # Half, Quarter turn in units

140 '''(INTERNAL) Helper for functions C{atan2} and C{atan2d}.

141 '''

142 if fabs(y) > fabs(x) > 0:

143 if y < 0: # q = 3

144 r = _2u(_atan2(x, -y)) - Q

145 else: # q = 2

146 r = Q - _2u(_atan2(x, y))

147 elif isnan(x) or isnan(y):

148 return NAN

149 elif y:

150 if x > 0: # q = 0

151 r = _2u(_atan2(y, x))

152 elif x < 0: # q = 1

153 r = _copysign(H, y) - _2u(_atan2(y, -x))

154 else: # x == 0

155 r = _copysign(Q, y)

156 else: # preserve signBit(y) like Python's math.atan2

157 r = _copysign(H, y) if x < 0 else _0_0

158 return r

159

160

161 def _azireversed(azimuth): # in .rhumbBase

162 '''(INTERNAL) Return the I{reverse} B{C{azimuth}} in degrees M{[-180..+180]}.

163 '''

164 return azimuth + (_N_180_0 if azimuth > 0 else _180_0)

165

166

167 def chain2m(chains):

168 '''Convert I{UK} chains to meter.

169

170 @arg chains: Value in chains (C{scalar}).

171

172 @return: Value in C{meter} (C{float}).

173

174 @raise ValueError: Invalid B{C{chains}}.

175 '''

176 return Meter(Float(chains=chains) * _M_CHAIN)

177

178

179 def circle4(earth, lat):

180 '''Get the equatorial or a parallel I{circle of latitude}.

181

182 @arg earth: The earth radius (C{meter}), ellipsoid or datum

183 (L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple}).

184 @arg lat: Geodetic latitude (C{degrees90}, C{str}).

185

186 @return: A L{Circle4Tuple}C{(radius, height, lat, beta)}.

187

188 @raise RangeError: Latitude B{C{lat}} outside valid range and

189 L{rangerrors<pygeodesy.rangerrors>} is C{True}.

190

191 @raise TypeError: Invalid B{C{earth}}.

192

193 @raise ValueError: B{C{earth}} or B{C{lat}}.

194 '''

195 E = _MODS.datums._earth_ellipsoid(earth)

196 return E.circle4(lat)

197

198

199 def cot(rad, **raiser_kwds):

200 '''Return the C{cotangent} of an angle in C{radians}.

201

202 @arg rad: Angle (C{radians}).

203 @kwarg raiser_kwds: Use C{B{raiser}=False} to avoid

204 ValueErrors or optionally, additional

205 ValueError keyword argments.

206

207 @return: C{cot(B{rad})}.

208

209 @raise Error: If L{pygeodesy.isnear0}C{(sin(B{rad})}.

210 '''

211 try:

212 return _cotu(*sincos2(rad), **raiser_kwds)

213 except ZeroDivisionError:

214 raise _valueError(cot, rad, **raiser_kwds)

215

216

217 def cot_(*rads, **raiser_kwds):

218 '''Return the C{cotangent} of angle(s) in C{radians}.

219

220 @arg rads: One or more angles (each in C{radians}).

221

222 @return: Yield the C{cot(B{rad})} for each angle.

223

224 @see: Function L{pygeodesy.cot} for further details.

225 '''

226 try:

227 for r in rads:

228 yield _cotu(*sincos2(r), **raiser_kwds)

229 except ZeroDivisionError:

230 raise _valueError(cot_, r, **raiser_kwds)

231

232

233 def cotd(deg, **raiser_kwds):

234 '''Return the C{cotangent} of an angle in C{degrees}.

235

236 @arg deg: Angle (C{degrees}).

237 @kwarg raiser_kwds: Use C{B{raiser}=False} to avoid

238 ValueErrors or optionally, additional

239 ValueError keyword argments.

240

241 @return: C{cot(B{deg})}.

242

243 @raise Error: If L{pygeodesy.isnear0}C{(sin(B{deg})}.

244 '''

245 try:

246 return _cotu(*sincos2d(deg), **raiser_kwds)

247 except ZeroDivisionError:

248 raise _valueError(cotd, deg, **raiser_kwds)

249

250

251 def cotd_(*degs, **raiser_kwds):

252 '''Return the C{cotangent} of angle(s) in C{degrees}.

253

254 @arg degs: One or more angles (each in C{degrees}).

255

256 @return: Yield the C{cot(B{deg})} for each angle.

257

258 @see: Function L{pygeodesy.cotd} for further details.

259 '''

260 try:

261 for d in degs:

262 yield _cotu(*sincos2d(d), **raiser_kwds)

263 except ZeroDivisionError:

264 raise _valueError(cotd_, d, **raiser_kwds)

265

266

267 def _cotu(s, c, **raiser_kwds):

268 '''(INTERNAL) Helper for functions C{cot}, C{cotd}, C{cot_} and C{cotd_}.

269 '''

270 return _tanu(c, s, **raiser_kwds)

271

272

273 def degrees90(rad):

274 '''Convert radians to degrees and wrap M{[-90..+90)}.

275

276 @arg rad: Angle (C{radians}).

277

278 @return: Angle, wrapped (C{degrees90}).

279 '''

280 return wrap90(degrees(rad))

281

282

283 def degrees180(rad):

284 '''Convert radians to degrees and wrap M{[-180..+180)}.

285

286 @arg rad: Angle (C{radians}).

287

288 @return: Angle, wrapped (C{degrees180}).

289 '''

290 return wrap180(degrees(rad))

291

292

293 def degrees360(rad):

294 '''Convert radians to degrees and wrap M{[0..+360)}.

295

296 @arg rad: Angle (C{radians}).

297

298 @return: Angle, wrapped (C{degrees360}).

299 '''

300 return _umod_360(degrees(rad))

301

302

303 def degrees2grades(deg):

304 '''Convert degrees to I{grades} (aka I{gons} or I{gradians}).

305

306 @arg deg: Angle (C{degrees}).

307

308 @return: Angle (C{grades}).

309 '''

310 return Float(grades=Degrees(deg) * _G_DEG)

311

312

313 def degrees2m(deg, radius=R_M, lat=0):

314 '''Convert an angle to a distance along the equator or along a parallel

315 at (geodetic) latitude.

316

317 @arg deg: The angle (C{degrees}).

318 @kwarg radius: Mean earth radius (C{meter}), an ellipsoid or datum

319 (L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple}).

320 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

321

322 @return: Distance (C{meter}, same units as B{C{radius}} or polar and

323 equatorial radii) or C{0.0} for near-polar B{C{lat}}.

324

325 @raise RangeError: Latitude B{C{lat}} outside valid range and

326 L{rangerrors<pygeodesy.rangerrors>} is C{True}.

327

328 @raise TypeError: Invalid B{C{radius}}.

329

330 @raise ValueError: Invalid B{C{deg}}, B{C{radius}} or B{C{lat}}.

331

332 @see: Function L{radians2m} and L{m2degrees}.

333 '''

334 return _Radians2m(Lamd(deg=deg, clip=0), radius, lat)

335

336

337 def fathom2m(fathoms):

338 '''Convert I{Imperial} fathom to meter.

339

340 @arg fathoms: Value in fathoms (C{scalar}).

341

342 @return: Value in C{meter} (C{float}).

343

344 @raise ValueError: Invalid B{C{fathoms}}.

345

346 @see: Function L{toise2m}, U{Fathom<https://WikiPedia.org/wiki/Fathom>}

347 and U{Klafter<https://WikiPedia.org/wiki/Klafter>}.

348 '''

349 return Meter(Float(fathoms=fathoms) * _M_FATHOM)

350

351

352 def ft2m(feet, usurvey=False, pied=False, fuss=False):

353 '''Convert I{International}, I{US Survey}, I{French} or I{German}

354 B{C{feet}} to C{meter}.

355

356 @arg feet: Value in feet (C{scalar}).

357 @kwarg usurvey: If C{True}, convert I{US Survey} foot else ...

358 @kwarg pied: If C{True}, convert French I{pied-du-Roi} else ...

359 @kwarg fuss: If C{True}, convert German I{Fuss}, otherwise

360 I{International} foot to C{meter}.

361

362 @return: Value in C{meter} (C{float}).

363

364 @raise ValueError: Invalid B{C{feet}}.

365 '''

366 return Meter(Feet(feet) * (_M_FOOT_US if usurvey else

367 (_M_FOOT_FR if pied else

368 (_M_FOOT_GE if fuss else _M_FOOT))))

369

370

371 def furlong2m(furlongs):

372 '''Convert a furlong to meter.

373

374 @arg furlongs: Value in furlongs (C{scalar}).

375

376 @return: Value in C{meter} (C{float}).

377

378 @raise ValueError: Invalid B{C{furlongs}}.

379 '''

380 return Meter(Float(furlongs=furlongs) * _M_FURLONG)

381

382

383 def grades(rad):

384 '''Convert radians to I{grades} (aka I{gons} or I{gradians}).

385

386 @arg rad: Angle (C{radians}).

387

388 @return: Angle (C{grades}).

389 '''

390 return Float(grades=Radians(rad) * _G_RAD)

391

392

393 def grades400(rad):

394 '''Convert radians to I{grades} (aka I{gons} or I{gradians}) and wrap M{[0..+400)}.

395

396 @arg rad: Angle (C{radians}).

397

398 @return: Angle, wrapped (C{grades}).

399 '''

400 return Float(grades400=wrapPI2(rad) * _G_RAD)

401

402

403 def grades2degrees(gon):

404 '''Convert I{grades} (aka I{gons} or I{gradians}) to C{degrees}.

405

406 @arg gon: Angle (C{grades}).

407

408 @return: Angle (C{degrees}).

409 '''

410 return Degrees(Float(gon=gon) / _G_DEG)

411

412

413 def grades2radians(gon):

414 '''Convert I{grades} (aka I{gons} or I{gradians}) to C{radians}.

415

416 @arg gon: Angle (C{grades}).

417

418 @return: Angle (C{radians}).

419 '''

420 return Radians(Float(gon=gon) / _G_RAD)

421

422

423 def km2m(km):

424 '''Convert kilo meter to meter (m).

425

426 @arg km: Value in kilo meter (C{scalar}).

427

428 @return: Value in meter (C{float}).

429

430 @raise ValueError: Invalid B{C{km}}.

431 '''

432 return Meter(Float(km=km) * _M_KM)

433

434

435 def _loneg(lon):

436 '''(INTERNAL) "Complement" of C{lon}.

437 '''

438 return _180_0 - lon

439

440

441 def m2chain(meter):

442 '''Convert meter to I{UK} chains.

443

444 @arg meter: Value in meter (C{scalar}).

445

446 @return: Value in C{chains} (C{float}).

447

448 @raise ValueError: Invalid B{C{meter}}.

449 '''

450 return Float(chain=Meter(meter) / _M_CHAIN) # * 0.049709695378986715

451

452

453 def m2degrees(distance, radius=R_M, lat=0):

454 '''Convert a distance to an angle along the equator or along a parallel

455 at (geodetic) latitude.

456

457 @arg distance: Distance (C{meter}, same units as B{C{radius}}).

458 @kwarg radius: Mean earth radius (C{meter}), an ellipsoid or datum

459 (L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple}).

460 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

461

462 @return: Angle (C{degrees}) or C{INF} for near-polar B{C{lat}}.

463

464 @raise RangeError: Latitude B{C{lat}} outside valid range and

465 L{rangerrors<pygeodesy.rangerrors>} is C{True}.

466

467 @raise TypeError: Invalid B{C{radius}}.

468

469 @raise ValueError: Invalid B{C{distance}}, B{C{radius}} or B{C{lat}}.

470

471 @see: Function L{m2radians} and L{degrees2m}.

472 '''

473 return degrees(m2radians(distance, radius=radius, lat=lat))

474

475

476 def m2fathom(meter):

477 '''Convert meter to I{Imperial} fathoms.

478

479 @arg meter: Value in meter (C{scalar}).

480

481 @return: Value in C{fathoms} (C{float}).

482

483 @raise ValueError: Invalid B{C{meter}}.

484

485 @see: Function L{m2toise}, U{Fathom<https://WikiPedia.org/wiki/Fathom>}

486 and U{Klafter<https://WikiPedia.org/wiki/Klafter>}.

487 '''

488 return Float(fathom=Meter(meter) / _M_FATHOM) # * 0.546806649

489

490

491 def m2ft(meter, usurvey=False, pied=False, fuss=False):

492 '''Convert meter to I{International}, I{US Survey}, I{French} or

493 or I{German} feet (C{ft}).

494

495 @arg meter: Value in meter (C{scalar}).

496 @kwarg usurvey: If C{True}, convert to I{US Survey} foot else ...

497 @kwarg pied: If C{True}, convert to French I{pied-du-Roi} else ...

498 @kwarg fuss: If C{True}, convert to German I{Fuss}, otherwise to

499 I{International} foot.

500

501 @return: Value in C{feet} (C{float}).

502

503 @raise ValueError: Invalid B{C{meter}}.

504 '''

505 # * 3.2808333333333333, US Survey 3937 / 1200

506 # * 3.2808398950131235, Int'l 10_000 / (254 * 12)

507 return Float(feet=Meter(meter) / (_M_FOOT_US if usurvey else

508 (_M_FOOT_FR if pied else

509 (_M_FOOT_GE if fuss else _M_FOOT))))

510

511

512 def m2furlong(meter):

513 '''Convert meter to furlongs.

514

515 @arg meter: Value in meter (C{scalar}).

516

517 @return: Value in C{furlongs} (C{float}).

518

519 @raise ValueError: Invalid B{C{meter}}.

520 '''

521 return Float(furlong=Meter(meter) / _M_FURLONG) # * 0.00497096954

522

523

524 def m2km(meter):

525 '''Convert meter to kilo meter (Km).

526

527 @arg meter: Value in meter (C{scalar}).

528

529 @return: Value in Km (C{float}).

530

531 @raise ValueError: Invalid B{C{meter}}.

532 '''

533 return Float(km=Meter(meter) / _M_KM)

534

535

536 def m2NM(meter):

537 '''Convert meter to nautical miles (NM).

538

539 @arg meter: Value in meter (C{scalar}).

540

541 @return: Value in C{NM} (C{float}).

542

543 @raise ValueError: Invalid B{C{meter}}.

544 '''

545 return Float(NM=Meter(meter) / _M_NM) # * 5.39956804e-4

546

547

548 def m2radians(distance, radius=R_M, lat=0):

549 '''Convert a distance to an angle along the equator or along a parallel

550 at (geodetic) latitude.

551

552 @arg distance: Distance (C{meter}, same units as B{C{radius}}).

553 @kwarg radius: Mean earth radius (C{meter}, an ellipsoid or datum

554 (L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple}).

555 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

556

557 @return: Angle (C{radians}) or C{INF} for near-polar B{C{lat}}.

558

559 @raise RangeError: Latitude B{C{lat}} outside valid range and

560 L{rangerrors<pygeodesy.rangerrors>} is C{True}.

561

562 @raise TypeError: Invalid B{C{radius}}.

563

564 @raise ValueError: Invalid B{C{distance}}, B{C{radius}} or B{C{lat}}.

565

566 @see: Function L{m2degrees} and L{radians2m}.

567 '''

568 m = circle4(radius, lat).radius

569 return INF if m < EPS0 else Radians(Float(distance=distance) / m)

570

571

572 def m2SM(meter):

573 '''Convert meter to statute miles (SM).

574

575 @arg meter: Value in meter (C{scalar}).

576

577 @return: Value in C{SM} (C{float}).

578

579 @raise ValueError: Invalid B{C{meter}}.

580 '''

581 return Float(SM=Meter(meter) / _M_SM) # * 6.21369949e-4 == 1 / 1609.344

582

583

584 def m2toise(meter):

585 '''Convert meter to French U{toises<https://WikiPedia.org/wiki/Toise>}.

586

587 @arg meter: Value in meter (C{scalar}).

588

589 @return: Value in C{toises} (C{float}).

590

591 @raise ValueError: Invalid B{C{meter}}.

592

593 @see: Function L{m2fathom}.

594 '''

595 return Float(toise=Meter(meter) / _M_TOISE) # * 0.513083632632119

596

597

598 def m2yard(meter):

599 '''Convert meter to I{UK} yards.

600

601 @arg meter: Value in meter (C{scalar}).

602

603 @return: Value in C{yards} (C{float}).

604

605 @raise ValueError: Invalid B{C{meter}}.

606 '''

607 return Float(yard=Meter(meter) / _M_YARD_UK) # * 1.0936132983377078

608

609

610 def NM2m(nm):

611 '''Convert nautical miles to meter (m).

612

613 @arg nm: Value in nautical miles (C{scalar}).

614

615 @return: Value in meter (C{float}).

616

617 @raise ValueError: Invalid B{C{nm}}.

618 '''

619 return Meter(Float(nm=nm) * _M_NM)

620

621

622 def radians2m(rad, radius=R_M, lat=0):

623 '''Convert an angle to a distance along the equator or along a parallel

624 at (geodetic) latitude.

625

626 @arg rad: The angle (C{radians}).

627 @kwarg radius: Mean earth radius (C{meter}) or an ellipsoid or datum

628 (L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple}).

629 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

630

631 @return: Distance (C{meter}, same units as B{C{radius}} or polar and

632 equatorial radii) or C{0.0} for near-polar B{C{lat}}.

633

634 @raise RangeError: Latitude B{C{lat}} outside valid range and

635 L{rangerrors<pygeodesy.rangerrors>} is C{True}.

636

637 @raise TypeError: Invalid B{C{radius}}.

638

639 @raise ValueError: Invalid B{C{rad}}, B{C{radius}} or B{C{lat}}.

640

641 @see: Function L{degrees2m} and L{m2radians}.

642 '''

643 return _Radians2m(Lam(rad=rad, clip=0), radius, lat)

644

645

646 def _Radians2m(rad, radius, lat):

647 '''(INTERNAL) Helper for C{degrees2m} and C{radians2m}.

648 '''

649 m = circle4(radius, lat).radius

650 return _0_0 if m < EPS0 else (rad * m)

651

652

653 def radiansPI(deg):

654 '''Convert and wrap degrees to radians M{[-PI..+PI]}.

655

656 @arg deg: Angle (C{degrees}).

657

658 @return: Radians, wrapped (C{radiansPI})

659 '''

660 return wrapPI(radians(deg))

661

662

663 def radiansPI2(deg):

664 '''Convert and wrap degrees to radians M{[0..+2PI)}.

665

666 @arg deg: Angle (C{degrees}).

667

668 @return: Radians, wrapped (C{radiansPI2})

669 '''

670 return _umod_PI2(radians(deg))

671

672

673 def radiansPI_2(deg):

674 '''Convert and wrap degrees to radians M{[-3PI/2..+PI/2]}.

675

676 @arg deg: Angle (C{degrees}).

677

678 @return: Radians, wrapped (C{radiansPI_2})

679 '''

680 return wrapPI_2(radians(deg))

681

682

683 def _sin0cos2(q, r, sign):

684 '''(INTERNAL) 2-tuple (C{sin(r), cos(r)}) in quadrant C{0 <= B{q} <= 3}

685 and C{sin} zero I{signed} with B{C{sign}}.

686 '''

687 if r < PI_2:

688 s, c = sin(r), cos(r)

689 t = s, c, -s, -c, s

690 else: # r == PI_2

691 t = _1_0, _0_0, _N_1_0, _0_0, _1_0

692# else: # r == 0, testUtility failures

693# t = _0_0, _1_0, _0_0, _N_1_0, _0_0

694# q &= 3

695 s = t[q] or _copysign_0_0(sign)

696 c = t[q + 1] or _0_0

697 return s, c

698

699

700 def SinCos2(x):

701 '''Get C{sin} and C{cos} of I{typed} angle.

702

703 @arg x: Angle (L{Degrees}, L{Radians} or scalar C{radians}).

704

705 @return: 2-Tuple (C{sin(B{x})}, C{cos(B{x})}).

706 '''

707 return sincos2d(x) if isinstanceof(x, Degrees, Degrees_) else (

708 sincos2(x) if isinstanceof(x, Radians, Radians_) else

709 sincos2(float(x))) # assume C{radians}

710

711

712 def sincos2(rad):

713 '''Return the C{sine} and C{cosine} of an angle in C{radians}.

714

715 @arg rad: Angle (C{radians}).

716

717 @return: 2-Tuple (C{sin(B{rad})}, C{cos(B{rad})}).

718

719 @see: U{GeographicLib<https://GeographicLib.SourceForge.io/C++/doc/

720 classGeographicLib_1_1Math.html#sincosd>} function U{sincosd

721 <https://SourceForge.net/p/geographiclib/code/ci/release/tree/

722 python/geographiclib/geomath.py#l155>} and C++ U{sincosd

723 <https://SourceForge.net/p/geographiclib/code/ci/release/tree/

724 include/GeographicLib/Math.hpp#l558>}.

725 '''

726 if _isfinite(rad):

727 q = int(rad * _2__PI) # int(math.floor)

728 if q < 0:

729 q -= 1

730 t = _sin0cos2(q & 3, rad - q * PI_2, rad)

731 else:

732 t = NAN, NAN

733 return t

734

735

736 def sincos2_(*rads):

737 '''Return the C{sine} and C{cosine} of angle(s) in C{radians}.

738

739 @arg rads: One or more angles (C{radians}).

740

741 @return: Yield the C{sin(B{rad})} and C{cos(B{rad})} for each angle.

742

743 @see: function L{sincos2}.

744 '''

745 for r in rads:

746 s, c = sincos2(r)

747 yield s

748 yield c

749

750

751 def sincos2d(deg, **adeg):

752 '''Return the C{sine} and C{cosine} of an angle in C{degrees}.

753

754 @arg deg: Angle (C{degrees}).

755 @kwarg adeg: Optional correction (C{degrees}).

756

757 @return: 2-Tuple (C{sin(B{deg_})}, C{cos(B{deg_})}, C{B{deg_} =

758 B{deg} + B{adeg}}).

759

760 @see: U{GeographicLib<https://GeographicLib.SourceForge.io/C++/doc/

761 classGeographicLib_1_1Math.html#sincosd>} function U{sincosd

762 <https://SourceForge.net/p/geographiclib/code/ci/release/tree/

763 python/geographiclib/geomath.py#l155>} and C++ U{sincosd

764 <https://SourceForge.net/p/geographiclib/code/ci/release/tree/

765 include/GeographicLib/Math.hpp#l558>}.

766 '''

767 if _isfinite(deg):

768 q = int(deg * _1__90) # int(math.floor)

769 if q < 0:

770 q -= 1

771 d = deg - q * _90_0

772 if adeg:

773 t = _xkwds_get1(adeg, adeg=_0_0)

774 d = _MODS.karney._around(d + t)

775 t = _sin0cos2(q & 3, radians(d), deg)

776 else:

777 t = NAN, NAN

778 return t

779

780

781 def sincos2d_(*degs):

782 '''Return the C{sine} and C{cosine} of angle(s) in C{degrees}.

783

784 @arg degs: One or more angles (C{degrees}).

785

786 @return: Yield the C{sin(B{deg})} and C{cos(B{deg})} for each angle.

787

788 @see: Function L{sincos2d}.

789 '''

790 for d in degs:

791 s, c = sincos2d(d)

792 yield s

793 yield c

794

795

796 def sincostan3(rad):

797 '''Return the C{sine}, C{cosine} and C{tangent} of an angle in C{radians}.

798

799 @arg rad: Angle (C{radians}).

800

801 @return: 3-Tuple (C{sin(B{rad})}, C{cos(B{rad})}, C{tan(B{rad})}).

802

803 @see: Function L{sincos2}.

804 '''

805 s, c = sincos2(float(rad))

806 t = NAN if s is NAN else (_over(s, c) if s else neg(s, neg0=c < 0))

807 return s, c, t

808

809

810 def SM2m(sm):

811 '''Convert statute miles to meter (m).

812

813 @arg sm: Value in statute miles (C{scalar}).

814

815 @return: Value in meter (C{float}).

816

817 @raise ValueError: Invalid B{C{sm}}.

818 '''

819 return Meter(Float(sm=sm) * _M_SM)

820

821

822 def tan_2(rad, **semi): # edge=1

823 '''Compute the tangent of half angle.

824

825 @arg rad: Angle (C{radians}).

826 @kwarg semi: Angle or edge name and index

827 for semi-circular error.

828

829 @return: M{tan(rad / 2)} (C{float}).

830

831 @raise ValueError: If B{C{rad}} is semi-circular

832 and B{C{semi}} is given.

833 '''

834 # .formy.excessKarney_, .sphericalTrigonometry.areaOf

835 if semi and isnear0(fabs(rad) - PI):

836 for n, v in semi.items():

837 break

838 n = _SPACE_(n, _radians_) if not isint(v) else \

839 _SPACE_(_MODS.streprs.Fmt.SQUARE(**semi), _edge_)

840 raise _ValueError(n, rad, txt=_semi_circular_)

841

842 return _tan(rad * _0_5) if _isfinite(rad) else NAN

843

844

845 def tan(rad, **raiser_kwds):

846 '''Return the C{tangent} of an angle in C{radians}.

847

848 @arg rad: Angle (C{radians}).

849 @kwarg raiser_kwds: Use C{B{raiser}=False} to avoid

850 ValueErrors or optionally, additional

851 ValueError keyword argments.

852

853 @return: C{tan(B{rad})}.

854

855 @raise Error: If L{pygeodesy.isnear0}C{(cos(B{rad})}.

856 '''

857 try:

858 return _tanu(*sincos2(rad), **raiser_kwds)

859 except ZeroDivisionError:

860 raise _valueError(tan, rad, **raiser_kwds)

861

862

863 def tan_(*rads, **raiser_kwds):

864 '''Return the C{tangent} of angle(s) in C{radians}.

865

866 @arg rads: One or more angles (each in C{radians}).

867

868 @return: Yield the C{tan(B{rad})} for each angle.

869

870 @see: Function L{pygeodesy.tan} for futher details.

871 '''

872 try:

873 for r in rads:

874 yield _tanu(*sincos2(r), **raiser_kwds)

875 except ZeroDivisionError:

876 raise _valueError(tan_, r, **raiser_kwds)

877

878

879 def tand(deg, **raiser_kwds):

880 '''Return the C{tangent} of an angle in C{degrees}.

881

882 @arg deg: Angle (C{degrees}).

883 @kwarg raiser_kwds: Use C{B{raiser}=False} to avoid

884 ValueErrors or optionally, additional

885 ValueError keyword argments.

886

887 @return: C{tan(B{deg})}.

888

889 @raise Error: If L{pygeodesy.isnear0}C{(cos(B{deg})}.

890 '''

891 try:

892 return _tanu(*sincos2d(deg), **raiser_kwds)

893 except ZeroDivisionError:

894 raise _valueError(tand, deg, **raiser_kwds)

895

896

897 def tand_(*degs, **raiser_kwds):

898 '''Return the C{tangent} of angle(s) in C{degrees}.

899

900 @arg degs: One or more angles (each in C{degrees}).

901

902 @return: Yield the C{tan(B{deg})} for each angle.

903

904 @see: Function L{pygeodesy.tand} for futher details.

905 '''

906 try:

907 for d in degs:

908 yield _tanu(*sincos2d(d), **raiser_kwds)

909 except ZeroDivisionError:

910 raise _valueError(tand_, d, **raiser_kwds)

911

912

913 def tanPI_2_2(rad):

914 '''Compute the tangent of half angle, 90 degrees rotated.

915

916 @arg rad: Angle (C{radians}).

917

918 @return: M{tan((rad + PI/2) / 2)} (C{float}).

919 '''

920 return _tan((rad + PI_2) * _0_5) if _isfinite(rad) else (

921 NAN if isnan(rad) else _copysign(_90_0, rad))

922

923

924 def _tanu(s, c, raiser=True, **unused):

925 '''(INTERNAL) Helper for functions C{_cotu}, C{tan}, C{tan_}, C{tand} and C{tand_}.

926 '''

927 if s:

928 if raiser and isnear0(c):

929 raise ZeroDivisionError()

930 s = _over(s, c)

931 elif c < 0:

932 s = -s # negate-0

933 return s

934

935

936 def toise2m(toises):

937 '''Convert French U{toises<https://WikiPedia.org/wiki/Toise>} to meter.

938

939 @arg toises: Value in toises (C{scalar}).

940

941 @return: Value in C{meter} (C{float}).

942

943 @raise ValueError: Invalid B{C{toises}}.

944

945 @see: Function L{fathom2m}.

946 '''

947 return Meter(Float(toises=toises) * _M_TOISE)

948

949

950 def truncate(x, ndigits=None):

951 '''Truncate to the given number of digits.

952

953 @arg x: Value to truncate (C{scalar}).

954 @kwarg ndigits: Number of digits (C{int}),

955 aka I{precision}.

956

957 @return: Truncated B{C{x}} (C{float}).

958

959 @see: Python function C{round}.

960 '''

961 if isint(ndigits):

962 p = _10_0**ndigits

963 x = int(x * p) / p

964 return x

965

966

967 def unroll180(lon1, lon2, wrap=True):

968 '''Unroll longitudinal delta and wrap longitude in degrees.

969

970 @arg lon1: Start longitude (C{degrees}).

971 @arg lon2: End longitude (C{degrees}).

972 @kwarg wrap: If C{True}, wrap and unroll to the M{(-180..+180]}

973 range (C{bool}).

974

975 @return: 2-Tuple C{(B{lon2}-B{lon1}, B{lon2})} unrolled (C{degrees},

976 C{degrees}).

977

978 @see: Capability C{LONG_UNROLL} in U{GeographicLib

979 <https://GeographicLib.SourceForge.io/html/python/interface.html#outmask>}.

980 '''

981 d = lon2 - lon1

982 if wrap:

983 u = wrap180(d)

984 if u != d:

985 return u, (lon1 + u)

986 return d, lon2

987

988

989 def _unrollon(p1, p2, wrap=False): # unroll180 == .karney._unroll2

990 '''(INTERNAL) Wrap/normalize, unroll and replace longitude.

991 '''

992 lat, lon = p2.lat, p2.lon

993 if wrap and _Wrap.normal:

994 lat, lon = _Wrap.latlon(lat, lon)

995 _, lon = unroll180(p1.lon, lon, wrap=True)

996 if lat != p2.lat or fabs(lon - p2.lon) > EPS:

997 p2 = p2.dup(lat=lat, lon=wrap180(lon))

998 # p2 = p2.copy(); p2.latlon = lat, wrap180(lon)

999 return p2

1000

1001

1002 def _unrollon3(p1, p2, p3, wrap=False):

1003 '''(INTERNAL) Wrap/normalize, unroll 2 points.

1004 '''

1005 w = wrap

1006 if w:

1007 w = _Wrap.normal

1008 p2 = _unrollon(p1, p2, wrap=w)

1009 p3 = _unrollon(p1, p3, wrap=w)

1010 p2 = _unrollon(p2, p3)

1011 return p2, p3, w # was wrapped?

1012

1013

1014 def unrollPI(rad1, rad2, wrap=True):

1015 '''Unroll longitudinal delta and wrap longitude in radians.

1016

1017 @arg rad1: Start longitude (C{radians}).

1018 @arg rad2: End longitude (C{radians}).

1019 @kwarg wrap: If C{True}, wrap and unroll to the M{(-PI..+PI]}

1020 range (C{bool}).

1021

1022 @return: 2-Tuple C{(B{rad2}-B{rad1}, B{rad2})} unrolled

1023 (C{radians}, C{radians}).

1024

1025 @see: Capability C{LONG_UNROLL} in U{GeographicLib

1026 <https://GeographicLib.SourceForge.io/html/python/interface.html#outmask>}.

1027 '''

1028 r = rad2 - rad1

1029 if wrap:

1030 u = wrapPI(r)

1031 if u != r:

1032 return u, (rad1 + u)

1033 return r, rad2

1034

1035

1036 def _valueError(where, x, raiser=True, **kwds):

1037 '''(INTERNAL) Return a C{_ValueError} or C{None}.

1038 '''

1039 t = _MODS.streprs.Fmt.PAREN(where.__name__, x)

1040 return _ValueError(t, **kwds) if raiser else None

1041

1042

1043 class _Wrap(object):

1044

1045 _normal = False # default

1046

1047 @property

1048 def normal(self):

1049 '''Get the current L{normal} setting (C{True},

1050 C{False} or C{None}).

1051 '''

1052 return self._normal

1053

1054 @normal.setter # PYCHOK setter!

1055 def normal(self, setting):

1056 '''Set L{normal} to C{True}, C{False} or C{None}.

1057 '''

1058 m = _MODS.formy

1059 t = {True: (m.normal, m.normal_),

1060 False: (self.wraplatlon, self.wraphilam),

1061 None: (_passargs, _passargs)}.get(setting, ())

1062 if t:

1063 self.latlon, self.philam = t

1064 self._normal = setting

1065

1066 def latlonDMS2(self, lat, lon, **DMS2_kwds):

1067 if isstr(lat) or isstr(lon):

1068 kwds = _xkwds(DMS2_kwds, clipLon=0, clipLat=0)

1069 lat, lon = _MODS.dms.parseDMS2(lat, lon, **kwds)

1070 return self.latlon(lat, lon)

1071

1072# def normalatlon(self, *latlon):

1073# return _MODS.formy.normal(*latlon)

1074

1075# def normalamphi(self, *philam):

1076# return _MODS.formy.normal_(*philam)

1077

1078 def wraplatlon(self, lat, lon):

1079 return wrap90(lat), wrap180(lon)

1080

1081 latlon = wraplatlon # default

1082

1083 def latlon3(self, lon1, lat2, lon2, wrap):

1084 if wrap:

1085 lat2, lon2 = self.latlon(lat2, lon2)

1086 lon21, lon2 = unroll180(lon1, lon2)

1087 else:

1088 lon21 = lon2 - lon1

1089 return lon21, lat2, lon2

1090

1091 def _latlonop(self, wrap):

1092 if wrap and self._normal is not None:

1093 return self.latlon

1094 else:

1095 return _passargs

1096

1097 def wraphilam(self, phi, lam):

1098 return wrapPI_2(phi), wrapPI(lam)

1099

1100 philam = wraphilam # default

1101

1102 def philam3(self, lam1, phi2, lam2, wrap):

1103 if wrap:

1104 phi2, lam2 = self.philam(phi2, lam2)

1105 lam21, lam2 = unrollPI(lam1, lam2)

1106 else:

1107 lam21 = lam2 - lam1

1108 return lam21, phi2, lam2

1109

1110 def _philamop(self, wrap):

1111 if wrap and self._normal is not None:

1112 return self.philam

1113 else:

1114 return _passargs

1115

1116 def point(self, ll, wrap=True): # in .points._fractional, -.PointsIter.iterate, ...

1117 '''Return C{ll} or a copy, I{normalized} or I{wrap}'d.

1118 '''

1119 if wrap and self._normal is not None:

1120 lat, lon = ll.latlon

1121 if fabs(lon) > 180 or fabs(lat) > 90:

1122 _n = self.latlon

1123 ll = ll.copy(name=_n.__name__)

1124 ll.latlon = _n(lat, lon)

1125 return ll

1126

1127_Wrap = _Wrap() # PYCHOK singleton

1128

1129

1130# def _wrap(angle, wrap, modulo):

1131# '''(INTERNAL) Angle wrapper M{((wrap-modulo)..+wrap]}.

1132#

1133# @arg angle: Angle (C{degrees}, C{radians} or C{grades}).

1134# @arg wrap: Range (C{degrees}, C{radians} or C{grades}).

1135# @arg modulo: Upper limit (360 C{degrees}, PI2 C{radians} or 400 C{grades}).

1136#

1137# @return: The B{C{angle}}, wrapped (C{degrees}, C{radians} or C{grades}).

1138# '''

1139# a = float(angle)

1140# if not (wrap - modulo) <= a < wrap:

1141# # math.fmod(-1.5, 3.14) == -1.5, but -1.5 % 3.14 == 1.64

1142# # math.fmod(-1.5, 360) == -1.5, but -1.5 % 360 == 358.5

1143# a %= modulo

1144# if a > wrap:

1145# a -= modulo

1146# return a

1147

1148

1149 def wrap90(deg):

1150 '''Wrap degrees to M{[-90..+90]}.

1151

1152 @arg deg: Angle (C{degrees}).

1153

1154 @return: Degrees, wrapped (C{degrees90}).

1155 '''

1156 return _wrapu(wrap180(deg), _180_0, _90_0)

1157

1158

1159 def wrap180(deg):

1160 '''Wrap degrees to M{[-180..+180]}.

1161

1162 @arg deg: Angle (C{degrees}).

1163

1164 @return: Degrees, wrapped (C{degrees180}).

1165 '''

1166 d = float(deg)

1167 w = _umod_360(d)

1168 if w > _180_0:

1169 w -= _360_0

1170 elif d < 0 and w == _180_0:

1171 w = -w

1172 return w

1173

1174

1175 def wrap360(deg): # see .streprs._umod_360

1176 '''Wrap degrees to M{[0..+360)}.

1177

1178 @arg deg: Angle (C{degrees}).

1179

1180 @return: Degrees, wrapped (C{degrees360}).

1181 '''

1182 return _umod_360(float(deg))

1183

1184

1185 def wrapPI(rad):

1186 '''Wrap radians to M{[-PI..+PI]}.

1187

1188 @arg rad: Angle (C{radians}).

1189

1190 @return: Radians, wrapped (C{radiansPI}).

1191 '''

1192 r = float(rad)

1193 w = _umod_PI2(r)

1194 if w > PI:

1195 w -= PI2

1196 elif r < 0 and w == PI:

1197 w = -PI

1198 return w

1199

1200

1201 def wrapPI2(rad):

1202 '''Wrap radians to M{[0..+2PI)}.

1203

1204 @arg rad: Angle (C{radians}).

1205

1206 @return: Radians, wrapped (C{radiansPI2}).

1207 '''

1208 return _umod_PI2(float(rad))

1209

1210

1211 def wrapPI_2(rad):

1212 '''Wrap radians to M{[-PI/2..+PI/2]}.

1213

1214 @arg rad: Angle (C{radians}).

1215

1216 @return: Radians, wrapped (C{radiansPI_2}).

1217 '''

1218 return _wrapu(wrapPI(rad), PI, PI_2)

1219

1220

1221# def wraplatlon(lat, lon):

1222# '''Both C{wrap90(B{lat})} and C{wrap180(B{lon})}.

1223# '''

1224# return wrap90(lat), wrap180(lon)

1225

1226

1227 def wrap_normal(*normal):

1228 '''Define the operation for the keyword argument C{B{wrap}=True},

1229 across L{pygeodesy}: I{wrap}, I{normalize} or I{no-op}. For

1230 backward compatibility, the default is I{wrap}.

1231

1232 @arg normal: If C{True}, I{normalize} lat- and longitude using

1233 L{normal} or L{normal_}, if C{False}, I{wrap} the

1234 lat- and longitude individually by L{wrap90} or

1235 L{wrapPI_2} respectively L{wrap180}, L{wrapPI} or

1236 if C{None}, leave lat- and longitude I{unchanged}.

1237 Do not supply any value to get the current setting.

1238

1239 @return: The previous L{wrap_normal} setting (C{bool} or C{None}).

1240 '''

1241 t = _Wrap.normal

1242 if normal:

1243 _Wrap.normal = normal[0]

1244 return t

1245

1246

1247# def wraphilam(phi, lam,):

1248# '''Both C{wrapPI_2(B{phi})} and C{wrapPI(B{lam})}.

1249# '''

1250# return wrapPI_2(phi), wrapPI(lam)

1251

1252

1253 def _wrapu(w, H, Q):

1254 '''(INTERNAL) Helper for functions C{wrap180} and C{wrapPI}.

1255 '''

1256 return (w - H) if w > Q else ((w + H) if w < (-Q) else w)

1257

1258

1259 def yard2m(yards):

1260 '''Convert I{UK} yards to meter.

1261

1262 @arg yards: Value in yards (C{scalar}).

1263

1264 @return: Value in C{meter} (C{float}).

1265

1266 @raise ValueError: Invalid B{C{yards}}.

1267 '''

1268 return Float(yards=yards) * _M_YARD_UK

1269

1270# **) MIT License

1271#

1272# Copyright (C) 2016-2025 -- mrJean1 at Gmail -- All Rights Reserved.

1273#

1274# Permission is hereby granted, free of charge, to any person obtaining a

1275# copy of this software and associated documentation files (the "Software"),

1276# to deal in the Software without restriction, including without limitation

1277# the rights to use, copy, modify, merge, publish, distribute, sublicense,

1278# and/or sell copies of the Software, and to permit persons to whom the

1279# Software is furnished to do so, subject to the following conditions:

1280#

1281# The above copyright notice and this permission notice shall be included

1282# in all copies or substantial portions of the Software.

1283#

1284# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

1285# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

1286# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

1287# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

1288# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

1289# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

1290# OTHER DEALINGS IN THE SOFTWARE.