<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META content=text/html;charset=iso-8859-1 http-equiv=Content-Type>
<META name=GENERATOR content="MSHTML 8.00.7600.16535"></HEAD>
<BODY style="PADDING-LEFT: 10px; PADDING-RIGHT: 10px; PADDING-TOP: 15px"
id=MailContainerBody leftMargin=0 topMargin=0 CanvasTabStop="true"
name="Compose message area">
<DIV><FONT color=#a70401 size=2>Martin,</FONT></DIV>
<DIV><FONT color=#a70401 size=2></FONT> </DIV>
<DIV><FONT color=#a70401 size=2>Now I need some help with the inverse
projection.</FONT></DIV>
<DIV><FONT color=#a70401 size=2></FONT> </DIV>
<DIV><FONT color=#a70401 size=2>After many days messing with the Albers I have
found where the problems are but they do not make any sense. In the
projection for the ellipsoids, the following line caused the returned longitude
to always be either +/-90 which is what it is supposed to do. My problem
is: why is it there at all?</FONT></DIV>
<DIV><FONT color=#a70401 size=2></FONT> </DIV>
<DIV><FONT color=#a70401 size=2>// line 1 lpphi = lpphi < 0. ?
-ProjectionMath.HALFPI : ProjectionMath.HALFPI;</FONT><BR></DIV>
<DIV><FONT color=#a70401 size=2>When I commented it out the elliptical code gave
me very good results over all 64K tests.</FONT></DIV>
<DIV><FONT color=#a70401 size=2></FONT> </DIV>
<DIV><FONT color=#a70401 size=2>Just to keep things interesting, when working
with the sphere, the following line of code returned a NaN for all latitudes
between -89 and 16 degrees. From 17 thru 89 degrees it returned values
from 1.4964... through 0.2351....</FONT></DIV>
<DIV><FONT color=#a70401 size=2></FONT> </DIV>
<DIV><FONT color=#a70401 size=2> lpphi
= Math.asin(lpphi); // line 2 returns NaN<BR></FONT></DIV>
<DIV><FONT color=#a70401 size=2>When I replaced that line with the code used for
the ellipsoids:</FONT></DIV>
<DIV><FONT color=#a70401 size=2></FONT> </DIV>
<DIV><FONT color=#a70401 size=2> lpphi = (c - lpphi *
lpphi) / n;<BR> lpphi = phi1_(lpphi, e, one_es));</FONT></DIV>
<DIV><FONT color=#a70401 size=2></FONT> </DIV>
<DIV><FONT color=#a70401 size=2>the code for the sphere returned very accurate
results.</FONT></DIV>
<DIV><FONT color=#a70401 size=2></FONT> </DIV>
<DIV><FONT color=#a70401 size=2>Although I am getting good results, this does
not seem to be the way to do it. This just not make any sense. Over
to you for resolution.</FONT></DIV>
<DIV><FONT color=#a70401 size=2></FONT> </DIV>
<DIV><FONT color=#a70401 size=2>Fred</DIV></FONT>
<DIV><FONT color=#a70401 size=2></FONT> </DIV>
<DIV><FONT color=#a70401 size=2>/*<BR>Copyright 2006 Jerry Huxtable</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2>Licensed under the Apache License, Version 2.0
(the "License");<BR>you may not use this file except in compliance with the
License.<BR>You may obtain a copy of the License at</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> <A
title="http://www.apache.org/licenses/LICENSE-2.0 CTRL + Click to follow link"
href="http://www.apache.org/licenses/LICENSE-2.0">http://www.apache.org/licenses/LICENSE-2.0</A></FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2>Unless required by applicable law or agreed to
in writing, software<BR>distributed under the License is distributed on an "AS
IS" BASIS,<BR>WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
implied.<BR>See the License for the specific language governing permissions
and<BR>limitations under the License.<BR>*/</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2>/*<BR> * This file was semi-automatically
converted from the public-domain USGS PROJ source.<BR> */<BR>package
org.osgeo.proj4j.proj;</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2>import java.awt.geom.*;</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2>import
org.osgeo.proj4j.ProjectionMath;<BR>import
org.osgeo.proj4j.Projection;<BR>import
org.osgeo.proj4j.ProjectionException;</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2>public class AlbersProjection extends Projection
{</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> private final static double EPS10 =
1.e-10;<BR> private final static double TOL7 = 1.e-7;<BR> private
double ec;<BR> private double n;<BR> private double
c;<BR> private double dd;<BR> private double n2;<BR> private
double rho0;<BR> private double phi1;<BR> private double
phi2;<BR> private double[] en;</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> private final static int N_ITER =
15;<BR> private final static double EPSILON = 1.0e-7;<BR> private
final static double TOL = 1.0e-10;</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> //protected double projectionLatitude1 =
MapMath.degToRad(45.5);<BR> //protected double projectionLatitude2 =
MapMath.degToRad(29.5);</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> public AlbersProjection()
{<BR> minLatitude = Math.toRadians(0);<BR> maxLatitude =
Math.toRadians(80);<BR> projectionLatitude1 =
ProjectionMath.degToRad(45.5);<BR> projectionLatitude2 =
ProjectionMath.degToRad(29.5);<BR> initialize();<BR> }</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> private static double phi1_(double qs,
double Te, double Tone_es) {<BR> int i;<BR> double Phi,
sinpi, cospi, con, com, dphi;</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> Phi = Math.asin (.5 *
qs);<BR> if (Te < EPSILON)<BR> return( Phi
);<BR> i = N_ITER;<BR> do {<BR> sinpi =
Math.sin (Phi);<BR> cospi = Math.cos
(Phi);<BR> con = Te * sinpi;<BR> com = 1. -
con * con;<BR> dphi = .5 * com * com / cospi * (qs / Tone_es
-<BR> sinpi / com + .5 / Te * Math.log ((1. - con)
/<BR> (1. + con)));<BR> Phi +=
dphi;<BR> } while (Math.abs(dphi) > TOL && --i !=
0);<BR> return( i != 0 ? Phi : Double.MAX_VALUE
);<BR> }</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> public Point2D.Double project(double
lplam, double lpphi, Point2D.Double out) {<BR> double
rho;<BR> if ((rho = c - (!spherical ? n *
ProjectionMath.qsfn(Math.sin(lpphi), e, one_es) : n2 * Math.sin(lpphi))) <
0.)<BR> throw new ProjectionException("F");<BR> rho
= dd * Math.sqrt(rho);<BR> out.x = rho * Math.sin( lplam *= n
);<BR> out.y = rho0 - rho * Math.cos(lplam);<BR> return
out;<BR> }</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> public Point2D.Double
projectInverse(double xyx, double xyy, Point2D.Double out)<BR>
{<BR> double rho;<BR> if ((rho =
ProjectionMath.distance(xyx, xyy = rho0 - xyy)) !=
0)<BR> {<BR> double lpphi,
lplam;<BR> if (n <
0.)<BR>
{<BR> rho = -rho;<BR> xyx =
-xyx;<BR> xyy = -xyy;<BR>
}<BR> lpphi = rho / dd;<BR> if
(!spherical)<BR>
{<BR> lpphi = (c - lpphi * lpphi) /
n;<BR> if (Math.abs(ec - Math.abs(lpphi)) >
TOL7)<BR>
{<BR> if ((lpphi = phi1_(lpphi, e, one_es)) ==
Double.MAX_VALUE)<BR> throw new
ProjectionException("I");<BR>
}<BR>
else<BR> lpphi = lpphi < 0. ?
-ProjectionMath.HALFPI : ProjectionMath.HALFPI;<BR> } //
end if (!spherical)<BR> else if (Math.abs(out.y =
(c - lpphi * lpphi) / n2) <= 1.)<BR>
lpphi = Math.asin(lpphi); // line 2
returns NaN<BR> else<BR>// line 1 lpphi = lpphi < 0. ?
-ProjectionMath.HALFPI : ProjectionMath.HALFPI;<BR>
lplam = Math.atan2(xyx, xyy) / n;<BR> out.x =
lplam;<BR> out.y = lpphi;<BR> } // end
if ((rho = ProjectionMath<BR>
else<BR> {<BR> out.x =
0.;<BR> out.y = n > 0. ? ProjectionMath.HALFPI : -
ProjectionMath.HALFPI;<BR> }<BR> return
out;<BR> } // end projectInverse</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> public void initialize()
{<BR> super.initialize();<BR> double cosphi,
sinphi;<BR> boolean secant;</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> phi1 =
projectionLatitude1;<BR> phi2 = projectionLatitude2;</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> if (Math.abs(phi1 + phi2) <
EPS10)<BR> throw new
IllegalArgumentException("-21");<BR> n = sinphi =
Math.sin(phi1);<BR> cosphi = Math.cos(phi1);<BR> secant =
Math.abs(phi1 - phi2) >= EPS10;<BR> //spherical = es >
0.0;<BR> if (!spherical) {<BR> double ml1,
m1;</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> if ((en =
ProjectionMath.enfn(es)) == null)<BR> throw new
IllegalArgumentException("0");<BR> m1 =
ProjectionMath.msfn(sinphi, cosphi, es);<BR> ml1 =
ProjectionMath.qsfn(sinphi, e, one_es);<BR> if (secant) { /*
secant cone */<BR> double ml2, m2;</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> sinphi =
Math.sin(phi2);<BR> cosphi =
Math.cos(phi2);<BR> m2 = ProjectionMath.msfn(sinphi,
cosphi, es);<BR> ml2 = ProjectionMath.qsfn(sinphi, e,
one_es);<BR> n = (m1 * m1 - m2 * m2) / (ml2 -
ml1);<BR> }<BR> ec = 1. - .5 * one_es *
Math.log((1. - e) /<BR> (1. + e)) /
e;<BR> c = m1 * m1 + n * ml1;<BR> dd = 1. /
n;<BR> rho0 = dd * Math.sqrt(c - n *
ProjectionMath.qsfn(Math.sin(projectionLatitude),<BR> e,
one_es));<BR> } else {<BR> if (secant) n = .5 * (n +
Math.sin(phi2));<BR> n2 = n + n;<BR> c =
cosphi * cosphi + n2 * sinphi;<BR> dd = 1. /
n;<BR> rho0 = dd * Math.sqrt(c - n2 *
Math.sin(projectionLatitude));<BR> }<BR> }</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> /**<BR> * Returns true if this
projection is equal area<BR> */<BR> public boolean isEqualArea()
{<BR> return true;<BR> }</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> public boolean hasInverse()
{<BR> return true;<BR> }</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> /**<BR> * Returns the ESPG code for
this projection, or 0 if unknown.<BR> */<BR> public int getEPSGCode()
{<BR> return 9822;<BR> }</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2> public String toString()
{<BR> return "Albers Equal Area";<BR> }</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2>/*<BR> public String toString()
{<BR> return "Lambert Equal Area
Conic";<BR> }<BR>*/<BR>}</FONT></DIV>
<DIV> </DIV>
<DIV><FONT color=#a70401 size=2></FONT> </DIV></BODY></HTML>