1212
1313#define EPS 1.0e-12
1414
15- int pj_factors (PJ_LP lp, PJ *toplevel, const PJ *internal, double h,
16- struct FACTORS *fac) {
15+ static int pj_factors (PJ_LP lp, PJ *toplevel, const PJ *internal, double h,
16+ struct FACTORS *fac) {
1717 double cosphi, t, n, r;
1818 int err;
1919 PJ_COORD coo = {{0 , 0 , 0 , 0 }};
@@ -55,9 +55,17 @@ int pj_factors(PJ_LP lp, PJ *toplevel, const PJ *internal, double h,
5555 lp.phi = lp.phi < 0 . ? -(M_HALFPI - h) : (M_HALFPI - h);
5656
5757 /* Longitudinal distance from central meridian */
58- lp.lam -= internal->lam0 ;
59- if (!internal->over )
60- lp.lam = adjlon (lp.lam );
58+ /* Only apply that correction for a non-pipeline transformation, otherwise
59+ * it will be applied twice (here and by the fwd_prepare() method of the
60+ * first step)
61+ */
62+ const bool bIsPipeline =
63+ internal->short_name && strcmp (internal->short_name , " pipeline" 0 ;
64+ if (!bIsPipeline) {
65+ lp.lam -= internal->lam0 ;
66+ if (!internal->over )
67+ lp.lam = adjlon (lp.lam );
68+ }
6169
6270 /* Derivatives */
6371 if (pj_deriv (lp, h, internal, &(fac->der ))) {
@@ -66,6 +74,19 @@ int pj_factors(PJ_LP lp, PJ *toplevel, const PJ *internal, double h,
6674 return 1 ;
6775 }
6876
77+ /* This code path is intended to be taken for a derived projected CRS,
78+ * with first step being a map projection, and second step some other
79+ * transformation. We must undo the scaling from the unit sphere to the
80+ * actual semi major radius done in the fwd_finalize() method of the
81+ * first step, otherwise most factors will be wrong.
82+ */
83+ if (bIsPipeline) {
84+ fac->der .x_l *= internal->ra * internal->fr_meter ;
85+ fac->der .x_p *= internal->ra * internal->fr_meter ;
86+ fac->der .y_l *= internal->ra * internal->fr_meter ;
87+ fac->der .y_p *= internal->ra * internal->fr_meter ;
88+ }
89+
6990 /* Scale factors */
7091 cosphi = cos (lp.phi );
7192 fac->h = hypot (fac->der .x_p , fac->der .y_p );
@@ -137,21 +158,9 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {
137158 type = proj_get_type (pj);
138159 }
139160
140- if (type == PJ_TYPE_PROJECTED_CRS) {
141- // If it is a projected CRS, then compute the factors on the
142- // conversion associated to it. We need to start from a temporary
143- // geographic CRS using the same datum as the one of the projected
144- // CRS, and with input coordinates being in longitude, latitude
145- // order in radian, to be consistent with the expectations of the lp
146- // input parameter. We also need to create a modified projected CRS
147- // with a normalized easting/northing axis order in metre, so the
148- // resulting operation is just a single step pipeline with no
149- // axisswap or unitconvert steps.
150-
151- auto ctx = pj->ctx ;
152- auto geodetic_crs = proj_get_source_crs (ctx, pj);
153- assert (geodetic_crs);
154- auto pm = proj_get_prime_meridian (ctx, geodetic_crs);
161+ const auto createGeogCRSLonLatRad = [](PJ *base_geodetic_crs) {
162+ auto ctx = base_geodetic_crs->ctx ;
163+ auto pm = proj_get_prime_meridian (ctx, base_geodetic_crs);
155164 double pm_longitude = 0 ;
156165 proj_prime_meridian_get_parameters (ctx, pm, &pm_longitude, nullptr ,
157166 nullptr );
@@ -160,7 +169,7 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {
160169 auto cs = proj_create_ellipsoidal_2D_cs (
161170 ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, " Radian" 1.0 );
162171 if (pm_longitude != 0 ) {
163- auto ellipsoid = proj_get_ellipsoid (ctx, geodetic_crs );
172+ auto ellipsoid = proj_get_ellipsoid (ctx, base_geodetic_crs );
164173 double semi_major_metre = 0 ;
165174 double inv_flattening = 0 ;
166175 proj_ellipsoid_get_parameters (ctx, ellipsoid, &semi_major_metre,
@@ -172,15 +181,34 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {
172181 " reference prime meridian" 0 , nullptr , 0 , cs);
173182 proj_destroy (ellipsoid);
174183 } else {
175- auto datum = proj_crs_get_datum (ctx, geodetic_crs );
184+ auto datum = proj_crs_get_datum (ctx, base_geodetic_crs );
176185 auto datum_ensemble =
177- proj_crs_get_datum_ensemble (ctx, geodetic_crs );
186+ proj_crs_get_datum_ensemble (ctx, base_geodetic_crs );
178187 geogCRSNormalized = proj_create_geographic_crs_from_datum (
179188 ctx, " unnamed crs"
180189 proj_destroy (datum);
181190 proj_destroy (datum_ensemble);
182191 }
183192 proj_destroy (cs);
193+ return geogCRSNormalized;
194+ };
195+
196+ double semiMajorToDivide = 1 ;
197+ if (type == PJ_TYPE_PROJECTED_CRS) {
198+ // If it is a projected CRS, then compute the factors on the
199+ // conversion associated to it. We need to start from a temporary
200+ // geographic CRS using the same datum as the one of the projected
201+ // CRS, and with input coordinates being in longitude, latitude
202+ // order in radian, to be consistent with the expectations of the lp
203+ // input parameter. We also need to create a modified projected CRS
204+ // with a normalized easting/northing axis order in metre, so the
205+ // resulting operation is just a single step pipeline with no
206+ // axisswap or unitconvert steps.
207+
208+ auto ctx = pj->ctx ;
209+ auto geodetic_crs = proj_get_source_crs (ctx, pj);
210+ assert (geodetic_crs);
211+ auto geogCRSNormalized = createGeogCRSLonLatRad (geodetic_crs);
184212 auto conversion = proj_crs_get_coordoperation (ctx, pj);
185213 auto projCS = proj_create_cartesian_2D_cs (
186214 ctx, PJ_CART2D_EASTING_NORTHING, " metre" 1.0 );
@@ -201,6 +229,56 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {
201229
202230 P->cached_op_for_proj_factors = newOp;
203231 pj = newOp;
232+
233+ } else if (type == PJ_TYPE_DERIVED_PROJECTED_CRS) {
234+
235+ auto ctx = pj->ctx ;
236+ auto base_proj_crs = proj_get_source_crs (ctx, pj);
237+ assert (base_proj_crs);
238+ auto geodetic_crs = proj_get_source_crs (ctx, base_proj_crs);
239+ proj_destroy (base_proj_crs);
240+ assert (geodetic_crs);
241+ auto ellipsoid = proj_get_ellipsoid (ctx, geodetic_crs);
242+ proj_ellipsoid_get_parameters (ctx, ellipsoid, &semiMajorToDivide,
243+ nullptr , nullptr , nullptr );
244+ proj_destroy (ellipsoid);
245+ auto geogCRSNormalized = createGeogCRSLonLatRad (geodetic_crs);
246+ auto newOp = proj_create_crs_to_crs_from_pj (ctx, geogCRSNormalized,
247+ pj, nullptr , nullptr );
248+ proj_destroy (geodetic_crs);
249+ proj_destroy (geogCRSNormalized);
250+ assert (newOp);
251+
252+ PJ *cs = proj_crs_get_coordinate_system (ctx, pj);
253+ double cs_unit_conv_factor = 1.0 ;
254+ proj_cs_get_axis_info (ctx, cs, 0 , nullptr , nullptr , nullptr ,
255+ &cs_unit_conv_factor, nullptr , nullptr ,
256+ nullptr );
257+ proj_destroy (cs);
258+
259+ // Remove trailing unit conversion and axis swap
260+ std::string osProjStr =
261+ proj_as_proj_string (ctx, newOp, PJ_PROJ_5, nullptr );
262+ bool resized = false ;
263+ for (const char *suffix :
264+ {" +step +proj=unitconvert" " +step +proj=axiswap"
265+ auto nPos = osProjStr.rfind (suffix);
266+ if (nPos != std::string::npos) {
267+ resized = true ;
268+ osProjStr.resize (nPos);
269+ }
270+ }
271+ if (resized) {
272+ proj_destroy (newOp);
273+ newOp = proj_create (ctx, osProjStr.c_str ());
274+ assert (newOp);
275+ }
276+
277+ newOp->fr_meter = cs_unit_conv_factor;
278+
279+ P->cached_op_for_proj_factors = newOp;
280+ pj = newOp;
281+
204282 } else if (type != PJ_TYPE_CONVERSION &&
205283 type != PJ_TYPE_TRANSFORMATION &&
206284 type != PJ_TYPE_CONCATENATED_OPERATION &&
0 commit comments