source: josm/trunk/src/org/openstreetmap/josm/data/projection/LambertCC9Zones.java@ 2499

Last change on this file since 2499 was 2304, checked in by stoecker, 14 years ago

apply #3737 - patch by Pieren - new Lambert zones for France, still needs some rework
File size: 8.2 KB
Line 
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
Note: See TracBrowser for help on using the repository browser.