Skip to content

Commit

Permalink
Merge pull request #4289 from rouault/proj_factors_perf_enhancements
Browse files Browse the repository at this point in the history 
proj_factors(): enhance speed when called repeatedly on same compound or projected CRS
  • Loading branch information
@rouault
rouault authored Oct 28, 2024
2 parents 252b13c + 7e945f8 commit 36328ee
Show file tree
Hide file tree
Showing 6 changed files with 133 additions and 110 deletions.
4 changes: 4 additions & 0 deletions docs/source/development/reference/functions.rst
View file
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
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
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
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
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
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

0 comments on commit 36328ee

Please sign in to comment.