1

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

3

4 u'''Various utility functions.

5

6 After I{(C) Chris Veness 2011-2015} 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, _passargs

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

15 _float as _F, _isfinite, isnan, isnear0, _over, \

16 _umod_360, _umod_PI2, _M_KM, _M_NM, _M_SM, _0_0, \

17 _1__90, _0_5, _1_0, _N_1_0, _2__PI, _10_0, _90_0, \

18 _N_90_0, _180_0, _N_180_0, _360_0, _400_0

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

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

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

22 from pygeodesy.units import Degrees, Degrees_, Feet, Float, Lam, Lam_, \

23 Meter, Meter2, Radians, Radians_

24

25 from math import acos, asin, atan2, cos, degrees, fabs, radians, sin, tan # pow

26

27__all__ = _ALL_LAZY.utily

28__version__ = '24.01.12'

29

30# read constant name "_M_Unit" as "meter per Unit"

31_M_CHAIN = _F( 20.1168) # yard2m(1) * 22

32_M_FATHOM = _F( 1.8288) # yard2m(1) * 2 or _M_NM * 1e-3

33_M_FOOT = _F( 0.3048) # Int'l (1 / 3.2808398950131 = 10_000 / (254 * 12))

34_M_FOOT_GE = _F( 0.31608) # German Fuss (1 / 3.1637560111364)

35_M_FOOT_FR = _F( 0.3248406) # French Pied-du-Roi or pied (1 / 3.0784329298739)

36_M_FOOT_USVY = _F( 0.3048006096012192) # US Survey (1200 / 3937)

37_M_FURLONG = _F( 201.168) # 220 * yard2m(1) == 10 * m2chain(1)

38# _M_KM = _F(1000.0) # kilo meter

39# _M_NM = _F(1852.0) # nautical mile

40# _M_SM = _F(1609.344) # statute mile

41_M_TOISE = _F( 1.9490436) # French toise, 6 pieds (6 / 3.0784329298739)

42_M_YARD_UK = _F( 0.9144) # 254 * 12 * 3 / 10_000 == 3 * ft2m(1) Int'l

43

44

45 def _abs1nan(x):

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

47 '''

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

49

50

51 def acos1(x):

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

53 '''

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

55

56

57 def acre2ha(acres):

58 '''Convert acres to hectare.

59

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

61

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

63

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

65 '''

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

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

68

69

70 def acre2m2(acres):

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

72

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

74

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

76

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

78 '''

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

80 return Meter2(Float(acres) * 4046.8564224)

81

82

83 def asin1(x):

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

85 '''

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

87

88

89 def atan1(y, x=_1_0):

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

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

92 '''

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

94

95

96 def atan1d(y, x=_1_0):

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

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

99

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

101 '''

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

103

104

105 def atan2b(y, x):

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

107

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

109 '''

110 b = atan2d(y, x)

111 if b < 0:

112 b += _360_0

113 return b

114

115

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

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

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

119

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

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

122 '''

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

124 if y < 0: # q = 3

125 d = degrees(atan2(x, -y)) - _90_0

126 else: # q = 2

127 d = _90_0 - degrees(atan2(x, y))

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

129 return NAN

130 elif y:

131 if x > 0: # q = 0

132 d = degrees(atan2(y, x))

133 elif x < 0: # q = 1

134 d = _copysign(_180_0, y) - degrees(atan2(y, -x))

135 else: # x == 0

136 d = _copysign(_90_0, y)

137 else:

138 d = _180_0 if x < 0 else _0_0

139 return _azireversed(d) if reverse else d

140

141

142 def _azireversed(azimuth): # in .rhumbBase

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

144 '''

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

146

147

148 def chain2m(chains):

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

150

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

152

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

154

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

156 '''

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

158

159

160 def circle4(earth, lat):

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

162

163 @arg earth: The earth radius, ellipsoid or datum

164 (C{meter}, L{Ellipsoid}, L{Ellipsoid2},

165 L{Datum} or L{a_f2Tuple}).

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

167

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

169 instance.

170

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

172 L{pygeodesy.rangerrors} set to C{True}.

173

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

175

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

177 '''

178 E = _MODS.datums._earth_ellipsoid(earth)

179 return E.circle4(lat)

180

181

182 def cot(rad, **error_kwds):

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

184

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

