Skip to content

Commit f89e5a7

Browse files
JoostGevaertJoostGevaert
committed
Do additional tests for ellips - ortho transformations
1 parent 351526a commit f89e5a7

File tree

1 file changed

+226
-56
lines changed

1 file changed

+226
-56
lines changed

‎sandbox/pyproj_3d_transformations.py‎

Lines changed: 226 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
import marimo
22

3-
__generated_with = "0.13.11"
3+
__generated_with = "0.15.2"
44
app = marimo.App(width="medium")
55

66

@@ -9,91 +9,261 @@ def _():
99
import marimo as mo
1010
from pyproj import CRS, Transformer
1111
from pyproj.crs.crs import CompoundCRS
12-
return CRS, CompoundCRS, Transformer
12+
return CRS, CompoundCRS, Transformer, mo
1313

1414

1515
@app.cell
16-
def _(CRS, Transformer):
17-
# https://pyproj4.github.io/pyproj/stable/advanced_examples.html#promote-crs-to-3d
18-
wgs84 = CRS("EPSG:4326")
19-
swiss_proj = CRS("EPSG:2056")
20-
transformer = Transformer.from_crs(wgs84, swiss_proj, always_xy=True)
21-
# 2D Transformation
22-
print(f"2D transform of point {(8.37909, 47.01987, 1000)} from {wgs84} to {swiss_proj} gives:")
23-
print(transformer.transform(8.37909, 47.01987, 1000))
16+
def _(CRS, CompoundCRS):
17+
wgs84 = CRS(4326)
18+
wgs84_3d = wgs84.to_3d() # = CRS(4979)
19+
20+
egm2008_3855 = CRS(3855)
21+
wgs84_egm2008_9518 = CRS(9518)
22+
23+
world_mercator_3395 = CRS(3395)
24+
world_mercator_3d = world_mercator_3395.to_3d()
25+
world_mercator_egm2008_6893 = CRS(6893)
26+
27+
web_mercator_3857 = CRS(3857)
28+
web_mercator_3d = web_mercator_3857.to_3d()
2429

