[GRASS5] Projection programs

[Paul Kelly]

Hello Bob

On Mon, 3 Feb 2003, Bob Covill wrote:

> Paul,
>
> I have latest CVS GRASS and have found a couple of projection programs
> are not quite working. I discovered the problem when using m.ll2db which
> converts to and from GRASS database and Lat. Long. coordinates. When I
> run the program I get the following error ...

BTW I couldn't find a man page for this program; could you give me an
example of a typical command line you might use with it?

> cannot initialize pj
> cause: major axis or radius = 0 or not given

Oops. It looks like we need to specify the ellipsoid for the input
projection; it doesn't just assume it is the same as the output if not
specified. I had been testing this behaviour using different combinations
of arguments given to cs2cs, and it seemed happy with a simple +proj=latlong,
but I realise now it must have been silently adding a default ellipsoid to
the input projection definition, to get round this problem.

> The problem is with how the "ll" projection is initialized, namely ...
> parms_in[0] = '\0';
> pj_get_string(&info_in, parms_in);
>
> If I change this to a better qualified initialization such as ...
>
> sprintf(ll_buf, "+proj=ll +ellsp=wgs84"); /* +datum also works with
> current setting */
> pj_get_string(&info_in, ll_buf);

In this case, and for most of the programs that use the proj library,
there is no datum transformation, so the input ellipsoid should be the
same as the output ellipsoid? Following that logic, could you please test
the following patch for me and see if it gives you your expected results.

Index: src/misc/m.ll2db/main.c
===================================================================
RCS file: /grassrepository/grass/src/misc/m.ll2db/main.c,v
retrieving revision 1.3
diff -u -r1.3 main.c
--- src/misc/m.ll2db/main.c	26 Jan 2003 17:47:01 -0000	1.3
+++ src/misc/m.ll2db/main.c	3 Feb 2003 16:21:35 -0000
@@ -148,11 +148,6 @@
 exit(1);
 }

-/* In Info */
-parms_in[0] = '\0';
-pj_get_string(&info_in, parms_in);
-
-
 /* Out Info */
 out_proj_keys = G_get_projinfo();
 out_unit_keys = G_get_projunits();
@@ -160,6 +155,12 @@
 exit (0);
 }

+/* In Info */
+if( G_find_key_value("ellps", out_proj_keys) != NULL )
+    sprintf(parms_in, "proj=ll ellps=%s", G_find_key_value("ellps", out_proj_keys) );
+else
+    sprintf(parms_in, "proj=ll ellps=wgs84");
+pj_get_string(&info_in, parms_in);

     if (isatty(0))
     {

If it works I'll add something similar to the other modules that use proj
with no datum transformation.

Paul