[GRASS5] datum transforms...

[Eric_G_Miller]

After some prodding by Markus, I began looking at the datum transform
problem again. And, I see pj_transform() seems to already do what we
need (inverse project, datum shift, forward project) with any part
possibly being not applicable. Does anyone know of a good reason not to
use it?

After a little hacking, I have a modified s.proj that seems to work
well (not heavily tested).

Changes required:

pj_datums.c: s/ntv1 can.dat/ntv1_can.dat/
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
get_proj.c:
   pj_get_kv()
     * don't just return 1 on strncmp(proj_in,"Lat",3) == 0 ||
                            strcmp(info->proj,"ll") == 0
     * when in_proj_keys->key[i] = "proj"
            && in_proj_keys->value[i] = "ll"
      set buffa = "proj=latlong"
       (libproj uses "latlong" or "longlat" vs. GRASS's "ll")

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

That was about all I needed to do to get s.proj2 working with
pj_transform (it took me the longest to figure out the space vs.
underscore in pj_datums.c was the cause of a failure...).

Since pj_transform operates on arrays and has height values (required),
a slightly different interface from do_proj might be appropriate (to
take advantage of transforming arrays in one "pass").

Anyway, comments?

--
"...the plural of anecdote is [not?] data." - attrib. to George Stigler

[Paul_Kelly]

Hello Eric

On Fri, 29 Nov 2002, Eric G. Miller wrote:

After some prodding by Markus, I began looking at the datum transform

I've had a look at this also and seem to have it got it working too.
Only realised how important it was when I downloaded an ASTER 30m DEM.
It was in UTM / WGS84 and I needed to re-project it to Irish Grid, when
I found it was off by 30 to 50 metres everywhere. Hadn't realised how
vital datum transform functionality is to the [rsv].proj programs before
that. As well as getting datum transformation working I have also
implemented support for the Ireland 1965 datum in GRASS.

problem again. And, I see pj_transform() seems to already do what we
need (inverse project, datum shift, forward project) with any part
possibly being not applicable. Does anyone know of a good reason not to
use it?

Seems obvious to me also.

After a little hacking, I have a modified s.proj that seems to work
well (not heavily tested).

I have got [rsv].proj all working; I didn't need to make any changes to
those programs at all, only to files in the src/libes/proj directory (a
diff file is attached). Here is a summary:

src/libes/gis/datum.table, src/libes/proj/pj_datums.c:
Added a line for Ireland 1965 datum. gis/datum.table only has
3-parameter datum shift parameters, however if a 7-parameter
transformation is available it is more accurate and supported by the
proj library. I added the 3-parameter shift to this file and the
7-parameter (as supplied by Ordnance Survey of Northern Ireland and
Ordnance Survey Ireland) to the proj list in pj_datums.c.

src/libes/proj/ellipse.table:
Not sure what this file is used for but corrected missing underscores in
some of the ellipse names anyway.

src/libes/proj/pj_datums.c, pj_ellps.c:
Corrected for more missing underscores in datum and ellipse names and
filenames. In the standalone proj library versions of these files all
the descriptive names have underscores instead of spaces: it looks like
someone substituted all the underscores for spaces without checking. I
suppose the problem didn't manifest itself before as GRASS passed raw
ellipsoid parameters instead of names to proj wherever possible, but I
have changed this behaviour now (see below, get_proj.c)
Also added a definition for sphere to pj_ellps.c; this is also in the
standalone proj library distribution now.

src/libes/proj/get_proj.c:
Implemented Eric's suggestions to enable pj_transform() to be used:

   pj_get_kv()
     * don't just return 1 on strncmp(proj_in,"Lat",3) == 0 ||
                            strcmp(info->proj,"ll") == 0
     * when in_proj_keys->key[i] = "proj"
            && in_proj_keys->value[i] = "ll"
          set buffa = "proj=latlong"
       (libproj uses "latlong" or "longlat" vs. GRASS's "ll")

Also some gratuitous changes:
1. Don't pass the 'f' ellipsoid parameter from the PROJ_INFO file onto
proj as I don't know what it is and proj doesn't seem to need it

2. Change the logic for passing (a) datum names, (b) ellipsoid names and
(c) ellipsoid parameters.
(a) implies (b) implies (c), so just pass on the highest level that we
have available, i.e. if there is an ellipsoid name in the PROJ_INFO file
don't pass on the parameters, and if there is a datum name don't pass on
the ellipsoid name. This means the values in pj_datums.c and pj_ellps.c
will be used and are the definitive ones.

3. If a datum is specified don't pass on the dx dy dz transformation
parameters; proj will get these from pj_datums.c If a datum name isn't
specified but dx dy dz are, form these into a proper proj argument
'towgs84=dx,dy,dz' and pass that on.

src/libes/proj/do_proj.c:
Implemented Eric's ideas:

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:

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

I have put that in but you will see in the diff file it is currently
commented out. 'setenv' doesn't appear to be a function on my (IRIX 6.2)
system; I'm not sure what library might be needed or if there's a better
way of doing this? Perhaps the path to the conversion tables could be
added to the nadgrids line in pj_datums.c? I don't really have any
American data so I can't really test this without a lot of effort.
However the NAD27 and NAD83 datum conversions are currently still
handled via the set_datumshift() / pj_do_proj_nad() hack so this still
works.

And that's it! After recompiling [rsv].proj to link in the new library
they now all appear to be doing datum conversions properly. It seems
amazing to me that all the code to do WGS84 datum conversions was
already in GRASS and the pj_transform() function was just waiting to be
called.

Since pj_transform operates on arrays and has height values (required),

(I had to use a dummy variable for the height value.)

a slightly different interface from do_proj might be appropriate (to
take advantage of transforming arrays in one "pass").

Yes, and perhaps the height conversion could be useful too, maybe for
taking in GPS data in ellipsoidal height format and converting it to
orthometric height. But there could well be more to it than that; I
don't really understand it but it might be useful for s.proj at least.
I wouldn't mind looking at implementing one-pass array transforms as
well when I have time.

Anyway, comments?

Likewise. I feel that this is nearly ready to go into GRASS now but
needs a bit of polishing and some checking from other people. But really
I would like to remove the current NAD27 / NAD83 hack and have the
American datum transformations handled the same as all the other ones,
by pj_transform() alone. Need some help from somebody who can do that
and has some test data though.

Paul

(attachments)

proj.diff (16.5 KB)

[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@2ztr< eryyvZ .T pveR" | rot13 | reverse

[Morten_Hulden3]

On Sun, 12 Jan 2003, Paul Kelly wrote:

...

src/libes/proj/pj_datums.c, pj_ellps.c:
Corrected for more missing underscores in datum and ellipse names and
filenames. In the standalone proj library versions of these files all
the descriptive names have underscores instead of spaces: it looks like
someone substituted all the underscores for spaces without checking. I
suppose the problem didn't manifest itself before as GRASS passed raw
ellipsoid parameters instead of names to proj wherever possible, but I
have changed this behaviour now (see below, get_proj.c)

if the behaviour is changed to grass passing ellipsoid names to proj
instead of parameters we have to take care that grass' ellipsoid table
contains entries that proj can understand. last i looked there were a lot
of differences, not just in case (wgs84/WGS84), and abbreviations , but
many variants too.

and what about all ellipsoid names that people have in their
existing PROJ_INFO files? they have to be changed too to match names proj
uses.

Morten Hulden

[Paul_Kelly]

On Mon, 13 Jan 2003, Morten Hulden wrote:

if the behaviour is changed to grass passing ellipsoid names to proj
instead of parameters we have to take care that grass' ellipsoid table
contains entries that proj can understand. last i looked there were a lot
of differences, not just in case (wgs84/WGS84), and abbreviations , but
many variants too.

and what about all ellipsoid names that people have in their
existing PROJ_INFO files? they have to be changed too to match names proj
uses.

As far as I can see the two lists (datum.table / pj_datums.c and
ellipse.table / pj_ellps.c) *are* in sync, apart from the missing
underscores problem I mentioned previously. I didn't check all the numbers
meticulously but a random sample all corresponded between the two files.
The names in the GRASS version of the proj library are different from in the
standalone proj, and appear to have been changed to match GRASS's own names.
I just presumed a decision had already been made about this.

[Paul_Kelly]

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

[Ojaswa_Sharma]

hi all!
we create a dem from a sites file have x, y and spot heights using
s.surf.rest ( or s.surf.idw).
but the dem thus obtained has bumps at the points where the heights are
originally known.
is there any way to smooth the dem thus obtained or any other way to
create a dem from a sites file?

further, what is delaunay's triangulatuion, and is it useful to me since i
want a contour map from my sites file.

thanks .

long live GRASS

-ojaswa sharma
IIT Roorkee

[Frank_Warmerdam]

Paul Kelly wrote:

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:

Paul / Eric,

I would prefer to see GRASS use the pj_set_finder() function to install a
function to find PROJ support files. This finder function would form the
path to the PROJ files within the GRASS tree using methods similar to what
you use to set the PROJ_LIB now.

This is the preferred manner for applications to control the searching for
PROJ.4 support files. Let me know if you would like details on how to
accomplish this.

I would prefer to keep GRASS changes to PROJ.4 to a minimum so that upgrades
are easy. As always, I am willing to consider upstream changes to reduce the
need for GRASS to directly modify the PROJ.4 source.

Best regards,

--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, warmerdam@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent

[Alessandro_Frigeri4]

Hello,

Nice to hear datum transforms checks in GRASS.

Will it be possible to create a 'user defined' ellipse/datum set,
i.e. during startup of a new LOCATION, so that the 'user defined'
parameters (a,b||f and dx,dy,dz) could then be used in proj calculations?

In this way ellips.table and pj_ellips.c , datum.table and pj_datum.c
could be maintained im sync with 'default' values and users can define
specific ones.

I'm suggesting this just because I'm using GRASS for studying Mars and
,as here on the Earth, we deal with different ellipsoidal parameters on
different data-set. Since just the Mars case will be easily implementable
by listing the different parameters pj_datum.c,ellips.table,pj_ellips.c
and datum.table, a general 'user defined' parameter-set will make the
trick in order to make GRASS ready-to-go for planetary studies on all
(spherical/ellipsoidal!) non-Earth systems.

Alessandro Frigeri

--
Alessandro Frigeri <geoalf[at]libero[dot]it>

[Paul_Kelly]

Hello

On Mon, 13 Jan 2003, Alessandro Frigeri wrote:

Hello,

Nice to hear datum transforms checks in GRASS.

Will it be possible to create a 'user defined' ellipse/datum set,
i.e. during startup of a new LOCATION, so that the 'user defined'
parameters (a,b||f and dx,dy,dz) could then be used in proj calculations?

At present you can do this by manually creating or editing the PROJ_INFO
file in the GRASS database. It should be possible to change g.setproj to
allow a custom ellipsoid option, in which case it would prompt you for the
further parameters:

To extend it to do that you would need to know which combinations /
permutations of ellipsoid parameters were necessary and sufficient for
defining an ellipse. (I mean I'm sure it isn't very complicated; just that
I don't know ). E.g. 'a' and 'es' uniquely define an ellipse; perhaps
also 'a' and 'b' do. Do you know what 'f' is? proj doesn't seem to need it
for transformations.

As far as I understand it so far dx, dy and dz are only necessary for
transforming from one ellipsoid to another (e.g. from a local country's
ellipsoid to wgs84- I think this is their specific meaning in GRASS); they
aren't really part of the ellipsoid definition and are only needed for
datum transformations.

Paul Kelly

[Paul_Kelly]

Hello Frank

On Mon, 13 Jan 2003, Frank Warmerdam wrote:

Paul / Eric,

I would prefer to see GRASS use the pj_set_finder() function to install a
function to find PROJ support files. This finder function would form the
path to the PROJ files within the GRASS tree using methods similar to what
you use to set the PROJ_LIB now.

At last this sounds like an elegant solution! As I see it the finder
function would be very similar to Eric's set_proj_lib_env and would go in
do_proj.c

So before calling pj_transform(), the transform function ( pj_do_proj() )
would have a line something like this:

pj_set_finder( "set_proj_lib" );

and the finder function itself could look like:

#define GRIDDIR "/etc/nad"

const char * set_proj_lib (const char *name)
{
        const char *gisbase = G_gisbase();
        char *buf = NULL;
        size_t len = strlen(gisbase) + sizeof(GRIDDIR)
                     + strlen(name) + 1;
        buf = G_malloc(len);

        sprintf (buf, "%s%s/%s", gisbase, GRIDDIR, name);

        return buf;
}

How does this look?

Paul

[Eric_G_Miller]

On Mon, Jan 13, 2003 at 12:16:28PM +0000, Paul Kelly wrote:

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.

Sounds like a fine suggestion.

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

Hmm, I thought it was offsets into the arrays. You could be right....

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.

A braino on my part...

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.

Well, I was calling it at the beginning of the s.proj. I don't know
that modifying the environment is that much overhead. Anyway, I think
we wanted to get away from modifying the library code so we could just
use libproj as distributed rather than maintaining a modified copy (more
work).

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?

How you going to hardcode the name when you don't know the final install
directory? I don't like hardcoded paths.

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

[Paul_Kelly]

Hello everyone
Well I've had a bit more time to look at this now and with the
suggestions from Eric and Frank I think I've come up with something
usable. See below for my comments on what I've changed and the attached
diff file for the changes.

grass/src/libes/gis/datum.table
grass/src/libes/proj/pj_datums.c
grass/src/libes/proj/pj_ellps.c
Add support for Ireland 65 Datum, correct missing underscores in
ellipsoid and datum names. Add support for sphere in pj_ellps.c

grass/src/libes/proj/do_proj.c
Add two new functions written by Eric Miller, pj_do_proj, a direct
replacement for the existing pj_do_proj which uses the PROJ.4
pj_transform() function instead of pj_fwd() and pj_inv(), which takes
advantage of datum transformation. Second new function pj_do_transform()
is for future use and allows a height value to be tranformed as well
(from ellipsoidal to orthometric height?) and arrays of co-ordinates to
be transformed in one pass.

grass/src/libes/proj/get_proj.c

get_proj.c contains three functions for setting up co-ordinate systems:
pj_get_kv() (which uses PROJ.4 parameters read from a PROJ_INFO file),
pj_get_string() (which uses PROJ.4 parameters passed to it in a string)
and pj_zero_proj() (which initialises a pj_info struct to 'no data').
'latlong' is now treated the same as a projected co-ordinate system and
a parameter 'proj=latlong' must be passed to PROJ.4. pj_get_kv and
pj_get_string have been changed to do this.

The parameter-passing logic in pj_get_kv has been changed to treat
pj_ellps.c and pj_datums.c as the definitive source of datum and
ellipsoid parameters, i.e. ellipsoid parameters are only passed on if
there is no ellipsoid name, and ellipsoid name is only passed on if
there is no datum name. When using pj_get_string() the programmer can
override this behaviour.

To use pj_transform, there must be an instance of the PJ struct (which
holds the parameters passed to PROJ.4) within the pj_info struct (which
GRASS programs use to communicate with these wrapper functions). At its
simplest, the PJ struct will contain one parameter, 'proj=latlong', but
this still must be set up now.

A PJ struct is created when pj_init is called. The first two functions
(pj_get_kv and pj_get_string) call pj_init and create a PJ struct within
the pj_info struct; the last one (pj_zero_proj) doesn't so it should not
be used any more.

Some programs set up a latlong projection system as follows:

  pj_zero_proj(&oproj);
  sprintf(oproj.proj,"%s","ll");
This won't work anymore as pj_zero_proj() does not create a PJ struct;
some programs used the following syntax which will work:

  pj_zero_proj(&info_in);
  parms_in[0] = '\0';
  pj_get_string(&info_in, parms_in);
as pj_get_string creates calls pj_init and creates a PJ struct (passing
a null string creates a latlong projection). However the call to
pj_zero_proj() is superfluous as the 3 lines in that function are also
contained in pj_get_string. So I removed that call from these programs.
And the other programs that used the first method to set up a basic
latlong projection system have been changed to use the second method
(see details below). pj_get_string can also take NULL for the second
argument; is this the same as a string with its first character '\0'? As
I wasn't able to test most of these changes I kept the syntax exactly
the same as in other GRASS programs so that there should be no problems.

I would suggest removing pj_zero_proj() from the GRASS API now as it
doesn't seem to have any valid purpose and I removed all calls to it in
the user programs.

Also in get_proj.c, a new function added:
set_proj_lib() is a finder function for use with the PROJ.4 function
pj_set_finder(), and sets the path name where PROJ.4 should look for the
datum conversion tables. pj_set_finder is called with this function as
an argument, in pj_get_kv and pj_get_string, before pj_init is called.

(Deleted:
grass/src/libes/proj/datum_shift.c
grass/src/libes/proj/do_proj_nad.c)
These contained the set_datumshift() and do_proj_nad() functions which
were only necessary for the previous implementation of NAD/NAD datum
transformation, which used lower level PROJ.4 functions than
pj_transform().

grass/src/include/projects.h
Add prototypes for new functions and remove deleted ones
grass/src/libes/proj/Gmakefile
Remove datum_shift.o and do_proj_nad.o

Other GRASS programs which use PROJ.4 library:
1. grass/src/general/g.region/cmd/printwindow.c
Change from
  pj_zero_proj(&oproj);
  sprintf(oproj.proj,"%s","ll");
  style initialisation to
  parms_in[0] = '\0';
  pj_get_string(&info_in, parms_in);
  style.
2. grass/src/mapdev/v.in.gshhs/main.c
No problems; already used pj_get_string as above. Removed superfluous
call to pj_zero_proj().
3. grass/src/mapdev/v.proj/local_proto.h
grass/src/mapdev/v.proj/main.c
Revert program from using separate NAD/NAD datum conversion
4. grass/src/misc/m.ll2db/main.c
No problems; already used pj_get_string. Removed superfluous call to
pj_zero_proj().
5. grass/src/raster/r.proj/cmd/bordwalk.c
grass/src/raster/r.proj/cmd/main.c
Revert program from using separate NAD/NAD datum conversion
6. grass/src/raster/r.sun/main.c
Same as 1.
7. grass/src/raster/r.sunmask/g_solposition.c
Same as 1.
8. grass/src/sites/s.datum.shift/main.c
Same as 1.
9. grass/src/sites/s.proj/main.c
Revert program from using separate NAD/NAD datum conversion
10. grass/src/mapdev/v.digit/map_init_new.c
No problems; already used pj_get_string.
11. grass/src/libes/image3/convert_ll.c
12. grass/src/mapdev/v.mkquads/convert.c
13. grass/src/raster/r.in.gdal/main.c
No problems; all already used pj_get_kv(). Only v.mkquads and r.in.gdal
appear to read the PROJ_INFO file for two separate locations, and will
thus be able to take advantage of datum conversion straight away, as
will [rsv].proj.

14. grass/src/misc/m.proj/main.c
15. grass/src/misc/m.proj2/main.c
inter and cmd versions of this program will have no problems with the
new pj_do_transform() but from what I can see it looks like they won't
be able to take advantage of datum transformation either as they
manually form the projection parameters for PROJ.4 and don't pass on the
datum, even if it existed in the PROJ_INFO file. Would need to be looked
into more though.

More general notes:

PROJ.4 will perform datum conversion if a datum= parameter is passed to
it, or a towgs84= or nadgrids= parameter, for BOTH the input AND output
projections. All the datum= parameter really does is expand the
corresponding definition in the pj_datums.c file, which includes the
correct towgs84= or nadgrids= value, together with the ellipsoid name.

If the PROJ_INFO file contains a datum: entry, and the program passes
all the entries in this file onto PROJ.4 for both projections (as
[rsv].proj do and r.in.gdal and v.mkquads appear to do, as far as I can
see), then datum conversion will automatically happen. If only one
projection has a datum specified then datum conversion won't take place.

For a datum that isn't included in the pj_datums.c list, you could
manually add the necessary parameters to your PROJ_INFO file. Other
parts of GRASS would probably ignore them but I can't be sure. E.g.
  dx:
  dy:
  dz:
or
  towgs84: dx,dy,dz
should do the same thing,
or
  nadgrids: alaska,hawaii
The last example might work and would be useful if your location is one
of the NAD conversion areas not covered by the conus and ntv1_can.dat
files.

For more details on how the changes work look at the diff file.

If nobody has any objections I'd really like to have this put into the
CVS for testing. I wonder is the attached file enough for that; the CVS
server seems to have an access problem (I get a Connection refused error
upon 'cvs login') so I took the latest CVS snapshot from the website and
diffed against that. I hope it will do.

Paul

(attachments)

proj.diff (37.4 KB)

[Frank_Warmerdam]

Paul Kelly wrote:

+}
diff -ru grass50_exp_2003_01_10/src/libes/proj/pj_datums.c grass/src/libes/proj/pj_datums.c
--- grass50_exp_2003_01_10/src/libes/proj/pj_datums.c 2002-05-31 04:12:41.000000000 +0100
+++ grass/src/libes/proj/pj_datums.c 2003-01-15 12:34:16.549358000 +0000
@@ -56,7 +56,7 @@
/* id definition ellipse comments */
/* -- ---------- ------- -------- */
"ggrs87", "towgs84=-199.87,74.79,246.62", "grs80", "Greek Geodetic Reference System 1987",
-"nad27", "nadgrids=conus,ntv1 can.dat", "clark66", "North American Datum 1927",
+"nad27", "nadgrids=conus,ntv1_can.dat", "clark66", "North American Datum 1927",
"wgs84", "towgs84=0.0,0.0,0.0", "wgs84", "World Geodetic System 1984",
"wgs72", "towgs84=0.0,0.0,5.0", "wgs72", "World Geodetic System 1972",
"nad83", "towgs84=0.0,0.0,0.0", "grs80", "North American 1983",
@@ -79,5 +79,6 @@
"potsdam", "towgs84=606.0,23.0,413.0", "bessel", "Potsdam Rauenberg 1950 DHDN",
"carthage", "towgs84=-263.0,6.0,431.0", "clark80", "Carthage 1934 Tunisia",
"hermannskogel", "towgs84=653.0,-212.0,449.0", "bessel", "Hermannskogel",
+"ire65", "towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15", "modif_airy", "Ireland 1965",
NULL, NULL, NULL, NULL ,
};

Paul,

In general your changes look good to me. I am curious about a few things.
Do you have any idea how the underscores got stripped out of stuff? They are
fine in the PROJ.4 originals. I would note that even the underscores in the
comment fields for the datums are intended to be significant. The comment
fields are intended to hold the full EPSG/OGC name for the datums for ease
of translation.

I have tentatively added the potsdam, carthage, hermannskogel and ire65
datums to the pj_datums.c in PROJ.4, even though I am hesitant to extend the
datums list untill such time as I have a meaningful naming convention for
the datums. I am sticking with the upper case names for WGS84, GGRS87,
NAD83 and NAD27. Why does GRASS seem to force these to lower case? I hope
the comparison logic would be case insensitive.

For me, the key things are to try and minimize the differences in the core
of PROJ.4 between GRASS and the public PROJ.4 release. In GRASS 5.1 and
beyond, as GRASS becomes more shared library oriented I would hope that
GRASS would be able to work with pre-installed PROJ.4 libraries instead of
including a static copy in libgis. If this is ever to work we have to ensure
that the GRASS specific logic for PROJ.4 use is embedded in the GRASS
wrapper functions (such as pj_do_transform()), not in PROJ.4 code itself.

Best regards,

--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, warmerdam@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent

[Paul_Kelly]

Hello Frank

On Wed, 15 Jan 2003, Frank Warmerdam wrote:

I have tentatively added the potsdam, carthage, hermannskogel and ire65
datums to the pj_datums.c in PROJ.4, even though I am hesitant to extend the

There are lots more datums in the GRASS CVS version of pj_datums.c than
you can see in the diff file, if you are interested.

datums list untill such time as I have a meaningful naming convention for
the datums. I am sticking with the upper case names for WGS84, GGRS87,
NAD83 and NAD27. Why does GRASS seem to force these to lower case? I hope
the comparison logic would be case insensitive.

Some of the names are different (not just in case) from the PROJ.4
originals as well, especially in the ellipse table. But having these
already consistent with the internal GRASS names made the changes I made
much easier to do. I've no idea why it was changed like that but it was
definitely the simplest way to get it working.

A slightly related point: I had a problem with r.in.gdal automatically
creating a location for me (from an Aster DEM) using the ellipse name
'WGS84'. I couldn't reproject it and the problem was the capital letters
but it was a long time after I'd re-created the projection information
manually using g.setproj and got r.proj working that I realised what the
problem was.

For me, the key things are to try and minimize the differences in the core
of PROJ.4 between GRASS and the public PROJ.4 release. In GRASS 5.1 and
beyond, as GRASS becomes more shared library oriented I would hope that

I don't like relying on shared libraries. It is really tedious to get them
all installed and configured and can take hours or a day at least for each
one. It is not very convenient for users of minority operating systems on
which the software has more often than not not been tested before release. I
have not even been able to get the standard PROJ.4 release to compile on
my IRIX system (libtool problem I think) but it is OK as I don't need it
on its own.

If there is one version compiled as part of GRASS that has all the
features GRASS programs need, then it can stay fixed like that for ever
and there will be no problems with new releases of the library introducing
bugs unrelated to GRASS (I'm not talking about PROJ.4, just generally).

GRASS would be able to work with pre-installed PROJ.4 libraries instead of
including a static copy in libgis. If this is ever to work we have to ensure
that the GRASS specific logic for PROJ.4 use is embedded in the GRASS
wrapper functions (such as pj_do_transform()), not in PROJ.4 code itself.

That sounds like a good idea, what we are already doing for ll -> latlong.
We can't change existing PROJ_INFO files in GRASS installations all over
the world as Morten Hulden already said, so the GRASS internal names have
to stay consistent.

Paul

[Paul_Kelly]

Hello everyone
I've committed to CVS the changes described in my e-mail of 15th January
(See http://grass.itc.it/pipermail/grass5/2003-January/004562.html ) and
would appreciate testing and bug reports on the new datum transformation
capabilities in GRASS.

I intend to further improve on this, not changing the functionality, but
implementing it more cleanly. This will avoid the need for the GRASS and
PROJ.4 datum lists to be in sync (as at present) and should enable GRASS
to work with a shared PROJ.4 library by only passing the necessary
underlying datum and ellipsoid parameters to it. The GRASS datum table
format and functions need to be improved to do this and I will post my
suggestions and ask for comments in due course.

Hoping the improvements are useful for someone,

Paul

[Eric_G_Miller]

On Sun, Jan 26, 2003 at 06:02:13PM +0000, Paul Kelly wrote:

Hello everyone
I've committed to CVS the changes described in my e-mail of 15th January
(See http://grass.itc.it/pipermail/grass5/2003-January/004562.html ) and
would appreciate testing and bug reports on the new datum transformation
capabilities in GRASS.

I intend to further improve on this, not changing the functionality, but
implementing it more cleanly. This will avoid the need for the GRASS and
PROJ.4 datum lists to be in sync (as at present) and should enable GRASS
to work with a shared PROJ.4 library by only passing the necessary
underlying datum and ellipsoid parameters to it. The GRASS datum table
format and functions need to be improved to do this and I will post my
suggestions and ask for comments in due course.

Hoping the improvements are useful for someone,

Groovy. I'm sure all GRASS users will benefit. This has been a long
sought after functionality...

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

[neteler]

On Sun, Jan 26, 2003 at 01:44:43PM -0800, Eric G. Miller wrote:

On Sun, Jan 26, 2003 at 06:02:13PM +0000, Paul Kelly wrote:
> Hello everyone
> I've committed to CVS the changes described in my e-mail of 15th January
> (See http://grass.itc.it/pipermail/grass5/2003-January/004562.html ) and
> would appreciate testing and bug reports on the new datum transformation
> capabilities in GRASS.
>
> I intend to further improve on this, not changing the functionality, but
> implementing it more cleanly. This will avoid the need for the GRASS and
> PROJ.4 datum lists to be in sync (as at present) and should enable GRASS
> to work with a shared PROJ.4 library by only passing the necessary
> underlying datum and ellipsoid parameters to it. The GRASS datum table
> format and functions need to be improved to do this and I will post my
> suggestions and ask for comments in due course.
>
> Hoping the improvements are useful for someone,

Groovy. I'm sure all GRASS users will benefit. This has been a long
sought after functionality...

Great, Paul! This was really long awaited (some years only

Thanks,

Markus

[neteler]

Hello Paul,

(cc to 'grass5' as it is of general interest)

On Mon, Jan 27, 2003 at 10:40:34AM +0000, Paul Kelly wrote:

Hello Markus
>From what I can see m.proj does the same thing as the proj or cs2cs
programs included with the PROJ.4 library. Does it pre-date its general
availability or is there something else it does that I have missed?

The m.proj is there for a long time, maybe as a sort of convenient
tool for 'proj'.
The problem is that results calculated with m.proj/m.proj2 (this is the cmd
version) are shifted by around 100m.

Do you have any examples of what people use it for?

Yes: E.G. they have a paper map in a different projection (for
example in Germany UTM and Gauss-Krueger are used in parallel).
to check points etc some users use m.proj (maybe because they
don't know that nowadays the nice 'cs2cs' exists).

Then people may also use it to reproject point lists from GPS devices etc.
before importing them with s.proj.

It may be simpler to re-write it using most of the code from cs2cs, but
that seems a bit pointless.

Sounds right. But then we should think about removing m.proj[2]
from CVS to avoid that some people continue to use it and receive
shifted point data.

If m.proj doesn't actually operate on files in the GRASS database but uses
external text files, should it be there at all.....

This might be discussed on 'grass5' if we want to support that
in general in future. The problem is that many users take results for
truth when working with GRASS tools. So I feel that it is dangerous
when the results from m.proj[2] and cs2cs differ.

An example: UTM Zone 32-> Gauss-Krueger (Transverse Mercator)
# proj/nad/epsg
# DHDN / Germany zone 3
# datum from:
# http://www.colorado.Edu/geography/gcraft/notes/datum/dlist.html
cs2cs +init=epsg:32632 +to +init=epsg:31493 +towgs84=606,23,413
486000 5912000 0
3486066.49 5913928.42 6.189

(this result is rather identical to other proprietary software, it differs
by some tens of centimeters)

m.proj delivers:

Initializing INPUT PROJECTION:
Please specify projection name

utm

Initializing INPUT projection ELLIPSOID:

wgs84

    Enter INPUT Projection Zone [zone] : 32
INPUT Projection: Would you like to use South Hemisphere (y/[n]) ? n
Enter plural for INPUT units [meters]:
Parameters used in INPUT projection:
         +proj=utm +a=6378137.0000000000 +es=0.0066943800 +zone=32 +unfact=1.0000000000

Initializing OUTPUT PROJECTION:
Please specify projection name

tmerc

Initializing OUTPUT projection ELLIPSOID:
Please specify ellipsoid name

bessel

  Enter OUTPUT Central Parallel [lat_0] (23N) : 0N
  Enter OUTPUT Central Meridian [lon_0] (96W) : 9E
    Enter OUTPUT Scale Factor at the Central Meridian [k_0] (1.000000) :
    Enter OUTPUT False Easting [x_0] (0.000000) : 3500000
    Enter OUTPUT False Northing [y_0] (0.000000) :
Enter plural for OUTPUT units [meters]:
Parameters used in OUTPUT projection:
         +proj=tmerc +a=6377397.1550000003 +es=0.0066743722
+lat_0=0.0000000000 +lon_0=9.0000000000 +k=1.0000000000
+x_0=3500000.0000000000 +y_0=0.0000000000 +unfact=1.0000000000

    Universe Transverse Mercator -> Transverse Mercator Conversion:
        X_in (meters) Y_in (meters) X_out (meters) Y_out (meters)
        ------------ ----------- ------------- -------------
486000.00 5912000.00 3485996.11 5913755.52

-----------
Differences cs2cs/m.proj:
3486066.49-3485996.11=70.38
5913928.42-5913755.52=172.9

So the results of m.proj/m.proj2 are shifted as expected.

As conclusion:
Is it right that either m.proj/m.proj2 should ask for the
datum (and pass to proj) or both commands be removed?

Markus

[Helena_Mitasova]

I would vote for m.proj to be removed because it does not work with
GRASS data,
but it may be useful to put somewhere information about cs2cs (which I
did not know
about). Maybe keep m.proj in the manual for some time with a note that
it was
retired and people should use cs2cs (where do you get it?)
The best place of course for this (and other essential non-GRASS tools)
would be in an on-line GRASS tutorial if it ever happens.
But this is just my opinion (I have never used m.proj except for the
book)

Helena

Markus Neteler wrote:

Hello Paul,

(cc to 'grass5' as it is of general interest)

On Mon, Jan 27, 2003 at 10:40:34AM +0000, Paul Kelly wrote:
> Hello Markus
> >From what I can see m.proj does the same thing as the proj or cs2cs
> programs included with the PROJ.4 library. Does it pre-date its general
> availability or is there something else it does that I have missed?

The m.proj is there for a long time, maybe as a sort of convenient
tool for 'proj'.
The problem is that results calculated with m.proj/m.proj2 (this is the cmd
version) are shifted by around 100m.

> Do you have any examples of what people use it for?

Yes: E.G. they have a paper map in a different projection (for
example in Germany UTM and Gauss-Krueger are used in parallel).
to check points etc some users use m.proj (maybe because they
don't know that nowadays the nice 'cs2cs' exists).

Then people may also use it to reproject point lists from GPS devices etc.
before importing them with s.proj.

> It may be simpler to re-write it using most of the code from cs2cs, but
> that seems a bit pointless.

Sounds right. But then we should think about removing m.proj[2]
from CVS to avoid that some people continue to use it and receive
shifted point data.

> If m.proj doesn't actually operate on files in the GRASS database but uses
> external text files, should it be there at all.....

This might be discussed on 'grass5' if we want to support that
in general in future. The problem is that many users take results for
truth when working with GRASS tools. So I feel that it is dangerous
when the results from m.proj[2] and cs2cs differ.

An example: UTM Zone 32-> Gauss-Krueger (Transverse Mercator)
# proj/nad/epsg
# DHDN / Germany zone 3
# datum from:
# http://www.colorado.Edu/geography/gcraft/notes/datum/dlist.html
cs2cs +init=epsg:32632 +to +init=epsg:31493 +towgs84=606,23,413
486000 5912000 0
3486066.49 5913928.42 6.189

(this result is rather identical to other proprietary software, it differs
by some tens of centimeters)

m.proj delivers:

Initializing INPUT PROJECTION:
Please specify projection name
>utm
Initializing INPUT projection ELLIPSOID:
>wgs84
    Enter INPUT Projection Zone [zone] : 32
INPUT Projection: Would you like to use South Hemisphere (y/[n]) ? n
Enter plural for INPUT units [meters]:
Parameters used in INPUT projection:
         +proj=utm +a=6378137.0000000000 +es=0.0066943800 +zone=32 +unfact=1.0000000000

Initializing OUTPUT PROJECTION:
Please specify projection name
>tmerc
Initializing OUTPUT projection ELLIPSOID:
Please specify ellipsoid name
>bessel
  Enter OUTPUT Central Parallel [lat_0] (23N) : 0N
  Enter OUTPUT Central Meridian [lon_0] (96W) : 9E
    Enter OUTPUT Scale Factor at the Central Meridian [k_0] (1.000000) :
    Enter OUTPUT False Easting [x_0] (0.000000) : 3500000
    Enter OUTPUT False Northing [y_0] (0.000000) :
Enter plural for OUTPUT units [meters]:
Parameters used in OUTPUT projection:
         +proj=tmerc +a=6377397.1550000003 +es=0.0066743722
+lat_0=0.0000000000 +lon_0=9.0000000000 +k=1.0000000000
+x_0=3500000.0000000000 +y_0=0.0000000000 +unfact=1.0000000000

    Universe Transverse Mercator -> Transverse Mercator Conversion:
        X_in (meters) Y_in (meters) X_out (meters) Y_out (meters)
        ------------ ----------- ------------- -------------
486000.00 5912000.00 3485996.11 5913755.52

-----------
Differences cs2cs/m.proj:
3486066.49-3485996.11=70.38
5913928.42-5913755.52=172.9

So the results of m.proj/m.proj2 are shifted as expected.

As conclusion:
Is it right that either m.proj/m.proj2 should ask for the
datum (and pass to proj) or both commands be removed?

Markus

_______________________________________________
grass5 mailing list
grass5@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass5