[GRASS5] problems with [rsv].proj

[Glynn Clements]

Glynn Clements wrote:

> > > When computing the region of interest, cells which fail to project
> > > should just be ignored. When performing the actual projection, cells
> > > which fail to project should be NULL.
> > 
> > Should be fairly easy to fix: Find the places in the code where 
> > [src].proj exits with 'Error in do_proj' and change action according to 
> > above.
> > 
> > (I don't have a GRASS installation to work with right now, but maybe 
> > someone else is willing to take on this task)
> 
> OK.

OK, done.

> > The way r.proj works (omitting here the initial testing/trimming 
> > done in bordwalk.c) is:
> > 
> > For each cell in the destination region
> > 	-project the cell center into source map 
> > 	-find the category value of the nearest cell 
> > 		(alternatively using cubic etc interpolation)
> > 	-assign that value to the cell in the destination region
> > 
> > For a pseudo-cylindrical destination there will a number of cells outside 
> > the ellipse that should become NULL. (The real 'borders' in this case is 
> > actually the ellipse, not the window borders given in DEFAULT_WIND). 
> > Projecting a co-ordinate outside the ellipse should result in proj 
> > returning * (inf?) for either lon or lat or both.
> 
> That clearly isn't the case. If it were, the cells outside of the
> ellipse would be garbage, but they are actually "repeats" of the data
> inside the ellipse (albeit more distorted), which seems to suggest
> periodicity.
> 
> In this test case (Wagner I), the repeat happens in both latitude and
> longitude, but the latitudinal repeat is a mirror image in both
> directions (i.e. 180° rotation).
> 
> Anyway, I'll look into it further.

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.

Also, I had to disable the remaining call to bordwalk() in r.proj,
otherwise it cropped the result. Realistically, I think that we need a
flag to allow the bordwalk() call(s) to be disabled for the cases
where they don't work correctly.

In the longer term, we should re-think the architecture of the PROJ
library. There are many practical applications where it would be
useful to have more information about the overall "nature" of a
projection than just a pair of functions which project individual
points.

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