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

Last change on this file since 2509 was 2507, checked in by stoecker, 14 years ago

apply #2989 - patch by pieren - improve lambert projection
File size: 10.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.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 ntf_rgf93Grid = new NTV2GridShiftFile();
100 ntf_rgf93Grid.loadGridShiftFile(is, false);
101 //System.out.println("NTF<->RGF93 grid loaded.");
102 } catch (Exception e) {
103 e.printStackTrace();
104 }
105 }
106 }
107
108 /**
109 * @param p WGS84 lat/lon (ellipsoid GRS80) (in degree)
110 * @return eastnorth projection in Lambert Zone (ellipsoid Clark)
111 * @throws IOException
112 */
113 public EastNorth latlon2eastNorth(LatLon p) {
114 // translate ellipsoid GRS80 (WGS83) => Clark
115 LatLon geo = WGS84_to_NTF(p);
116 double lt = Math.toRadians(geo.lat()); // in radian
117 double lg = Math.toRadians(geo.lon());
118
119 // check if longitude and latitude are inside the French Lambert zones
120 if (lt >= cMinLatZone1Radian && lt <= cMaxLatZone1Radian && lg >= cMinLonZonesRadian && lg <= cMaxLonZonesRadian)
121 return ConicProjection(lt, lg, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]);
122 return ConicProjection(lt, lg, Xs[0], Ys[0], c[0], n[0]);
123 }
124
125 public LatLon eastNorth2latlon(EastNorth p) {
126 LatLon geo;
127 geo = Geographic(p, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]);
128 // translate geodetic system from NTF (ellipsoid Clark) to RGF93/WGS84 (ellipsoid GRS80)
129 return NTF_to_WGS84(geo);
130 }
131
132 @Override public String toString() {
133 return tr("Lambert 4 Zones (France)");
134 }
135
136 public String toCode() {
137 return "EPSG:"+(27561+layoutZone);
138 }
139
140 public String getCacheDirectoryName() {
141 return "lambert";
142 }
143
144 /**
145 * Initializes from geographic coordinates. Note that reference ellipsoid
146 * used by Lambert is the Clark ellipsoid.
147 *
148 * @param lat latitude in grad
149 * @param lon longitude in grad
150 * @param Xs false east (coordinate system origin) in meters
151 * @param Ys false north (coordinate system origin) in meters
152 * @param c projection constant
153 * @param n projection exponent
154 * @return EastNorth projected coordinates in meter
155 */
156 private EastNorth ConicProjection(double lat, double lon, double Xs, double Ys, double c, double n) {
157 double eslt = Ellipsoid.clarke.e * Math.sin(lat);
158 double l = Math.log(Math.tan(Math.PI / 4.0 + (lat / 2.0))
159 * Math.pow((1.0 - eslt) / (1.0 + eslt), Ellipsoid.clarke.e / 2.0));
160 double east = Xs + c * Math.exp(-n * l) * Math.sin(n * (lon - lg0));
161 double north = Ys - c * Math.exp(-n * l) * Math.cos(n * (lon - lg0));
162 return new EastNorth(east, north);
163 }
164
165 /**
166 * Initializes from projected coordinates (conic projection). Note that
167 * reference ellipsoid used by Lambert is Clark
168 *
169 * @param eastNorth projected coordinates pair in meters
170 * @param Xs false east (coordinate system origin) in meters
171 * @param Ys false north (coordinate system origin) in meters
172 * @param c projection constant
173 * @param n projection exponent
174 * @return LatLon in degree
175 */
176 private LatLon Geographic(EastNorth eastNorth, double Xs, double Ys, double c, double n) {
177 double dx = eastNorth.east() - Xs;
178 double dy = Ys - eastNorth.north();
179 double R = Math.sqrt(dx * dx + dy * dy);
180 double gamma = Math.atan(dx / dy);
181 double l = -1.0 / n * Math.log(Math.abs(R / c));
182 l = Math.exp(l);
183 double lon = lg0 + gamma / n;
184 double lat = 2.0 * Math.atan(l) - Math.PI / 2.0;
185 double delta = 1.0;
186 while (delta > epsilon) {
187 double eslt = Ellipsoid.clarke.e * Math.sin(lat);
188 double nlt = 2.0 * Math.atan(Math.pow((1.0 + eslt) / (1.0 - eslt), Ellipsoid.clarke.e / 2.0) * l) - Math.PI
189 / 2.0;
190 delta = Math.abs(nlt - lat);
191 lat = nlt;
192 }
193 return new LatLon(Math.toDegrees(lat), Math.toDegrees(lon)); // in radian
194 }
195
196 /**
197 * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to Lambert
198 * geographic, (ellipsoid Clark)
199 * @throws IOException
200 */
201 private LatLon WGS84_to_NTF(LatLon wgs) {
202 NTV2GridShift gs = new NTV2GridShift(wgs);
203 if (ntf_rgf93Grid != null) {
204 ntf_rgf93Grid.gridShiftReverse(gs);
205 return new LatLon(wgs.lat()+gs.getLatShiftDegrees(), wgs.lon()+gs.getLonShiftPositiveEastDegrees());
206 }
207 return new LatLon(0,0);
208 }
209
210 private LatLon NTF_to_WGS84(LatLon ntf) {
211 NTV2GridShift gs = new NTV2GridShift(ntf);
212 if (ntf_rgf93Grid != null) {
213 ntf_rgf93Grid.gridShiftForward(gs);
214 return new LatLon(ntf.lat()+gs.getLatShiftDegrees(), ntf.lon()+gs.getLonShiftPositiveEastDegrees());
215 }
216 return new LatLon(0,0);
217 }
218
219 public Bounds getWorldBoundsLatLon()
220 {
221 Bounds b= new Bounds(
222 new LatLon(zoneLimitsDegree[layoutZone][1] - cMaxOverlappingZonesDegree, -4.9074074074074059),
223 new LatLon(zoneLimitsDegree[layoutZone][0] + cMaxOverlappingZonesDegree, 10.2));
224 return b;
225 }
226
227 /**
228 * Returns the default zoom scale in pixel per degree ({@see #NavigatableComponent#scale}))
229 */
230 public double getDefaultZoomInPPD() {
231 // this will set the map scaler to about 1000 m (in default scale, 1 pixel will be 10 meters)
232 return 10.0;
233 }
234
235 public int getLayoutZone() {
236 return layoutZone;
237 }
238
239 public static String[] lambert4zones = {
240 tr("{0} ({1} to {2} degrees)", 1,"51.30","48.15"),
241 tr("{0} ({1} to {2} degrees)", 1,"48.15","45.45"),
242 tr("{0} ({1} to {2} degrees)", 1,"45.45","42.76"),
243 tr("{0} (Corsica)", 4)
244 };
245
246 public void setupPreferencePanel(JPanel p) {
247 JComboBox prefcb = new JComboBox(lambert4zones);
248
249 prefcb.setSelectedIndex(layoutZone);
250 p.setLayout(new GridBagLayout());
251 p.add(new JLabel(tr("Lambert CC Zone")), GBC.std().insets(5,5,0,5));
252 p.add(GBC.glue(1, 0), GBC.std().fill(GBC.HORIZONTAL));
253 /* Note: we use component position 2 below to find this again */
254 p.add(prefcb, GBC.eop().fill(GBC.HORIZONTAL));
255 p.add(GBC.glue(1, 1), GBC.eol().fill(GBC.BOTH));
256 }
257
258 public Collection<String> getPreferences(JPanel p) {
259 Object prefcb = p.getComponent(2);
260 if(!(prefcb instanceof JComboBox))
261 return null;
262 layoutZone = ((JComboBox)prefcb).getSelectedIndex();
263 return Collections.singleton(Integer.toString(layoutZone+1));
264 }
265
266 public void setPreferences(Collection<String> args) {
267 layoutZone = DEFAULT_ZONE;
268 if (args != null) {
269 try {
270 for(String s : args)
271 {
272 layoutZone = Integer.parseInt(s)-1;
273 if(layoutZone < 0 || layoutZone > 3) {
274 layoutZone = DEFAULT_ZONE;
275 }
276 break;
277 }
278 } catch(NumberFormatException e) {}
279 }
280 }
281
282 public Collection<String> getPreferencesFromCode(String code) {
283 if (code.startsWith("EPSG:2756") && code.length() == 9) {
284 try {
285 String zonestring = code.substring(9);
286 int zoneval = Integer.parseInt(zonestring);
287 if(zoneval >= 1 && zoneval <= 4)
288 return Collections.singleton(zonestring);
289 } catch(NumberFormatException e) {}
290 }
291 return null;
292 }
293
294}
Note: See TracBrowser for help on using the repository browser.