# [GRASS5] datum transforms...

Eric G. Miller egm2 at jps.net
Sun Jan 12 20:55:40 EST 2003

```On Sun, Jan 12, 2003 at 09:03:31PM +0000, Paul Kelly wrote:

Thanks for working on this.  Perhaps the following will help.

> > do_proj.c:
> >    * LL<->LL is no longer a no-op.
> >    * change all variations of pj_fwd and pj_inv to one call to
> >      pj_transform
>
> Had a bit of a problem with this part:

Here's the function I replaced do_proj with.  At the moment, I don't
remember why x & y aren't converted to radians for the ll<->ll case.
However, that case needs to be handle properly for strictly doing
datum translations.

static int do_transform(x,y,info_in,info_out)
double *x,*y;
struct pj_info *info_in,*info_out;
{
int ok;
double u,v;
double h = 0.0;

METERS_in = info_in->meters;
METERS_out = info_out->meters;

if (strncmp(info_in->proj,"ll",2) == 0 ) {
if (strncmp(info_out->proj,"ll",2) == 0 )
ok = pj_transform (info_in->pj, info_out->pj, 1, 0, x, y, &h);
else  {
ok = pj_transform(info_in->pj,info_out->pj, 1, 0, &u, &v, &h);
*x = u / METERS_out;
*y = v / METERS_out;
}
}
else {
if (strncmp(info_out->proj,"ll",2) == 0 ) {
u = *x * METERS_in;
v = *y * METERS_in;
ok = pj_transform(info_in->pj,info_out->pj, 1, 0, &u, &v, &h);
}
else {
u = *x * METERS_in;
v = *y * METERS_in;
ok = pj_transform(info_in->pj,info_out->pj, 1, 0, &u, &v, &h);
*x = u / METERS_out;
*y = v / METERS_out;
}
}
return ok;
}

For the future, a function like the following would be better.

/* a couple defines to simplify reading the function */
#define MULTIPLY_LOOP(x,y,c,m) \
do {\
int i; \
for (i = 0; i < c; ++i) {\
x[i] *= m; \
y[i] *= m; \
}\
} while (0)

#define DIVIDE_LOOP(x,y,c,m) \
do {\
int i; \
for (i = 0; i < c; ++i) {\
x[i] /= m; \
y[i] /= m; \
}\
} while (0)

int do_transform (int count, double *x, double *y, double *h,
struct pj_info *info_in, struct pj_info *info_out);
{
int ok;
int has_h = 1;

METERS_in = info_in->meters;
METERS_out = info_out->meters;

if (h == NULL) {
int i;
h = G_malloc (sizeof *h * count);
/* they say memset is only guaranteed for chars ;-( */
for (i = 0; i < count; ++i)
h[i] = 0.0;
has_h = 0;
}
if (strncmp(info_in->proj,"ll",2) == 0 ) {
if (strncmp(info_out->proj,"ll",2) == 0 )
ok = pj_transform (info_in->pj, info_out->pj, count, 0, x, y, h);
else  {
ok = pj_transform(info_in->pj,info_out->pj, count, 0, x, y, h);
DIVIDE_LOOP(x,y,count,METERS_out);
}
}
else {
if (strncmp(info_out->proj,"ll",2) == 0 ) {
MULTIPLY_LOOP(x,y,count,METERS_in);
ok = pj_transform(info_in->pj,info_out->pj, count, 0, x, y, h);
}
else {
MULTIPLY_LOOP(x,y,count,METERS_in);
ok = pj_transform(info_in->pj,info_out->pj, count, 0, x, y, h);
DIVIDE_LOOP(x,y,count,METERS_out);
}
}
if (!has_h)  free (h);

return ok;
}

> > Before using pj_transform(), set the environment variable PROJ_LIB to
> >
> >
> > static void set_proj_lib_env (void)
> > {
> >     const char *gisbase = G_gisbase();
> >     char *value = NULL;
> >     size_t len = strlen(gisbase) * sizeof(GRIDDIR) + 1;
> >     value = G_malloc (len);
> >     sprintf (value, "%s%s", gisbase, GRIDDIR);
> >     if (setenv ("PROJ_LIB", value, 1) == 0)
> >         fprintf(stderr, "set PROJ_LIB='%s'\n", value);
> >     else
> >         fprintf(stderr, "failed to set PROJ_LIB='%s'\n", value);
> >
> >     free (value);
> > }

Hmm, maybe s/setenv/putenv/? The function will have to be rewritten
to something like:

int set_proj_lib_env (void)
{
int ok;
const char *gisbase = G_gisbase();
char *buf = NULL;
size_t len = strlen(gisbase) + sizeof(GRIDDIR)
+ sizeof("PROJ_LIB=") + 1;
buf = G_malloc(len);
sprintf (buf, "PROJ_LIB=%s/%s", gisbase, GRIDDIR);
if ((ok = putenv (buf)) == -1)
G_warning ("%s:putenv(%s) failed!", "set_proj_lib_env", buf);
free (buf);
return ok;
}

--
echo ">gra.fcw at 2ztr< eryyvZ .T pveR" | rot13 | reverse

```