Ticket #12181: NetherlandsRD.patch

File NetherlandsRD.patch, 7.9 KB (added by vholten, 10 years ago)

Patch

  • data/projection/epsg

     
    110110<27563> +proj=lcc +lat_0=44.1 +lat_1=43d11'57.449" +lat_2=44d59'45.938" +lon_0=2d20'14.025" +x_0=600000 +y_0=200000 +a=6378249.2 +b=6356515 +nadgrids=ntf_r93_b.gsb +units=m +bounds=-4.416666666666665,41.49,9.18,46.95  <>
    111111# Lambert 4 Zones France (Corsica)
    112112<27564> +proj=lcc +lat_0=42.165 +lat_1=41d33'37.396" +lat_2=42d46'3.588" +lon_0=2d20'14.025" +x_0=234.358 +y_0=185861.369 +a=6378249.2 +b=6356515 +nadgrids=ntf_r93_b.gsb +units=m +bounds=-4.416666666666665,41.49,9.18,44.267667  <>
     113# Rijksdriehoekscoördinaten (Netherlands)
     114<28992> +proj=sterea +lat_0=52d9'22.178" +lon_0=5d23'15.5" +k_0=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957,0.343988,-1.877402,4.0725 +units=m +bounds=3.1,50.56,7.6,53.63  <>
    113115# Belgian Lambert 1972
    114116<31370> +proj=lcc +lat_0=90 +lat_1=49d50'0.00204" +lat_2=51d10'0.00204" +lon_0=4d22'2.952" +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-99.06,53.32,-112.49,0.419,-0.83,1.885,-1.0 +units=m +bounds=2.54,49.51,6.4,51.5  <>
    115117# Gauß-Krüger Zone 2

  • src/org/openstreetmap/josm/data/projection/Projections.java

     
    2626import org.openstreetmap.josm.data.projection.datum.NTV2GridShiftFileWrapper;
    2727import org.openstreetmap.josm.data.projection.datum.WGS84Datum;
    2828import org.openstreetmap.josm.data.projection.proj.ClassProjFactory;
     29import org.openstreetmap.josm.data.projection.proj.DoubleStereographic;
    2930import org.openstreetmap.josm.data.projection.proj.LambertConformalConic;
    3031import org.openstreetmap.josm.data.projection.proj.LonLat;
    3132import org.openstreetmap.josm.data.projection.proj.Mercator;
    … …  
    7677        registerBaseProjection("lcc", LambertConformalConic.class, "core");
    7778        registerBaseProjection("somerc", SwissObliqueMercator.class, "core");
    7879        registerBaseProjection("tmerc", TransverseMercator.class, "core");
     80        registerBaseProjection("sterea", DoubleStereographic.class, "core");
    7981
    8082        ellipsoids.put("airy", Ellipsoid.Airy);
    8183        ellipsoids.put("mod_airy", Ellipsoid.AiryMod);

  • src/org/openstreetmap/josm/data/projection/proj/DoubleStereographic.java

     
     1// License: GPL. For details, see LICENSE file.
     2package org.openstreetmap.josm.data.projection.proj;
     3
     4import static java.lang.Math.PI;
     5import static java.lang.Math.abs;
     6import static java.lang.Math.asin;
     7import static java.lang.Math.atan;
     8import static java.lang.Math.cos;
     9import static java.lang.Math.exp;
     10import static java.lang.Math.log;
     11import static java.lang.Math.pow;
     12import static java.lang.Math.sin;
     13import static java.lang.Math.sqrt;
     14import static java.lang.Math.tan;
     15import static java.lang.Math.toRadians;
     16import static org.openstreetmap.josm.tools.I18n.tr;
     17
     18import org.openstreetmap.josm.data.projection.Ellipsoid;
     19import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
     20
     21/**
     22 * Implementation of the stereographic double projection,
     23 * also known as Oblique Stereographic and the Schreiber double stereographic projection.
     24 *
     25 * @author vholten
     26 *
     27 * Source: IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2,
     28 * Sec. 1.3.7.1 Oblique and Equatorial Stereographic, http://www.epsg.org/GuidanceNotes
     29 */
     30
     31public class DoubleStereographic implements Proj {
     32
     33    private Ellipsoid ellps;
     34    private double e;
     35    private double n;
     36    private double c;
     37    private double chi0;
     38    private double R;
     39
     40    private static final double EPSILON = 1e-12;
     41
     42    @Override
     43    public String getName() {
     44        return tr("Double Stereographic");
     45    }
     46
     47    @Override
     48    public String getProj4Id() {
     49        return "sterea";
     50    }
     51
     52    @Override
     53    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
     54        if (params.lat0 == null)
     55            throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
     56        ellps = params.ellps;
     57        e = ellps.e;
     58        initialize(params.lat0);
     59    }
     60
     61    private void initialize(double lat_0) {
     62        double phi0 = toRadians(lat_0);
     63        double e2 = ellps.e2;
     64        R = sqrt(1-e2) / (1 - e2*pow(sin(phi0), 2));
     65        n = sqrt(1 + ellps.eb2 * pow(cos(phi0),4));
     66        double S1 = (1 + sin(phi0)) / (1 - sin(phi0));
     67        double S2 = (1 - e * sin(phi0)) / (1 + e * sin(phi0));
     68        double w1 = pow(S1 * pow(S2, e), n);
     69        double sinchi00 = (w1 - 1) / (w1 + 1);
     70        c = (n + sin(phi0)) * (1 - sinchi00) / ((n - sin(phi0)) * (1 + sinchi00));
     71        double w2 = c * w1;
     72        chi0 = asin((w2 - 1) / (w2 + 1));
     73    }
     74
     75    @Override
     76    public double[] project(double phi, double lambda) {
     77        double Lambda = n * lambda;
     78        double Sa = (1 + sin(phi)) / (1 - sin(phi));
     79        double Sb = (1 - e * sin(phi)) / (1 + e * sin(phi));
     80        double w = c * pow(Sa * pow(Sb, e), n);
     81        double chi = asin((w - 1) / (w + 1));
     82        double B = 1 + sin(chi) * sin(chi0) + cos(chi) * cos(chi0) * cos(Lambda);
     83        double x = 2 * R * cos(chi) * sin(Lambda) / B;
     84        double y = 2 * R * (sin(chi) * cos(chi0) - cos(chi) * sin(chi0) * cos(Lambda)) / B;
     85        return new double[] {x, y};
     86    }
     87
     88    @Override
     89    public double[] invproject(double x, double y) {
     90        double e2 = ellps.e2;
     91        double g = 2 * R * tan(PI/4 - chi0/2);
     92        double h = 4 * R * tan(chi0) + g;
     93        double i = atan(x/(h + y));
     94        double j = atan(x/(g - y)) - i;
     95        double chi = chi0 + 2 * atan((y - x * tan(j/2)) / (2 * R));
     96        double Lambda = j + 2*i;
     97        double lambda = Lambda / n;
     98        double psi = 0.5 * log((1 + sin(chi)) / (c*(1 - sin(chi)))) / n;
     99        double phiprev = -1000;
     100        int iteration = 0;
     101        double phi = 2 * atan(exp(psi)) - PI/2;
     102        while (abs(phi - phiprev) > EPSILON) {
     103            if (++iteration > 10)
     104                throw new RuntimeException("Too many iterations");
     105            phiprev = phi;
     106            double psii = log(tan(phi/2 + PI/4) * pow((1 - e * sin(phi)) / (1 + e * sin(phi)), e/2));
     107            phi = phi - (psii - psi) * cos(phi) * (1 - e2 * pow(sin(phi), 2)) / (1 - e2);
     108        }
     109        return new double[] {phi, lambda};
     110    }
     111}

  • src/org/openstreetmap/josm/gui/preferences/projection/ProjectionPreference.java 

    Property changes on: src/org/openstreetmap/josm/data/projection/proj/DoubleStereographic.java
___________________________________________________________________
Added: svn:executable
## -0,0 +1 ##
+*
\ No newline at end of property
     
    193193        registerProjectionChoice(tr("LKS-92 (Latvia TM)"), "core:tmerclv", 3059);                   // LV
    194194
    195195        /**
     196         * Netherlands RD projection
     197         *
     198         * @author vholten
     199         */
     200        registerProjectionChoice(tr("Rijksdriehoekscoördinaten (Netherlands)"), "core:dutchrd", 28992); // NL
     201
     202        /**
    196203         * PUWG 1992 and 2000 are the official cordinate systems in Poland.
    197204         *
    198205         * They use the same math as UTM only with different constants.