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

Jump to bottom

proj_factors(): enhance speed when called repeatedly on same compound or projected CRS #4289

Merged
merged 1 commit into from
Oct 28, 2024
Merged

proj_factors(): enhance speed when called repeatedly on same compound or projected CRS #4289

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
4 changes: 4 additions & 0 deletions docs/source/development/reference/functions.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -826,6 +826,10 @@ Various
instantiated from a EPSG CRS code. The factors computed will be those of the
map projection implied by the transformation from the base geographic CRS of
the projected CRS to the projected CRS.
Starting with PROJ 9.6, to improve performance on repeated calls on a
projected CRS object, the above steps will modify the internal state of the
provided P object, and thus calling this function concurrently from multiple
threads on the same P object will no longer be supported.

The input geodetic coordinate lp should be such that lp.lam is the longitude
in radian, and lp.phi the latitude in radian (thus independently of the
Expand Down
167 changes: 88 additions & 79 deletions src/4D_api.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2882,90 +2882,99 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {
if (nullptr == P)
return factors;

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;
auto pj = P;
auto type = proj_get_type(pj);

PJ *horiz = nullptr;
if (pj->cached_op_for_proj_factors) {
pj = pj->cached_op_for_proj_factors;
} else {
if (type == PJ_TYPE_COMPOUND_CRS) {
horiz = proj_crs_get_sub_crs(pj->ctx, pj, 0);
pj = horiz;
type = proj_get_type(pj);
}
}

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
// using the same datum as the one of the projected CRS, and with
// input coordinates being in longitude, latitude order in radian,
// to be consistent with the expectations of the lp input parameter.
// We also need to create a modified projected CRS with a normalized
// easting/northing axis order in metre, so the resulting operation is
// just a single step pipeline with no axisswap or unitconvert steps.

auto ctx = P->ctx;
auto geodetic_crs = proj_get_source_crs(ctx, P);
assert(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);
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);
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 using the same datum as the one of the projected
// CRS, and with input coordinates being in longitude, latitude
// order in radian, to be consistent with the expectations of the lp
// input parameter. We also need to create a modified projected CRS
// with a normalized easting/northing axis order in metre, so the
// resulting operation is just a single step pipeline with no
// axisswap or unitconvert steps.

auto ctx = pj->ctx;
auto geodetic_crs = proj_get_source_crs(ctx, pj);
assert(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);
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, pj);
auto projCS = proj_create_cartesian_2D_cs(
ctx, PJ_CART2D_EASTING_NORTHING, "metre", 1.0);
auto projCRSNormalized = proj_create_projected_crs(
ctx, nullptr, geodetic_crs, conversion, projCS);
assert(projCRSNormalized);
proj_destroy(geodetic_crs);
proj_destroy(conversion);
proj_destroy(projCS);
auto newOp = proj_create_crs_to_crs_from_pj(
ctx, geogCRSNormalized, projCRSNormalized, nullptr, nullptr);
proj_destroy(geogCRSNormalized);
proj_destroy(projCRSNormalized);
assert(newOp);
// For debugging:
// printf("%s\n", proj_as_proj_string(ctx, newOp, PJ_PROJ_5,
// nullptr));

P->cached_op_for_proj_factors = newOp;
pj = newOp;
} else if (type != PJ_TYPE_CONVERSION &&
type != PJ_TYPE_TRANSFORMATION &&
type != PJ_TYPE_CONCATENATED_OPERATION &&
type != PJ_TYPE_OTHER_COORDINATE_OPERATION) {
proj_log_error(P, _("Invalid type for P object"));
proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
if (horiz)
proj_destroy(horiz);
return factors;
}
proj_destroy(cs);
auto conversion = proj_crs_get_coordoperation(ctx, P);
auto projCS = proj_create_cartesian_2D_cs(
ctx, PJ_CART2D_EASTING_NORTHING, "metre", 1.0);
auto projCRSNormalized = proj_create_projected_crs(
ctx, nullptr, geodetic_crs, conversion, projCS);
assert(projCRSNormalized);
proj_destroy(geodetic_crs);
proj_destroy(conversion);
proj_destroy(projCS);
auto newOp = proj_create_crs_to_crs_from_pj(
ctx, geogCRSNormalized, projCRSNormalized, nullptr, nullptr);
proj_destroy(geogCRSNormalized);
proj_destroy(projCRSNormalized);
assert(newOp);
// For debugging:
// printf("%s\n", proj_as_proj_string(ctx, newOp, PJ_PROJ_5, nullptr));
auto ret = proj_factors(newOp, lp);
proj_destroy(newOp);
return ret;
}

