[Proj] cs2cs the fast way

[Noel Zinn]

Gerald, I agree.  You were right to abandon geodesy in favor of cartography.
-Noel

-----Original Message-----
From: proj-bounces at lists.maptools.org
[mailto:proj-bounces at lists.maptools.org] On Behalf Of Gerald I. Evenden
Sent: Friday, May 29, 2009 9:54 PM
To: PROJ.4 and general Projections Discussions
Subject: [Proj] cs2cs the fast way

I resurrected the Chebyshev code, cleaned out some error and built a quickie

prototype to show the fast way to do datum conversions for large data sets.

I did not do a true datum conversion but merely converting x-y from one 
ellipsoid to another---more on this later.  The method is outlined in the 
scaffold program below:

#include "projects.h"
#include <stdio.h>
#include <stdlib.h>
/* a structure to transport the proj constants through to the
 * func procedure 
 */
typedef struct {
	PROJ *P1, *P2;
} PP;
static PP Ps;
/* procedure called by the Chebychev routine in generating
 * coefficients.
 */
	PROJ_ANY
func(PROJ_ANY arg, PP *P) {
	PROJ_ANY v;
/* NOTE: here we use the "old" calls denigrated by the infamous email some 
time ago. Of course using the "modern call" suggested would have been a 
coding nightmare */
	v.xy = proj_fwd(proj_inv(arg.xy, P->P1), P->P2);  <<---- NOTE this
line.
	v.xy.x -= arg.xy.x; // return correction
	v.xy.y -= arg.xy.y;
	return v;
}
/* this binds it all together */
int main(int argc, char **argv) {
	PROJ_UV a={-300000.,4000000.}, b={300000.,4500000.}, resid;
	PROJ_Tseries *t;
	double r=0.001;

/* open up a projection constants for the Clark '66 ellipsoid and
 * the WGS84 ellipsoid */
	if (((Ps.P1 = proj_initstr("proj=tmerc ellps=clrk66 lon_0=0")) ==
NULL) ||
		((Ps.P2 = proj_initstr("proj=tmerc ellps=wgs84 lon_0=0")) ==
NULL)) {
		fprintf(stderr, "failed to open projections\n");
		exit(1);
	}
/*  call there coefficient generator with the data range (a, b in meters)
 *  r= desired precision, residual precision, function from above and
 *  projection constants, maximum coefficient array sizes, and "0" for
 *  Chebychev coefficients (rather than polynomial) */
	if ((t = proj_mk_cheby(a, b, r, &resid, func, &Ps, 14, 14, 0)) ==
NULL) {
		fprintf(stderr, "failed to open approximation\n");
		exit(1);
	}
/* print out max error and coefficient arrays */
	printf("resid: %g %g\n", resid.u, resid.v);
	proj_prnt_tseries(t, stdout, "%.4e");
}

The results of the above are:

resid: 0.000429018 0.000316773
u: 4
1 3 -1.5140e+01 -4.2727e-01 -1.9885e-03
3 1 2.2328e-03
v: 3
0 4 8.2507e+02 1.0405e+01 -5.3504e-01 -1.6595e-03
2 2 -2.5693e-01 -2.3924e-03

Those coefficients are all that have to be evaluated in doing the ellipsoid 
recalculation.

The transform error is kept under a meter and I would think that evaluating 
the above coefficients might be a bit faster than churning through tmerc 
twice.

Again, note the line emphasized in the "func" procedure.  If we are doing a 
true datum shift with a function prototype:

PROJ_LP datum_shift(PROJ_LP, void *Pconst);

the line would look something like:

v.xy = proj_fwd(
                 datum_shift(
                          proj_inv(arg.xy, P->P1),    P->Pconst),    P->P2);

The datum shift may add a few coefficients but we do not have to evaluated a

possibly complex datum routine.

Obviously, this is not yet ready for prime time.

-- 
The whole religious complexion of the modern world is due
to the absence from Jerusalem of a lunatic asylum.
-- Havelock Ellis (1859-1939) British psychologist
_______________________________________________
Proj mailing list
Proj at lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj