[Proj] Roussilhe Projection

[Oscar van Vlijmen]

Following a lead provided by Maciek Sieczka, I stumbled upon testpoints and
formulae for the Polish quasi-stereographic Roussilhe projection (Uklad
1965).
Testpoints:
Bugzilla Bug 818 report
<http://bugzilla.remotesensing.org/show_bug.cgi?id=818>
Formulae (in Polish):
Geodeta Magazine
<http://www.atomnet.pl/~geodeta/2000/59text1.htm>
It appears that the formulae are derived from the official publication
"Wytyczne techniczne G-1.10 ".

As with every other publication with formulae, you will find errors and
vagueries in the text. I think I've found them plus the solutions.

Forward projection:
* Seems to be correct, but it is not immediately clear what xG0 is. The text
describes it in Polish. I meditated on it a bit and found that xG0 is the
xG, computed for the central point B0,L0.
* x0 is the false northing, y0 is the false easting. There are more
countries (e.g. Switzerland) with the meaning of x and y interchanged.
Inverse projection:
* error 1: yG=yG-xq2*yQ/(4*R02)+... should be: yG=yQ-xq2*yQ/(4*R02)+...
* error 2: chsa=cos(h)*sin(alfa); fi=atan(chsa/(1-sqr(chsa)))
should be: fi=atan(chsa/sqrt(1-sqr(chsa)))

Testpoint for zone 1 (Strefa I) (units in degrees, meters):
lat=52d; lon=19d;
B0=50d 37m 30s; L0=21d 5m;
x0=5467000; y0=4637000; m0=0.9998;
I computed the xG0 as: 5610467.57697895
My result:
x65 = 5621984.46761389, y65 = 4493939.53052073
The Transpol program mentioned in the bug report gave:
5621984.46761, 4493939.53053
The inverse with my x65,y65 values gave:
lat = 51.9999999991681, lon =18.9999999227301

The text says (if I'm reading it correctly) that with a complete
forward-inverse-forward projection the error stays within 0.1 mm within the
Polish territory.

Interestingly, my corrected libproj-derived rouss gave:
x=4493939.52614385, y=5621984.47469338
take  y0=5467000; x0=4637000; m0=0.9998;
lat0=50d 37m 30s; lon0=21d 5m;
lat=52d; lon=19d;
Krasowski ellipsoid.
Note the x,y interchange.
Difference in easting with the Polish formulae: 4 mm, in northing: 7 mm.

Also interesting, the double stereographic projection (yes, sterea, at least
my translation of it) gives:
x=4493939.53118226, y=5621984.47176215
Difference in easting: 0.7 mm, in northing: 4 mm

Conclusion: now it should be possible to add - as soon as the developers and
maintainers can find the time of course - the (Polish) Roussilhe projection.
If you have difficulties understanding the article: I don't read Polish too
well, but I can email you some raw JavaScript code that works.

On reply: please try to delete some parts of this relatively long email!