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

Last change on this file since 4067 was 3779, checked in by Upliner, 13 years ago

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