Re: libproj4 stmerc = French Gauss-Laborde projection

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Re: libproj4 stmerc = French Gauss-Laborde projection

strebe

I'm not lost. I've analyzed the method and programmed it. It works fine.

Yes, p is the parametric co-latitude in the first formula. It's wrong to use the same variable in the later formulae because p refers to something else there -- in fact, it's a complex variable in the later instances.

The "trick" is this: Wallis uses the polar stereographic because it's the simplest way to get a conformal mapping to the plane. Once the ellipsoid is mapped, he treats the plane as the complex plane and looks for a complex "co-latitude" which can be used with the polar stereographic, but this time treating the polar stereographic as function of a complex variable. The reason he does this is (a) to preserve conformality; and (b) so that the central meridian (in which the imaginary axis is 0) maps back to the parametric colatitude. At this point the ellipsoid is mapped conformally in such a way that leaves the central meridian effectively unmapped.

Leaving the complex plane aside, using the colatitude as the parameter to the elliptic integral of the second kind gives the distance from the pole to the colatitude. Since this odd mapping Wallis contrived effectively leaves the central meridian unmapped, and since any analytic function applied to a conformal mapping results in a conformal mapping, and since the elliptic integral has an analytic form, all that is left is to push the mapping through the complex form of the elliptic integral of the second kind. This "straightens out" the central meridian to its true differential lengths whilst dragging the whole complex plane with it in a conformal fashion. The result must be the transverse Mercator, since the central meridian is projected with correct scale and since a conformal projection is unique except with respect to scale and rotation.

Regards,
-- daan Strebe



In a message dated 6/15/06 01:18:47, [hidden email] writes:


On Wed, 2006-06-14 at 13:50 -0400, [hidden email] wrote:
> Hm. I didn't know about that web page. Obviously it's wrong -- for some
> reason "p" appears in several different roles. I tend to think that's
> an error in conversion to a web page. (I see that the entire blurb is a
> single graphic, not HTML mark-up.) Certainly he's been pedantic and
> precise in all his communications with me.
>
> The p/2 exponent should read (e/2), where e is the eccentricity.

Yes, I agree.

> Use some other variable (perhaps p') in place of p in "Then, the
> complex variable tan (p/2) can be obtained..." and "...yields the
> argument p..."


Actually the argument p is simply the (ellipsoidal) co-latitude
90d - phi.

The common expression in u and v corresponds to exp(psi), where psi is
the _isometric latitude_, i.e., essentially the "northing" in a
traditional (non-transverse) Mercator map plane.

Isometric latitude and longitude (psi, lambda) together as (x,y)
co-ordinates in a plane define a conformal mapping from the curved
Earth's surface. Using (psi, lambda) directly as rectangular
co-ordinates produces classical Mercator. Using

u + iv = exp(psi + i * lambda)

i.e., polar co-ordinates, produces the stereographic projection. This is
very much what Dr Wallis's formula looks like. Apparently for him it is
only a trick leading somewhere... but then I also get lost.

Regards Martin V

PS you may want to look at

http://users.tkk.fi/~mvermeer/geom.pdf

pp 99-100 and around p. 90. Sorry it's in Fenno-ugrian formulese...




_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: libproj4 stmerc = French Gauss-Laborde projection

Martin Vermeer
On Thu, 2006-06-15 at 05:05 -0400, [hidden email] wrote:

>
> I'm not lost. I've analyzed the method and programmed it. It works
> fine.
>
> Yes, p is the parametric co-latitude in the first formula. It's wrong
> to use the same variable in the later formulae because p refers to
> something else there -- in fact, it's a complex variable in the later
> instances.
>
> The "trick" is this: Wallis uses the polar stereographic because it's
> the simplest way to get a conformal mapping to the plane. Once the
> ellipsoid is mapped, he treats the plane as the complex plane and
> looks for a complex "co-latitude" which can be used with the polar
> stereographic, but this time treating the polar stereographic as
> function of a complex variable. The reason he does this is (a) to
> preserve conformality; and (b) so that the central meridian (in which
> the imaginary axis is 0) maps back to the parametric colatitude. At
> this point the ellipsoid is mapped conformally in such a way that
> leaves the central meridian effectively unmapped.
>
> Leaving the complex plane aside, using the colatitude as the parameter
> to the elliptic integral of the second kind gives the distance from
> the pole to the colatitude. Since this odd mapping Wallis contrived
> effectively leaves the central meridian unmapped, and since any
> analytic function applied to a conformal mapping results in a
> conformal mapping, and since the elliptic integral has an analytic
> form, all that is left is to push the mapping through the complex form
> of the elliptic integral of the second kind. This "straightens out"
> the central meridian to its true differential lengths whilst dragging
> the whole complex plane with it in a conformal fashion. The result
> must be the transverse Mercator, since the central meridian is
> projected with correct scale and since a conformal projection is
> unique except with respect to scale and rotation.
>
> Regards,
> -- daan Strebe
>
Ah! But that's precisely what I have been doing numerically, using a
polynomial expansion rather than an elliptic integral! (And starting
from classical Mercator rather than stereographic, so it will run into
problems at high latitudes.) It's essentially solving a boundary value
problem, with the set of PDEs being the Cauchy-Riemann conformity
conditions and the central meridian the boundary.

I suppose I have to get it written up in english ;-)