if (type != PJ_TYPE_CONVERSION && type != PJ_TYPE_TRANSFORMATION &&
type != PJ_TYPE_CONCATENATED_OPERATION &&
type != PJ_TYPE_OTHER_COORDINATE_OPERATION) {
proj_log_error(P, _("Invalid type for P object"));
proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
return factors;
}

if (pj_factors(lp.lp, P, 0.0, &f))
const int ret = pj_factors(lp.lp, P, pj, 0.0, &f);
if (horiz)
proj_destroy(horiz);
if (ret)
return factors;

factors.meridional_scale = f.h;
Expand Down
50 changes: 22 additions & 28 deletions src/factors.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,38 +12,31 @@

#define EPS 1.0e-12

int pj_factors(PJ_LP lp, const PJ *P, double h, struct FACTORS *fac) {
int pj_factors(PJ_LP lp, PJ *toplevel, const PJ *internal, double h,
struct FACTORS *fac) {
double cosphi, t, n, r;
int err;
PJ_COORD coo = {{0, 0, 0, 0}};
coo.lp = lp;

/* Failing the 3 initial checks will most likely be due to */
/* earlier errors, so we leave errno alone */
if (nullptr == fac)
return 1;

if (nullptr == P)
return 1;

if (HUGE_VAL == lp.lam)
return 1;

/* But from here, we're ready to make our own mistakes */
err = proj_errno_reset(P);
err = proj_errno_reset(toplevel);

/* Indicate that all factors are numerical approximations */
fac->code = 0;

/* Check for latitude or longitude overange */
if ((fabs(lp.phi) - M_HALFPI) > EPS) {
proj_log_error(P, _("Invalid latitude"));
proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
proj_log_error(toplevel, _("Invalid latitude"));
proj_errno_set(toplevel, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
return 1;
}
if (fabs(lp.lam) > 10.) {
proj_log_error(P, _("Invalid longitude"));
proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
proj_log_error(toplevel, _("Invalid longitude"));
proj_errno_set(toplevel, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
return 1;
}

Expand All @@ -53,23 +46,23 @@ int pj_factors(PJ_LP lp, const PJ *P, double h, struct FACTORS *fac) {
h = DEFAULT_H;

/* If input latitudes are geocentric, convert to geographic */
if (P->geoc)
lp = pj_geocentric_latitude(P, PJ_INV, coo).lp;
if (internal->geoc)
lp = pj_geocentric_latitude(internal, PJ_INV, coo).lp;

/* If latitude + one step overshoots the pole, move it slightly inside, */
/* so the numerical derivative still exists */
if (fabs(lp.phi) > (M_HALFPI - h))
lp.phi = lp.phi < 0. ? -(M_HALFPI - h) : (M_HALFPI - h);

/* Longitudinal distance from central meridian */
lp.lam -= P->lam0;
if (!P->over)
lp.lam -= internal->lam0;
if (!internal->over)
lp.lam = adjlon(lp.lam);

/* Derivatives */
if (pj_deriv(lp, h, P, &(fac->der))) {
proj_log_error(P, _("Invalid latitude or longitude"));
proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
if (pj_deriv(lp, h, internal, &(fac->der))) {
proj_log_error(toplevel, _("Invalid latitude or longitude"));
proj_errno_set(toplevel, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
return 1;
}

Expand All @@ -78,13 +71,13 @@ int pj_factors(PJ_LP lp, const PJ *P, double h, struct FACTORS *fac) {
fac->h = hypot(fac->der.x_p, fac->der.y_p);
fac->k = hypot(fac->der.x_l, fac->der.y_l) / cosphi;

if (P->es != 0.0) {
if (internal->es != 0.0) {
t = sin(lp.phi);
t = 1. - P->es * t * t;
t = 1. - internal->es * t * t;
n = sqrt(t);
fac->h *= t * n / P->one_es;
fac->h *= t * n / internal->one_es;
fac->k *= n;
r = t * t / P->one_es;
r = t * t / internal->one_es;
} else
r = 1.;

Expand All @@ -96,7 +89,7 @@ int pj_factors(PJ_LP lp, const PJ *P, double h, struct FACTORS *fac) {
cosphi;

/* Meridian-parallel angle (theta prime) */
fac->thetap = aasin(P->ctx, fac->s / (fac->h * fac->k));
fac->thetap = aasin(internal->ctx, fac->s / (fac->h * fac->k));

/* Tissot ellipse axis */
t = fac->k * fac->k + fac->h * fac->h;
Expand All @@ -107,8 +100,9 @@ int pj_factors(PJ_LP lp, const PJ *P, double h, struct FACTORS *fac) {
fac->a = 0.5 * (fac->a + t);

/* Angular distortion */
fac->omega = 2. * aasin(P->ctx, (fac->a - fac->b) / (fac->a + fac->b));
fac->omega =
2. * aasin(internal->ctx, (fac->a - fac->b) / (fac->a + fac->b));

proj_errno_restore(P, err);
proj_errno_restore(toplevel, err);
return 0;
}
2 changes: 2 additions & 0 deletions src/malloc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -166,6 +166,8 @@ PJ *pj_default_destructor(PJ *P, int errlev) { /* Destructor */
proj_destroy(P->hgridshift);
proj_destroy(P->vgridshift);

proj_destroy(P->cached_op_for_proj_factors);

free(static_cast<struct pj_opaque *>(P->opaque));
delete P;
return nullptr;
Expand Down
6 changes: 5 additions & 1 deletion src/proj_internal.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -675,6 +675,9 @@ struct PJconsts {
true; /* to remove in PROJ 10? */
bool skipNonInstantiable = true;

// Used internally by proj_factors()
PJ *cached_op_for_proj_factors = nullptr;

/*************************************************************************************

E N D O F G E N E R A L P A R A M E T E R S T R U C T
Expand Down Expand Up @@ -916,7 +919,8 @@ COMPLEX pj_zpoly1(COMPLEX, const COMPLEX *, int);
COMPLEX pj_zpolyd1(COMPLEX, const COMPLEX *, int, COMPLEX *);

int pj_deriv(PJ_LP, double, const PJ *, struct DERIVS *);
int pj_factors(PJ_LP, const PJ *, double, struct FACTORS *);
int pj_factors(PJ_LP, PJ *toplevel, const PJ *internal, double,
struct FACTORS *);

void *proj_mdist_ini(double);
double proj_mdist(double, double, double, const void *);
Expand Down
14 changes: 12 additions & 2 deletions test/unit/gie_self_tests.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -611,9 +611,19 @@ TEST(gie, info_functions) {
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);
{
const auto factors2 = proj_factors(P, c);
EXPECT_NEAR(factors2.meridional_scale, 1 - 28.54587730 * 1e-5,
1e-8);
}

// Try again to test caching of internal objects
{
const auto factors2 = proj_factors(P, c);
EXPECT_NEAR(factors2.meridional_scale, 1 - 28.54587730 * 1e-5,
1e-8);
}

proj_destroy(P);
}
Expand Down