[GRASS5] problems with [rsv].proj

[Glynn Clements]

Morten Hulden wrote:

> > Done. The following patch produces the "desired" result (i.e. the way
> > that pseudo-cylindrical world maps normally look in illustrations),
> > but is it the right thing to do?
> > 
> > diff -u -r1.2 PJ_urmfps.c
> > --- src/libes/proj/PJ_urmfps.c	20 Apr 2002 19:13:44 -0000	1.2
> > +++ src/libes/proj/PJ_urmfps.c	30 May 2002 04:30:08 -0000
> > @@ -17,8 +17,12 @@
> >  }
> >  INVERSE(s_inverse); /* sphere */
> >  	xy.y /= P->C_y;
> > +	if (fabs(xy.y) > HALFPI)
> > +		I_ERROR
> >  	lp.phi = aasin(sin(xy.y) / P->n);
> >  	lp.lam = xy.x / (C_x * cos(xy.y));
> > +	if (fabs(lp.lam) > PI)
> > +		I_ERROR
> >  	return (lp);
> >  }
> >  FREEUP; if (P) pj_dalloc(P); }
> > 
> > If it is, instinct suggests that the same issue will apply to many
> > (most?) of the projections, not just Wagner I.
> 
> But these are changes to proj, not to grass. You should take them up on 
> the proj mailing list and try to get them accepted there, otherwise we'll 
> soon be out of sync with the proj code.
> 
> I was thinking of a more general solution: catch the x & y return values 
> from pj_do_proj ín main.c just before r.proj goes into the 'readcell'
> routine, and just loop if x or y are not correct.
> 
> In the command line application proj returns a star symbol (*) for 
> lon or lat if you try to project a point outside the bounding ellipse in 
> a pseudo-cylindical projections. Should be possible to catch that 
> situation in [rsv].proj.

This may be true for some projections, but certainly not for all of
them. E.g. Wagner I (this just happens to be the first one I tried):

[glynn at cerise grass]$ /usr/src/proj-4.4.5/src/proj -I +proj=wag1 +a=6378137.0
20000000 20000000
77d46'10.176"E	52d37'14.123"N
[glynn at cerise grass]$ /usr/src/proj-4.4.5/src/proj +proj=wag1 +a=6378137.0
77d45'56.247"E	52d37'21.503"N
5510652.47	6371061.63

Look at the code; s_inverse doesn't care about the magnitude of the y
coordinate; it just gets passed to sin(), hence the behaviour is bound
to be periodic.

In this case, the longitudinal wrap could be trapped within GRASS, but
the latitudinal wrap can't be trapped without knowing the specifics of
the projection; and that is supposed to be PROJ's job.

-- 
Glynn Clements <glynn.clements at virgin.net>