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

Last change on this file since 3715 was 3480, checked in by bastiK, 14 years ago

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