-
Notifications
You must be signed in to change notification settings - Fork 1
/
faq_geo.py
74 lines (62 loc) · 2.59 KB
/
faq_geo.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
# coding: utf-8
import math
from typing import Tuple
class constants93:
GRS80E = 0.081819191042816
LONG_0 = 3
XS = 700000
YS = 12655612.0499
n = 0.7256077650532670
C = 11754255.4261
def lambert93_to_WGPS(lambertE: float, lambertN: float) -> Tuple[float, float]:
"""
Converts coordinates given in
`Lambert 93 <https://fr.wikipedia.org/wiki/Projection_conique_conforme_de_Lambert>`_
system, this system is used by `IGN <https://www.ign.fr/>`_
and their :epkg:`GEOFLA` file format.
:param lambertE: east
:param lambertN: north
:return: longitude, latitude
The function is inspired from
`lam93toLatLon.py <https://gist.github.com/flire/0a305eeec77bc84a73af8ddc8f9ec043>`_.
.. faqref::
:tag: geo
:title: Les fichiers GEOFLA ne contiennent pas de longitude, latitude ?
Les coordonnées contenues dans les fichiers :epkg:`GEOFLA`
ne sont pas toujours des longitudes, latitudes mais des coordonnées exprimées dans un système
de projection conique `Lambert 93 <https://fr.wikipedia.org/wiki/Projection_conique_conforme_de_Lambert>`_.
Il faut convertir les coordonnées avant de pouvoir tracer la carte ou changer la projection
utilisée par :epkg:`cartopy` :
`Lambert Conformal Projection <https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html#lambertconformal>`_.
"""
delX = lambertE - constants93.XS
delY = lambertN - constants93.YS
gamma = math.atan(-delX / delY)
R = math.sqrt(delX * delX + delY * delY)
latiso = math.log(constants93.C / R) / constants93.n
sinPhiit0 = math.tanh(
latiso + constants93.GRS80E * math.atanh(constants93.GRS80E * math.sin(1))
)
sinPhiit1 = math.tanh(
latiso + constants93.GRS80E * math.atanh(constants93.GRS80E * sinPhiit0)
)
sinPhiit2 = math.tanh(
latiso + constants93.GRS80E * math.atanh(constants93.GRS80E * sinPhiit1)
)
sinPhiit3 = math.tanh(
latiso + constants93.GRS80E * math.atanh(constants93.GRS80E * sinPhiit2)
)
sinPhiit4 = math.tanh(
latiso + constants93.GRS80E * math.atanh(constants93.GRS80E * sinPhiit3)
)
sinPhiit5 = math.tanh(
latiso + constants93.GRS80E * math.atanh(constants93.GRS80E * sinPhiit4)
)
sinPhiit6 = math.tanh(
latiso + constants93.GRS80E * math.atanh(constants93.GRS80E * sinPhiit5)
)
longRad = math.asin(sinPhiit6)
latRad = gamma / constants93.n + constants93.LONG_0 / 180 * math.pi
longitude = latRad / math.pi * 180
latitude = longRad / math.pi * 180
return longitude, latitude