186 @kwarg error_kwds: Error to raise (C{ValueError}).

187

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

189

190 @raise ValueError: L{pygeodesy.isnear0}C{(sin(B{rad})}.

191 '''

192 s, c = sincos2(rad)

193 if c:

194 if isnear0(s):

195 raise _valueError(cot, rad, **error_kwds)

196 c = c / s # /= chokes PyChecker

197 return c

198

199

200 def cot_(*rads, **error_kwds):

201 '''Return the C{cotangent} of angle(s) in C{radiansresection}.

202

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

204 @kwarg error_kwds: Error to raise (C{ValueError}).

205

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

207

208 @raise ValueError: See L{pygeodesy.cot}.

209 '''

210 try:

211 for r in rads:

212 yield cot(r)

213 except ValueError:

214 raise _valueError(cot_, r, **error_kwds)

215

216

217 def cotd(deg, **error_kwds):

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

219

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

221 @kwarg error_kwds: Error to raise (C{ValueError}).

222

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

224

225 @raise ValueError: L{pygeodesy.isnear0}C{(sin(B{deg})}.

226 '''

227 s, c = sincos2d(deg)

228 if c:

229 if isnear0(s):

230 raise _valueError(cotd, deg, **error_kwds)

231 c = c / s # /= chokes PyChecker

232 elif s < 0:

233 c = -c # negate-0

234 return c

235

236

237 def cotd_(*degs, **error_kwds):

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

239

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

241 @kwarg error_kwds: Error to raise (C{ValueError}).

242

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

244

245 @raise ValueError: See L{pygeodesy.cotd}.

246 '''

247 try:

248 for d in degs:

249 yield cotd(d)

250 except ValueError:

251 raise _valueError(cotd_, d, **error_kwds)

252

253

254 def degrees90(rad):

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

256

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

258

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

260 '''

261 return wrap90(degrees(rad))

262

263

264 def degrees180(rad):

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

266

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

268

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

270 '''

271 return wrap180(degrees(rad))

272

273

274 def degrees360(rad):

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

276

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

278

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

280 '''

281 return _umod_360(degrees(rad))

282

283

284 def degrees2grades(deg):

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

286

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

288

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

290 '''

291 return Degrees(deg) * _400_0 / _360_0

292

293

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

295 '''Convert an angle to a distance along the equator or

296 along the parallel at an other (geodetic) latitude.

297

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

299 @kwarg radius: Mean earth radius, ellipsoid or datum

300 (C{meter}, L{Ellipsoid}, L{Ellipsoid2},

301 L{Datum} or L{a_f2Tuple}).

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

303

304 @return: Distance (C{meter}, same units as B{C{radius}}

305 or equatorial and polar radii) or C{0.0} for

306 near-polar B{C{lat}}.

307

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

309 L{pygeodesy.rangerrors} set to C{True}.

310

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

312

313 @raise ValueError: Invalid B{C{deg}}, B{C{radius}} or

314 B{C{lat}}.

315

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

317 '''

318 return _Radians2m(Lam_(deg=deg, clip=0), radius, lat)

319

320

321 def fathom2m(fathoms):

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

323

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

325

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

327

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

329

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

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

332 '''

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

334

335

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

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

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

339

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

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

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

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

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

345

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

347

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

349 '''

350 return Meter(Feet(feet) * (_M_FOOT_USVY if usurvey else

351 (_M_FOOT_FR if pied else

352 (_M_FOOT_GE if fuss else _M_FOOT))))

353

354

355 def furlong2m(furlongs):

356 '''Convert a furlong to meter.

357

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

359

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

361

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

363 '''

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

365

366

367 def grades(rad):

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

369

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

371

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

373 '''

374 return Float(grades=Float(rad=rad) * _400_0 / PI2)

375

376

377 def grades400(rad):

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

379

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

381

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

383 '''

384 return Float(grades400=(grades(rad) % _400_0) or _0_0) # _umod_400

385

386

387 def grades2degrees(gon):

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

389

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

391

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

393 '''

394 return Degrees(Float(gon=gon) * _360_0 / _400_0)

395

396

397 def grades2radians(gon):

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

399

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

401

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

403 '''

404 return Radians(Float(gon=gon) * PI2 / _400_0)

405

406

407 def km2m(km):

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

409

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

411

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

413

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

415 '''

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

417

418

419 def _loneg(lon):

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

421 '''

422 return _180_0 - lon

423

424

425 def m2chain(meter):

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

427

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

429

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

431

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

433 '''

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

435

436

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

438 '''Convert a distance to an angle along the equator or

439 along the parallel at an other (geodetic) latitude.

440

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

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

443 an L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or

444 L{a_f2Tuple}).

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

446

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

448

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

450 L{pygeodesy.rangerrors} set to C{True}.

451

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

453

454 @raise ValueError: Invalid B{C{distance}}, B{C{radius}}

455 or B{C{lat}}.

456

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

458 '''

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

460

461

462 def m2fathom(meter):

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

464

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

466

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

468

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

470

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

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

473 '''

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

475

476

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

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

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

480

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

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

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

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

485 I{International} foot.

486

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

488

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

490 '''

491 # * 3.2808333333333333, US Survey 3937 / 1200

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

493 return Float(feet=Meter(meter) / (_M_FOOT_USVY if usurvey else

494 (_M_FOOT_FR if pied else

495 (_M_FOOT_GE if fuss else _M_FOOT))))

496

497

498 def m2furlong(meter):

499 '''Convert meter to furlongs.

500

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

502

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

504

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

506 '''

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

508

509

510 def m2km(meter):

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

512

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

514

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

516

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

518 '''

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

520

521

522 def m2NM(meter):

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

524

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

526

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

528

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

530 '''

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

532

533

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

535 '''Convert a distance to an angle along the equator or along the

536 parallel at an other (geodetic) latitude.

537

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

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

540 an L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or

541 L{a_f2Tuple}).

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

543

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

545

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

547 L{pygeodesy.rangerrors} set to C{True}.

548

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

550

551 @raise ValueError: Invalid B{C{distance}}, B{C{radius}}

552 or B{C{lat}}.

553

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

555 '''

556 m = circle4(radius, lat).radius

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

558

559

560 def m2SM(meter):

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

562

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

564

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

566

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

568 '''

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

570

571

572 def m2toise(meter):

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

574

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

576

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

578

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

580

581 @see: Function L{m2fathom}.

582 '''

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

584

585

586 def m2yard(meter):

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

588

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

590

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

592

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

594 '''

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

596

597

598 def NM2m(nm):

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

600

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

602

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

604

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

606 '''

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

608

609

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

611 '''Convert an angle to a distance along the equator or

612 along the parallel at an other (geodetic) latitude.

613

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

615 @kwarg radius: Mean earth radius, ellipsoid or datum

616 (C{meter}, L{Ellipsoid}, L{Ellipsoid2},

617 L{Datum} or L{a_f2Tuple}).

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

619

620 @return: Distance (C{meter}, same units as B{C{radius}}

621 or equatorial and polar radii) or C{0.0} for

622 near-polar B{C{lat}}.

623

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

625 L{pygeodesy.rangerrors} set to C{True}.

626

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

628

629 @raise ValueError: Invalid B{C{rad}}, B{C{radius}} or

630 B{C{lat}}.

631

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

633 '''

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

635

636

637 def _Radians2m(rad, radius, lat):

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

639 '''

640 m = circle4(radius, lat).radius

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

642

643

644 def radiansPI(deg):

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

646

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

648

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

650 '''

651 return wrapPI(radians(deg))

652

653

654 def radiansPI2(deg):

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

656

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

658

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

660 '''

661 return _umod_PI2(radians(deg))

662

663

664 def radiansPI_2(deg):

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

666

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

668

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

670 '''

671 return wrapPI_2(radians(deg))

672

673

674 def _sin0cos2(q, r, sign):

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

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

677 '''

678 if r < PI_2:

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

680 t = s, c, -s, -c, s

681 else: # r == PI_2

682 t = _1_0, _0_0, _N_1_0, _0_0, _1_0

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

684# t = _0_0, _1_0, _0_0, _N_1_0, _0_0

685# q &= 3

686 s = t[q] or (NEG0 if sign < 0 else _0_0)

687 c = t[q + 1] or _0_0

688 return s, c

689

690

691 def SinCos2(x):

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

693

694 @arg x: Angle (L{Degrees}, L{Degrees_}, L{Radians}, L{Radians_}

695 or scalar C{radians}).

696

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

698 '''

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

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

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

702

703

704 def sincos2(rad):

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

706

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

708

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

710

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

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

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

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

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

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

717 '''

718 if _isfinite(rad):

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

720 if q < 0:

721 q -= 1

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

723 else:

724 t = NAN, NAN

725 return t

726

727

728 def sincos2_(*rads):

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

730

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

732

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

734

735 @see: function L{sincos2}.

736 '''

737 for r in rads:

738 s, c = sincos2(r)

739 yield s

740 yield c

741

742

743 def sincos2d(deg, **adeg):

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

745

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

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

748

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

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

751

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

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

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

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

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

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

758 '''

759 if _isfinite(deg):

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

761 if q < 0:

762 q -= 1

763 d = deg - q * _90_0

764 if adeg:

765 t = _xkwds_get(adeg, adeg=_0_0)

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

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

768 else:

769 t = NAN, NAN

770 return t

771

772

773 def sincos2d_(*degs):

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

775

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

777

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

779

780 @see: Function L{sincos2d}.

781 '''

782 for d in degs:

783 s, c = sincos2d(d)

784 yield s

785 yield c

786

787

788 def sincostan3(rad):

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

790

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

792

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

794

795 @see: Function L{sincos2}.

796 '''

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

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

799 return s, c, t

800

801

802 def SM2m(sm):

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

804

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

806

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

808

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

810 '''

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

812

813

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

815 '''Compute the tangent of half angle.

816

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

818 @kwarg semi: Angle or edge name and index

819 for semi-circular error.

820

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

822

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

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

825 '''

826 # .formy.excessKarney_, .sphericalTrigonometry.areaOf

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

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

829 break

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

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

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

833

834 return tan(rad * _0_5) if _isfinite(rad) else NAN

835

836

837 def tand(deg, **error_kwds):

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

839

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

841 @kwarg error_kwds: Error to raise (C{ValueError}).

842

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

844

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

846 '''

847 s, c = sincos2d(deg)

848 if s:

849 if isnear0(c):

850 raise _valueError(tand, deg, **error_kwds)

851 s = s / c # /= chokes PyChecker

852 elif c < 0:

853 s = -s # negate-0

854 return s

855

856

857 def tand_(*degs, **error_kwds):

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

859

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

861 @kwarg error_kwds: Error to raise (C{ValueError}).

862

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

864

865 @raise ValueError: See L{pygeodesy.tand}.

866 '''

867 for d in degs:

868 yield tand(d, **error_kwds)

869

870

871 def tanPI_2_2(rad):

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

873

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

875

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

877 '''

878 return tan((rad + PI_2) * _0_5) if _isfinite(rad) else (

879 NAN if isnan(rad) else (_N_90_0 if rad < 0 else _90_0))

880

881

882 def toise2m(toises):

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

884

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

886

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

888

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

890

891 @see: Function L{fathom2m}.

892 '''

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

894

895

896 def truncate(x, ndigits=None):

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

898

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

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

901 aka I{precision}.

902

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

904

905 @see: Python function C{round}.

906 '''

907 if isint(ndigits):

908 p = _10_0**ndigits

909 x = int(x * p) / p

910 return x

911

912

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

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

915

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

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

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

919 range (C{bool}).

920

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

922 C{degrees}).

923

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

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

926 '''

927 d = lon2 - lon1

928 if wrap:

929 u = wrap180(d)

930 if u != d:

931 return u, (lon1 + u)

932 return d, lon2

933

934

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

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

937 '''

938 lat, lon = p2.lat, p2.lon

939 if wrap and _Wrap.normal:

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

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

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

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

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

945 return p2

946

947

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

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

950 '''

951 w = wrap

952 if w:

953 w = _Wrap.normal

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

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

956 p2 = _unrollon(p2, p3)

957 return p2, p3, w # was wrapped?

958

959

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

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

962

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

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

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

966 range (C{bool}).

967

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

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

970

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

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

973 '''

974 r = rad2 - rad1

975 if wrap:

976 u = wrapPI(r)

977 if u != r:

978 return u, (rad1 + u)

979 return r, rad2

980

981

982 def _valueError(where, x, **kwds):

983 '''(INTERNAL) Return a C{_ValueError}.

984 '''

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

986 return _ValueError(t, **kwds)

987

988

989 class _Wrap(object):

990

991 _normal = False # default

992

993 @property

994 def normal(self):

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

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

997 '''

998 return self._normal

999

1000 @normal.setter # PYCHOK setter!

1001 def normal(self, setting):

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

1003 '''

1004 t = {True: (_MODS.formy.normal, _MODS.formy.normal_),

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

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

1007 if t:

1008 self.latlon, self.philam = t

1009 self._normal = setting

1010

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

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

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

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

1015 return self.latlon(lat, lon)

1016

1017# def normalatlon(self, *latlon):

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

1019

1020# def normalamphi(self, *philam):

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

1022

1023 def wraplatlon(self, lat, lon):

1024 return wrap90(lat), wrap180(lon)

1025

1026 latlon = wraplatlon # default

1027

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

1029 if wrap:

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

1031 lon21, lon2 = unroll180(lon1, lon2)

1032 else:

1033 lon21 = lon2 - lon1

1034 return lon21, lat2, lon2

1035

1036 def _latlonop(self, wrap):

1037 if wrap and self._normal is not None:

1038 return self.latlon

1039 else:

1040 return _passargs

1041

1042 def wraphilam(self, phi, lam):

1043 return wrapPI_2(phi), wrapPI(lam)

1044

1045 philam = wraphilam # default

1046

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

1048 if wrap:

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

1050 lam21, lam2 = unrollPI(lam1, lam2)

1051 else:

1052 lam21 = lam2 - lam1

1053 return lam21, phi2, lam2

1054

1055 def _philamop(self, wrap):

1056 if wrap and self._normal is not None:

1057 return self.philam

1058 else:

1059 return _passargs

1060

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

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

1063 '''

1064 if wrap and self._normal is not None:

1065 lat, lon = ll.latlon

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

1067 _n = self.latlon

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

1069 ll.latlon = _n(lat, lon)

1070 return ll

1071

1072_Wrap = _Wrap() # PYCHOK singleton

1073

1074

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

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

1077#

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

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

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

1081#

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

1083# '''

1084# a = float(angle)

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

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

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

1088# a %= modulo

1089# if a > wrap:

1090# a -= modulo

1091# return a

1092

1093

1094 def wrap90(deg):

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

1096

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

1098

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

1100 '''

1101 w = wrap180(deg)

1102 return (w - _180_0) if w > 90 else ((w + _180_0) if w < -90 else w)

1103

1104

1105 def wrap180(deg):

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

1107

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

1109

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

1111 '''

1112 d = float(deg)

1113 w = _umod_360(d)

1114 if w > 180:

1115 w -= _360_0

1116 elif d < 0 and w == 180:

1117 w = -w

1118 return w

1119

1120

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

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

1123

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

1125

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

1127 '''

1128 return _umod_360(float(deg))

1129

1130

1131 def wrapPI(rad):

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

1133

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

1135

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

1137 '''

1138 r = float(rad)

1139 w = _umod_PI2(r)

1140 if w > PI:

1141 w -= PI2

1142 elif r < 0 and w == PI:

1143 w = -PI

1144 return w

1145

1146

1147 def wrapPI2(rad):

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

1149

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

1151

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

1153 '''

1154 return _umod_PI2(float(rad))

1155

1156

1157 def wrapPI_2(rad):

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

1159

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

1161

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

1163 '''

1164 w = wrapPI(rad)

1165 return (w - PI) if w > PI_2 else ((w + PI) if w < (-PI_2) else w)

1166

1167

1168# def wraplatlon(lat, lon):

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

1170# '''

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

1172

1173

1174 def wrap_normal(*normal):

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

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

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

1178

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

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

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

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

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

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

1185

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

1187 '''

1188 t = _Wrap.normal

1189 if normal:

1190 _Wrap.normal = normal[0]

1191 return t

1192

1193

1194# def wraphilam(phi, lam,):

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

1196# '''

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

1198

1199

1200 def yard2m(yards):

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

1202

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

1204

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

1206

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

1208 '''

1209 return Float(yards=yards) * _M_YARD_UK

1210

1211# **) MIT License

1212#

1213# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

1214#

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

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

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

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

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

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

1221#

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

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

1224#

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

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

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

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

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

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

1231# OTHER DEALINGS IN THE SOFTWARE.