Index: trunk/src/org/openstreetmap/josm/data/projection/proj/AbstractProj.java =================================================================== --- trunk/src/org/openstreetmap/josm/data/projection/proj/AbstractProj.java (revision 9112) +++ trunk/src/org/openstreetmap/josm/data/projection/proj/AbstractProj.java (revision 9112) @@ -0,0 +1,123 @@ +// License: GPL. For details, see LICENSE file. +package org.openstreetmap.josm.data.projection.proj; + +import org.openstreetmap.josm.data.projection.ProjectionConfigurationException; + +/** + * Abstract base class providing utilities for implementations of the Proj + * interface. + * + * This class has been derived from the implementation of the Geotools project; + * git 8cbf52d, org.geotools.referencing.operation.projection.MapProjection + * at the time of migration. + *

+ * + * @author André Gosselin + * @author Martin Desruisseaux (PMO, IRD) + * @author Rueben Schulz +*/ +public abstract class AbstractProj implements Proj { + + /** + * Maximum number of iterations for iterative computations. + */ + private static final int MAXIMUM_ITERATIONS = 15; + + /** + * Relative iteration precision used in the mlfn method + */ + private static final double MLFN_TOL = 1E-11; + + /** + * Constants used to calculate {@link #en0}, {@link #en1}, + * {@link #en2}, {@link #en3}, {@link #en4}. + */ + private static final double C00= 1.0, + C02= 0.25, + C04= 0.046875, + C06= 0.01953125, + C08= 0.01068115234375, + C22= 0.75, + C44= 0.46875, + C46= 0.01302083333333333333, + C48= 0.00712076822916666666, + C66= 0.36458333333333333333, + C68= 0.00569661458333333333, + C88= 0.3076171875; + + /** + * Constant needed for the mlfn method. + * Setup at construction time. + */ + protected double en0,en1,en2,en3,en4; + + /** + * The square of excentricity: e² = (a²-b²)/a² where + * e is the {@linkplain #excentricity excentricity}, + * a is the {@linkplain #semiMajor semi major} axis length and + * b is the {@linkplain #semiMinor semi minor} axis length. + */ + protected double e2; + + @Override + public void initialize(ProjParameters params) throws ProjectionConfigurationException { + e2 = params.ellps.e2; + // Compute constants for the mlfn + double t; + en0 = C00 - e2 * (C02 + e2 * + (C04 + e2 * (C06 + e2 * C08))); + en1 = e2 * (C22 - e2 * + (C04 + e2 * (C06 + e2 * C08))); + en2 = (t = e2 * e2) * + (C44 - e2 * (C46 + e2 * C48)); + en3 = (t *= e2) * (C66 - e2 * C68); + en4 = t * e2 * C88; + } + + /** + * Calculates the meridian distance. This is the distance along the central + * meridian from the equator to {@code phi}. Accurate to < 1e-5 meters + * when used in conjuction with typical major axis values. + * + * @param phi latitude to calculate meridian distance for. + * @param sphi sin(phi). + * @param cphi cos(phi). + * @return meridian distance for the given latitude. + */ + protected final double mlfn(final double phi, double sphi, double cphi) { + cphi *= sphi; + sphi *= sphi; + return en0 * phi - cphi * + (en1 + sphi * + (en2 + sphi * + (en3 + sphi * + (en4)))); + } + + /** + * Calculates the latitude ({@code phi}) from a meridian distance. + * Determines phi to TOL (1e-11) radians, about 1e-6 seconds. + * + * @param arg meridian distance to calulate latitude for. + * @return the latitude of the meridian distance. + * @throws RuntimeException if the itteration does not converge. + */ + protected final double inv_mlfn(double arg) { + double s, t, phi, k = 1.0/(1.0 - e2); + int i; + phi = arg; + for (i=MAXIMUM_ITERATIONS; true;) { // rarely goes over 5 iterations + if (--i < 0) { + throw new RuntimeException("Too many iterations"); + } + s = Math.sin(phi); + t = 1.0 - e2 * s * s; + t = (mlfn(phi, s, Math.cos(phi)) - arg) * (t * Math.sqrt(t)) * k; + phi -= t; + if (Math.abs(t) < MLFN_TOL) { + return phi; + } + } + } + +} Index: trunk/src/org/openstreetmap/josm/data/projection/proj/TransverseMercator.java =================================================================== --- trunk/src/org/openstreetmap/josm/data/projection/proj/TransverseMercator.java (revision 9108) +++ trunk/src/org/openstreetmap/josm/data/projection/proj/TransverseMercator.java (revision 9112) @@ -2,9 +2,4 @@ package org.openstreetmap.josm.data.projection.proj; -import static java.lang.Math.cos; -import static java.lang.Math.pow; -import static java.lang.Math.sin; -import static java.lang.Math.sqrt; -import static java.lang.Math.tan; import static org.openstreetmap.josm.tools.I18n.tr; @@ -12,13 +7,101 @@ /** - * Transverse Mercator projection. + * Transverse Mercator Projection (EPSG code 9807). This + * is a cylindrical projection, in which the cylinder has been rotated 90°. + * Instead of being tangent to the equator (or to an other standard latitude), + * it is tangent to a central meridian. Deformation are more important as we + * are going futher from the central meridian. The Transverse Mercator + * projection is appropriate for region wich have a greater extent north-south + * than east-west. + *

* - * @author Dirk Stöcker - * code based on JavaScript from Chuck Taylor + * The elliptical equations used here are series approximations, and their accuracy + * decreases as points move farther from the central meridian of the projection. + * The forward equations here are accurate to a less than a mm ±10 degrees from + * the central meridian, a few mm ±15 degrees from the + * central meridian and a few cm ±20 degrees from the central meridian. + * The spherical equations are not approximations and should always give the + * correct values. + *

* + * There are a number of versions of the transverse mercator projection + * including the Universal (UTM) and Modified (MTM) Transverses Mercator + * projections. In these cases the earth is divided into zones. For the UTM + * the zones are 6 degrees wide, numbered from 1 to 60 proceeding east from + * 180 degrees longitude, and between lats 84 degrees North and 80 + * degrees South. The central meridian is taken as the center of the zone + * and the latitude of origin is the equator. A scale factor of 0.9996 and + * false easting of 500000m is used for all zones and a false northing of 10000000m + * is used for zones in the southern hemisphere. + *

+ * + * NOTE: formulas used below are not those of Snyder, but rather those + * from the {@code proj4} package of the USGS survey, which + * have been reproduced verbatim. USGS work is acknowledged here. + *

+ * + * This class has been derived from the implementation of the Geotools project; + * git 8cbf52d, org.geotools.referencing.operation.projection.TransverseMercator + * at the time of migration. + *

+ * + * References: + *

+ * + * @see Transverse Mercator projection on MathWorld + * @see "Transverse_Mercator" on RemoteSensing.org + * + * @author André Gosselin + * @author Martin Desruisseaux (PMO, IRD) + * @author Rueben Schulz */ -public class TransverseMercator implements Proj { +public class TransverseMercator extends AbstractProj { - protected double a, b; + /** + * Contants used for the forward and inverse transform for the eliptical + * case of the Transverse Mercator. + */ + private static final double FC1= 1.00000000000000000000000, // 1/1 + FC2= 0.50000000000000000000000, // 1/2 + FC3= 0.16666666666666666666666, // 1/6 + FC4= 0.08333333333333333333333, // 1/12 + FC5= 0.05000000000000000000000, // 1/20 + FC6= 0.03333333333333333333333, // 1/30 + FC7= 0.02380952380952380952380, // 1/42 + FC8= 0.01785714285714285714285; // 1/56 + + /** + * Maximum difference allowed when comparing real numbers. + */ + private static final double EPSILON = 1E-6; + + /** + * A derived quantity of excentricity, computed by e'² = (a²-b²)/b² = es/(1-es) + * where a is the semi-major axis length and b is the semi-minor axis + * length. + */ + private double eb2; + + /** + * Latitude of origin in radians. Default value is 0, the equator. + * This is called 'phi0' in Snyder. + *

+ * Consider this field as final. It is not final only + * because some classes need to modify it at construction time. + */ + protected double latitudeOfOrigin; + + /** + * Meridian distance at the {@code latitudeOfOrigin}. + * Used for calculations for the ellipsoid. + */ + private double ml0; @Override @@ -34,253 +117,66 @@ @Override public void initialize(ProjParameters params) throws ProjectionConfigurationException { - this.a = params.ellps.a; - this.b = params.ellps.b; + super.initialize(params); + eb2 = params.ellps.eb2; + latitudeOfOrigin = params.lat0 == null ? 0 : Math.toRadians(params.lat0); + ml0 = mlfn(latitudeOfOrigin, Math.sin(latitudeOfOrigin), Math.cos(latitudeOfOrigin)); } - /** - * Converts a latitude/longitude pair to x and y coordinates in the - * Transverse Mercator projection. Note that Transverse Mercator is not - * the same as UTM; a scale factor is required to convert between them. - * - * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J., - * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994. - * - * @param phi Latitude of the point, in radians - * @param lambda Longitude of the point, in radians - * @return A 2-element array containing the x and y coordinates - * of the computed point - */ @Override - public double[] project(double phi, double lambda) { + public double[] project(double y, double x) { + double sinphi = Math.sin(y); + double cosphi = Math.cos(y); - /* Precalculate ep2 */ - double ep2 = (pow(a, 2.0) - pow(b, 2.0)) / pow(b, 2.0); + double t = (Math.abs(cosphi) > EPSILON) ? sinphi/cosphi : 0; + t *= t; + double al = cosphi*x; + double als = al*al; + al /= Math.sqrt(1.0 - e2 * sinphi*sinphi); + double n = eb2 * cosphi*cosphi; - /* Precalculate nu2 */ - double nu2 = ep2 * pow(cos(phi), 2.0); + /* NOTE: meridinal distance at latitudeOfOrigin is always 0 */ + y = (mlfn(y, sinphi, cosphi) - ml0 + + sinphi * al * x * + FC2 * ( 1.0 + + FC4 * als * (5.0 - t + n*(9.0 + 4.0*n) + + FC6 * als * (61.0 + t * (t - 58.0) + n*(270.0 - 330.0*t) + + FC8 * als * (1385.0 + t * ( t*(543.0 - t) - 3111.0)))))); - /* Precalculate N / a */ - double N_a = a / (b * sqrt(1 + nu2)); + x = al*(FC1 + FC3 * als*(1.0 - t + n + + FC5 * als * (5.0 + t*(t - 18.0) + n*(14.0 - 58.0*t) + + FC7 * als * (61.0+ t*(t*(179.0 - t) - 479.0 ))))); - /* Precalculate t */ - double t = tan(phi); - double t2 = t * t; - - /* Precalculate l */ - double l = lambda; - - /* Precalculate coefficients for l**n in the equations below - so a normal human being can read the expressions for easting - and northing - -- l**1 and l**2 have coefficients of 1.0 */ - double l3coef = 1.0 - t2 + nu2; - - double l4coef = 5.0 - t2 + 9 * nu2 + 4.0 * (nu2 * nu2); - - double l5coef = 5.0 - 18.0 * t2 + (t2 * t2) + 14.0 * nu2 - - 58.0 * t2 * nu2; - - double l6coef = 61.0 - 58.0 * t2 + (t2 * t2) + 270.0 * nu2 - - 330.0 * t2 * nu2; - - double l7coef = 61.0 - 479.0 * t2 + 179.0 * (t2 * t2) - (t2 * t2 * t2); - - double l8coef = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2); - - return new double[] { - /* Calculate easting (x) */ - N_a * cos(phi) * l - + (N_a / 6.0 * pow(cos(phi), 3.0) * l3coef * pow(l, 3.0)) - + (N_a / 120.0 * pow(cos(phi), 5.0) * l5coef * pow(l, 5.0)) - + (N_a / 5040.0 * pow(cos(phi), 7.0) * l7coef * pow(l, 7.0)), - /* Calculate northing (y) */ - ArcLengthOfMeridian(phi) / a - + (t / 2.0 * N_a * pow(cos(phi), 2.0) * pow(l, 2.0)) - + (t / 24.0 * N_a * pow(cos(phi), 4.0) * l4coef * pow(l, 4.0)) - + (t / 720.0 * N_a * pow(cos(phi), 6.0) * l6coef * pow(l, 6.0)) - + (t / 40320.0 * N_a * pow(cos(phi), 8.0) * l8coef * pow(l, 8.0)) }; + return new double[] { x, y }; } - /** - * Converts x and y coordinates in the Transverse Mercator projection to - * a latitude/longitude pair. Note that Transverse Mercator is not - * the same as UTM; a scale factor is required to convert between them. - * - * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J., - * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994. - * - * Remarks: - * The local variables Nf, nuf2, tf, and tf2 serve the same purpose as - * N, nu2, t, and t2 in MapLatLonToXY, but they are computed with respect - * to the footpoint latitude phif. - * - * x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and - * to optimize computations. - * - * @param x The easting of the point, in meters, divided by the semi major axis of the ellipsoid - * @param y The northing of the point, in meters, divided by the semi major axis of the ellipsoid - * @return A 2-element containing the latitude and longitude - * in radians - */ @Override public double[] invproject(double x, double y) { - /* Get the value of phif, the footpoint latitude. */ - double phif = footpointLatitude(y); + double phi = inv_mlfn(ml0 + y); - /* Precalculate ep2 */ - double ep2 = (a*a - b*b) - / (b*b); + if (Math.abs(phi) >= Math.PI/2) { + y = y<0.0 ? -(Math.PI/2) : (Math.PI/2); + x = 0.0; + } else { + double sinphi = Math.sin(phi); + double cosphi = Math.cos(phi); + double t = (Math.abs(cosphi) > EPSILON) ? sinphi/cosphi : 0.0; + double n = eb2 * cosphi*cosphi; + double con = 1.0 - e2 * sinphi*sinphi; + double d = x * Math.sqrt(con); + con *= t; + t *= t; + double ds = d*d; - /* Precalculate cos(phif) */ - double cf = cos(phif); + y = phi - (con*ds / (1.0 - e2)) * + FC2 * (1.0 - ds * + FC4 * (5.0 + t*(3.0 - 9.0*n) + n*(1.0 - 4*n) - ds * + FC6 * (61.0 + t*(90.0 - 252.0*n + 45.0*t) + 46.0*n - ds * + FC8 * (1385.0 + t*(3633.0 + t*(4095.0 + 1574.0*t)))))); - /* Precalculate nuf2 */ - double nuf2 = ep2 * pow(cf, 2.0); - - /* Precalculate Nf / a and initialize Nfpow */ - double Nf_a = a / (b * sqrt(1 + nuf2)); - double Nfpow = Nf_a; - - /* Precalculate tf */ - double tf = tan(phif); - double tf2 = tf * tf; - double tf4 = tf2 * tf2; - - /* Precalculate fractional coefficients for x**n in the equations - below to simplify the expressions for latitude and longitude. */ - double x1frac = 1.0 / (Nfpow * cf); - - Nfpow *= Nf_a; /* now equals Nf**2) */ - double x2frac = tf / (2.0 * Nfpow); - - Nfpow *= Nf_a; /* now equals Nf**3) */ - double x3frac = 1.0 / (6.0 * Nfpow * cf); - - Nfpow *= Nf_a; /* now equals Nf**4) */ - double x4frac = tf / (24.0 * Nfpow); - - Nfpow *= Nf_a; /* now equals Nf**5) */ - double x5frac = 1.0 / (120.0 * Nfpow * cf); - - Nfpow *= Nf_a; /* now equals Nf**6) */ - double x6frac = tf / (720.0 * Nfpow); - - Nfpow *= Nf_a; /* now equals Nf**7) */ - double x7frac = 1.0 / (5040.0 * Nfpow * cf); - - Nfpow *= Nf_a; /* now equals Nf**8) */ - double x8frac = tf / (40320.0 * Nfpow); - - /* Precalculate polynomial coefficients for x**n. - -- x**1 does not have a polynomial coefficient. */ - double x2poly = -1.0 - nuf2; - double x3poly = -1.0 - 2 * tf2 - nuf2; - double x4poly = 5.0 + 3.0 * tf2 + 6.0 * nuf2 - 6.0 * tf2 * nuf2 - 3.0 * (nuf2 *nuf2) - 9.0 * tf2 * (nuf2 * nuf2); - double x5poly = 5.0 + 28.0 * tf2 + 24.0 * tf4 + 6.0 * nuf2 + 8.0 * tf2 * nuf2; - double x6poly = -61.0 - 90.0 * tf2 - 45.0 * tf4 - 107.0 * nuf2 + 162.0 * tf2 * nuf2; - double x7poly = -61.0 - 662.0 * tf2 - 1320.0 * tf4 - 720.0 * (tf4 * tf2); - double x8poly = 1385.0 + 3633.0 * tf2 + 4095.0 * tf4 + 1575 * (tf4 * tf2); - - return new double[] { - /* Calculate latitude */ - phif + x2frac * x2poly * (x * x) - + x4frac * x4poly * pow(x, 4.0) - + x6frac * x6poly * pow(x, 6.0) - + x8frac * x8poly * pow(x, 8.0), - /* Calculate longitude */ - x1frac * x - + x3frac * x3poly * pow(x, 3.0) - + x5frac * x5poly * pow(x, 5.0) - + x7frac * x7poly * pow(x, 7.0) }; - } - - /** - * ArcLengthOfMeridian - * - * Computes the ellipsoidal distance from the equator to a point at a - * given latitude. - * - * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J., - * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994. - * - * @param phi Latitude of the point, in radians - * @return The ellipsoidal distance of the point from the equator - * (in meters, divided by the semi major axis of the ellipsoid) - */ - private double ArcLengthOfMeridian(double phi) { - /* Precalculate n */ - double n = (a - b) / (a + b); - - /* Precalculate alpha */ - double alpha = ((a + b) / 2.0) - * (1.0 + (pow(n, 2.0) / 4.0) + (pow(n, 4.0) / 64.0)); - - /* Precalculate beta */ - double beta = (-3.0 * n / 2.0) + (9.0 * pow(n, 3.0) / 16.0) - + (-3.0 * pow(n, 5.0) / 32.0); - - /* Precalculate gamma */ - double gamma = (15.0 * pow(n, 2.0) / 16.0) - + (-15.0 * pow(n, 4.0) / 32.0); - - /* Precalculate delta */ - double delta = (-35.0 * pow(n, 3.0) / 48.0) - + (105.0 * pow(n, 5.0) / 256.0); - - /* Precalculate epsilon */ - double epsilon = 315.0 * pow(n, 4.0) / 512.0; - - /* Now calculate the sum of the series and return */ - return alpha - * (phi + (beta * sin(2.0 * phi)) - + (gamma * sin(4.0 * phi)) - + (delta * sin(6.0 * phi)) - + (epsilon * sin(8.0 * phi))); - } - - /** - * FootpointLatitude - * - * Computes the footpoint latitude for use in converting transverse - * Mercator coordinates to ellipsoidal coordinates. - * - * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J., - * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994. - * - * @param y northing coordinate, in meters, divided by the semi major axis of the ellipsoid - * @return The footpoint latitude, in radians - */ - private double footpointLatitude(double y) { - /* Precalculate n (Eq. 10.18) */ - double n = (a - b) / (a + b); - - /* Precalculate alpha_ (Eq. 10.22) */ - /* (Same as alpha in Eq. 10.17) */ - double alpha_ = ((a + b) / 2.0) - * (1 + (pow(n, 2.0) / 4) + (pow(n, 4.0) / 64)); - - /* Precalculate y_ (Eq. 10.23) */ - double y_ = y / alpha_ * a; - - /* Precalculate beta_ (Eq. 10.22) */ - double beta_ = (3.0 * n / 2.0) + (-27.0 * pow(n, 3.0) / 32.0) - + (269.0 * pow(n, 5.0) / 512.0); - - /* Precalculate gamma_ (Eq. 10.22) */ - double gamma_ = (21.0 * pow(n, 2.0) / 16.0) - + (-55.0 * pow(n, 4.0) / 32.0); - - /* Precalculate delta_ (Eq. 10.22) */ - double delta_ = (151.0 * pow(n, 3.0) / 96.0) - + (-417.0 * pow(n, 5.0) / 128.0); - - /* Precalculate epsilon_ (Eq. 10.22) */ - double epsilon_ = 1097.0 * pow(n, 4.0) / 512.0; - - /* Now calculate the sum of the series (Eq. 10.21) */ - return y_ + (beta_ * sin(2.0 * y_)) - + (gamma_ * sin(4.0 * y_)) - + (delta_ * sin(6.0 * y_)) - + (epsilon_ * sin(8.0 * y_)); + x = d*(FC1 - ds * FC3 * (1.0 + 2.0*t + n - + ds*FC5*(5.0 + t*(28.0 + 24* t + 8.0*n) + 6.0*n - + ds*FC7*(61.0 + t*(662.0 + t*(1320.0 + 720.0*t))))))/cosphi; + } + return new double[] { y, x }; } }