001 package edu.nrao.sss.geom;
002
003 import java.awt.geom.AffineTransform;
004 import java.awt.geom.NoninvertibleTransformException;
005 import java.awt.geom.Point2D;
006 import java.math.BigDecimal;
007
008 import edu.nrao.sss.measure.ArcUnits;
009 import edu.nrao.sss.measure.Latitude;
010 import edu.nrao.sss.measure.Longitude;
011
012 /**
013 * A coordinate transformer that maps a point on a sphere to a point on a
014 * plane using an azimuthal equidistant projection. This projection is
015 * such that the relative distances from a point at the center of the
016 * plane to other points elsewhere on the plane are accurate in their
017 * ratios to one another. This makes it a useful projection for comparing
018 * relative distances from a fixed point.
019 * <p>
020 * The algorithms used in this class were taken from:
021 * <blockquote>
022 * <a href="http://mathworld.wolfram.com/about/author.html">
023 * Weisstein, Eric W</a>,
024 * <a href="http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html">
025 * "Azimuthal Equidistant Projection"</a>, from
026 * <a href="http://mathworld.wolfram.com/">
027 * <i>MathWorld</i></a>--A Wolfram Web Resource.
028 * </blockquote></p>
029 * <p>
030 * <b>Version Info:</b>
031 * <table style="margin-left:2em">
032 * <tr><td>$Revision: 1239 $</td></tr>
033 * <tr><td>$Date: 2008-04-25 10:34:57 -0600 (Fri, 25 Apr 2008) $</td></tr>
034 * <tr><td>$Author: dharland $ (last person to modify)</td></tr>
035 * </table></p>
036 *
037 * @author David M. Harland
038 * @since 2007-06-07
039 */
040 public class AzimuthalEquidistantProjector
041 {
042 private SphericalPosition centerPoint;
043 private AffineTransform xyTransformer;
044
045 private double latCtr, lonCtr;
046 private double cosLatCtr, sinLatCtr;
047
048 /** Creates a new instance. */
049 public AzimuthalEquidistantProjector()
050 {
051 setCenter(new EarthPosition());
052 xyTransformer = new AffineTransform();
053 }
054
055 /**
056 * Sets the center point, in longitude and latitude, for this projector.
057 * @param centerPosition
058 * the center point, in longitude and latitude, for this projector.
059 */
060 public final void setCenter(SphericalPosition centerPosition)
061 {
062 centerPoint = centerPosition;
063
064 latCtr = centerPoint.getLatitude().toUnits(ArcUnits.RADIAN).doubleValue();
065 lonCtr = centerPoint.getLongitude().toUnits(ArcUnits.RADIAN).doubleValue();
066
067 cosLatCtr = Math.cos(latCtr);
068 sinLatCtr = Math.sin(latCtr);
069 }
070
071 /**
072 * Returns the 2D-to-2D planar transformer used by this projector.
073 * This transformer is initially created as an identity transformer.
074 * The returned transformer is a reference to the one held by this
075 * projector, so changes made to it are reflected herein.
076 * <p>
077 * This projector turns a longitude or latitude of ninety degrees
078 * into an x or y value of one-half pi. By properly configuring
079 * the returned transformer, you can change this projector's center
080 * coordinate and x and y ranges to other values.</p>
081 *
082 * @return the affine transform used by this projector to transform the
083 * x,y coordinates it produces into x',y' coordinates.
084 */
085 public AffineTransform getXyTransformer()
086 {
087 return xyTransformer;
088 }
089
090 /**
091 * Projects the longitude, latitude position on a sphere into an x,y
092 * position on a plane. The returned value is first run through
093 * an azimuthal equidistant projection, and then through an
094 * {@link #getXyTransformer() affine transform} that may be configured
095 * by clients.
096 *
097 * @param sphericalPosition
098 * a longitude, latitude coordinate pair on a sphere.
099 * (The distance attribute is ignored.)
100 *
101 * @return an x,y point on a plane created by an azimuthal equidistant
102 * projection from the given position on a sphere.
103 */
104 public Point2D getXyFor(SphericalPosition sphericalPosition)
105 {
106 double latPos = sphericalPosition.getLatitude().toUnits(ArcUnits.RADIAN).doubleValue();
107 double lonPos = sphericalPosition.getLongitude().toUnits(ArcUnits.RADIAN).doubleValue();
108
109 double cosLatPos = Math.cos(latPos);
110 double sinLatPos = Math.sin(latPos);
111
112 double deltaLon = lonPos - lonCtr;
113 double cosDeltaLon = Math.cos(deltaLon);
114
115 double c =
116 Math.acos(sinLatCtr * sinLatPos + cosLatCtr * cosLatPos * cosDeltaLon);
117
118 double k = (c == 0.0) ? 0.0 : c / Math.sin(c);
119
120 double x = k * cosLatPos * Math.sin(deltaLon);
121
122 double y = k * (cosLatCtr*sinLatPos - sinLatCtr*cosLatPos*cosDeltaLon);
123
124 Point2D rawPoint = new Point2D.Double(x, y);
125
126 return
127 xyTransformer.isIdentity() ? rawPoint
128 : xyTransformer.transform(rawPoint, rawPoint);
129 }
130
131 /**
132 * Performs the inverse projection of an x,y point on a plane to a
133 * longitude, latitude position on a sphere. The x,y point is first run
134 * through the inverse of an
135 * {@link #getXyTransformer() affine transform} that may be configured
136 * by clients. It is then run through the inverse of an azimuthal
137 * equidistant projection.
138 *
139 * @param xyPosition an x,y position on a plane.
140 *
141 * @return a longitude, latitude position on a sphere created by an
142 * inverse azimuthal equidistant projection from the given position
143 * on a plane.
144 */
145 public SphericalPosition getLatLonFor(Point2D xyPosition)
146 throws NoninvertibleTransformException
147 {
148 double x = xyPosition.getX();
149 double y = xyPosition.getY();
150
151 //Apply inverse transformation
152 if (!xyTransformer.isIdentity())
153 {
154 Point2D raw = new Point2D.Double(x, y);
155 Point2D xf = xyTransformer.inverseTransform(raw, null);
156
157 x = xf.getX();
158 y = xf.getY();
159 }
160
161 double c = Math.sqrt(x*x + y*y);
162
163 double cosC = Math.cos(c);
164 double sinC = Math.sin(c);
165
166 double sinLat = cosC * sinLatCtr + y * sinC * cosLatCtr / c;
167
168 double lat = Math.asin(sinLat);
169 double lon = lonCtr;
170
171 if (sinLat >= +1)
172 {
173 lon += Math.atan(-x/y);
174 }
175 else if (sinLat <= -1)
176 {
177 lon += Math.atan(+x/y);
178 }
179 else
180 {
181 lon += Math.atan(x * sinC /
182 (c * cosLatCtr * cosC - y * sinLatCtr * sinC));
183 }
184
185 Latitude latitude = new Latitude (new BigDecimal(lat), ArcUnits.RADIAN);
186 Longitude longitude = new Longitude(new BigDecimal(lon), ArcUnits.RADIAN);
187 EarthPosition ep = new EarthPosition(longitude, latitude, null);
188
189 return ep;
190 }
191
192 //============================================================================
193 //
194 //============================================================================
195 /*
196 public static void main(String[] args) throws Exception
197 {
198 AzimuthalEquidistantProjector projector = new AzimuthalEquidistantProjector();
199
200 AffineTransform xform = projector.getXyTransformer();
201 double scale = 100.0 / (Math.PI / 2.0);
202 //xform.translate(100, 100);
203 xform.scale(scale, -scale);
204
205 SphericalPosition center =
206 new edu.nrao.sss.astronomy.SimpleSkyPosition(); //Lat=0, Lon=0
207
208 projector.setCenter(center);
209
210 System.out.println(" 0d, 0d => " + projector.getXyFor(center));
211
212 System.out.println();
213
214 SphericalPosition pos = new edu.nrao.sss.astronomy.SimpleSkyPosition();
215 pos.getLongitude().set(90.0, ArcUnits.DEGREE);
216 System.out.println(" 90d, 0d => " + projector.getXyFor(pos));
217
218 pos.getLongitude().set(0.0, ArcUnits.DEGREE);
219 pos.getLatitude().set(90.0, ArcUnits.DEGREE);
220 System.out.println(" 0d, 90d => " + projector.getXyFor(pos));
221
222 pos.getLongitude().set(-90.0, ArcUnits.DEGREE);
223 pos.getLatitude().set(0.0, ArcUnits.DEGREE);
224 System.out.println(" -90d, 0d => " + projector.getXyFor(pos));
225
226 pos.getLongitude().set(0.0, ArcUnits.DEGREE);
227 pos.getLatitude().set(-90.0, ArcUnits.DEGREE);
228 System.out.println(" 0d, -90d => " + projector.getXyFor(pos));
229
230 System.out.println();
231
232 pos.getLongitude().set(180.0, ArcUnits.DEGREE);
233 pos.getLatitude().set(0.0, ArcUnits.DEGREE);
234 System.out.println(" 180d, 0d => " + projector.getXyFor(pos));
235
236 pos.getLongitude().set(-180.0, ArcUnits.DEGREE);
237 pos.getLatitude().set(0.0, ArcUnits.DEGREE);
238 System.out.println("-180d, 0d => " + projector.getXyFor(pos));
239
240 pos.getLongitude().set(360.0, ArcUnits.DEGREE);
241 pos.getLatitude().set(0.0, ArcUnits.DEGREE);
242 System.out.println(" 360d, 0d => " + projector.getXyFor(pos));
243
244 System.out.println();
245
246 pos.getLongitude().set(180.0, ArcUnits.DEGREE);
247 pos.getLatitude().set(90.0, ArcUnits.DEGREE);
248 System.out.println(" 180d, 90d => " + projector.getXyFor(pos));
249
250 pos.getLongitude().set(360.0, ArcUnits.DEGREE);
251 pos.getLatitude().set(90.0, ArcUnits.DEGREE);
252 System.out.println(" 360d, 90d => " + projector.getXyFor(pos));
253
254 System.out.println();
255
256 pos.getLongitude().set(0.0, ArcUnits.DEGREE);
257 pos.getLatitude().set(90.0, ArcUnits.DEGREE);
258 Point2D xy = projector.getXyFor(pos);
259 System.out.print(" 0d, 90d => " + xy);
260 pos = projector.getLatLonFor(xy);
261 System.out.println("; INVERSE => " + pos.getLongitude().toUnits(ArcUnits.DEGREE) + ", " + pos.getLatitude().toUnits(ArcUnits.DEGREE));
262
263 pos.getLongitude().set(90.0, ArcUnits.DEGREE);
264 pos.getLatitude().set(0.0, ArcUnits.DEGREE);
265 xy = projector.getXyFor(pos);
266 System.out.print(" 90d, 0d => " + xy);
267 pos = projector.getLatLonFor(xy);
268 System.out.println("; INVERSE => " + pos.getLongitude().toUnits(ArcUnits.DEGREE) + ", " + pos.getLatitude().toUnits(ArcUnits.DEGREE));
269
270 pos.getLongitude().set(30.0, ArcUnits.DEGREE);
271 pos.getLatitude().set(60.0, ArcUnits.DEGREE);
272 xy = projector.getXyFor(pos);
273 System.out.print(" 30d, 60d => " + xy);
274 pos = projector.getLatLonFor(xy);
275 System.out.println("; INVERSE => " + pos.getLongitude().toUnits(ArcUnits.DEGREE) + ", " + pos.getLatitude().toUnits(ArcUnits.DEGREE));
276
277 pos.getLongitude().set(45.0, ArcUnits.DEGREE);
278 pos.getLatitude().set(45.0, ArcUnits.DEGREE);
279 xy = projector.getXyFor(pos);
280 System.out.print(" 45d, 45d => " + xy);
281 pos = projector.getLatLonFor(xy);
282 System.out.println("; INVERSE => " + pos.getLongitude().toUnits(ArcUnits.DEGREE) + ", " + pos.getLatitude().toUnits(ArcUnits.DEGREE));
283
284 pos.getLongitude().set(60.0, ArcUnits.DEGREE);
285 pos.getLatitude().set(30.0, ArcUnits.DEGREE);
286 xy = projector.getXyFor(pos);
287 System.out.print(" 60d, 30d => " + xy);
288 pos = projector.getLatLonFor(xy);
289 System.out.println("; INVERSE => " + pos.getLongitude().toUnits(ArcUnits.DEGREE) + ", " + pos.getLatitude().toUnits(ArcUnits.DEGREE));
290 }
291 */
292 }