svn-gvsig-desktop / tags / v1_1_Build_1001 / libraries / libCq_CMS_praster / src / org / cresques / geo / Mercator.java @ 11984
History | View | Annotate | Download (12 KB)
1 | 8026 | nacho | /*
|
---|---|---|---|
2 | * Cresques Mapping Suite. Graphic Library for constructing mapping applications.
|
||
3 | *
|
||
4 | * Copyright (C) 2004-5.
|
||
5 | *
|
||
6 | * This program is free software; you can redistribute it and/or
|
||
7 | * modify it under the terms of the GNU General Public License
|
||
8 | * as published by the Free Software Foundation; either version 2
|
||
9 | * of the License, or (at your option) any later version.
|
||
10 | *
|
||
11 | * This program is distributed in the hope that it will be useful,
|
||
12 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||
13 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||
14 | * GNU General Public License for more details.
|
||
15 | *
|
||
16 | * You should have received a copy of the GNU General Public License
|
||
17 | * along with this program; if not, write to the Free Software
|
||
18 | * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,USA.
|
||
19 | *
|
||
20 | * For more information, contact:
|
||
21 | *
|
||
22 | * cresques@gmail.com
|
||
23 | */
|
||
24 | package org.cresques.geo; |
||
25 | |||
26 | import org.cresques.cts.ICoordTrans; |
||
27 | import org.cresques.cts.IDatum; |
||
28 | import org.cresques.cts.IProjection; |
||
29 | |||
30 | import org.cresques.px.Extent; |
||
31 | |||
32 | import java.awt.FontMetrics; |
||
33 | import java.awt.Graphics2D; |
||
34 | import java.awt.geom.AffineTransform; |
||
35 | import java.awt.geom.Point2D; |
||
36 | 10645 | nacho | import java.awt.geom.Rectangle2D; |
37 | 8026 | nacho | |
38 | import java.util.TreeMap; |
||
39 | |||
40 | |||
41 | /**
|
||
42 | * Proyeccion Mercator
|
||
43 | * @author "Luis W. Sevilla" <sevilla_lui@gva.es>* @author administrador
|
||
44 | */
|
||
45 | public class Mercator extends Projection { |
||
46 | static String name = "Mercator"; |
||
47 | static String abrev = "Mer"; |
||
48 | private static TreeMap projPool = new TreeMap(); |
||
49 | public static final Mercator hayford = new Mercator(Ellipsoid.hayford); |
||
50 | public static final Mercator wgs84 = new Mercator(Ellipsoid.wgs84); |
||
51 | private double a; |
||
52 | private double f; |
||
53 | private double b; |
||
54 | private double Eps2; |
||
55 | private double EE2; |
||
56 | private double EE3; |
||
57 | private double Epps2; |
||
58 | |||
59 | public Mercator(Ellipsoid eli) {
|
||
60 | super(eli);
|
||
61 | grid = new Graticule(this); |
||
62 | |||
63 | double[] p = eli.getParam(); |
||
64 | a = p[1];
|
||
65 | f = 1 / p[2]; |
||
66 | b = p[3];
|
||
67 | |||
68 | Eps2 = p[5];
|
||
69 | EE2 = Eps2 * Eps2; |
||
70 | EE3 = EE2 * Eps2; |
||
71 | Epps2 = p[7];
|
||
72 | } |
||
73 | |||
74 | public String getAbrev() { |
||
75 | return abrev;
|
||
76 | } |
||
77 | |||
78 | public static Mercator getProjection(Ellipsoid eli) { |
||
79 | Mercator ret = null;
|
||
80 | |||
81 | if (projPool.containsKey(eli.getName())) {
|
||
82 | ret = (Mercator) Mercator.projPool.get(eli.getName()); |
||
83 | } else {
|
||
84 | if (eli == Ellipsoid.hayford) {
|
||
85 | ret = hayford; |
||
86 | } else if (eli == Ellipsoid.wgs84) { |
||
87 | ret = wgs84; |
||
88 | } else {
|
||
89 | ret = new Mercator(eli);
|
||
90 | } |
||
91 | |||
92 | projPool.put(eli.getName(), ret); |
||
93 | } |
||
94 | |||
95 | return ret;
|
||
96 | } |
||
97 | |||
98 | /**
|
||
99 | *
|
||
100 | */
|
||
101 | public static IProjection getProjectionByName(IDatum eli, String name) { |
||
102 | if (name.indexOf("ME") < 0) { |
||
103 | return null; |
||
104 | } |
||
105 | |||
106 | return getProjection((Ellipsoid) eli);
|
||
107 | } |
||
108 | |||
109 | /**
|
||
110 | *
|
||
111 | */
|
||
112 | public Point2D createPoint(double x, double y) { |
||
113 | return new ProjPoint(this, x, y); |
||
114 | } |
||
115 | |||
116 | /**
|
||
117 | *
|
||
118 | * @param uPt
|
||
119 | * @return
|
||
120 | */
|
||
121 | public Point2D toGeo(Point2D mPt) { |
||
122 | GeoPoint gPt = new GeoPoint();
|
||
123 | |||
124 | return toGeo((ProjPoint) mPt, gPt);
|
||
125 | } |
||
126 | |||
127 | /**
|
||
128 | *
|
||
129 | * @param mPt
|
||
130 | * @param gPt
|
||
131 | * @return
|
||
132 | */
|
||
133 | public GeoPoint toGeo(ProjPoint mPt, GeoPoint gPt) {
|
||
134 | double t = Math.pow(Math.E, (-mPt.getY() / a)); |
||
135 | |||
136 | double x1;
|
||
137 | double x = ((Math.PI / 2) - (2 * Math.atan(t))); |
||
138 | |||
139 | do {
|
||
140 | x1 = x; |
||
141 | x = (Math.PI / 2) - |
||
142 | (2 * Math.atan(t * (Math.pow((1 - |
||
143 | (Math.sqrt(Eps2) * Math.sin(x))) / (1 + |
||
144 | (Math.sqrt(Eps2) * Math.sin(x))), |
||
145 | (Math.sqrt(Eps2) / 2))))); |
||
146 | } while ((x - x1) > 0.0000000001); |
||
147 | |||
148 | double lat = (Math.PI / 2) - |
||
149 | (2 * Math.atan(t * (Math.pow((1 - |
||
150 | (Math.sqrt(Eps2) * Math.sin(x))) / (1 + |
||
151 | (Math.sqrt(Eps2) * Math.sin(x))), |
||
152 | (Math.sqrt(Eps2) / 2))))); |
||
153 | |||
154 | double lng = mPt.getX() / a;
|
||
155 | gPt.setLocation((lng * 180.0) / Math.PI, (lat * 180.0) / Math.PI); |
||
156 | gPt.proj = Geodetic.getProjection(((Projection) mPt.proj).eli); |
||
157 | |||
158 | return gPt;
|
||
159 | } |
||
160 | |||
161 | /**
|
||
162 | *
|
||
163 | * @param gPt
|
||
164 | * @param uPt
|
||
165 | * @return
|
||
166 | */
|
||
167 | public Point2D fromGeo(Point2D gPt, Point2D mPt) { |
||
168 | double sl = Math.sin(((GeoPoint) gPt).Latitude.ToRadians()); |
||
169 | double cl = Math.cos(((GeoPoint) gPt).Latitude.ToRadians()); |
||
170 | double tl = (1 + sl) / (1 - sl); |
||
171 | |||
172 | // Calcula Easting
|
||
173 | double x = a * ((GeoPoint) gPt).Longitude.ToRadians();
|
||
174 | |||
175 | // Calcula Northing
|
||
176 | double y = Math.pow(((1 - (Math.sqrt(Eps2) * sl)) / (1 + |
||
177 | (Math.sqrt(Eps2) * sl))), (Math.sqrt(Eps2))); |
||
178 | y = a / 2 * (Math.log(tl * y)); |
||
179 | ((ProjPoint) mPt).setLocation(x, y); |
||
180 | ((ProjPoint) mPt).proj = this;
|
||
181 | |||
182 | return mPt;
|
||
183 | } |
||
184 | |||
185 | // Calcula el step en funci?n del zoom
|
||
186 | private void generateGrid(Graphics2D g, Extent extent, AffineTransform mat) { |
||
187 | // calculo del step en funci?n del zoom
|
||
188 | Point2D pt1 = extent.getMin();
|
||
189 | |||
190 | double step = 3.0; |
||
191 | double x = pt1.getX();
|
||
192 | double dist = 0.0; |
||
193 | ProjPoint ppt1; |
||
194 | ProjPoint ppt2; |
||
195 | GeoPoint gp1; |
||
196 | GeoPoint gp2; |
||
197 | ppt1 = (ProjPoint) createPoint(x, pt1.getY()); |
||
198 | ppt2 = (ProjPoint) createPoint(x + 100, pt1.getY() - 100); |
||
199 | gp1 = (GeoPoint) ppt1.toGeo(); |
||
200 | gp2 = (GeoPoint) ppt2.toGeo(); |
||
201 | |||
202 | /* GeoPoint gp1, gp2;
|
||
203 | gp1 = (GeoPoint) createPoint( x, (int) pt1.getY());
|
||
204 | mat.transform(gp1, gp1);
|
||
205 | gp2 = (GeoPoint) createPoint(gp1.getX()+100, gp1.getY()-100);
|
||
206 | try {
|
||
207 | mat.inverseTransform(gp2, gp2);
|
||
208 | } catch (NoninvertibleTransformException e) {
|
||
209 | // TODO Auto-generated catch block
|
||
210 | e.printStackTrace();
|
||
211 | }
|
||
212 | dist = (gp2.getX()-x);
|
||
213 | System.err.println("distX = " + dist);
|
||
214 | |||
215 | if (dist > 30.0) { step = 30.0;
|
||
216 | } else if (dist > 18.0) { step = 18.0;
|
||
217 | } else if (dist > 12.0) { step = 12.0;
|
||
218 | } else if (dist > 6.0) { step = 6.0;
|
||
219 | } else if (dist > 3.0) { step = 3.0;
|
||
220 | } else if (dist > 2.0) { step = 2.0;
|
||
221 | } else if (dist > 1.0) { step = 1.0;
|
||
222 | } else if (dist > .5) { step =.5;
|
||
223 | } else if (dist > .25) { step =.25;
|
||
224 | } else if (dist > 1.0/60*5.0) { step = 1.0/60*5.0;
|
||
225 | } else { step = 1.0/60*2.0;
|
||
226 | }
|
||
227 | //step = 1.0;
|
||
228 | */
|
||
229 | generateGrid(g, extent, mat, step); |
||
230 | } |
||
231 | |||
232 | private void generateGrid(Graphics2D g, Extent extent, AffineTransform mat, |
||
233 | double step) {
|
||
234 | grid = new Graticule(this); |
||
235 | |||
236 | Point2D pt1 = extent.getMin();
|
||
237 | Point2D pt2 = extent.getMax();
|
||
238 | Point2D.Double ptx = new Point2D.Double(0.0, 0.0); |
||
239 | GeoPoint gp1; |
||
240 | GeoPoint gp2; |
||
241 | ProjPoint up1 = (ProjPoint) createPoint(0, 0); |
||
242 | ProjPoint up2 = (ProjPoint) createPoint(0, 0); |
||
243 | Geodetic geoProj = Geodetic.getProjection((Ellipsoid) getDatum()); |
||
244 | double xAxis;
|
||
245 | double yAxis;
|
||
246 | |||
247 | // Calculos para el texto
|
||
248 | FontMetrics fm = g.getFontMetrics();
|
||
249 | int fmWidth = 0; |
||
250 | int fmHeight = fm.getAscent();
|
||
251 | String tit = ""; |
||
252 | String fmt = "%G?%N"; |
||
253 | |||
254 | if (step < 1.0) { |
||
255 | fmt = "%G?%M'%N";
|
||
256 | } |
||
257 | |||
258 | // Lineas Horzontales
|
||
259 | gp1 = (GeoPoint) toGeo(new ProjPoint(pt1));
|
||
260 | gp2 = (GeoPoint) toGeo(new ProjPoint(pt2));
|
||
261 | xAxis = gp1.getX(); |
||
262 | yAxis = gp2.getY(); |
||
263 | System.err.println(name + ": ViewPort Extent = (" + gp1 + "," + gp2 + |
||
264 | ")");
|
||
265 | |||
266 | double xMin = (int) gp1.getX() - 1.0; |
||
267 | xMin -= (xMin % step); |
||
268 | |||
269 | double xMax = (int) gp2.getX() + 1.0; |
||
270 | double yMin = (int) gp1.getY() - 1.0; |
||
271 | yMin -= (yMin % step); |
||
272 | |||
273 | double yMax = (int) gp2.getY() + 1.0; |
||
274 | |||
275 | if (xMin < -180.0) { |
||
276 | xMin = -180.0;
|
||
277 | } |
||
278 | |||
279 | if (xMax > 180.0) { |
||
280 | xMax = 180.0;
|
||
281 | } |
||
282 | |||
283 | if (yMin < -80.0) { |
||
284 | yMin = -80.0;
|
||
285 | } |
||
286 | |||
287 | if (yMax > 80.0) { |
||
288 | yMax = 80.0;
|
||
289 | } |
||
290 | |||
291 | if (xAxis < -180.0) { |
||
292 | xAxis = -180.0;
|
||
293 | } |
||
294 | |||
295 | if (yAxis > 80.0) { |
||
296 | yAxis = 80.0;
|
||
297 | } |
||
298 | |||
299 | for (double y = yMin; y <= yMax; y += step) { |
||
300 | gp1 = (GeoPoint) geoProj.createPoint(xAxis, y); |
||
301 | gp2 = (GeoPoint) geoProj.createPoint(xMax, y); |
||
302 | fromGeo(gp1, up1); |
||
303 | fromGeo(gp2, up2); |
||
304 | mat.transform(up1, up1); |
||
305 | mat.transform(up2, up2); |
||
306 | grid.addLine(up1, up2); |
||
307 | |||
308 | tit = coordToString(y, fmt, true);
|
||
309 | |||
310 | //fmWidth = fm.stringWidth(tit);
|
||
311 | ptx.setLocation(up1.getX() + 3, up1.getY() - 2); |
||
312 | grid.addText(tit, ptx); |
||
313 | } |
||
314 | |||
315 | // Lineas Verticales
|
||
316 | for (double x = xMin; x <= xMax; x += step) { |
||
317 | gp1 = (GeoPoint) geoProj.createPoint(x, yMin); |
||
318 | gp2 = (GeoPoint) geoProj.createPoint(x, yAxis); |
||
319 | fromGeo(gp1, up1); |
||
320 | fromGeo(gp2, up2); |
||
321 | mat.transform(up1, up1); |
||
322 | mat.transform(up2, up2); |
||
323 | grid.addLine(up1, up2); |
||
324 | |||
325 | tit = coordToString(x, fmt, false);
|
||
326 | |||
327 | //fmWidth = fm.stringWidth(tit);
|
||
328 | ptx.setLocation(up2.getX() + 3, up2.getY() + fmHeight);
|
||
329 | grid.addText(tit, ptx); |
||
330 | } |
||
331 | } |
||
332 | |||
333 | public void drawGrid(Graphics2D g, ViewPortData vp) { |
||
334 | generateGrid(g, vp.getExtent(), vp.getMat()); |
||
335 | grid.setColor(gridColor); |
||
336 | grid.draw(g, vp); |
||
337 | } |
||
338 | |||
339 | /* (non-Javadoc)
|
||
340 | * @see org.cresques.cts.IProjection#getScale(double, double, double, double)
|
||
341 | */
|
||
342 | public double getScale(double minX, double maxX, double width, double dpi) { |
||
343 | Projection prj = Geodetic.getProjection((Ellipsoid) getDatum()); |
||
344 | GeoPoint pt1 = (GeoPoint) prj.createPoint(1.0, 0.0); |
||
345 | GeoPoint pt2 = (GeoPoint) prj.createPoint(2.0, 0.0); |
||
346 | ProjPoint ppt1 = (ProjPoint) createPoint(0.0, 0.0); |
||
347 | ProjPoint ppt2 = (ProjPoint) createPoint(0.0, 0.0); |
||
348 | fromGeo(pt1, ppt1); |
||
349 | fromGeo(pt2, ppt2); |
||
350 | |||
351 | //scale = ppt2.getX()-ppt1.getX();
|
||
352 | 10645 | nacho | double scale = (((maxX - minX) / (ppt2.getX() - ppt1.getX())) *
|
353 | 8026 | nacho | //scale = ((extent.maxX()-extent.minX())/ getWidth());// *
|
354 | (dpi / 2.54 * 100.0 * 1852.0 * 60.0)) / width; |
||
355 | |||
356 | return scale;
|
||
357 | } |
||
358 | |||
359 | public ICoordTrans getCT(IProjection dest) {
|
||
360 | // TODO Auto-generated method stub
|
||
361 | return null; |
||
362 | } |
||
363 | 10645 | nacho | |
364 | public Rectangle2D getExtent(Rectangle2D extent, double scale, double wImage, double hImage, double changeUnits, double dpi) { |
||
365 | Projection prj = Geodetic.getProjection((Ellipsoid) getDatum()); |
||
366 | GeoPoint pt1 = (GeoPoint) prj.createPoint(1.0, 0.0); |
||
367 | GeoPoint pt2 = (GeoPoint) prj.createPoint(2.0, 0.0); |
||
368 | ProjPoint ppt1 = (ProjPoint) createPoint(0.0, 0.0); |
||
369 | ProjPoint ppt2 = (ProjPoint) createPoint(0.0, 0.0); |
||
370 | fromGeo(pt1, ppt1); |
||
371 | fromGeo(pt2, ppt2); |
||
372 | double w =0; |
||
373 | double h =0; |
||
374 | double wExtent =0; |
||
375 | double hExtent =0; |
||
376 | w = ((wImage / dpi) * 2.54);
|
||
377 | h = ((hImage / dpi) * 2.54);
|
||
378 | wExtent =((w*scale)/ (ppt2.getX() - ppt1.getX()))/ (changeUnits*1852.0*60.0); |
||
379 | hExtent =((h*scale)/ (ppt2.getX() - ppt1.getX()))/ (changeUnits*1852.0*60.0); |
||
380 | double xExtent = extent.getCenterX() - wExtent/2; |
||
381 | double yExtent = extent.getCenterY() - hExtent/2; |
||
382 | Rectangle2D rec=new Rectangle2D.Double(xExtent,yExtent,wExtent,hExtent); |
||
383 | return rec;
|
||
384 | } |
||
385 | 11345 | dguerrero | |
386 | public String getFullCode() { |
||
387 | return getAbrev();
|
||
388 | } |
||
389 | 8026 | nacho | } |