Skip to content

Commit

Permalink
Merge pull request #3950 from OSGeo/backport-3949-to-9.3
Browse files Browse the repository at this point in the history 
[Backport 9.3] proj_factor: fix when input is a compound CRS of a projected CRS, and…
  • Loading branch information
@rouault
rouault committed Nov 14, 2023
2 parents ab34bf5 + dc01f9c commit d331569
Show file tree
Hide file tree
Showing 5 changed files with 123 additions and 15 deletions.
42 changes: 36 additions & 6 deletions src/4D_api.cpp
View file
Expand Up @@ -2886,6 +2886,16 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {

const auto type = proj_get_type(P);

if (type == PJ_TYPE_COMPOUND_CRS) {
auto ctx = P->ctx;
auto horiz = proj_crs_get_sub_crs(ctx, P, 0);
if (horiz) {
auto ret = proj_factors(horiz, lp);
proj_destroy(horiz);
return ret;
}
}

if (type == PJ_TYPE_PROJECTED_CRS) {
// If it is a projected CRS, then compute the factors on the conversion
// associated to it. We need to start from a temporary geographic CRS
Expand All @@ -2899,14 +2909,34 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {
auto ctx = P->ctx;
auto geodetic_crs = proj_get_source_crs(ctx, P);
assert(geodetic_crs);
auto datum = proj_crs_get_datum(ctx, geodetic_crs);
auto datum_ensemble = proj_crs_get_datum_ensemble(ctx, geodetic_crs);
auto pm = proj_get_prime_meridian(ctx, geodetic_crs);
double pm_longitude = 0;
proj_prime_meridian_get_parameters(ctx, pm, &pm_longitude, nullptr,
nullptr);
proj_destroy(pm);
PJ *geogCRSNormalized;
auto cs = proj_create_ellipsoidal_2D_cs(
ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, "Radian", 1.0);
auto geogCRSNormalized = proj_create_geographic_crs_from_datum(
ctx, "unnamed crs", datum ? datum : datum_ensemble, cs);
proj_destroy(datum);
proj_destroy(datum_ensemble);
if (pm_longitude != 0) {
auto ellipsoid = proj_get_ellipsoid(ctx, geodetic_crs);
double semi_major_metre = 0;
double inv_flattening = 0;
proj_ellipsoid_get_parameters(ctx, ellipsoid, &semi_major_metre,
nullptr, nullptr, &inv_flattening);
geogCRSNormalized = proj_create_geographic_crs(
ctx, "unname crs", "unnamed datum", proj_get_name(ellipsoid),
semi_major_metre, inv_flattening, "reference prime meridian", 0,
nullptr, 0, cs);
proj_destroy(ellipsoid);
} else {
auto datum = proj_crs_get_datum(ctx, geodetic_crs);
auto datum_ensemble =
proj_crs_get_datum_ensemble(ctx, geodetic_crs);
geogCRSNormalized = proj_create_geographic_crs_from_datum(
ctx, "unnamed crs", datum ? datum : datum_ensemble, cs);
proj_destroy(datum);
proj_destroy(datum_ensemble);
}
proj_destroy(cs);
auto conversion = proj_crs_get_coordoperation(ctx, P);
auto projCS = proj_create_cartesian_2D_cs(
Expand Down
51 changes: 42 additions & 9 deletions src/apps/proj.cpp
View file
Expand Up @@ -537,7 +537,20 @@ int main(int argc, char **argv) {
eargc--;
// logic copied from proj_factors function
if (PJ *P = proj_create(nullptr, ocrs.c_str())) {
const auto type = proj_get_type(P);
auto type = proj_get_type(P);
auto ctx = P->ctx;
if (type == PJ_TYPE_COMPOUND_CRS) {
auto horiz = proj_crs_get_sub_crs(ctx, P, 0);
if (horiz) {
if (proj_get_type(horiz) == PJ_TYPE_PROJECTED_CRS) {
proj_destroy(P);
P = horiz;
type = proj_get_type(P);
} else {
proj_destroy(horiz);
}
}
}
if (type == PJ_TYPE_PROJECTED_CRS) {
try {
auto crs = dynamic_cast<const NS_PROJ::crs::ProjectedCRS *>(
Expand All @@ -548,18 +561,38 @@ int main(int argc, char **argv) {
dir == NS_PROJ::cs::AxisDirection::SOUTH;
} catch (...) {
}
auto ctx = P->ctx;
auto geodetic_crs = proj_get_source_crs(ctx, P);
assert(geodetic_crs);
auto datum = proj_crs_get_datum(ctx, geodetic_crs);
auto datum_ensemble =
proj_crs_get_datum_ensemble(ctx, geodetic_crs);
auto pm = proj_get_prime_meridian(ctx, geodetic_crs);
double pm_longitude = 0;
proj_prime_meridian_get_parameters(ctx, pm, &pm_longitude,
nullptr, nullptr);
proj_destroy(pm);
PJ *geogCRSNormalized;
auto cs = proj_create_ellipsoidal_2D_cs(
ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, "Radian", 1.0);
auto geogCRSNormalized = proj_create_geographic_crs_from_datum(
ctx, "unnamed crs", datum ? datum : datum_ensemble, cs);
proj_destroy(datum);
proj_destroy(datum_ensemble);
if (pm_longitude != 0) {
auto ellipsoid = proj_get_ellipsoid(ctx, geodetic_crs);
double semi_major_metre = 0;
double inv_flattening = 0;
proj_ellipsoid_get_parameters(ctx, ellipsoid,
&semi_major_metre, nullptr,
nullptr, &inv_flattening);
geogCRSNormalized = proj_create_geographic_crs(
ctx, "unname crs", "unnamed datum",
proj_get_name(ellipsoid), semi_major_metre,
inv_flattening, "reference prime meridian", 0, nullptr,
0, cs);
proj_destroy(ellipsoid);
} else {
auto datum = proj_crs_get_datum(ctx, geodetic_crs);
auto datum_ensemble =
proj_crs_get_datum_ensemble(ctx, geodetic_crs);
geogCRSNormalized = proj_create_geographic_crs_from_datum(
ctx, "unnamed crs", datum ? datum : datum_ensemble, cs);
proj_destroy(datum);
proj_destroy(datum_ensemble);
}
proj_destroy(cs);
Proj = proj_create_crs_to_crs_from_pj(ctx, geogCRSNormalized, P,
nullptr, nullptr);
Expand Down
15 changes: 15 additions & 0 deletions test/cli/testproj
View file
Expand Up @@ -41,6 +41,21 @@ echo "Test CRS option"
$EXE EPSG:32620 -S >>${OUT} <<EOF
-63 0
EOF
echo "Test projection factors on projected CRS with non-Greenwhich prime meridian"
#
$EXE EPSG:27571 -S >>${OUT} <<EOF
2.33722917 49.5
EOF
echo "Test projection factors on compound CRS with a projected CRS"
#
$EXE EPSG:5972 -S >>${OUT} <<EOF
9 0
EOF

# On some architectures the angular distortion (omega) of EPSG:27571 is
# not exactly 0, but 8.53878e-07
cat ${OUT} | sed "s/8.53878e-07/0/" > ${OUT}.tmp
mv ${OUT}.tmp ${OUT}

#
# do 'diff' with distribution results
Expand Down
2 changes: 2 additions & 0 deletions test/cli/testproj_out.dist
View file
@@ -1,2 +1,4 @@
2 49 2.000 49.000
500000.00 0.00 <0.9996 0.9996 0.9992 0 0.9996 0.9996>
600000.00 1200000.00 <0.999877 0.999877 0.999755 0 0.999877 0.999877>
500000.00 0.00 <0.9996 0.9996 0.9992 0 0.9996 0.9996>
28 changes: 28 additions & 0 deletions test/unit/gie_self_tests.cpp
View file
Expand Up @@ -590,6 +590,34 @@ TEST(gie, info_functions) {
proj_destroy(P);
}

// Test with a projected CRS whose datum has a non-Greenwich prime meridian
{
P = proj_create(PJ_DEFAULT_CTX, "EPSG:27571");

PJ_COORD c;
c.lp.lam = proj_torad(0.0689738);
c.lp.phi = proj_torad(49.508567);
const auto factors2 = proj_factors(P, c);

EXPECT_NEAR(factors2.meridional_scale, 1 - 12.26478760 * 1e-5, 1e-8);

proj_destroy(P);
}

// Test with a compound CRS of a projected CRS
{
P = proj_create(PJ_DEFAULT_CTX, "EPSG:5972");

PJ_COORD c;
c.lp.lam = proj_torad(10.729030600);
c.lp.phi = proj_torad(59.916494500);
const auto factors2 = proj_factors(P, c);

EXPECT_NEAR(factors2.meridional_scale, 1 - 28.54587730 * 1e-5, 1e-8);

proj_destroy(P);
}

// Test with a geographic CRS --> error
{
P = proj_create(PJ_DEFAULT_CTX, "EPSG:4326");
Expand Down

0 comments on commit d331569

Please sign in to comment.