Changeset 1929 in josm for trunk/src

Timestamp:

2009-08-08T11:35:46+02:00 (15 years ago)

Author:

Gubaer

Message:

applied #3214: patch by pieren: new projections for the French cadastre

Location:

trunk/src/org/openstreetmap/josm/data/projection

Files:

• Ellipsoid.java (modified) (1 diff)
• GaussLaborde_Reunion.java (added)
• Projection.java (modified) (2 diffs)
• UTM_20N_France_DOM.java (added)
• UTM_20N_Guadeloupe_Fort_Marigot.java (added)
• UTM_20N_Guadeloupe_Ste_Anne.java (added)
• UTM_20N_Martinique_Fort_Desaix.java (added)

trunk/src/org/openstreetmap/josm/data/projection/Ellipsoid.java

@@ -12 +12 @@
  */
 class Ellipsoid {
-  /**
-   * Clarke's elipsoid (NTF system)
-   */
-  public static final Ellipsoid clarke = new Ellipsoid(6378249.2, 6356515.0);
-  /**
-   * Hayford's ellipsoid (ED50 system)
-   */
-  public static final Ellipsoid hayford =
-    new Ellipsoid(6378388.0, 6356911.9461);
-  /**
-   * WGS84 elipsoid
-   */
-  public static final Ellipsoid GRS80 = new Ellipsoid(6378137.0, 6356752.314);
+    /**
+     * Clarke's ellipsoid (NTF system)
+     */
+    public static final Ellipsoid clarke = new Ellipsoid(6378249.2, 6356515.0);
+    /**
+     * Hayford's ellipsoid (ED50 system)
+     */
+    public static final Ellipsoid hayford =
+        new Ellipsoid(6378388.0, 6356911.9461);
+    /**
+     * WGS84 ellipsoid
+     */
+    public static final Ellipsoid GRS80 = new Ellipsoid(6378137.0, 6356752.314);
 
-  /**
-   * half long axis
-   */
-  public final double a;
-  /**
-   * half short axis
-   */
-  public final double b;
-  /**
-   * first excentricity
-   */
-  public final double e;
-  /**
-   * first excentricity squared
-   */
-  public final double e2;
+    /**
+     * half long axis
+     */
+    public final double a;
+    /**
+     * half short axis
+     */
+    public final double b;
+    /**
+     * first eccentricity
+     */
+    public final double e;
+    /**
+     * first eccentricity squared
+     */
+    public final double e2;
 
-  /**
-   * create a new ellipsoid and precompute its parameters
-   *
-   * @param a ellipsoid long axis (in meters)
-   * @param b ellipsoid short axis (in meters)
-   */
-  public Ellipsoid(double a, double b) {
-    this.a = a;
-    this.b = b;
-    e2 = (a*a - b*b) / (a*a);
-    e = Math.sqrt(e2);
-  }
+    /**
+     * square of the second eccentricity
+     */
+    public final double eb2;
+
+    /**
+     * create a new ellipsoid and precompute its parameters
+     *
+     * @param a ellipsoid long axis (in meters)
+     * @param b ellipsoid short axis (in meters)
+     */
+    public Ellipsoid(double a, double b) {
+        this.a = a;
+        this.b = b;
+        e2 = (a*a - b*b) / (a*a);
+        e = Math.sqrt(e2);
+        eb2 = e2 / (1.0 - e2);
+    }
+
+    /**
+     * Returns the <i>radius of curvature in the prime vertical</i>
+     * for this reference ellipsoid at the specified latitude.
+     *
+     * @param phi The local latitude (radians).
+     * @return The radius of curvature in the prime vertical (meters).
+     */
+    public double verticalRadiusOfCurvature(final double phi) {
+        return a / Math.sqrt(1.0 - (e2 * sqr(Math.sin(phi))));
+    }
+
+    private static double sqr(final double x) {
+        return x * x;
+    }
+
+    /**
+     *  Returns the meridional arc, the true meridional distance on the
+     * ellipsoid from the equator to the specified latitude, in meters.
+     *
+     * @param phi   The local latitude (in radians).
+     * @return  The meridional arc (in meters).
+     */
+    public double meridionalArc(final double phi) {
+        final double sin2Phi = Math.sin(2.0 * phi);
+        final double sin4Phi = Math.sin(4.0 * phi);
+        final double sin6Phi = Math.sin(6.0 * phi);
+        final double sin8Phi = Math.sin(8.0 * phi);
+        // TODO . calculate 'f'
+        //double f = 1.0 / 298.257222101; // GRS80
+        double f = 1.0 / 298.257223563; // WGS84
+        final double n = f / (2.0 - f);
+        final double n2 = n * n;
+        final double n3 = n2 * n;
+        final double n4 = n3 * n;
+        final double n5 = n4 * n;
+        final double n1n2 = n - n2;
+        final double n2n3 = n2 - n3;
+        final double n3n4 = n3 - n4;
+        final double n4n5 = n4 - n5;
+        final double ap = a * (1.0 - n + (5.0 / 4.0) * (n2n3) + (81.0 / 64.0) * (n4n5));
+        final double bp = (3.0 / 2.0) * a * (n1n2 + (7.0 / 8.0) * (n3n4) + (55.0 / 64.0) * n5);
+        final double cp = (15.0 / 16.0) * a * (n2n3 + (3.0 / 4.0) * (n4n5));
+        final double dp = (35.0 / 48.0) * a * (n3n4 + (11.0 / 16.0) * n5);
+        final double ep = (315.0 / 512.0) * a * (n4n5);
+        return ap * phi - bp * sin2Phi + cp * sin4Phi - dp * sin6Phi + ep * sin8Phi;
+    }
+
+    /**
+     *  Returns the <i>radius of curvature in the meridian<i>
+     *  for this reference ellipsoid at the specified latitude.
+     *
+     * @param phi The local latitude (in radians).
+     * @return  The radius of curvature in the meridian (in meters).
+     */
+    public double meridionalRadiusOfCurvature(final double phi) {
+        return verticalRadiusOfCurvature(phi)
+        / (1.0 + eb2 * sqr(Math.cos(phi)));
+    }
+
+    /**
+     * Returns isometric latitude of phi on given first eccentricity (e)
+     * @param phi The local latitude (radians).
+     * @return isometric latitude of phi on first eccentricity (e)
+     */
+    public double latitudeIsometric(double phi, double e) {
+        double v1 = 1-e*Math.sin(phi);
+        double v2 = 1+e*Math.sin(phi);
+        return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2,e/2));
+    }
+
+    /**
+     * Returns isometric latitude of phi on first eccentricity (e)
+     * @param phi The local latitude (radians).
+     * @return isometric latitude of phi on first eccentricity (e)
+     */
+    public double latitudeIsometric(double phi) {
+        double v1 = 1-e*Math.sin(phi);
+        double v2 = 1+e*Math.sin(phi);
+        return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2,e/2));
+    }
+
+    /*
+     * Returns geographic latitude of isometric latitude of first eccentricity (e)
+     * and epsilon precision
+     */
+    public double latitude(double latIso, double e, double epsilon) {
+        double lat0 = 2*Math.atan(Math.exp(latIso))-Math.PI/2;
+        double lati = lat0;
+        double lati1 = 1.0; // random value to start the iterative processus
+        while(Math.abs(lati1-lati)>=epsilon) {
+            lati = lati1;
+            double v1 = 1+e*Math.sin(lati);
+            double v2 = 1-e*Math.sin(lati);
+            lati1 = 2*Math.atan(Math.pow(v1/v2,e/2)*Math.exp(latIso))-Math.PI/2;
+        }
+        return lati1;
+    }
 }
 
 
+

trunk/src/org/openstreetmap/josm/data/projection/Projection.java

@@ -5 +5 @@
 import org.openstreetmap.josm.data.coor.LatLon;
 import org.openstreetmap.josm.data.Bounds;
-import org.openstreetmap.josm.data.ProjectionBounds;
 
 /**
@@ -25 +24 @@
         new Epsg4326(),
         new Mercator(),
+        new LambertEST(),
         new Lambert(),
-        new LambertEST(),
         new SwissGrid(),
-        new UTM()
+        new UTM(),
+        new UTM_20N_Guadeloupe_Ste_Anne(),
+        new UTM_20N_Guadeloupe_Fort_Marigot(),
+        new UTM_20N_Martinique_Fort_Desaix(),
+        new GaussLaborde_Reunion()
     };