source: josm/trunk/src/org/openstreetmap/josm/data/projection/Lambert.java@ 3204

Last change on this file since 3204 was 3083, checked in by bastiK, 14 years ago

added svn:eol-style=native to source files
  • Property svn:eol-style set to native
File size: 11.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 java.awt.GridBagLayout;
7import java.io.IOException;
8import java.io.InputStream;
9import java.util.Collection;
10import java.util.Collections;
11
12import javax.swing.JComboBox;
13import javax.swing.JLabel;
14import javax.swing.JPanel;
15
16import org.openstreetmap.josm.Main;
17import org.openstreetmap.josm.data.Bounds;
18import org.openstreetmap.josm.data.coor.EastNorth;
19import org.openstreetmap.josm.data.coor.LatLon;
20import org.openstreetmap.josm.tools.GBC;
21
22/**
23 * This class provides the two methods <code>latlon2eastNorth</code> and <code>eastNorth2latlon</code>
24 * converting the JOSM LatLon coordinates in WGS84 system (GPS) to and from East North values in
25 * the projection Lambert conic conform 4 zones using the French geodetic system NTF.
26 * This newer version uses the grid translation NTF<->RGF93 provided by IGN for a submillimetric accuracy.
27 * (RGF93 is the French geodetic system similar to WGS84 but not mathematically equal)
28 * @author Pieren
29 */
30public class Lambert implements Projection, ProjectionSubPrefs {
31 /**
32 * Lambert I, II, III, and IV projection exponents
33 */
34 public static final double n[] = { 0.7604059656, 0.7289686274, 0.6959127966, 0.6712679322 };
35
36 /**
37 * Lambert I, II, III, and IV projection constants
38 */
39 public static final double c[] = { 11603796.98, 11745793.39, 11947992.52, 12136281.99 };
40
41 /**
42 * Lambert I, II, III, and IV false east
43 */
44 public static final double Xs[] = { 600000.0, 600000.0, 600000.0, 234.358 };
45
46 /**
47 * Lambert I, II, III, and IV false north
48 */
49 public static final double Ys[] = { 5657616.674, 6199695.768, 6791905.085, 7239161.542 };
50
51 /**
52 * Lambert I, II, III, and IV longitudinal offset to Greenwich meridian
53 */
54 public static final double lg0 = 0.04079234433198; // 2deg20'14.025"
55
56 /**
57 * precision in iterative schema
58 */
59
60 public static final double epsilon = 1e-11;
61
62 /**
63 * France is divided in 4 Lambert projection zones (1,2,3 + 4th for Corsica)
64 */
65 public static final double cMaxLatZone1Radian = Math.toRadians(57 * 0.9);
66 public static final double cMinLatZone1Radian = Math.toRadians(46.1 * 0.9);// lowest latitude of Zone 4 (South Corsica)
67
68 public static final double zoneLimitsDegree[][] = {
69 {Math.toDegrees(cMaxLatZone1Radian), (53.5 * 0.9)}, // Zone 1 (reference values in grad *0.9)
70 {(53.5 * 0.9), (50.5 * 0.9)}, // Zone 2
71 {(50.5 * 0.9), (47.0 * 0.9)}, // Zone 3
72 {(47.51963 * 0.9), Math.toDegrees(cMinLatZone1Radian)} // Zone 4
73 };
74
75 public static final double cMinLonZonesRadian = Math.toRadians(-4.9074074074074059 * 0.9);
76
77 public static final double cMaxLonZonesRadian = Math.toRadians(10.2 * 0.9);
78
79 /**
80 * Allow some extension beyond the theoretical limits
81 */
82 public static final double cMaxOverlappingZonesDegree = 1.5;
83
84 public static final int DEFAULT_ZONE = 0;
85
86 private static int layoutZone = DEFAULT_ZONE;
87
88 private static NTV2GridShiftFile ntf_rgf93Grid = null;
89
90 public static NTV2GridShiftFile getNtf_rgf93Grid() {
91 return ntf_rgf93Grid;
92 }
93
94 public Lambert() {
95 if (ntf_rgf93Grid == null) {
96 try {
97 String gridFileName = "ntf_r93_b.gsb";
98 InputStream is = Main.class.getResourceAsStream("/data/"+gridFileName);
99 if (is == null) {
100 System.err.println(tr("Warning: failed to open input stream for resource ''/data/{0}''. Cannot load NTF<->RGF93 grid", gridFileName));
101 return;
102 }
103 ntf_rgf93Grid = new NTV2GridShiftFile();
104 ntf_rgf93Grid.loadGridShiftFile(is, false);
105 //System.out.println("NTF<->RGF93 grid loaded.");
106 } catch (Exception e) {
107 e.printStackTrace();
108 }
109 }
110 }
111
112 /**
113 * @param p WGS84 lat/lon (ellipsoid GRS80) (in degree)
114 * @return eastnorth projection in Lambert Zone (ellipsoid Clark)
115 * @throws IOException
116 */
117 public EastNorth latlon2eastNorth(LatLon p) {
118 // translate ellipsoid GRS80 (WGS83) => Clark
119 LatLon geo = WGS84_to_NTF(p);
120 double lt = Math.toRadians(geo.lat()); // in radian
121 double lg = Math.toRadians(geo.lon());
122
123 // check if longitude and latitude are inside the French Lambert zones
124 if (lt >= cMinLatZone1Radian && lt <= cMaxLatZone1Radian && lg >= cMinLonZonesRadian && lg <= cMaxLonZonesRadian)
125 return ConicProjection(lt, lg, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]);
126 return ConicProjection(lt, lg, Xs[0], Ys[0], c[0], n[0]);
127 }
128
129 public LatLon eastNorth2latlon(EastNorth p) {
130 LatLon geo;
131 geo = Geographic(p, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]);
132 // translate geodetic system from NTF (ellipsoid Clark) to RGF93/WGS84 (ellipsoid GRS80)
133 return NTF_to_WGS84(geo);
134 }
135
136 @Override public String toString() {
137 return tr("Lambert 4 Zones (France)");
138 }
139
140 public String toCode() {
141 return "EPSG:"+(27561+layoutZone);
142 }
143
144 @Override
145 public int hashCode() {
146 return getClass().getName().hashCode()+layoutZone; // our only real variable
147 }
148
149 public String getCacheDirectoryName() {
150 return "lambert";
151 }
152
153 /**
154 * Initializes from geographic coordinates. Note that reference ellipsoid
155 * used by Lambert is the Clark ellipsoid.
156 *
157 * @param lat latitude in grad
158 * @param lon longitude in grad
159 * @param Xs false east (coordinate system origin) in meters
160 * @param Ys false north (coordinate system origin) in meters
161 * @param c projection constant
162 * @param n projection exponent
163 * @return EastNorth projected coordinates in meter
164 */
165 private EastNorth ConicProjection(double lat, double lon, double Xs, double Ys, double c, double n) {
166 double eslt = Ellipsoid.clarke.e * Math.sin(lat);
167 double l = Math.log(Math.tan(Math.PI / 4.0 + (lat / 2.0))
168 * Math.pow((1.0 - eslt) / (1.0 + eslt), Ellipsoid.clarke.e / 2.0));
169 double east = Xs + c * Math.exp(-n * l) * Math.sin(n * (lon - lg0));
170 double north = Ys - c * Math.exp(-n * l) * Math.cos(n * (lon - lg0));
171 return new EastNorth(east, north);
172 }
173
174 /**
175 * Initializes from projected coordinates (conic projection). Note that
176 * reference ellipsoid used by Lambert is Clark
177 *
178 * @param eastNorth projected coordinates pair in meters
179 * @param Xs false east (coordinate system origin) in meters
180 * @param Ys false north (coordinate system origin) in meters
181 * @param c projection constant
182 * @param n projection exponent
183 * @return LatLon in degree
184 */
185 private LatLon Geographic(EastNorth eastNorth, double Xs, double Ys, double c, double n) {
186 double dx = eastNorth.east() - Xs;
187 double dy = Ys - eastNorth.north();
188 double R = Math.sqrt(dx * dx + dy * dy);
189 double gamma = Math.atan(dx / dy);
190 double l = -1.0 / n * Math.log(Math.abs(R / c));
191 l = Math.exp(l);
192 double lon = lg0 + gamma / n;
193 double lat = 2.0 * Math.atan(l) - Math.PI / 2.0;
194 double delta = 1.0;
195 while (delta > epsilon) {
196 double eslt = Ellipsoid.clarke.e * Math.sin(lat);
197 double nlt = 2.0 * Math.atan(Math.pow((1.0 + eslt) / (1.0 - eslt), Ellipsoid.clarke.e / 2.0) * l) - Math.PI
198 / 2.0;
199 delta = Math.abs(nlt - lat);
200 lat = nlt;
201 }
202 return new LatLon(Math.toDegrees(lat), Math.toDegrees(lon)); // in radian
203 }
204
205 /**
206 * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to Lambert
207 * geographic, (ellipsoid Clark)
208 * @throws IOException
209 */
210 private LatLon WGS84_to_NTF(LatLon wgs) {
211 NTV2GridShift gs = new NTV2GridShift(wgs);
212 if (ntf_rgf93Grid != null) {
213 ntf_rgf93Grid.gridShiftReverse(gs);
214 return new LatLon(wgs.lat()+gs.getLatShiftDegrees(), wgs.lon()+gs.getLonShiftPositiveEastDegrees());
215 }
216 return new LatLon(0,0);
217 }
218
219 private LatLon NTF_to_WGS84(LatLon ntf) {
220 NTV2GridShift gs = new NTV2GridShift(ntf);
221 if (ntf_rgf93Grid != null) {
222 ntf_rgf93Grid.gridShiftForward(gs);
223 return new LatLon(ntf.lat()+gs.getLatShiftDegrees(), ntf.lon()+gs.getLonShiftPositiveEastDegrees());
224 }
225 return new LatLon(0,0);
226 }
227
228 public Bounds getWorldBoundsLatLon()
229 {
230 Bounds b= new Bounds(
231 new LatLon(zoneLimitsDegree[layoutZone][1] - cMaxOverlappingZonesDegree, -4.9074074074074059),
232 new LatLon(zoneLimitsDegree[layoutZone][0] + cMaxOverlappingZonesDegree, 10.2));
233 return b;
234 }
235
236 /**
237 * Returns the default zoom scale in pixel per degree ({@see #NavigatableComponent#scale}))
238 */
239 public double getDefaultZoomInPPD() {
240 // this will set the map scaler to about 1000 m (in default scale, 1 pixel will be 10 meters)
241 return 10.0;
242 }
243
244 public int getLayoutZone() {
245 return layoutZone;
246 }
247
248 public static String[] lambert4zones = {
249 tr("{0} ({1} to {2} degrees)", 1,"51.30","48.15"),
250 tr("{0} ({1} to {2} degrees)", 2,"48.15","45.45"),
251 tr("{0} ({1} to {2} degrees)", 3,"45.45","42.76"),
252 tr("{0} (Corsica)", 4)
253 };
254
255 public void setupPreferencePanel(JPanel p) {
256 JComboBox prefcb = new JComboBox(lambert4zones);
257
258 prefcb.setSelectedIndex(layoutZone);
259 p.setLayout(new GridBagLayout());
260 p.add(new JLabel(tr("Lambert CC Zone")), GBC.std().insets(5,5,0,5));
261 p.add(GBC.glue(1, 0), GBC.std().fill(GBC.HORIZONTAL));
262 /* Note: we use component position 2 below to find this again */
263 p.add(prefcb, GBC.eop().fill(GBC.HORIZONTAL));
264 p.add(GBC.glue(1, 1), GBC.eol().fill(GBC.BOTH));
265 }
266
267 public Collection<String> getPreferences(JPanel p) {
268 Object prefcb = p.getComponent(2);
269 if(!(prefcb instanceof JComboBox))
270 return null;
271 layoutZone = ((JComboBox)prefcb).getSelectedIndex();
272 return Collections.singleton(Integer.toString(layoutZone+1));
273 }
274
275 public void setPreferences(Collection<String> args) {
276 layoutZone = DEFAULT_ZONE;
277 if (args != null) {
278 try {
279 for(String s : args)
280 {
281 layoutZone = Integer.parseInt(s)-1;
282 if(layoutZone < 0 || layoutZone > 3) {
283 layoutZone = DEFAULT_ZONE;
284 }
285 break;
286 }
287 } catch(NumberFormatException e) {}
288 }
289 }
290
291 public Collection<String> getPreferencesFromCode(String code) {
292 if (code.startsWith("EPSG:2756") && code.length() == 9) {
293 try {
294 String zonestring = code.substring(9);
295 int zoneval = Integer.parseInt(zonestring);
296 if(zoneval >= 1 && zoneval <= 4)
297 return Collections.singleton(zonestring);
298 } catch(NumberFormatException e) {}
299 }
300 return null;
301 }
302
303}
Note: See TracBrowser for help on using the repository browser.