Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Backport 9.3] proj_factor: fix when input is a compound CRS of a projected CRS, and… #3950

Merged
merged 2 commits into from
Nov 14, 2023
Merged

[Backport 9.3] proj_factor: fix when input is a compound CRS of a projected CRS, and… #3950

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e04c6a6
proj_factor: fix when input is a compound CRS of a projected CRS, and…
rouault Nov 10, 2023
dc01f9c
proj: projection factors: fix when input is a compound CRS of a proje…
rouault Nov 10, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 36 additions & 6 deletions src/4D_api.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
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
Open in desktop
Original file line number Diff line number Diff line change
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
Open in desktop
Original file line number Diff line number Diff line change
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
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Open in desktop
Original file line number Diff line number Diff line change
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