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
|
import pytest
from numpy.testing import assert_almost_equal
from pyproj import Proj, transform
# illustrates the use of the transform function to
# perform coordinate transformations with datum shifts.
#
# This example is from Roberto Vidmar
#
# Test point is Trieste, Molo Sartorio
#
# This data come from the Istituto Geografico Militare (IGM), as well as
# the 7 parameters to transform from Gauss-Boaga (our reference frame)
# to WGS84
#
# WGS84 Lat: 45d38'49.879" (45.647188611)
# WGS84 Lon: 13d45'34.397" (13.759554722)
# WGS84 z: 52.80;
# UTM 33: 403340.97 5055597.17
# GB: 2423346.99 5055619.87
UTM_x = 403340.9672367854
UTM_y = 5055597.175553089
GB_x = 2423346.99
GB_y = 5055619.87
WGS84_lat = 45.647188611 # Degrees
WGS84_lon = 13.759554722 # Degrees
UTM_z = WGS84_z = 52.8 # Ellipsoidical height in meters
WGS84_PROJ = Proj(proj="latlong", datum="WGS84")
UTM_33_PROJ = Proj(proj="utm", zone="33")
with pytest.warns(FutureWarning):
GAUSSSB_PROJ = Proj(
init="epsg:3004", towgs84="-122.74,-34.27,-22.83,-1.884,-3.400,-3.030,-15.62"
)
def test_shift_wgs84_to_utm33():
with pytest.warns(FutureWarning):
xutm33, yutm33, zutm33 = transform(
WGS84_PROJ, UTM_33_PROJ, WGS84_lon, WGS84_lat, WGS84_z
)
assert_almost_equal((xutm33, yutm33, zutm33), (UTM_x, UTM_y, UTM_z))
def test_shift_utm33_to_wgs84():
with pytest.warns(FutureWarning):
back_lon, back_lat, back_z = transform(
UTM_33_PROJ, WGS84_PROJ, UTM_x, UTM_y, UTM_z
)
assert_almost_equal((back_lon, back_lat, back_z), (WGS84_lon, WGS84_lat, WGS84_z))
def test_shift_wgs84_to_gaussb_no_ellipsoidal_height():
with pytest.warns(FutureWarning):
xgb, ygb, zgb = transform(WGS84_PROJ, GAUSSSB_PROJ, WGS84_lon, WGS84_lat, 0)
assert_almost_equal((xgb, ygb, zgb), (GB_x, 5055619.899, 0), decimal=2)
def test_shift_gaussb_to_wgs84_no_ellipsoidal_height():
with pytest.warns(FutureWarning):
back_lon, back_lat, back_z = transform(GAUSSSB_PROJ, WGS84_PROJ, GB_x, GB_y, 0)
assert_almost_equal(
(back_lon, back_lat, back_z), (WGS84_lon, WGS84_lat, 0), decimal=3
)
|