Thanks

Martin V

>
> In a message dated 6/15/06 01:18:47, [hidden email] writes:
>
>
> > On Wed, 2006-06-14 at 13:50 -0400, [hidden email] wrote:
> > > Hm. I didn't know about that web page. Obviously it's wrong -- for
> > some
> > > reason "p" appears in several different roles. I tend to think
> > that's
> > > an error in conversion to a web page. (I see that the entire blurb
> > is a
> > > single graphic, not HTML mark-up.) Certainly he's been pedantic
> > and
> > > precise in all his communications with me.
> > >
> > > The p/2 exponent should read (e/2), where e is the eccentricity.
> >
> > Yes, I agree.
> >
> > > Use some other variable (perhaps p') in place of p in "Then, the
> > > complex variable tan (p/2) can be obtained..." and "...yields the
> > > argument p..."
> >
> >
> > Actually the argument p is simply the (ellipsoidal) co-latitude
> > 90d - phi.
> >
> > The common expression in u and v corresponds to exp(psi), where psi
> > is
> > the _isometric latitude_, i.e., essentially the "northing" in a
> > traditional (non-transverse) Mercator map plane.
> >
> > Isometric latitude and longitude (psi, lambda) together as (x,y)
> > co-ordinates in a plane define a conformal mapping from the curved
> > Earth's surface. Using (psi, lambda) directly as rectangular
> > co-ordinates produces classical Mercator. Using
> >
> > u + iv = exp(psi + i * lambda)
> >
> > i.e., polar co-ordinates, produces the stereographic projection.
> > This is
> > very much what Dr Wallis's formula looks like. Apparently for him it
> > is
> > only a trick leading somewhere... but then I also get lost.
> >
> > Regards Martin V
> >
> > PS you may want to look at
> >
> > http://users.tkk.fi/~mvermeer/geom.pdf
> >
> > pp 99-100 and around p. 90. Sorry it's in Fenno-ugrian formulese...
> >
>
>

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj

signature.asc (198 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Re: libproj4 stmerc = French Gauss-Laborde projection

Martin Vermeer
On Thu, 2006-06-15 at 12:27 +0300, Martin Vermeer wrote:

> On Thu, 2006-06-15 at 05:05 -0400, [hidden email] wrote:
> >
> > I'm not lost. I've analyzed the method and programmed it. It works
> > fine.
> >
> > Yes, p is the parametric co-latitude in the first formula. It's wrong
> > to use the same variable in the later formulae because p refers to
> > something else there -- in fact, it's a complex variable in the later
> > instances.
> >
> > The "trick" is this: Wallis uses the polar stereographic because it's
> > the simplest way to get a conformal mapping to the plane. Once the
> > ellipsoid is mapped, he treats the plane as the complex plane and
> > looks for a complex "co-latitude" which can be used with the polar
> > stereographic, but this time treating the polar stereographic as
> > function of a complex variable. The reason he does this is (a) to
> > preserve conformality; and (b) so that the central meridian (in which
> > the imaginary axis is 0) maps back to the parametric colatitude. At
> > this point the ellipsoid is mapped conformally in such a way that
> > leaves the central meridian effectively unmapped.
> >
> > Leaving the complex plane aside, using the colatitude as the parameter
> > to the elliptic integral of the second kind gives the distance from
> > the pole to the colatitude. Since this odd mapping Wallis contrived
> > effectively leaves the central meridian unmapped, and since any
> > analytic function applied to a conformal mapping results in a
> > conformal mapping, and since the elliptic integral has an analytic
> > form, all that is left is to push the mapping through the complex form
> > of the elliptic integral of the second kind. This "straightens out"
> > the central meridian to its true differential lengths whilst dragging
> > the whole complex plane with it in a conformal fashion. The result
> > must be the transverse Mercator, since the central meridian is
> > projected with correct scale and since a conformal projection is
> > unique except with respect to scale and rotation.
> >
> > Regards,
> > -- daan Strebe
> >
>
> Ah! But that's precisely what I have been doing numerically, using a
> polynomial expansion rather than an elliptic integral! (And starting
> from classical Mercator rather than stereographic, so it will run into
> problems at high latitudes.) It's essentially solving a boundary value
> problem, with the set of PDEs being the Cauchy-Riemann conformity
> conditions and the central meridian the boundary.
>
> I suppose I have to get it written up in english ;-)

Just found my Matlab routines for playing with this. Posted them here:

http://users.tkk.fi/~mvermeer/gk.zip

You run it as

>> gkmain(phi, lab, [phimin phimax delta maxk refphi]);

e.g.

>> gkmain(60, 5, [55 65 0.5 12 60])

(lab is the distance from the central meridian)

With this I have also obtained sub-millimeter precision (in typical
situations). Yes, it breaks down at the poles due to using Mercator as
the starting projection. Replacing this by stereographic as Dr Wallis
did, shouldn't be hard.

Using a polynomial expansion (instead of the elliptic integral /
Newton-Raphson thing) has the advantage of simplicity in implementation.
One disadvantage that I noticed is, that the normal matrix for
estimating the coefficients gets poorly conditioned for expansion degree
maxk = 15 or more.

I used this code in an exercise for my students; in a production
environment you would print out and hardwire the estimated coefficients
for the area of interest into your code.

The exercise instruction (in Finnish :) at

http://users.tkk.fi/~mvermeer/kotiharjB.pdf


- Martin

> >
> > In a message dated 6/15/06 01:18:47, [hidden email] writes:
> >
> >
> > > On Wed, 2006-06-14 at 13:50 -0400, [hidden email] wrote:
> > > > Hm. I didn't know about that web page. Obviously it's wrong -- for
> > > some
> > > > reason "p" appears in several different roles. I tend to think
> > > that's
> > > > an error in conversion to a web page. (I see that the entire blurb
> > > is a
> > > > single graphic, not HTML mark-up.) Certainly he's been pedantic
> > > and
> > > > precise in all his communications with me.
> > > >
> > > > The p/2 exponent should read (e/2), where e is the eccentricity.
> > >
> > > Yes, I agree.
> > >
> > > > Use some other variable (perhaps p') in place of p in "Then, the
> > > > complex variable tan (p/2) can be obtained..." and "...yields the
> > > > argument p..."
> > >
> > >
> > > Actually the argument p is simply the (ellipsoidal) co-latitude
> > > 90d - phi.
> > >
> > > The common expression in u and v corresponds to exp(psi), where psi
> > > is
> > > the _isometric latitude_, i.e., essentially the "northing" in a
> > > traditional (non-transverse) Mercator map plane.
> > >
> > > Isometric latitude and longitude (psi, lambda) together as (x,y)
> > > co-ordinates in a plane define a conformal mapping from the curved
> > > Earth's surface. Using (psi, lambda) directly as rectangular
> > > co-ordinates produces classical Mercator. Using
> > >
> > > u + iv = exp(psi + i * lambda)
> > >
> > > i.e., polar co-ordinates, produces the stereographic projection.
> > > This is
> > > very much what Dr Wallis's formula looks like. Apparently for him it
> > > is
> > > only a trick leading somewhere... but then I also get lost.
> > >
> > > Regards Martin V
> > >
> > > PS you may want to look at
> > >
> > > http://users.tkk.fi/~mvermeer/geom.pdf
> > >
> > > pp 99-100 and around p. 90. Sorry it's in Fenno-ugrian formulese...
> > >
> >
> >
> _______________________________________________
> Proj mailing list
> [hidden email]
> http://lists.maptools.org/mailman/listinfo/proj

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj

signature.asc (198 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: libproj4 stmerc = French Gauss-Laborde projection

strebe
In reply to this post by strebe

I heartily recommend the AGM method for computing the elliptic integral. It's fast and accurate everywhere, and can be implemented to arbitrary precision. Polynomials seem a poor choice for this.

See Abramowitz & Stegun, p. 598 (17.6). Bulirsch (1965) supplies an excellent implementation of the AGM, but you can use it as described in Abramowitz & Stegun with great efficiency.

Regards
-- daan Strebe


In a message dated 6/22/06 00:39:37, [hidden email] writes:


On Thu, 2006-06-15 at 12:27 +0300, Martin Vermeer wrote:
> On Thu, 2006-06-15 at 05:05 -0400, [hidden email] wrote:
> >
> > I'm not lost. I've analyzed the method and programmed it. It works
> > fine.
> >
> > Yes, p is the parametric co-latitude in the first formula. It's wrong
> > to use the same variable in the later formulae because p refers to
> > something else there -- in fact, it's a complex variable in the later
> > instances.
> >
> > The "trick" is this: Wallis uses the polar stereographic because it's
> > the simplest way to get a conformal mapping to the plane. Once the
> > ellipsoid is mapped, he treats the plane as the complex plane and
> > looks for a complex "co-latitude" which can be used with the polar
> > stereographic, but this time treating the polar stereographic as
> > function of a complex variable. The reason he does this is (a) to
> > preserve conformality; and (b) so that the central meridian (in which
> > the imaginary axis is 0) maps back to the parametric colatitude. At
> > this point the ellipsoid is mapped conformally in such a way that
> > leaves the central meridian effectively unmapped.
> >
> > Leaving the complex plane aside, using the colatitude as the parameter
> > to the elliptic integral of the second kind gives the distance from
> > the pole to the colatitude. Since this odd mapping Wallis contrived
> > effectively leaves the central meridian unmapped, and since any
> > analytic function applied to a conformal mapping results in a
> > conformal mapping, and since the elliptic integral has an analytic
> > form, all that is left is to push the mapping through the complex form
> > of the elliptic integral of the second kind. This "straightens out"
> > the central meridian to its true differential lengths whilst dragging
> > the whole complex plane with it in a conformal fashion. The result
> > must be the transverse Mercator, since the central meridian is
> > projected with correct scale and since a conformal projection is
> > unique except with respect to scale and rotation.
> >
> > Regards,
> > -- daan Strebe
> >
>
> Ah! But that's precisely what I have been doing numerically, using a
> polynomial expansion rather than an elliptic integral! (And starting
> from classical Mercator rather than stereographic, so it will run into
> problems at high latitudes.) It's essentially solving a boundary value
> problem, with the set of PDEs being the Cauchy-Riemann conformity
> conditions and the central meridian the boundary.
>
> I suppose I have to get it written up in english ;-)


Just found my Matlab routines for playing with this. Posted them here:

http://users.tkk.fi/~mvermeer/gk.zip

You run it as

>> gkmain(phi, lab, [phimin phimax delta maxk refphi]);

e.g.

>> gkmain(60, 5, [55 65 0.5 12 60])

(lab is the distance from the central meridian)

With this I have also obtained sub-millimeter precision (in typical
situations). Yes, it breaks down at the poles due to using Mercator as
the starting projection. Replacing this by stereographic as Dr Wallis
did, shouldn't be hard.

Using a polynomial expansion (instead of the elliptic integral /
Newton-Raphson thing) has the advantage of simplicity in implementation.
One disadvantage that I noticed is, that the normal matrix for
estimating the coefficients gets poorly conditioned for expansion degree
maxk = 15 or more.

I used this code in an exercise for my students; in a production
environment you would print out and hardwire the estimated coefficients
for the area of interest into your code.

The exercise instruction (in Finnish :) at

http://users.tkk.fi/~mvermeer/kotiharjB.pdf


- Martin



_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj