Ticket #3214: Projection_dir.patch

File Projection_dir.patch, 32.6 KB (added by pieren, 15 years ago)

The whole directory patch file

  • Ellipsoid.java

     
    1111 * the reference ellipsoids
    1212 */
    1313class Ellipsoid {
    14   /**
    15    * Clarke's elipsoid (NTF system)
    16    */
    17   public static final Ellipsoid clarke = new Ellipsoid(6378249.2, 6356515.0);
    18   /**
    19    * Hayford's ellipsoid (ED50 system)
    20    */
    21   public static final Ellipsoid hayford =
    22     new Ellipsoid(6378388.0, 6356911.9461);
    23   /**
    24    * WGS84 elipsoid
    25    */
    26   public static final Ellipsoid GRS80 = new Ellipsoid(6378137.0, 6356752.314);
     14    /**
     15     * Clarke's ellipsoid (NTF system)
     16     */
     17    public static final Ellipsoid clarke = new Ellipsoid(6378249.2, 6356515.0);
     18    /**
     19     * Hayford's ellipsoid (ED50 system)
     20     */
     21    public static final Ellipsoid hayford =
     22        new Ellipsoid(6378388.0, 6356911.9461);
     23    /**
     24     * WGS84 ellipsoid
     25     */
     26    public static final Ellipsoid GRS80 = new Ellipsoid(6378137.0, 6356752.314);
    2727
    28   /**
    29    * half long axis
    30    */
    31   public final double a;
    32   /**
    33    * half short axis
    34    */
    35   public final double b;
    36   /**
    37    * first excentricity
    38    */
    39   public final double e;
    40   /**
    41    * first excentricity squared
    42    */
    43   public final double e2;
     28    /**
     29     * half long axis
     30     */
     31    public final double a;
     32    /**
     33     * half short axis
     34     */
     35    public final double b;
     36    /**
     37     * first eccentricity
     38     */
     39    public final double e;
     40    /**
     41     * first eccentricity squared
     42     */
     43    public final double e2;
    4444
    45   /**
    46    * create a new ellipsoid and precompute its parameters
    47    *
    48    * @param a ellipsoid long axis (in meters)
    49    * @param b ellipsoid short axis (in meters)
    50    */
    51   public Ellipsoid(double a, double b) {
    52     this.a = a;
    53     this.b = b;
    54     e2 = (a*a - b*b) / (a*a);
    55     e = Math.sqrt(e2);
    56   }
     45    /**
     46     * square of the second eccentricity
     47     */
     48    public final double eb2;
     49
     50    /**
     51     * create a new ellipsoid and precompute its parameters
     52     *
     53     * @param a ellipsoid long axis (in meters)
     54     * @param b ellipsoid short axis (in meters)
     55     */
     56    public Ellipsoid(double a, double b) {
     57        this.a = a;
     58        this.b = b;
     59        e2 = (a*a - b*b) / (a*a);
     60        e = Math.sqrt(e2);
     61        eb2 = e2 / (1.0 - e2);
     62    }
     63
     64    /**
     65     * Returns the <i>radius of curvature in the prime vertical</i>
     66     * for this reference ellipsoid at the specified latitude.
     67     *
     68     * @param phi The local latitude (radians).
     69     * @return The radius of curvature in the prime vertical (meters).
     70     */
     71    public double verticalRadiusOfCurvature(final double phi) {
     72        return a / Math.sqrt(1.0 - (e2 * sqr(Math.sin(phi))));
     73    }
     74
     75    private static double sqr(final double x) {
     76        return x * x;
     77    }
     78
     79    /**
     80     *  Returns the meridional arc, the true meridional distance on the
     81     * ellipsoid from the equator to the specified latitude, in meters.
     82     *
     83     * @param phi   The local latitude (in radians).
     84     * @return  The meridional arc (in meters).
     85     */
     86    public double meridionalArc(final double phi) {
     87        final double sin2Phi = Math.sin(2.0 * phi);
     88        final double sin4Phi = Math.sin(4.0 * phi);
     89        final double sin6Phi = Math.sin(6.0 * phi);
     90        final double sin8Phi = Math.sin(8.0 * phi);
     91        // TODO . calculate 'f'
     92        //double f = 1.0 / 298.257222101; // GRS80
     93        double f = 1.0 / 298.257223563; // WGS84
     94        final double n = f / (2.0 - f);
     95        final double n2 = n * n;
     96        final double n3 = n2 * n;
     97        final double n4 = n3 * n;
     98        final double n5 = n4 * n;
     99        final double n1n2 = n - n2;
     100        final double n2n3 = n2 - n3;
     101        final double n3n4 = n3 - n4;
     102        final double n4n5 = n4 - n5;
     103        final double ap = a * (1.0 - n + (5.0 / 4.0) * (n2n3) + (81.0 / 64.0) * (n4n5));
     104        final double bp = (3.0 / 2.0) * a * (n1n2 + (7.0 / 8.0) * (n3n4) + (55.0 / 64.0) * n5);
     105        final double cp = (15.0 / 16.0) * a * (n2n3 + (3.0 / 4.0) * (n4n5));
     106        final double dp = (35.0 / 48.0) * a * (n3n4 + (11.0 / 16.0) * n5);
     107        final double ep = (315.0 / 512.0) * a * (n4n5);
     108        return ap * phi - bp * sin2Phi + cp * sin4Phi - dp * sin6Phi + ep * sin8Phi;
     109    }
     110
     111    /**
     112     *  Returns the <i>radius of curvature in the meridian<i>
     113     *  for this reference ellipsoid at the specified latitude.
     114     *
     115     * @param phi The local latitude (in radians).
     116     * @return  The radius of curvature in the meridian (in meters).
     117     */
     118    public double meridionalRadiusOfCurvature(final double phi) {
     119        return verticalRadiusOfCurvature(phi)
     120        / (1.0 + eb2 * sqr(Math.cos(phi)));
     121    }
     122
     123    /**
     124     * Returns isometric latitude of phi on given first eccentricity (e)
     125     * @param phi The local latitude (radians).
     126     * @return isometric latitude of phi on first eccentricity (e)
     127     */
     128    public double latitudeIsometric(double phi, double e) {
     129        double v1 = 1-e*Math.sin(phi);
     130        double v2 = 1+e*Math.sin(phi);
     131        return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2,e/2));
     132    }
     133
     134    /**
     135     * Returns isometric latitude of phi on first eccentricity (e)
     136     * @param phi The local latitude (radians).
     137     * @return isometric latitude of phi on first eccentricity (e)
     138     */
     139    public double latitudeIsometric(double phi) {
     140        double v1 = 1-e*Math.sin(phi);
     141        double v2 = 1+e*Math.sin(phi);
     142        return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2,e/2));
     143    }
     144
     145    /*
     146     * Returns geographic latitude of isometric latitude of first eccentricity (e)
     147     * and epsilon precision
     148     */
     149    public double latitude(double latIso, double e, double epsilon) {
     150        double lat0 = 2*Math.atan(Math.exp(latIso))-Math.PI/2;
     151        double lati = lat0;
     152        double lati1 = 1.0; // random value to start the iterative processus
     153        while(Math.abs(lati1-lati)>=epsilon) {
     154            lati = lati1;
     155            double v1 = 1+e*Math.sin(lati);
     156            double v2 = 1-e*Math.sin(lati);
     157            lati1 = 2*Math.atan(Math.pow(v1/v2,e/2)*Math.exp(latIso))-Math.PI/2;
     158        }
     159        return lati1;
     160    }
    57161}
    58162
    59163
     164

  • GaussLaborde_Reunion.java

     
     1// License: GPL. For details, see LICENSE file.
     2package org.openstreetmap.josm.data.projection;
     3
     4import static org.openstreetmap.josm.tools.I18n.tr;
     5
     6import org.openstreetmap.josm.data.Bounds;
     7import org.openstreetmap.josm.data.coor.EastNorth;
     8import org.openstreetmap.josm.data.coor.LatLon;
     9
     10public class GaussLaborde_Reunion implements Projection {
     11
     12    private static final double lon0 = Math.toRadians(55.53333333333333333333);
     13    private static final double lat0 = Math.toRadians(-21.11666666666666666666);
     14    private static final double x0 = 160000.0;
     15    private static final double y0 = 50000.0;
     16    private static final double k0 = 1;
     17
     18    private static double sinLat0 = Math.sin(lat0);
     19    private static double cosLat0 = Math.cos(lat0);
     20    private static double n1 = Math.sqrt(1 + cosLat0*cosLat0*cosLat0*cosLat0*Ellipsoid.hayford.e2/(1-Ellipsoid.hayford.e2));
     21    private static double phic = Math.asin(sinLat0/n1);
     22    private static double c = Ellipsoid.hayford.latitudeIsometric(phic, 0) - n1*Ellipsoid.hayford.latitudeIsometric(lat0, Ellipsoid.hayford.e);
     23    private static double n2 = k0*Ellipsoid.hayford.a*Math.sqrt(1-Ellipsoid.hayford.e2)/(1-Ellipsoid.hayford.e2*sinLat0*sinLat0);
     24    private static double xs = x0;
     25    private static double ys = y0 - n2*phic;
     26    private static final double epsilon = 1e-11;
     27    private static final double scaleDiff = -32.3241E-6;
     28    private static final double Tx = 789.524;
     29    private static final double Ty = -626.486;
     30    private static final double Tz = -89.904;
     31    private static final double rx = Math.toRadians(0.6006/3600);
     32    private static final double ry = Math.toRadians(76.7946/3600);
     33    private static final double rz = Math.toRadians(-10.5788/3600);
     34    private static final double rx2 = rx*rx;
     35    private static final double ry2 = ry*ry;
     36    private static final double rz2 = rz*rz;
     37
     38    public LatLon eastNorth2latlon(EastNorth p) {
     39        // plan Gauss-Laborde to geographic Piton-des-Neiges
     40        LatLon geo = Geographic(p);
     41
     42        // geographic Piton-des-Neiges/Hayford to geographic WGS84/GRS80
     43        LatLon wgs = PTN2GRS80(geo);
     44        return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon()));
     45    }
     46
     47    /*
     48     * Convert projected coordinates (Gauss-Laborde) to reference ellipsoid Hayforf geographic Piton-des-Neiges
     49     * @return geographic coordinates in radian
     50     */
     51    private LatLon Geographic(EastNorth eastNorth) {
     52        double dxn = (eastNorth.east()-xs)/n2;
     53        double dyn = (eastNorth.north()-ys)/n2;
     54        double lambda = Math.atan(sinh(dxn)/Math.cos(dyn));
     55        double latIso = Ellipsoid.hayford.latitudeIsometric(Math.asin(Math.sin(dyn)/cosh(dxn)), 0);
     56        double lon = lon0 + lambda/n1;
     57        double lat = Ellipsoid.hayford.latitude((latIso-c)/n1, Ellipsoid.hayford.e, 1E-12);
     58        return new LatLon(lat, lon);
     59    }
     60
     61    /**
     62     * Convert geographic Piton-des-Neiges ellipsoid Hayford to geographic WGS84/GRS80
     63     * @param PTN in radian
     64     * @return
     65     */
     66    private LatLon PTN2GRS80(LatLon PTN) {
     67        double lat = PTN.lat();
     68        double lon = PTN.lon();
     69        double N = Ellipsoid.hayford.a / (Math.sqrt(1.0 - Ellipsoid.hayford.e2 * Math.sin(lat) * Math.sin(lat)));
     70        double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
     71        double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
     72        double Z = (N * (1.0 - Ellipsoid.hayford.e2)/* + height */) * Math.sin(lat);
     73
     74        // translation
     75        double coord[] = sevenParametersTransformation(X, Y, Z);
     76
     77        // WGS84 cartesian => WGS84 geographic
     78        return cart2LatLon(coord[0], coord[1], coord[2], Ellipsoid.GRS80);
     79    }
     80
     81    /**
     82     * 7 parameters transformation
     83     * @param coord X, Y, Z in array
     84     * @return transformed X, Y, Z in array
     85     */
     86    private double[] sevenParametersTransformation(double Xa, double Ya, double Za){
     87        double Xb = Tx + Xa*(1+scaleDiff) + Za*ry - Ya*rz;
     88        double Yb = Ty + Ya*(1+scaleDiff) + Xa*rz - Za*rx;
     89        double Zb = Tz + Za*(1+scaleDiff) + Ya*rx - Xa*ry;
     90        return new double[]{Xb, Yb, Zb};
     91    }
     92
     93    public EastNorth latlon2eastNorth(LatLon wgs) {
     94        // translate ellipsoid GRS80 (WGS83) => reference ellipsoid geographic Réunion
     95        LatLon geo = GRS802Hayford(wgs);
     96
     97        // reference ellipsoid geographic => GaussLaborde plan
     98        return GaussLabordeProjection(geo);
     99    }
     100
     101    private LatLon GRS802Hayford(LatLon wgs) {
     102        double lat = Math.toRadians(wgs.lat()); // degree to radian
     103        double lon = Math.toRadians(wgs.lon());
     104        // WGS84 geographic => WGS84 cartesian
     105        double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(lat) * Math.sin(lat)));
     106        double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
     107        double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
     108        double Z = (N * (1.0 - Ellipsoid.GRS80.e2)/* + height */) * Math.sin(lat);
     109        // translation
     110        double coord[] = invSevenParametersTransformation(X, Y, Z);
     111        // Gauss cartesian => Gauss geographic
     112        return Geographic(coord[0], coord[1], coord[2], Ellipsoid.hayford);
     113    }
     114
     115    /**
     116     * 7 parameters inverse transformation
     117     * @param coord X, Y, Z in array
     118     * @return transformed X, Y, Z in array
     119     */
     120    private double[] invSevenParametersTransformation(double Vx, double Vy, double Vz){
     121        Vx = Vx - Tx;
     122        Vy = Vy - Ty;
     123        Vz = Vz - Tz;
     124        double e = 1 + scaleDiff;
     125        double e2 = e*e;
     126        double det = e*(e2 + rx2 + ry2 + rz2);
     127        double Ux = ((e2 + rx2)*Vx + (e*rz + rx*ry)*Vy + (rx*rz - e*ry)*Vz)/det;
     128        double Uy = ((-e*rz + rx*ry)*Vx + (e2 + ry2)*Vy + (e*rx + ry*rz)*Vz)/det;
     129        double Uz = ((e*ry + rx*rz)*Vx + (-e*rx + ry*rz)*Vy + (e2 + rz2)*Vz)/det;
     130        return new double[]{Ux, Uy, Uz};
     131    }
     132
     133    private LatLon Geographic(double X, double Y, double Z, Ellipsoid ell) {
     134        double norm = Math.sqrt(X * X + Y * Y);
     135        double lg = 2.0 * Math.atan(Y / (X + norm));
     136        double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
     137        double delta = 1.0;
     138        while (delta > epsilon) {
     139            double s2 = Math.sin(lt);
     140            s2 *= s2;
     141            double l = Math.atan((Z / norm)
     142                    / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
     143            delta = Math.abs(l - lt);
     144            lt = l;
     145        }
     146        double s2 = Math.sin(lt);
     147        s2 *= s2;
     148        // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
     149        return new LatLon(lt, lg);
     150    }
     151
     152    private EastNorth GaussLabordeProjection(LatLon geo) {
     153        double lambda = n1*(geo.lon()-lon0);
     154        double latIso = c + n1*Ellipsoid.hayford.latitudeIsometric(geo.lat());
     155        double x = xs + n2*Ellipsoid.hayford.latitudeIsometric(Math.asin(Math.sin(lambda)/((Math.exp(latIso) + Math.exp(-latIso))/2)),0);
     156        double y = ys + n2*Math.atan(((Math.exp(latIso) - Math.exp(-latIso))/2)/(Math.cos(lambda)));
     157        return new EastNorth(x, y);
     158    }
     159
     160    /**
     161     * initializes from cartesian coordinates
     162     *
     163     * @param X 1st coordinate in meters
     164     * @param Y 2nd coordinate in meters
     165     * @param Z 3rd coordinate in meters
     166     * @param ell reference ellipsoid
     167     */
     168    private LatLon cart2LatLon(double X, double Y, double Z, Ellipsoid ell) {
     169        double norm = Math.sqrt(X * X + Y * Y);
     170        double lg = 2.0 * Math.atan(Y / (X + norm));
     171        double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
     172        double delta = 1.0;
     173        while (delta > epsilon) {
     174            double s2 = Math.sin(lt);
     175            s2 *= s2;
     176            double l = Math.atan((Z / norm)
     177                    / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
     178            delta = Math.abs(l - lt);
     179            lt = l;
     180        }
     181        double s2 = Math.sin(lt);
     182        s2 *= s2;
     183        // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
     184        return new LatLon(lt, lg);
     185    }
     186
     187    /*
     188     * hyperbolic sine
     189     */
     190    public static final double sinh(double x) {
     191        return (Math.exp(x) - Math.exp(-x))/2;
     192    }
     193    /*
     194     * hyperbolic cosine
     195     */
     196    public static final double cosh(double x) {
     197        return (Math.exp(x) + Math.exp(-x))/2;
     198    }
     199
     200    public String getCacheDirectoryName() {
     201        return this.toString();
     202    }
     203
     204    public Bounds getWorldBoundsLatLon() {
     205        return new Bounds(
     206                new LatLon(-21.5, 55.14),
     207                new LatLon(-20.76, 55.94));
     208    }
     209
     210    public String toCode() {
     211        return "EPSG::3727";
     212    }
     213
     214    @Override public String toString() {
     215        return tr("Gauss-Laborde Réunion 1947");
     216    }
     217
     218}

  • Projection.java

     
    44import org.openstreetmap.josm.data.coor.EastNorth;
    55import org.openstreetmap.josm.data.coor.LatLon;
    66import org.openstreetmap.josm.data.Bounds;
    7 import org.openstreetmap.josm.data.ProjectionBounds;
    87
    98/**
    109 * Classes implementing this are able to convert lat/lon values to
    … …  
    2423    public static Projection[] allProjections = new Projection[]{
    2524        new Epsg4326(),
    2625        new Mercator(),
     26        new LambertEST(),
    2727        new Lambert(),
    28         new LambertEST(),
    2928        new SwissGrid(),
    30         new UTM()
     29        new UTM(),
     30        new UTM_20N_Guadeloupe_Ste_Anne(),
     31        new UTM_20N_Guadeloupe_Fort_Marigot(),
     32        new UTM_20N_Martinique_Fort_Desaix(),
     33        new GaussLaborde_Reunion()
    3134    };
    3235
    3336    /**

  • UTM_20N_France_DOM.java

     
     1// License: GPL. For details, see LICENSE file.
     2package org.openstreetmap.josm.data.projection;
     3
     4/**
     5 * This class is not a direct implementation of Projection. It collects all methods used
     6 * for the projection of the French departements in the Caribbean Sea UTM zone 20N
     7 * but using each time different local geodesic settings for the 7 parameters transformation algorithm.
     8 */
     9import org.openstreetmap.josm.data.coor.EastNorth;
     10import org.openstreetmap.josm.data.coor.LatLon;
     11
     12public class UTM_20N_France_DOM {
     13
     14    /**
     15     * false east in meters (constant)
     16     */
     17    private static final double Xs = 500000.0;
     18    /**
     19     * false north in meters (0 in northern hemisphere, 10000000 in southern
     20     * hemisphere)
     21     */
     22    private static double Ys = 0;
     23    /**
     24     * origin meridian longitude
     25     */
     26    protected double lg0;
     27    /**
     28     * UTM zone (from 1 to 60)
     29     */
     30    int zone = 20;
     31    /**
     32     * whether north or south hemisphere
     33     */
     34    private boolean isNorth = true;
     35    /**
     36     * 7 parameters transformation
     37     */
     38    double tx = 0.0;
     39    double ty = 0.0;
     40    double tz = 0.0;
     41    double rx = 0;
     42    double ry = 0;
     43    double rz = 0;
     44    double scaleDiff = 0;
     45    /**
     46     * precision in iterative schema
     47     */
     48    public static final double epsilon = 1e-11;
     49    public final static double DEG_TO_RAD = Math.PI / 180;
     50    public final static double RAD_TO_DEG = 180 / Math.PI;
     51
     52    public UTM_20N_France_DOM(double[] translation, double[] rotation, double scaleDiff) {
     53        this.tx = translation[0];
     54        this.ty = translation[1];
     55        this.tz = translation[2];
     56        this.rx = rotation[0]/206264.806247096355; // seconds to radian
     57        this.ry = rotation[1]/206264.806247096355;
     58        this.rz = rotation[2]/206264.806247096355;
     59        this.scaleDiff = scaleDiff;
     60    }
     61    public EastNorth latlon2eastNorth(LatLon p) {
     62        // translate ellipsoid GRS80 (WGS83) => reference ellipsoid geographic
     63        LatLon geo = GRS802Hayford(p);
     64
     65        // reference ellipsoid geographic => UTM projection
     66        return MTProjection(geo, Ellipsoid.hayford.a, Ellipsoid.hayford.e);
     67    }
     68
     69    /**
     70     * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to UTM
     71     * geographic, (ellipsoid Hayford)
     72     */
     73    private LatLon GRS802Hayford(LatLon wgs) {
     74        double lat = Math.toRadians(wgs.lat()); // degree to radian
     75        double lon = Math.toRadians(wgs.lon());
     76        // WGS84 geographic => WGS84 cartesian
     77        double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(lat) * Math.sin(lat)));
     78        double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
     79        double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
     80        double Z = (N * (1.0 - Ellipsoid.GRS80.e2)/* + height */) * Math.sin(lat);
     81        // translation
     82        double coord[] = invSevenParametersTransformation(X, Y, Z);
     83        // UTM cartesian => UTM geographic
     84        return Geographic(coord[0], coord[1], coord[2], Ellipsoid.hayford);
     85    }
     86
     87    /**
     88     * initializes from cartesian coordinates
     89     *
     90     * @param X
     91     *            1st coordinate in meters
     92     * @param Y
     93     *            2nd coordinate in meters
     94     * @param Z
     95     *            3rd coordinate in meters
     96     * @param ell
     97     *            reference ellipsoid
     98     */
     99    private LatLon Geographic(double X, double Y, double Z, Ellipsoid ell) {
     100        double norm = Math.sqrt(X * X + Y * Y);
     101        double lg = 2.0 * Math.atan(Y / (X + norm));
     102        double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
     103        double delta = 1.0;
     104        while (delta > epsilon) {
     105            double s2 = Math.sin(lt);
     106            s2 *= s2;
     107            double l = Math.atan((Z / norm)
     108                    / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
     109            delta = Math.abs(l - lt);
     110            lt = l;
     111        }
     112        double s2 = Math.sin(lt);
     113        s2 *= s2;
     114        // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
     115        return new LatLon(lt, lg);
     116    }
     117
     118    /**
     119     * initalizes from geographic coordinates
     120     *
     121     * @param coord geographic coordinates triplet
     122     * @param a reference ellipsoid long axis
     123     * @param e reference ellipsoid excentricity
     124     */
     125    private EastNorth MTProjection(LatLon coord, double a, double e) {
     126        double n = 0.9996 * a;
     127        Ys = (coord.lat() >= 0.0) ? 0.0 : 10000000.0;
     128        double r6d = Math.PI / 30.0;
     129        //zone = (int) Math.floor((coord.lon() + Math.PI) / r6d) + 1;
     130        lg0 = r6d * (zone - 0.5) - Math.PI;
     131        double e2 = e * e;
     132        double e4 = e2 * e2;
     133        double e6 = e4 * e2;
     134        double e8 = e4 * e4;
     135        double C[] = {
     136                1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0 - 175.0*e8/16384.0,
     137                e2/8.0 - e4/96.0 - 9.0*e6/1024.0 - 901.0*e8/184320.0,
     138                13.0*e4/768.0 + 17.0*e6/5120.0 - 311.0*e8/737280.0,
     139                61.0*e6/15360.0 + 899.0*e8/430080.0,
     140                49561.0*e8/41287680.0
     141        };
     142        double s = e * Math.sin(coord.lat());
     143        double l = Math.log(Math.tan(Math.PI/4.0 + coord.lat()/2.0) *
     144                Math.pow((1.0 - s) / (1.0 + s), e/2.0));
     145        double phi = Math.asin(Math.sin(coord.lon() - lg0) /
     146                ((Math.exp(l) + Math.exp(-l)) / 2.0));
     147        double ls = Math.log(Math.tan(Math.PI/4.0 + phi/2.0));
     148        double lambda = Math.atan(((Math.exp(l) - Math.exp(-l)) / 2.0) /
     149                Math.cos(coord.lon() - lg0));
     150
     151        double north = C[0] * lambda;
     152        double east = C[0] * ls;
     153        for(int k = 1; k < 5; k++) {
     154            double r = 2.0 * k * lambda;
     155            double m = 2.0 * k * ls;
     156            double em = Math.exp(m);
     157            double en = Math.exp(-m);
     158            double sr = Math.sin(r)/2.0 * (em + en);
     159            double sm = Math.cos(r)/2.0 * (em - en);
     160            north += C[k] * sr;
     161            east += C[k] * sm;
     162        }
     163        east *= n;
     164        east += Xs;
     165        north *= n;
     166        north += Ys;
     167        return new EastNorth(east, north);
     168    }
     169
     170    public LatLon eastNorth2latlon(EastNorth p) {
     171        MTProjection(p.east(), p.north(), zone, isNorth);
     172        LatLon geo = Geographic(p, Ellipsoid.hayford.a, Ellipsoid.hayford.e, 0.0 /* z */);
     173
     174        // reference ellipsoid geographic => reference ellipsoid cartesian
     175        //Hayford2GRS80(ellipsoid, geo);
     176        double N = Ellipsoid.hayford.a / (Math.sqrt(1.0 - Ellipsoid.hayford.e2 * Math.sin(geo.lat()) * Math.sin(geo.lat())));
     177        double X = (N /*+ h*/) * Math.cos(geo.lat()) * Math.cos(geo.lon());
     178        double Y = (N /*+ h*/) * Math.cos(geo.lat()) * Math.sin(geo.lon());
     179        double Z = (N * (1.0-Ellipsoid.hayford.e2) /*+ h*/) * Math.sin(geo.lat());
     180
     181        // translation
     182        double coord[] = sevenParametersTransformation(X, Y, Z);
     183
     184        // WGS84 cartesian => WGS84 geographic
     185        LatLon wgs = cart2LatLon(coord[0], coord[1], coord[2], Ellipsoid.GRS80);
     186        return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon()));
     187    }
     188
     189    /**
     190     * initializes new projection coordinates (in north hemisphere)
     191     *
     192     * @param east east from origin in meters
     193     * @param north north from origin in meters
     194     * @param zone zone number (from 1 to 60)
     195     * @param isNorth true in north hemisphere, false in south hemisphere
     196     */
     197    private void MTProjection(double east, double north, int zone, boolean isNorth) {
     198        Ys = isNorth ? 0.0 : 10000000.0;
     199        double r6d = Math.PI / 30.0;
     200        lg0 = r6d * (zone - 0.5) - Math.PI;
     201    }
     202
     203    public double scaleFactor() {
     204        return 1/Math.PI/2;
     205    }
     206
     207    /**
     208     * initalizes from projected coordinates (Mercator transverse projection)
     209     *
     210     * @param coord projected coordinates pair
     211     * @param e reference ellipsoid excentricity
     212     * @param a reference ellipsoid long axis
     213     * @param z altitude in meters
     214     */
     215    private LatLon Geographic(EastNorth coord, double a, double e, double z) {
     216        double n = 0.9996 * a;
     217        double e2 = e * e;
     218        double e4 = e2 * e2;
     219        double e6 = e4 * e2;
     220        double e8 = e4 * e4;
     221        double C[] = {
     222                1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0 - 175.0*e8/16384.0,
     223                e2/8.0 + e4/48.0 + 7.0*e6/2048.0 + e8/61440.0,
     224                e4/768.0 + 3.0*e6/1280.0 + 559.0*e8/368640.0,
     225                17.0*e6/30720.0 + 283.0*e8/430080.0,
     226                4397.0*e8/41287680.0
     227        };
     228        double l = (coord.north() - Ys) / (n * C[0]);
     229        double ls = (coord.east() - Xs) / (n * C[0]);
     230        double l0 = l;
     231        double ls0 = ls;
     232        for(int k = 1; k < 5; k++) {
     233            double r = 2.0 * k * l0;
     234            double m = 2.0 * k * ls0;
     235            double em = Math.exp(m);
     236            double en = Math.exp(-m);
     237            double sr = Math.sin(r)/2.0 * (em + en);
     238            double sm = Math.cos(r)/2.0 * (em - en);
     239            l -= C[k] * sr;
     240            ls -= C[k] * sm;
     241        }
     242        double lon = lg0 + Math.atan(((Math.exp(ls) - Math.exp(-ls)) / 2.0) /
     243                Math.cos(l));
     244        double phi = Math.asin(Math.sin(l) /
     245                ((Math.exp(ls) + Math.exp(-ls)) / 2.0));
     246        l = Math.log(Math.tan(Math.PI/4.0 + phi/2.0));
     247        double lat = 2.0 * Math.atan(Math.exp(l)) - Math.PI / 2.0;
     248        double lt0;
     249        do {
     250            lt0 = lat;
     251            double s = e * Math.sin(lat);
     252            lat = 2.0 * Math.atan(Math.pow((1.0 + s) / (1.0 - s), e/2.0) *
     253                    Math.exp(l)) - Math.PI / 2.0;
     254        }
     255        while(Math.abs(lat-lt0) >= epsilon);
     256        //h = z;
     257
     258        return new LatLon(lat, lon);
     259    }
     260
     261    /**
     262     * initializes from cartesian coordinates
     263     *
     264     * @param X 1st coordinate in meters
     265     * @param Y 2nd coordinate in meters
     266     * @param Z 3rd coordinate in meters
     267     * @param ell reference ellipsoid
     268     */
     269    private LatLon cart2LatLon(double X, double Y, double Z, Ellipsoid ell) {
     270        double norm = Math.sqrt(X * X + Y * Y);
     271        double lg = 2.0 * Math.atan(Y / (X + norm));
     272        double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
     273        double delta = 1.0;
     274        while (delta > epsilon) {
     275            double s2 = Math.sin(lt);
     276            s2 *= s2;
     277            double l = Math.atan((Z / norm)
     278                    / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
     279            delta = Math.abs(l - lt);
     280            lt = l;
     281        }
     282        double s2 = Math.sin(lt);
     283        s2 *= s2;
     284        // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
     285        return new LatLon(lt, lg);
     286    }
     287
     288    /**
     289     * 7 parameters transformation
     290     * @param coord X, Y, Z in array
     291     * @return transformed X, Y, Z in array
     292     */
     293    private double[] sevenParametersTransformation(double Xa, double Ya, double Za){
     294        double Xb = tx + Xa*(1+scaleDiff) + Za*ry - Ya*rz;
     295        double Yb = ty + Ya*(1+scaleDiff) + Xa*rz - Za*rx;
     296        double Zb = tz + Za*(1+scaleDiff) + Ya*rx - Xa*ry;
     297        return new double[]{Xb, Yb, Zb};
     298    }
     299
     300    /**
     301     * 7 parameters inverse transformation
     302     * @param coord X, Y, Z in array
     303     * @return transformed X, Y, Z in array
     304     */
     305    private double[] invSevenParametersTransformation(double Xa, double Ya, double Za){
     306        double Xb = (1-scaleDiff)*(-tx + Xa + ((-tz+Za)*(-ry) - (-ty+Ya)*(-rz)));
     307        double Yb = (1-scaleDiff)*(-ty + Ya + ((-tx+Xa)*(-rz) - (-tz+Za)*(-rx)));
     308        double Zb = (1-scaleDiff)*(-tz + Za + ((-ty+Ya)*(-rx) - (-tx+Xa)*(-ry)));
     309        return new double[]{Xb, Yb, Zb};
     310    }
     311
     312}

  • UTM_20N_Guadeloupe_Fort_Marigot.java

     
     1// License: GPL. For details, see LICENSE file.
     2package org.openstreetmap.josm.data.projection;
     3
     4import static org.openstreetmap.josm.tools.I18n.tr;
     5
     6import org.openstreetmap.josm.data.Bounds;
     7import org.openstreetmap.josm.data.coor.LatLon;
     8
     9/*
     10 * Local geodisic system with UTM zone 20N projection.
     11 * Apply to Guadeloupe, France - St Martin and St Barthelemy islands
     12 */
     13public class UTM_20N_Guadeloupe_Fort_Marigot extends UTM_20N_France_DOM implements Projection {
     14    public UTM_20N_Guadeloupe_Fort_Marigot() {
     15        super(new double[]{136.596, 248.148, -429.789},
     16                new double[]{0, 0, 0},
     17                0);
     18    }
     19
     20    public String getCacheDirectoryName() {
     21        return this.toString();
     22    }
     23
     24    public Bounds getWorldBoundsLatLon() {
     25        return new Bounds(
     26                new LatLon(17.6,-63.25),
     27                new LatLon(18.5,-62.5));
     28    }
     29
     30    public String toCode() {
     31        return "EPSG::2969";
     32    }
     33
     34    @Override public String toString() {
     35        return tr("UTM20N Guadeloupe Fort-Marigot 1949");
     36    }
     37
     38}

  • UTM_20N_Guadeloupe_Ste_Anne.java

     
     1// License: GPL. For details, see LICENSE file.
     2package org.openstreetmap.josm.data.projection;
     3
     4import static org.openstreetmap.josm.tools.I18n.tr;
     5
     6import org.openstreetmap.josm.data.Bounds;
     7import org.openstreetmap.josm.data.coor.LatLon;
     8
     9/*
     10 * Local geodisic system with UTM zone 20N projection.
     11 * Apply to Guadeloupe, France - Grande-Terre and surrounding islands.
     12 */
     13public class UTM_20N_Guadeloupe_Ste_Anne extends UTM_20N_France_DOM implements Projection {
     14    public UTM_20N_Guadeloupe_Ste_Anne() {
     15        super (new double[]{-472.29, -5.63, -304.12},
     16                new double[]{0.4362, -0.8374, 0.2563},
     17                1.8984E-6);
     18    }
     19
     20    public String getCacheDirectoryName() {
     21        return this.toString();
     22    }
     23
     24    public Bounds getWorldBoundsLatLon() {
     25        return new Bounds(
     26                new LatLon(15.8,-61.8),
     27                new LatLon(16.6,-60.9));
     28    }
     29
     30    public String toCode() {
     31        return "EPSG::2970";
     32    }
     33
     34    @Override public String toString() {
     35        return tr("UTM20N Guadeloupe Ste-Anne 1948");
     36    }
     37
     38}

  • UTM_20N_Martinique_Fort_Desaix.java

     
     1// License: GPL. For details, see LICENSE file.
     2package org.openstreetmap.josm.data.projection;
     3
     4import static org.openstreetmap.josm.tools.I18n.tr;
     5
     6import org.openstreetmap.josm.data.Bounds;
     7import org.openstreetmap.josm.data.coor.LatLon;
     8
     9/*
     10 * Local geodisic system with UTM zone 20N projection.
     11 * Apply to Martinique, France and surrounding islands
     12 */
     13public class UTM_20N_Martinique_Fort_Desaix extends UTM_20N_France_DOM implements Projection {
     14    public UTM_20N_Martinique_Fort_Desaix() {
     15        super(new double[]{126.926, 547.939, 130.409},
     16                new double[]{-2.78670, 5.16124,  -0.85844},
     17                13.82265E-6);
     18    }
     19
     20    public String getCacheDirectoryName() {
     21        return this.toString();
     22    }
     23
     24    public Bounds getWorldBoundsLatLon() {
     25        return new Bounds(
     26                new LatLon(14.25,-61.25),
     27                new LatLon(15.025,-60.725));
     28    }
     29
     30    public String toCode() {
     31        return "EPSG::2973";
     32    }
     33
     34    @Override public String toString() {
     35        return tr("UTM20N Martinique Fort Desaix 1952");
     36    }
     37}