[GRASS5] Projection programs

[Bob Covill]

Paul,

I tried the patch listed below, and it does the trick. The other two
places I know where the projection is initialized this way are in
m.proj(2), and v.in.gshhs. It may occur in other places as well.

Paul Kelly wrote:
> 
> 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?
> 

The m.ll2db program is derived from m.ll2u (and m.u2ll) and accepts
pretty much the same format. Instead of the user supplying projection
parameters, they are obtained from PROJ_* files. It accepts forward and
reverse projection to and from Lat. Long. coordinates. In other words it
should work in all GRASS Projections with the exception of XY and
LatLong (ll). Also, the program uses the PROJ library instead of the
coorcnv library.

For example a table of Lat. Long coordinates can be converted to the
users current projection with the following ....
m.ll2db in=coords_file.LL out=coords_file.UTM 
There are flags for reverse coordinate order, and inverse projection.
Data can also be piped through the program.

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


Thanks for your help.

-- 
Bob Covill

Tekmap Consulting
P.O. Box 2016
Fall River, N.S.
B2T 1K6
Canada

E-Mail: bcovill at tekmap.ns.ca
Phone: 902-860-1496
Fax: 902-860-1498