[OSRS-PROJ] Horrors of --- Krovak projection

[Ron Russell]

Gerald I. Evenden wrote:
> 
> I modified the title to be more explicit of the topic of this
> thread.
> 
> I found and am using as a bench FORTRAN code from:
> 
> http://pecny.asu.cas.cz/publ/KROVAK/krovak.html
> 
> The author is listed as J. Kostelecky but I wonder if it is
> your code.  Attribution is not one of the strong points
> of the web. :-(  This code seems to be generating the bench
> results.  At least the math in the front of the article
> matched the code---more than I can say for other sites.
> The FORTRAN code starts out with:
> 
> C TRANSFORMATION BETWEEN ELLIPSOIDAL AND PLANAR COORDINATES
> C IN CZECH AND SLOVAK NATIONAL PROJECTION "KROVAK"
>         PROGRAM KROTEST
> C INPUT DATA must be IN FILE 'KROTEST.DAT'
> C FIRST RAW : NO OF THE VARIANT: VARIANT 1 MEANS FI,
> C LAMBDA --> X,Y
> C VARIANT 2 MEANS X,Y --> FI, LAMBDA
> C SECOND AND OTHER LINES: NP, FI(DEG,MIN,SEC), LAMBDA(DEG,MIN,SEC)
>         ...
> Is that yours?

No, I am never so confident in my code to offer it outside the
building!

> 
> The thing I do not like so far, is I would like a general
> development of the projection such that it could be transfered
> to other regions.  Right now it is stuck with the specific area
> and the eqns to define the constants, pulled out of the air,
> are not presented.

These constants are used in the Conformal projection of the ellipsoid
to The Gausian Conformal Sphere. I think I tracked this down at
one time.

When I look back at my code, I have used your formulae to calculate
what you call B and I call the constant alpha (1.0005974983). This
is the Longitude factor.

I have a hand written formulae (from the Czecks) for the other 
Latitude constant k (0.9965924869) or 1/k (1.0034191640), but I 
do not seem to have used that; the constant is still in my code. 
Perhaps I couldn't get the formulae to give the same result or 
perhaps the code was finished by the time I had tracked down the data.

The formulae is (I hope):

k = top/bottom

top = (tan(phi0/2 + 45)*
        ((1.0 - e*sin(phi0))/(1.0 - e*sin(phi0)))^^(e/2))^^alpha

bottom = tan(U0/2 + 45)

phi0 is the latitude for the conformal projection to the sphere (49 30n)
in this case, and U0 is the conformal latitude of this, so perhaps
this needs iteration. 45 is 45 degrees.

I think I have got the brackets right, but don't sue me!

> 
> I did find someone's effort to translate it to PROJ form in C and
> I am using that as a base *after* major surgery.  It came from
> Thomas Flamming in thread:
> 
> http://remotesensing.org/lists/proj_archive/msg00306.html
> 
> As for odd usage of x-y, what's one to do.  At least *our*
> traditional form allows for making proper maps.
> 
> Thanks for the additional digits.
> 
> On Thu, 2003-02-13 at 04:42, Ron Russell wrote:
> > Gerald I. Evenden wrote:
> > >
> > > A few days ago I took up the seemingly minor task of installing
> > > the Czech projection Krovak Oblique Conformal Conic in my version
> > > of PROJ.4.  There seemed to be documentation and benchmarks on
>         ...
> >
> > 10 years ago I implemented this projection in Fortran in our (extended)
> > version of GCTPLIB. I was given a photocopy of the relevent chapter
> > of a Czeck publication (in Czeck!) by our clients, the Czeck Institute
> > for Surveying and Mapping and after a lot of interpretation of the
> > diagrams and figures (but with only the odd word translated)
> > produced the code.
> >
> > This paper used two constants for the Projection of the (Bessel)
> > ellipsoid onto A Gausian Conformal Sphere, (using a latitude of
> > 49 30N), what they called alpha and k.
> >
> > I finally tracked these down, alpha is what you call B and k is much
> > to complicated to write out in an email!
> >
> > As if the source document was not enough of a problem, the Czecks use
> > a left handed coordinate system (Westings and Southings) with the origin
> > at the apex of the cone. Our display systems assume a Easting/Northing
> > so in fact we produced a pseudo Krovak system with a large False Northing;
> > The numbers are correct but they all have the wrong sign, however the
> > maps come out the correct way on the paper!
> >
> > The Czecks also sometimes use Ferra as their prime meridian.
> >
> > My code gives exactly the same results for the points you have listed.
> >
> > The last two points are the test points given in the paper I used (with
> > a few more decimal places)
> >
> > Phi 50 12 32.44140 Lambda 34 30 59.17900 => 568990.9973W 1050538.6493S
> > Phi 50 15 30.70080 lambda 35 05 51.17920 => 527189.1948W 1049221.1899S
> >
> > (They assume Ferra is 17 40e of Greenwich.)
> >
> > Given the problems I had interpreting this document, if it is the
> > definitive source, I am not surprised it there are several intepretations
> > about!
> >
> > Ron
> 
> ----------------------------------------
> PROJ.4 Discussion List
> See http://www.remotesensing.org/proj for subscription, unsubscription
> and other information.

Ron
-- 


---------------------------------------------------------------------
Ron Russell
Senior Software Engineer, Laser-Scan Ltd
+44(0)1223 420414

http://www.laser-scan.com

----------------------------------------
PROJ.4 Discussion List
See http://www.remotesensing.org/proj for subscription, unsubscription
and other information.