Skip to content

Commit fd6e432

Browse files
rouaultrouault
committed
gie: various fixes around crs_src/crs_dst support
and bootstrap test/gie/epsg_grid.gie and test/gie/epsg_no_grid.gie Fixes #4740
1 parent 56e879b commit fd6e432

5 files changed

Lines changed: 226 additions & 18 deletions

File tree

‎src/apps/gie.cpp‎

Lines changed: 130 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -104,6 +104,8 @@ Thomas Knudsen, thokn@sdfe.dk, 2017-10-01/2017-10-08
104104
105105
***********************************************************************/
106106

107+
#define FROM_PROJ_CPP
108+
107109
#include <ctype.h>
108110
#include <errno.h>
109111
#include <math.h>
@@ -116,9 +118,13 @@ Thomas Knudsen, thokn@sdfe.dk, 2017-10-01/2017-10-08
116118
#include "proj.h"
117119
#include "proj_internal.h"
118120
#include "proj_strtod.h"
121+
122+
#include <algorithm>
119123
#include <cmath> /* for isnan */
120124
#include <math.h>
121125

126+
#include <proj/internal/internal.hpp>
127+
122128
#include "optargpm.h"
123129

124130
/* Package for flexible format I/O - ffio */
@@ -192,6 +198,7 @@ typedef struct {
192198
char crs_dst[MAX_OPERATION + 1];
193199
char crs_src[MAX_OPERATION + 1];
194200
PJ *P;
201+
bool crs_dst_is_lat_lon_or_y_x;
195202
PJ_COORD a, b, e;
196203
PJ_DIRECTION dir;
197204
int verbosity;
@@ -561,6 +568,21 @@ static int require_grid(const char *args) {
561568
const char *grid_filename = column(args, 1);
562569
grid_info = proj_grid_info(grid_filename);
563570
if (strlen(grid_info.filename) == 0) {
571+
572+
if (proj_context_is_network_enabled(nullptr)) {
573+
auto dbContext = NS_PROJ::io::DatabaseContext::create();
574+
std::string fullFilename, packageName, url;
575+
bool directDownload, openLicense, gridAvailable = false;
576+
if (dbContext->lookForGridInfo(
577+
grid_filename,
578+
/* considerKnownGridsAsAvailable = */ true, fullFilename,
579+
packageName, url, directDownload, openLicense,
580+
gridAvailable) &&
581+
gridAvailable) {
582+
return 0;
583+
}
584+
}
585+
564586
if (T.verbosity > 1) {
565587
fprintf(T.fout, "Test skipped because of missing grid %s\n",
566588
grid_filename);
@@ -644,6 +666,56 @@ static int operation(const char *args) {
644666
return 0;
645667
}
646668

669+
static bool isLatOrNorthingFirst(const PJ *crs) {
670+
auto crsType = proj_get_type(crs);
671+
672+
if (crsType == PJ_TYPE_COMPOUND_CRS) {
673+
auto horiz_crs = proj_crs_get_sub_crs(nullptr, crs, 0);
674+
assert(horiz_crs);
675+
bool ret = isLatOrNorthingFirst(horiz_crs);
676+
proj_destroy(horiz_crs);
677+
return ret;
678+
}
679+
680+
if (crsType != PJ_TYPE_GEOGRAPHIC_2D_CRS &&
681+
crsType != PJ_TYPE_GEOGRAPHIC_3D_CRS &&
682+
crsType != PJ_TYPE_GEOCENTRIC_CRS && crsType != PJ_TYPE_PROJECTED_CRS) {
683+
fprintf(stderr, "Bug in gie.cpp:%d: unsupported CRS type %d!\n",
684+
__LINE__, crsType);
685+
return false;
686+
}
687+
688+
auto cs = proj_crs_get_coordinate_system(nullptr, crs);
689+
assert(cs);
690+
691+
const char *axisName = "";
692+
const char *unitName = "";
693+
proj_cs_get_axis_info(nullptr, cs, 0,
694+
&axisName, // name,
695+
nullptr, // abbrev
696+
nullptr, // direction
697+
nullptr, // unit conv factor
698+
&unitName, // unit name
699+
nullptr, // unit authority
700+
nullptr // unit code
701+
);
702+
const bool isLatOrNorthingFirst =
703+
NS_PROJ::internal::ci_find(std::string(axisName), "latitude") !=
704+
std::string::npos ||
705+
NS_PROJ::internal::ci_find(std::string(axisName), "northing") !=
706+
std::string::npos;
707+
proj_destroy(cs);
708+
709+
if (strcmp(unitName, "degree") != 0 && strcmp(unitName, "metre") != 0) {
710+
// We should add normalization to this unit in expect()
711+
fprintf(stderr, "Bug in gie.cpp:%d: unsupported unit %s!\n", __LINE__,
712+
unitName);
713+
return false;
714+
}
715+
716+
return isLatOrNorthingFirst;
717+
}
718+
647719
static int crs_to_crs_operation() {
648720
T.op_id++;
649721
T.operation_lineno = F->lineno;
@@ -671,6 +743,14 @@ static int crs_to_crs_operation() {
671743
proj_errno_reset(nullptr);
672744
proj_context_use_proj4_init_rules(nullptr, T.use_proj4_init_rules);
673745

746+
T.crs_dst_is_lat_lon_or_y_x = false;
747+
PJ *pj_dst = proj_create(nullptr, T.crs_dst);
748+
if (!pj_dst) {
749+
fprintf(stderr, "Cannot instantiate crs_dst = %s\n", T.crs_dst);
750+
}
751+
T.crs_dst_is_lat_lon_or_y_x = isLatOrNorthingFirst(pj_dst);
752+
proj_destroy(pj_dst);
753+
674754
T.P = proj_create_crs_to_crs(nullptr, T.crs_src, T.crs_dst, nullptr);
675755

676756
strcpy(T.crs_src, "");
@@ -925,18 +1005,6 @@ static int expect_failure_with_errno_message(int expected, int got) {
9251005
return 1;
9261006
}
9271007

928-
/* For test purposes, we want to call a transformation of the same */
929-
/* dimensionality as the number of dimensions given in accept */
930-
static PJ_COORD expect_trans_n_dim(const PJ_COORD &ci) {
931-
if (4 == T.dimensions_given_at_last_accept)
932-
return proj_trans(T.P, T.dir, ci);
933-
934-
if (3 == T.dimensions_given_at_last_accept)
935-
return pj_approx_3D_trans(T.P, T.dir, ci);
936-
937-
return pj_approx_2D_trans(T.P, T.dir, ci);
938-
}
939-
9401008
/*****************************************************************************/
9411009
static int expect(const char *args) {
9421010
/*****************************************************************************
@@ -989,7 +1057,7 @@ static int expect(const char *args) {
9891057
/* Try to carry out the operation - and expect failure */
9901058
ci =
9911059
proj_angular_input(T.P, T.dir) ? torad_coord(T.P, T.dir, T.a) : T.a;
992-
co = expect_trans_n_dim(ci);
1060+
co = proj_trans(T.P, T.dir, ci);
9931061

9941062
if (expect_failure_with_errno) {
9951063
if (proj_errno(T.P) == expect_failure_with_errno)
@@ -1043,7 +1111,7 @@ static int expect(const char *args) {
10431111

10441112
/* do the transformation, but mask off dimensions not given in expect-ation
10451113
*/
1046-
co = expect_trans_n_dim(ci);
1114+
co = proj_trans(T.P, T.dir, ci);
10471115
if (T.dimensions_given < 4)
10481116
co.v[3] = 0;
10491117
if (T.dimensions_given < 3)
@@ -1056,23 +1124,68 @@ static int expect(const char *args) {
10561124
co.v[1], co.v[2], co.v[3]);
10571125

10581126
#if 0
1059-
/* We need to handle unusual axis orders - that'll be an item for version 5.1 */
1127+
/* We need to handle unusual axis orders when axisswap explicitly present.
1128+
* Otherwise when using dst_crs code below will take care of that.
1129+
*/
10601130
if (T.P->axisswap) {
10611131
ce = proj_trans (T.P->axisswap, T.dir, ce);
10621132
co = proj_trans (T.P->axisswap, T.dir, co);
10631133
}
10641134
#endif
1135+
10651136
if (std::isnan(co.v[0]) && std::isnan(ce.v[0])) {
10661137
d = 0.0;
1067-
} else if (proj_angular_output(T.P, T.dir)) {
1138+
} else if (proj_angular_output(T.P, T.dir)) { // angular = radians...
1139+
d = proj_lpz_dist(T.P, ce, co);
1140+
} else if (proj_degree_output(T.P, T.dir)) {
1141+
co.v[0] = proj_torad(co.v[0]);
1142+
co.v[1] = proj_torad(co.v[1]);
1143+
1144+
ce.v[0] = proj_torad(ce.v[0]);
1145+
ce.v[1] = proj_torad(ce.v[1]);
1146+
1147+
if (T.crs_dst_is_lat_lon_or_y_x) {
1148+
std::swap(co.v[0], co.v[1]);
1149+
std::swap(ce.v[0], ce.v[1]);
1150+
}
1151+
10681152
d = proj_lpz_dist(T.P, ce, co);
10691153
} else {
1070-
d = proj_xyz_dist(co, ce);
1154+
1155+
if (T.crs_dst_is_lat_lon_or_y_x) {
1156+
std::swap(co.v[0], co.v[1]);
1157+
std::swap(ce.v[0], ce.v[1]);
1158+
}
1159+
1160+
d = proj_xyz_dist(ce, co);
10711161
}
10721162

10731163
// Test written like that to handle NaN
10741164
if (!(d <= T.tolerance))
10751165
return expect_message(d, args);
1166+
1167+
// Somewhat arbitrary temporal theshold but should be fine for all intended
1168+
// purposes.
1169+
constexpr double TEMPORAL_THESHOLD_IN_YEAR = 1e-4;
1170+
if (T.dimensions_given == 4 &&
1171+
std::fabs(ce.v[3] - co.v[3]) > TEMPORAL_THESHOLD_IN_YEAR) {
1172+
another_failure();
1173+
1174+
if (T.verbosity < 0)
1175+
return 1;
1176+
if (0 == T.op_ko && T.verbosity < 2)
1177+
banner(T.operation);
1178+
fprintf(T.fout, "%s", T.op_ko ? " -----\n" : delim);
1179+
fprintf(T.fout, " FAILURE in %s(%d):\n",
1180+
opt_strip_path(T.curr_file), (int)F->lineno);
1181+
fprintf(T.fout, " expected: %s\n", args);
1182+
fprintf(T.fout, " got: %.12f %.12f %.9f %.9f\n",
1183+
T.b.v[0], T.b.v[1], T.b.v[2], T.b.v[3]);
1184+
fprintf(T.fout, " deviation: %.4f year, %.4f maximum allowed\n",
1185+
std::fabs(ce.v[3] - co.v[3]), TEMPORAL_THESHOLD_IN_YEAR);
1186+
return 1;
1187+
}
1188+
10761189
succs++;
10771190

10781191
another_success();

‎test/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -79,6 +79,7 @@ proj_add_gie_test("guyou" "gie/guyou.gie")
7979
proj_add_gie_test("peirce_q" "gie/peirce_q.gie")
8080
proj_add_gie_test("tinshift" "gie/tinshift.gie")
8181
proj_add_gie_test("spilhaus" "gie/spilhaus.gie")
82+
proj_add_gie_test("epsg_no_grid" "gie/epsg_no_grid.gie")
8283

8384
if(TIFF_ENABLED)
8485
proj_add_gie_test("Deformation" "gie/deformation.gie")
@@ -91,6 +92,7 @@ endif()
9192

9293
if(TIFF_ENABLED AND CURL_ENABLED AND RUN_NETWORK_DEPENDENT_TESTS)
9394
proj_add_gie_network_dependent_test("nkg" "gie/nkg.gie")
95+
proj_add_gie_network_dependent_test("epsg_grid" "gie/epsg_grid.gie")
9496
endif()
9597

9698
# GIGS tests. Uncommented tests are expected to fail due to issues with

‎test/gie/deformation.gie‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -152,7 +152,7 @@ expect -147.0 64.0 0.0 2011.0
152152
roundtrip 100
153153

154154
accept -147.0 64.0 0.0 2011.0
155-
expect -147.0 64.0 0.0 2020.0
155+
expect -147.0 64.0 0.0 2011.0
156156
roundtrip 100
157157

158158
-------------------------------------------------------------------------------

‎test/gie/epsg_grid.gie‎

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
### This file must contain only (crs_src, crs_dst) tuples where the best transformation
2+
### DOES involve the use a grid/json/etc. resource file
3+
### Otherwise use epsg_no_grid.gie
4+
5+
<gie-strict>
6+
7+
###################################
8+
#### Generic / worldwide scope ####
9+
###################################
10+
11+
# Test EGM2008 grid
12+
crs_src EPSG:4979 # WGS 84 geographic 3D
13+
crs_dst EPSG:9518 # WGS 84 + EGM2008
14+
tolerance 0.1 mm
15+
require_grid us_nga_egm08_25.tif
16+
17+
accept 55.0 12.0 10.0
18+
expect 55.0 12.0 -27.8257
19+
20+
################
21+
#### FRANCE ####
22+
################
23+
24+
crs_src EPSG:9785 # RGF93 v2b + NGF-IGN69 height
25+
crs_dst EPSG:9781 # RGF93 v2b geographic 3D
26+
tolerance 0.1 mm
27+
require_grid fr_ign_RAF20.tif
28+
29+
accept 43.6109 3.8761 0.0
30+
expect 43.6109 3.8761 49.6904
31+
32+
</gie-strict>

‎test/gie/epsg_no_grid.gie‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
### This file must contain only (crs_src, crs_dst) tuples where the best transformation
2+
### does NOT involve the use a grid/json/etc. resource file
3+
### Otherwise use epsg_grid.gie
4+
5+
<gie-strict>
6+
7+
###################################
8+
#### Generic / worldwide scope ####
9+
###################################
10+
11+
# ETRS89 -> ETRS89/UTM32
12+
crs_src EPSG:4258
13+
crs_dst EPSG:25832
14+
tolerance 0.1 mm
15+
accept 55.0 12.0
16+
expect 691875.6321 6098907.8250
17+
18+
# WGS 84 geographic 3D -> WGS 84 geocentric
19+
crs_src EPSG:4979
20+
crs_dst EPSG:4978
21+
tolerance 0.1 mm
22+
accept 55.0 12.0 10.0
23+
expect 3586475.2672 762328.8513 5201391.7147
24+
25+
# WGS 84 geocentric -> WGS 84 geographic 3D
26+
crs_src EPSG:4979
27+
crs_dst EPSG:4978
28+
tolerance 0.1 mm
29+
expect 3586475.2672 762328.8513 5201391.7147
30+
accept 55.0 12.0 10.0
31+
32+
###################
33+
#### AUSTRALIA ####
34+
###################
35+
36+
# GDA2020 -> ITRF2014, using EPSG:8049 "ITRF2014 to GDA2020 (1)" 15-parameter Helmert with 2020.0 central epoch
37+
crs_src EPSG:7843 # GDA2020
38+
crs_dst EPSG:7912 # ITRF2014
39+
tolerance 0.1 mm
40+
accept -33.8623 151.2077 0.0 2020.0
41+
expect -33.8623 151.2077 0.0 2020.0
42+
43+
crs_src EPSG:7843 # GDA2020
44+
crs_dst EPSG:7912 # ITRF2014
45+
tolerance 0.1 mm
46+
accept -33.8623 151.2077 0.0 2026.0
47+
expect -33.862297056 151.207701181 -0.00101 2026.0
48+
49+
###################
50+
###### FINLAND ####
51+
###################
52+
53+
# Test northing, easting output (actually more a gie test itself, than a EPSG one)
54+
crs_src EPSG:4123 # Finland KKJ
55+
crs_dst EPSG:2393 # Finland YKJ Northing, Easting
56+
57+
tolerance 0.1 mm
58+
accept 60.1699 24.9384
59+
expect 6674944.7742 3385559.8151
60+
61+
</gie-strict>

0 commit comments

Comments
 (0)