[GRASS5] datum transforms...

Paul Kelly paul-grass at stjohnspoint.co.uk
Mon Jan 13 07:16:28 EST 2003


On Sun, 12 Jan 2003, Eric G. Miller wrote:

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

Yes it was a good help, thank you.

> 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;
> {
...
>
> For the future, a function like the following would be better.
>
...
>
> int do_transform (int count, double *x, double *y, double *h,
> 		struct pj_info *info_in, struct pj_info *info_out);
> {

I suggest calling the first function pj_to_proj (i.e. a direct replacement
for the current pj_do_proj) as it takes the same parameters, the only
difference being that it internally uses pj_transform to take advantage of
datum conversion. And the new function could be called pj_do_transform and
put into do_proj.c as well for future use and as a starting point for
people working on this in the future.

...
>    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);
>       }

Should these calls to pj_transform not have 1 for the parameter after
count, so it will step through the array one double at a time? With 0 I
think it means it won't increment at all after reading each value
(acceptable behaviour in the first function as only one point is being
transformed at a time, but not here).

...
> 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);

There is a leading slash in GRIDDIR anyway so maybe it is not needed
between the two %s , but I could be wrong as I don't know what format
G_gisbase() returns.

>     if ((ok = putenv (buf)) == -1)
>         G_warning ("%s:putenv(%s) failed!", "set_proj_lib_env", buf);
>     free (buf);
>     return ok;
> }

Yes, the change from setenv() to putenv() made this compile OK for me. But
I still don't really like this way of setting the directory proj should
look in for the conversion tables. Is setting environment variables not a
big overhead? And should it be done inside the pj_do_proj() function every
time it is called or just done once in the user programs [rsv].proj? I
would prefer not to have to change them at all as all the datum
transformation can be done seamlessly within the proj library without
changing the function called from the user programs.

An alternative idea I had came from some hints in a note by Frank
Warmerdam (http://remotesensing.org/lists/proj_archive/msg00162.html) on
use of the proj library: modify pj_open_lib.c to tell it where to
search for the files. I feel this would be acceptable as GRASS is already
using a modified version of the proj library anyway.

Looking at that file I see as well as getting the location from an
environment variable PROJ_LIB, another option is to have PROJ_LIB
#define'd somewhere in the source code. But as it varies depending on the
installation I don't think that would be possible (or is it?). I propose a
small change something like the following:
--- ../../../../../grass/src/libes/proj/pj_open_lib.c	2002-07-18 11:28:44.000000000 +0100
+++ ./pj_open_lib.c	2003-01-13 12:05:19.520147099 +0000
@@ -8,6 +8,7 @@
 #include <string.h>
 #include <errno.h>

+#define GRIDDIR "/etc/nad"
 static const char *(*pj_finder)(const char *) = NULL;
 static char * proj_lib_name =
 #ifdef PROJ_LIB
@@ -35,8 +36,24 @@
 	char fname[MAX_PATH_FILENAME+1], *sysname;
 	FILE *fid;
 	int n = 0;
-
-	/* check if ~/name */
+        const char *gisbase = G_gisbase();
+        char *buf = NULL;
+        size_t len;
+
+        if( !proj_lib_name )
+        {
+                len = strlen(gisbase) + sizeof(GRIDDIR);
+                buf = G_malloc(len);
+                sprintf( buf, "%s%s", gisbase, GRIDDIR );
+	}
+        else
+        {
+                len = strlen(proj_lib_name);
+                buf = G_malloc(len);
+                sprintf( buf, proj_lib_name);
+	}
+
+        /* check if ~/name */
 	if (*name == '~' && name[1] == DIR_CHAR)
 		if (sysname = getenv("HOME")) {
 			(void)strcpy(fname, sysname);
@@ -57,7 +74,7 @@
             sysname = pj_finder( name );

 	/* or is environment PROJ_LIB defined */
-	else if ((sysname = getenv("PROJ_LIB")) || (sysname = proj_lib_name)) {
+	else if ((sysname = getenv("PROJ_LIB")) || (sysname = buf)) {
 		(void)strcpy(fname, sysname);
 		fname[n = strlen(fname)] = DIR_CHAR;
 		fname[++n] = '\0';

(proj_lib_name is set in the preamble as follows:
static char * proj_lib_name =
#ifdef PROJ_LIB
PROJ_LIB;
#else
0;
#endif)

Looking at the function pj_load_nadgrids in pj_apply_gridshift.c though,
I'm still not too sure that we couldn't just include the full pathnames in
the nadgrids line in pj_datums.c . Does anybody who knows about the NAD
gridshifting care to test it a bit?

Paul




More information about the grass-dev mailing list