forked from OSGeo/gdal
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtest_osr_util.py
155 lines (118 loc) · 6.76 KB
/
test_osr_util.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
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
#!/usr/bin/env pytest
# -*- coding: utf-8 -*-
###############################################################################
#
# Project: GDAL/OGR Test Suite
# Purpose: osgeo_utils.auxiliary.osr_utils (gdal-utils) testing
# Author: Idan Miara <[email protected]>
#
###############################################################################
# Copyright (c) 2021, Idan Miara <[email protected]>
#
# SPDX-License-Identifier: MIT
###############################################################################
import pytest
from osgeo import gdal, osr
# test that osgeo_utils and numpy are available, otherwise skip all tests
pytest.importorskip("osgeo_utils")
from osgeo_utils.auxiliary import array_util, osr_util
def test_gis_order():
pj4326_gis2 = osr_util.get_srs(4326) # tests the correct default
if gdal.GetConfigOption("OSR_DEFAULT_AXIS_MAPPING_STRATEGY") is None:
assert pj4326_gis2.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
pj4326_auth = osr_util.get_srs(4326, axis_order=osr.OAMS_AUTHORITY_COMPLIANT)
assert pj4326_auth.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
pj4326_gis = osr_util.get_srs(4326, axis_order=osr.OAMS_TRADITIONAL_GIS_ORDER)
assert pj4326_gis.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
assert not osr_util.are_srs_equivalent(pj4326_gis, pj4326_auth)
if gdal.GetConfigOption("OSR_DEFAULT_AXIS_MAPPING_STRATEGY") is None:
assert osr_util.are_srs_equivalent(pj4326_auth, 4326)
assert not osr_util.are_srs_equivalent(pj4326_gis, 4326)
pj4326_str = osr_util.get_srs_pj(pj4326_auth)
pj4326_str2 = osr_util.get_srs_pj(pj4326_gis)
# axis order is not reflected in proj strings
assert isinstance(pj4326_str, str) and pj4326_str == pj4326_str2
if gdal.GetConfigOption("OSR_DEFAULT_AXIS_MAPPING_STRATEGY") is None:
assert osr_util.are_srs_equivalent(pj4326_str, 4326)
assert osr_util.are_srs_equivalent(pj4326_auth, pj4326_str)
assert not osr_util.are_srs_equivalent(pj4326_gis, pj4326_str)
osr_util.set_default_axis_order(osr.OAMS_TRADITIONAL_GIS_ORDER) # sets gis order
srs = osr_util.get_srs(4326) # check the default was changed
if gdal.GetConfigOption("OSR_DEFAULT_AXIS_MAPPING_STRATEGY") is None:
assert srs.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
# check that srs object is not affected by default
srs = osr_util.get_srs(pj4326_auth)
assert srs.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
assert osr_util.are_srs_equivalent(srs, pj4326_auth)
# check that srs object is also affected if explicitly set
srs = osr_util.get_srs(pj4326_auth, axis_order=osr.OAMS_TRADITIONAL_GIS_ORDER)
assert srs.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
# check that the default does not effect explicit order
srs = osr_util.get_srs(4326, axis_order=osr.OAMS_TRADITIONAL_GIS_ORDER)
assert srs.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
srs = osr_util.get_srs(pj4326_str)
assert srs.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
srs = osr_util.get_srs(4326, axis_order=osr.OAMS_AUTHORITY_COMPLIANT)
assert srs.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
assert not osr_util.are_srs_equivalent(pj4326_gis, pj4326_auth)
assert not osr_util.are_srs_equivalent(pj4326_auth, 4326)
assert osr_util.are_srs_equivalent(pj4326_gis, 4326)
# restore the default and repeat some tests
osr_util.set_default_axis_order()
srs = osr_util.get_srs(4326) # check the default was restored
if gdal.GetConfigOption("OSR_DEFAULT_AXIS_MAPPING_STRATEGY") is None:
assert srs.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
srs = osr_util.get_srs(pj4326_str)
if gdal.GetConfigOption("OSR_DEFAULT_AXIS_MAPPING_STRATEGY") is None:
assert srs.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
# check that srs object is not affected by default
srs = osr_util.get_srs(pj4326_gis) # check that srs object is also affected
assert srs.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
# check that srs object is also affected if explicitly set
srs = osr_util.get_srs(pj4326_gis, axis_order=osr.OAMS_AUTHORITY_COMPLIANT)
assert srs.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
srs = osr_util.get_srs(4326, axis_order=osr.OAMS_TRADITIONAL_GIS_ORDER)
assert srs.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
def test_gis_order2():
pj4326_gis = osr_util.get_srs(4326, axis_order=osr.OAMS_TRADITIONAL_GIS_ORDER)
pj4326_str = osr_util.get_srs_pj(4326)
srs_from_epsg = osr.SpatialReference()
srs_from_epsg.ImportFromEPSG(4326)
if gdal.GetConfigOption("OSR_DEFAULT_AXIS_MAPPING_STRATEGY") is None:
assert srs_from_epsg.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
srs_from_str = osr.SpatialReference()
srs_from_str.ImportFromProj4(pj4326_str)
if gdal.GetConfigOption("OSR_DEFAULT_AXIS_MAPPING_STRATEGY") is None:
assert srs_from_str.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
assert srs_from_epsg.IsSame(srs_from_str)
# testing that explicitly setting OAMS_AUTHORITY_COMPLIANT does not effect equivalence
srs_from_epsg.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
srs_from_str.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
assert srs_from_epsg.IsSame(srs_from_str)
srs_from_epsg.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
srs_from_str.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
srs_from_epsg2 = osr.SpatialReference()
srs_from_epsg2.ImportFromEPSG(4326)
srs_from_epsg2.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
assert srs_from_epsg.IsSame(srs_from_epsg2)
test_gis_order_proj_str_vs_epsg = False
# explicitly setting OAMS_TRADITIONAL_GIS_ORDER triggers inequality between srs from proj-string and srs from epsg
# if this issue is resolved these tests can be enabled
if test_gis_order_proj_str_vs_epsg:
assert srs_from_epsg.IsSame(srs_from_str)
assert osr_util.are_srs_equivalent(pj4326_str, 4326)
assert osr_util.are_srs_equivalent(pj4326_gis, pj4326_str)
def test_transform():
pj_utm = osr_util.get_srs(32636)
utm_x = [690950.4640, 688927.6381]
utm_y = [3431318.8435, 3542183.4911]
for gis_order in (False, True):
axis_order = osr_util.get_axis_order_from_gis_order(gis_order)
pj4326 = osr_util.get_srs(4326, axis_order=axis_order)
ct = osr_util.get_transform(pj4326, pj_utm)
lon = [35, 35]
lat = [31, 32]
x, y = (lon, lat) if gis_order else (lat, lon)
osr_util.transform_points(ct, x, y)
d = array_util.array_dist(x, utm_x), array_util.array_dist(y, utm_y)
assert max(d) < 0.01