[GRASS5] datum transforms...

[Eric G. Miller]

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  {
            u = (*x) / RAD_TO_DEG;
            v = (*y) / RAD_TO_DEG;
            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);
            *x = u * RAD_TO_DEG;
            *y = v * RAD_TO_DEG;
          }
          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  {
	DIVIDE_LOOP(x,y,count,RAD_TO_DEG);
	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);
	MULTIPLY_LOOP(x,y,count,RAD_TO_DEG);
      }
      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
> > $GISBASE/etc/nad
> >
> > #define GRIDDIR "/etc/nad"
> >
> > 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:

#define GRIDDIR "/etc/nad"
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