25-
wgs84_3d = wgs84.to_3d()
26-
swiss_proj_ellipsoidal_3d = swiss_proj.to_3d()
27-
transformer_ellipsoidal_3d = Transformer.from_crs(
30+
nl_rdnew_28992 = CRS(28992)
31+
nl_3d = nl_rdnew_28992.to_3d()
32+
nl_rdnew_nap_7415 = CRS(7415)
33+
34+
uk_grid_27700 = CRS(27700)
35+
uk_3d = uk_grid_27700.to_3d()
36+
uk_grid_odn_7405 = CRS(7405)
37+
38+
swiss_2056 = CRS(2056)
39+
swiss_3d = swiss_2056.to_3d()
40+
# https://pyproj4.github.io/pyproj/stable/build_crs.html#compound-crs
41+
swiss_lhn95_height = CRS("EPSG:5729")
42+
swiss_compound = CompoundCRS(
43+
name="CH1903+ / LV95 + LHN95 height",
44+
components=[swiss_2056, swiss_lhn95_height]
45+
)
46+
return (
47+
swiss_2056,
48+
swiss_3d,
49+
swiss_compound,
50+
swiss_lhn95_height,
51+
wgs84,
2852
wgs84_3d,
29-
swiss_proj_ellipsoidal_3d,
30-
always_xy=True,
53+
wgs84_egm2008_9518,
54+
world_mercator_3395,
55+
world_mercator_3d,
56+
world_mercator_egm2008_6893,
3157
)
32-
# 3D Transformation
33-
print(f"3D transform of point {(8.37909, 47.01987, 1000)} from {wgs84_3d.to_string()} to {swiss_proj_ellipsoidal_3d} gives:")
34-
print(transformer_ellipsoidal_3d.transform(8.37909, 47.01987, 1000))
3558

36-
return swiss_proj, transformer_ellipsoidal_3d, wgs84_3d
59+
60+
@app.cell
61+
def _(world_mercator_3d):
62+
world_mercator_3d
63+
return
64+
65+
66+
@app.cell
67+
def _(mo):
68+
mo.md(
69+
r"""
70+
## WGS 84 to WGS84 + EGM 2008 orthometric height
71+
72+
- WGS 84 (EPSG:4979)
73+
- WGS 84 (EPSG:4326) + EGM2008 orthometric height (EPSG:3855) => (EPSG:9518)
74+
"""
75+
)
76+
return
3777

3878

3979
@app.cell
4080
def _(
4181
CRS,
42-
CompoundCRS,
4382
Transformer,
44-
swiss_proj,
45-
transformer_ellipsoidal_3d,
83+
print_srs,
84+
wgs84,
4685
wgs84_3d,
86+
wgs84_egm2008_9518,
87+
world_mercator_3d,
4788
):
48-
# https://pyproj4.github.io/pyproj/stable/build_crs.html#compound-crs
49-
swiss_lhn95_height = CRS("EPSG:5729")
50-
swiss_compound = CompoundCRS(
51-
name="CH1903+ / LV95 + LHN95 height",
52-
components=[swiss_proj, swiss_lhn95_height]
89+
# 2D Transformation
90+
null_transform_2d = Transformer.from_crs(
91+
wgs84, wgs84, always_xy=True
92+
)
93+
print(
94+
f"2D transform of point {(8.37909, 47.01987, 1000)} from {print_srs(wgs84)} to {print_srs(wgs84)} gives:"
95+
)
96+
print(null_transform_2d.transform(8.37909, 47.01987, 1000))
97+
98+
# 3D Transformations
99+
null_transform_3d = Transformer.from_crs(
100+
CRS(4979), wgs84_3d, always_xy=True
101+
)
102+
print(
103+
f"3D transform of point {(8.37909, 47.01987, 1000)} from {print_srs(wgs84_3d)} to {print_srs(world_mercator_3d)} gives:"
104+
)
105+
print(
106+
null_transform_3d.transform(
107+
8.37909, 47.01987, 1000
108+
)
53109
)
54-
transformer_wgs84_3d_to_swiss_compound = Transformer.from_crs(
110+
111+
transformer_ellips_to_egm2008_ortho = Transformer.from_crs(
112+
wgs84_3d, wgs84_egm2008_9518, always_xy=True
113+
)
114+
print(
115+
f"3D transform of point {(8.37909, 47.01987, 1000)} from {print_srs(wgs84_3d)} to {print_srs(wgs84_egm2008_9518)} gives:"
116+
)
117+
print(
118+
transformer_ellips_to_egm2008_ortho.transform(
119+
8.37909, 47.01987, 1000
120+
)
121+
)
122+
return
123+
124+
125+
@app.cell(hide_code=True)
126+
def _(mo):
127+
mo.md(
128+
r"""
129+
## WGS 84 to World Mercator + EGM 2008 orthometric height
130+
131+
- WGS 84 (EPSG:4979)
132+
- World Mercator (EPSG:3395) + EGM2008 orthometric height (EPSG:3855) => (EPSG:6893)
133+
"""
134+
)
135+
return
136+
137+
138+
@app.cell
139+
def _(
140+
Transformer,
141+
print_srs,
142+
wgs84,
143+
wgs84_3d,
144+
world_mercator_3395,
145+
world_mercator_3d,
146+
world_mercator_egm2008_6893,
147+
):
148+
# 2D Transformation
149+
transformer2d_wgs84_to_world_mercator = Transformer.from_crs(
150+
wgs84, world_mercator_3395, always_xy=True
151+
)
152+
print(
153+
f"2D transform of point {(8.37909, 47.01987, 1000)} from {print_srs(wgs84)} to {print_srs(world_mercator_3395)} gives:"
154+
)
155+
print(transformer2d_wgs84_to_world_mercator.transform(8.37909, 47.01987, 1000))
156+
157+
# 3D Transformations
158+
transformer3d_wgs84_to_world_mercator_ellips = Transformer.from_crs(
159+
wgs84_3d, world_mercator_3d, always_xy=True
160+
)
161+
print(
162+
f"3D transform of point {(8.37909, 47.01987, 1000)} from {print_srs(wgs84_3d)} to {print_srs(world_mercator_3d)} gives:"
163+
)
164+
print(
165+
transformer3d_wgs84_to_world_mercator_ellips.transform(
166+
8.37909, 47.01987, 1000
167+
)
168+
)
169+
170+
transformer3d_wgs84_to_world_mercator_egm2008 = Transformer.from_crs(
171+
wgs84_3d, world_mercator_egm2008_6893, always_xy=True
172+
)
173+
print(
174+
f"3D transform of point {(8.37909, 47.01987, 1000)} from {print_srs(wgs84_3d)} to {print_srs(world_mercator_egm2008_6893)} gives:"
175+
)
176+
print(
177+
transformer3d_wgs84_to_world_mercator_egm2008.transform(
178+
8.37909, 47.01987, 1000
179+
)
180+
)
181+
return
182+
183+
184+
@app.cell(hide_code=True)
185+
def _(mo):
186+
mo.md(
187+
r"""
188+
## WGS 84 to Swiss national grid
189+
190+
- WGS 84 (EPSG:4979)
191+
- Swiss national grid projected SRS: CH1903+ / LV95 (EPSG:2056) + LHN95 height (EPSG:5729)
192+
193+
"""
194+
)
195+
return
196+
197+
198+
@app.cell
199+
def _(
200+
Transformer,
201+
print_srs,
202+
swiss_2056,
203+
swiss_3d,
204+
swiss_compound,
205+
swiss_lhn95_height,
206+
wgs84,
207+
wgs84_3d,
208+
):
209+
# https://pyproj4.github.io/pyproj/stable/advanced_examples.html#promote-crs-to-3d
210+
# 2D Transformation
211+
transformer2d_wgs84_to_swiss = Transformer.from_crs(wgs84, swiss_2056, always_xy=True)
212+
print(
213+
f"2D transform of point {(8.37909, 47.01987, 1000)} from {print_srs(wgs84)} to {print_srs(swiss_2056)} gives:"
214+
)
215+
print(transformer2d_wgs84_to_swiss.transform(8.37909, 47.01987, 1000))
216+
217+
# 3D Transformation
218+
transformer3d_wgs84_to_swiss_ellips = Transformer.from_crs(
219+
wgs84_3d,
220+
swiss_3d,
221+
always_xy=True,
222+
)
223+
print(
224+
f"3D transform of point {(8.37909, 47.01987, 1000)} from {print_srs(wgs84_3d)} to {print_srs(swiss_3d)} gives:"
225+
)
226+
print(transformer3d_wgs84_to_swiss_ellips.transform(8.37909, 47.01987, 1000))
227+
transformer3d_wgs84_to_swiss_ortho = Transformer.from_crs(
55228
wgs84_3d,
56229
swiss_compound,
57230
always_xy=True,
58231
)
59-
print(transformer_ellipsoidal_3d.transform(8.37909, 47.01987, 1000))
232+
print(
233+
f"3D transform of point {(8.37909, 47.01987, 1000)} from {print_srs(wgs84_3d)} to {print_srs(swiss_compound)} gives:"
234+
)
235+
print(transformer3d_wgs84_to_swiss_ortho.transform(8.37909, 47.01987, 1000))
236+
237+
# wgs84_3d
238+
# swiss_2056
239+
swiss_lhn95_height
240+
# swiss_3d
241+
# swiss_compound
60242
return
61243

62244

63245
@app.cell
64246
def _(CRS):
65-
uk_grid_27700 = CRS(27700)
66-
uk_3d = uk_grid_27700.to_3d()
67-
swiss_2056 = CRS("EPSG:2056")
68-
swiss_3d = swiss_2056.to_3d()
69-
egm2008_3855 = CRS(3855)
70-
rdnew_nap_7415 = CRS(7415)
71-
wgs84_egm2008_9518 = CRS(9518)
72-
crs_components = extract_crs_components(swiss_3d)
73-
crs_components
74-
return
247+
def print_srs(srs: CRS) -> str:
248+
axes = [axis.name for axis in srs.axis_info]
249+
return f"{srs.name} {axes}"
250+
return (print_srs,)
75251

76252

77-
@app.function
78-
def extract_crs_components(compound_crs):
79-
"""Extract horizontal and vertical CRS from a compound CRS"""
80-
if not compound_crs.is_compound:
81-
return compound_crs, None
253+
@app.cell
254+
def _(mo):
255+
mo.md(
256+
r"""
257+
from https://bertt.wordpress.com/2023/08/24/vertical-coordinate-reprojection-from-geoid-to-ellipsoid/
82258
83-
horizontal_crs = None
84-
vertical_crs = None
259+
When doing surveys in the field, it’s critical to know the correct vertical elevation. But the earth has the shape of a lumpy potato, so it’s not easy to calculate the vertical elevation.
85260
86-
for sub_crs in compound_crs.sub_crs_list:
87-
if sub_crs.is_projected or sub_crs.is_geographic:
88-
print(f"Horizontal CRS {sub_crs.name} is a {sub_crs.type_name} and has EPSG:{sub_crs.to_epsg()}.")
89-
horizontal_crs = sub_crs
90-
elif sub_crs.is_vertical:
91-
print(f"Vertical CRS {sub_crs.name} has EPSG:{sub_crs.to_epsg()}.")
92-
vertical_crs = sub_crs
93-
else:
94-
print(f"This CRS is not horizontal (projected or geographic) nor vertical: {sub_crs.type_name}")
261+
Traditionally the vertical elevation is calculated relative to a local system: the geoid. For the Netherlands the geoid ‘NAP height’ (https://www.rijkswaterstaat.nl/zakelijk/open-data/normaal-amsterdams-peil) is used – https://spatialreference.org/ref/epsg/5709/. The elevation is between minus 6.78 meter (Nieuwerkerk aan den IJssel) and 322.38 meter (Vaals).
95262
96-
return horizontal_crs, vertical_crs
263+
For example the Dam square in Amsterdam has about 2.6 meter elevation:
264+
"""
265+
)
266+
return
97267

98268

99269
if __name__ == "__main__":

0 commit comments

Comments
 (0)