| 1 | //License: GPL. For details, see LICENSE file. |
| 2 | package org.openstreetmap.josm.data.projection; |
| 3 | |
| 4 | import static org.openstreetmap.josm.tools.I18n.tr; |
| 5 | |
| 6 | import javax.swing.JOptionPane; |
| 7 | |
| 8 | import org.openstreetmap.josm.Main; |
| 9 | import org.openstreetmap.josm.data.Bounds; |
| 10 | import org.openstreetmap.josm.data.coor.EastNorth; |
| 11 | import org.openstreetmap.josm.data.coor.LatLon; |
| 12 | import org.openstreetmap.josm.data.projection.Projection; |
| 13 | import 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 | */ |
| 21 | public 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 | |