[Proj] cs2cs the fast way
Gerald I. Evenden
geraldi.evenden at gmail.com
Fri May 29 21:54:18 EST 2009
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
More information about the Proj
mailing list