Ticket #3737: projection.patch

File projection.patch, 10.6 KB (added by pieren, 15 years ago)
  • Ellipsoid.java

     
    1010/**
    1111 * the reference ellipsoids
    1212 */
    13 class Ellipsoid {
     13public class Ellipsoid {
    1414    /**
    1515     * Clarke's ellipsoid (NTF system)
    1616     */
     
    2121    public static final Ellipsoid hayford =
    2222        new Ellipsoid(6378388.0, 6356911.9461);
    2323    /**
     24     * GRS80 ellipsoid
     25     */
     26    public static final Ellipsoid GRS80 = new Ellipsoid(6378137.0, 6356752.3141);
     27
     28    /**
    2429     * WGS84 ellipsoid
    2530     */
    26     public static final Ellipsoid GRS80 = new Ellipsoid(6378137.0, 6356752.314);
     31    public static final Ellipsoid WGS84 = new Ellipsoid(6378137.0, 6356752.3142);
    2732
    2833    /**
    2934     * half long axis
  • Lambert.java

     
    140140    }
    141141
    142142    @Override public String toString() {
    143         return tr("Lambert Zone (France)");
     143        return tr("Lambert 4 Zones (France)");
    144144    }
    145145
    146146    public String toCode() {
    147         return "EPSG:"+(27571+currentZone);
     147        return "EPSG:"+(27561+currentZone);
    148148    }
    149149
    150150    public String getCacheDirectoryName() {
  • LambertCC9Zones.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 javax.swing.JOptionPane;
     7
     8import org.openstreetmap.josm.Main;
     9import org.openstreetmap.josm.data.Bounds;
     10import org.openstreetmap.josm.data.coor.EastNorth;
     11import org.openstreetmap.josm.data.coor.LatLon;
     12import org.openstreetmap.josm.data.projection.Projection;
     13import org.openstreetmap.josm.data.projection.Ellipsoid;
     14
     15/**
     16 * This class implements the Lambert Conic Conform 9 Zones projection as specified by the IGN
     17 * in this document http://professionnels.ign.fr/DISPLAY/000/526/700/5267002/transformation.pdf
     18 * @author Pieren
     19 *
     20 */
     21public class LambertCC9Zones implements Projection {
     22
     23    /**
     24     * Lambert 9 zones projection exponents
     25     */
     26    public static final double n[] = { 0.6691500006885269, 0.682018118346418, 0.6946784863203991, 0.7071272481559119,
     27        0.7193606118567315, 0.7313748510399917, 0.7431663060711892, 0.7547313851789208, 0.7660665655489937};
     28
     29    /**
     30     * Lambert 9 zones projection constants
     31     */
     32    public static final double c[] = { 1.215363305807804E7, 1.2050261119223533E7, 1.195716926884592E7, 1.18737533925172E7,
     33        1.1799460698022118E7, 1.17337838820243E7, 1.16762559948139E7, 1.1626445901183508E7, 1.1583954251630554E7};
     34
     35    /**
     36     * Lambert 9 zones false east
     37     */
     38    public static final double Xs = 1700000;
     39
     40    /**
     41     * Lambert 9 zones false north
     42     */
     43    public static final double Ys[] = { 8293467.503439436, 9049604.665107645, 9814691.693461388, 1.0588107871787189E7,
     44        1.1369285637569271E7, 1.2157704903382052E7, 1.2952888086405803E7, 1.3754395745267643E7, 1.4561822739114787E7};
     45
     46    /**
     47     * Lambert I, II, III, and IV longitudinal offset to Greenwich meridian
     48     */
     49    public static final double lg0 = 0.04079234433198; // 2deg20'14.025"
     50
     51    /**
     52     * precision in iterative schema
     53     */
     54
     55    public static final double epsilon = 1e-12;
     56
     57    /**
     58     * France is divided in 9 Lambert projection zones, CC42 to CC50.
     59     */
     60    public static final double cMaxLatZones = Math.toRadians(51.1);
     61
     62    public static final double cMinLatZones = Math.toRadians(41.0);
     63
     64    public static final double cMinLonZones = Math.toRadians(-5.0);
     65
     66    public static final double cMaxLonZones = Math.toRadians(10.2);
     67
     68    public static final double lambda0 = Math.toRadians(3);
     69    public static final double e = Ellipsoid.GRS80.e; // but in doc=0.08181919112
     70    public static final double e2 =Ellipsoid.GRS80.e2;
     71    public static final double a = Ellipsoid.GRS80.a;
     72    /**
     73     *  Because josm cannot work correctly if two zones are displayed, we allow some overlapping
     74     */
     75    public static final double cMaxOverlappingZones = Math.toRadians(1.5);
     76
     77    public static int layoutZone = -1;
     78
     79    private double L(double phi, double e) {
     80        double sinphi = Math.sin(phi);
     81        return (0.5*Math.log((1+sinphi)/(1-sinphi))) - e/2*Math.log((1+e*sinphi)/(1-e*sinphi));
     82    }
     83
     84    /**
     85     * @param p  WGS84 lat/lon (ellipsoid GRS80) (in degree)
     86     * @return eastnorth projection in Lambert Zone (ellipsoid Clark)
     87     */
     88    public EastNorth latlon2eastNorth(LatLon p) {
     89        double lt = Math.toRadians(p.lat());
     90        double lg = Math.toRadians(p.lon());
     91        // check if longitude and latitude are inside the French Lambert zones and seek a zone number
     92        // if it is not already defined in layoutZone
     93        int possibleZone = 0;
     94        boolean outOfLambertZones = false;
     95        if (lt >= cMinLatZones && lt <= cMaxLatZones && lg >= cMinLonZones && lg <= cMaxLonZones) {
     96            /* with Lambert CC9 zones, each latitude is present in two zones. If the layout
     97               zone is not already set, we choose arbitrary the first possible zone */
     98            possibleZone = (int)p.lat()-42;
     99            if (possibleZone > 8) {
     100                possibleZone = 8;
     101            }
     102            if (possibleZone < 0) {
     103                possibleZone = 0;
     104            }
     105        } else {
     106            outOfLambertZones = true; // possible when MAX_LAT is used
     107        }
     108        if (!outOfLambertZones) {
     109            if (layoutZone == -1) {
     110                if (layoutZone != possibleZone) {
     111                    System.out.println("change Lambert zone from "+layoutZone+" to "+possibleZone);
     112                }
     113                layoutZone = possibleZone;
     114            } else if (Math.abs(layoutZone - possibleZone) > 1) {
     115                if (farawayFromLambertZoneFrance(lt, lg)) {
     116                    JOptionPane.showMessageDialog(Main.parent,
     117                            tr("IMPORTANT : data positioned far away from\n"
     118                                    + "the current Lambert zone limits.\n"
     119                                    + "Do not upload any data after this message.\n"
     120                                    + "Undo your last action, save your work\n"
     121                                    + "and start a new layer on the new zone."),
     122                                    tr("Warning"),
     123                                    JOptionPane.WARNING_MESSAGE);
     124                    layoutZone = -1;
     125                } else {
     126                    System.out.println("temporarily extend Lambert zone " + layoutZone + " projection at lat,lon:"
     127                            + lt + "," + lg);
     128                }
     129            }
     130        }
     131        if (layoutZone == -1)
     132            return ConicProjection(lt, lg, possibleZone);
     133        return ConicProjection(lt, lg, layoutZone);
     134    }
     135
     136    /**
     137     *
     138     * @param lat latitude in grad
     139     * @param lon longitude in grad
     140     * @param nz Lambert CC zone number (from 1 to 9) - 1 !
     141     * @return EastNorth projected coordinates in meter
     142     */
     143    private EastNorth ConicProjection(double lat, double lon, int nz) {
     144        double R = c[nz]*Math.exp(-n[nz]*L(lat,e));
     145        double gamma = n[nz]*(lon-lambda0);
     146        double X = Xs + R*Math.sin(gamma);
     147        double Y = Ys[nz] + -R*Math.cos(gamma);
     148        return new EastNorth(X, Y);
     149    }
     150
     151    public LatLon eastNorth2latlon(EastNorth p) {
     152        layoutZone  = north2ZoneNumber(p.north());
     153        return Geographic(p, layoutZone);
     154    }
     155
     156    private LatLon Geographic(EastNorth ea, int nz) {
     157        double R = Math.sqrt(Math.pow(ea.getX()-Xs,2)+Math.pow(ea.getY()-Ys[nz], 2));
     158        double gamma = Math.atan((ea.getX()-Xs)/(Ys[nz]-ea.getY()));
     159        double lon = lambda0+gamma/n[nz];
     160        double latIso = (-1/n[nz])*Math.log(Math.abs(R/c[nz]));
     161        double lat = Ellipsoid.GRS80.latitude(latIso, e, epsilon);
     162        return new LatLon(Math.toDegrees(lat), Math.toDegrees(lon));
     163    }
     164
     165    @Override public String toString() {
     166        return tr("Lambert CC9 Zone (France)");
     167    }
     168
     169    public static int north2ZoneNumber(double north) {
     170        int nz = (int)(north /1000000) - 1;
     171        if (nz < 0) return 0;
     172        else if (nz > 8) return 8;
     173        else return nz;
     174    }
     175
     176    public static boolean isInL9CCZones(LatLon p) {
     177        double lt = Math.toRadians(p.lat());
     178        double lg = Math.toRadians(p.lon());
     179        if (lg >= cMinLonZones && lg <= cMaxLonZones && lt >= cMinLatZones && lt <= cMaxLatZones)
     180            return true;
     181        return false;
     182    }
     183
     184    public String toCode() {
     185        if (layoutZone == -1)
     186            return "EPSG:"+(3942);
     187        return "EPSG:"+(3942+layoutZone); //CC42 is EPSG:3842 (up to EPSG:3950 for CC50)
     188    }
     189
     190    public String getCacheDirectoryName() {
     191        return "lambert";
     192    }
     193
     194    private boolean farawayFromLambertZoneFrance(double lat, double lon) {
     195        if (lat < (cMinLatZones - cMaxOverlappingZones) || (lat > (cMaxLatZones + cMaxOverlappingZones))
     196                || (lon < (cMinLonZones - cMaxOverlappingZones)) || (lon > (cMaxLonZones + cMaxOverlappingZones)))
     197            return true;
     198        return false;
     199    }
     200
     201    public double getDefaultZoomInPPD() {
     202        // TODO Auto-generated method stub
     203        return 0;
     204    }
     205
     206    public Bounds getWorldBoundsLatLon()
     207    {
     208        // These are not the Lambert Zone boundaries but we keep these values until coordinates outside the
     209        // projection boundaries are handled correctly.
     210        return new Bounds(
     211                new LatLon(-85.05112877980659, -180.0),
     212                new LatLon(85.05112877980659, 180.0));
     213        /*return new Bounds(
     214                new LatLon(45.0, -4.9074074074074059),
     215                new LatLon(57.0, 10.2));*/
     216    }
     217
     218}
     219
  • Projection.java

     
    3030        new UTM_20N_Guadeloupe_Ste_Anne(),
    3131        new UTM_20N_Guadeloupe_Fort_Marigot(),
    3232        new UTM_20N_Martinique_Fort_Desaix(),
    33         new GaussLaborde_Reunion()
     33        new GaussLaborde_Reunion(),
     34        new LambertCC9Zones()    // Still needs proper default zoom
    3435    };
    3536
    3637    /**