Polar极坐标投影(Java)

Polar极坐标地图投影，主要用于714数字雷达、多普勒天气雷达图像的显示、处理，以及屏幕坐标(像素点)与经纬度坐标的转换。

---

Polar.java

---

001 /**
002  *
003  *  Polar 投影(扫描方式，自正北方向顺时针)
004  *
005  *     PACKAGE: cma.common.projection
006  *    FILENAME: Polar.java
007  *    LANGUAGE: Java2 v1.4
008  *    ORIGINAL: 无
009  * DESCRIPTION: 极坐标投影(主要用于雷达图像处理)
010  *     RELATED: cma.common.projection.Lambert(兰勃特投影)
011  *      EDITOR: UltraEdit-32 v12.20a(Windows) NEdit(Linux)
012  *      CREATE: 2007-05-06 20:08:23
013  *      UPDATE: 2007-06-09
014  *      AUTHOR: 刘泽军 ([e-mail])
015  *              广西气象减灾研究所
016  *    Guangxi Institude of Meteorology and Disaster-reducing Research(GIMDR)
017  *
018  *    Compile : javac Polar.java
019  *
020  * How to use Polar class:
021  *
022  * Polar polar = new Polar(new Point(240, 240), 109.24, 24.35, 1.5);//构造函数
023  * polar.setScale(1.0);//设置比例尺，1公里对应1个像素点
024  * ...
025  *
026  */
027 
028 /**
029  *
030  *         扫描平面
031  *            /
032  *           /
033  *          /
034  *         /
035  *        /
036  *       /  仰角
037  *      -------------------- 0度平面
038  *
039  * 如图所示：
040  *          扫描平面=>0度平面，需要乘以cos(仰角)
041  *          0度平面=>扫描平面，需要除以cos(仰角)
042  *
043  * 注意，日常显示的雷达图是扫描平面上的图。本类所说的屏幕指扫描平面。
044  *
045  */
046 
047 /**
048  * 雷达扫描示意图
049  *
050  *                    359 0
051  *                        |     radius
052  *                        |       /
053  *                        |      /
054  *                        |angle/
055  *                        |    /
056  *                        | ^ /
057  *                        |  /
058  *                        | /
059  *                        |/
060  * 270 -----------------中心----------------- 90
061  *                        |
062  *                        |
063  *                        |
064  *                        |
065  *                        |
066  *                        |
067  *                        |
068  *                        |
069  *                        |
070  *                       180
071  */
072 
073 /**
074  *  Polar.java的函数列表
075  *  1、计算球面距离
076  *  double distanceOfSphere(double lon1, double lat1, double lon2, double lat2);
077  *  2、构造函数，不考虑仰角
078  *  Polar(Point pos, double lon, double lat);
079  *  3、构造函数，考虑仰角
080  *  Polar(Point pos, double lon, double lat, double elv);
081  *  4、获得极坐标中心点的屏幕位置
082  *  getCenterPosition();
083  *  5、设置极坐标中心点的屏幕位置
084  *  setCenterPosition(Point pos);
085  *  6、获得极坐标中心点的经度
086  *  getCenterLongitude();
087  *  7、获得极坐标中心点的纬度
088  *  getCenterLatitude();
089  *  8、设置极坐标中心点的经纬度
090  *  setCenterCoordinate(double lon, double lat);
091  *  9、获得仰角
092  *  getElevation();
093  *  10、设置仰角
094  *  setElevation(double elv);
095  *  11、获得缩放比例尺
096  *  getScale();
097  *  12、设置缩放比例尺
098  *  setScale(double value);
099  *  13、获得经纬度对应的屏幕坐标(扫描平面)，主要用于体扫数据显示、底图叠加等。
100  *  getPosition(double lon, double lat);
101  *  14、获得极坐标对应的屏幕坐标(扫描平面)，主要用于体扫数据显示、底图叠加等。
102  *  getXY(double radius, double angle);
103  *  15、根据屏幕坐标获得极半径
104  *  getRadius(int x, int y);
105  *  16、根据经纬度坐标获得
106  *  getRadius(double lon, double lat);
107  *  17、根据屏幕坐标获得极角
108  *  getAngle(int x, int y);
109  *  18、根据经纬度坐标获得极角
110  *  getAngle(double lon, double lat);
111  *  19、根据屏幕坐标获得对应的经度值，主要用于雷达产品的定位、底图叠加、转换为经纬度网格产品、拼图等。
112  *  getLongitude(int x, int y);
113  *  20、根据屏幕坐标获得对应的纬度值，主要用于雷达产品的定位、底图叠加、转换为经纬度网格产品、拼图等。
114  *  getLatitude(int x, int y);
115  *  21、获得屏幕坐标对应的经纬度
116  *  public Point2D.Double getCoordinate(int x, int y);
117  */
118 
119 package cma.common.projection;
120 
121 import java.awt.*;
122 import java.awt.geom.*;
123 import java.lang.Math.*;
124 
125 public class Polar {
126 
127     //静态常量，地球半径，来源：《大气科学常用公式》，P601，附录
128     public static double RADIUS         = 6371.004;//地球平均半径，单位：公里(Km)。
129     public static double RADIUS_POLAR   = 6356.755;//地球两极半径，单位：公里(Km)。
130     public static double RADIUS_EQUATOR = 6373.140;//地球赤道半径，单位：公里(Km)。
131 
132     //私有成员
133     private Point   centerPosition;             //中心对应的屏幕位置
134     private double  centerLongitude = 0.0;      //中心经度
135     private double  centerLatitude  = 0.0;      //中心纬度
136     private double  perKilometer    = 0.0;      //比例尺：一公里对应的像素点数(扫描平面)
137     private double  elevation       = 0.0;      //仰角
138     private double  cosineElevation = 1.0;      //仰角的余弦值
139     private double  kmPerDegreeX    = 0.0;      //1经度对应的距离(公里)，不同纬度数值不同
140     private double  kmPerDegreeY    = 0.0;      //1纬度对应的距离(公里)，不同纬度数值不同
141 
142 /**
143  * 功能：计算球面上两点间的距离(单位：公里)，原在edu.gimdr.Atmos.Meteorology类中写有，为避免import过多的类，故重写一份
144  * 参数：
145  *  lon1,lat1   - 第1点的位置(经纬度)
146  *  lon2,lat2   - 第2点的位置(经纬度)
147  * 返回值：
148  *      球面距离
149  */
150     public static double distanceOfSphere(double lon1, double lat1, double lon2, double lat2) {
151         /*  公式：
152             A(x,y)  B(a,b)
153             AB点的球面距离=R*{arccos[cos(b)*cos(y)*cos(a-x)+sin(b)*sin(y)]}, by Google
154         */
155 
156         double  rlon1   = Math.toRadians(lon1);
157         double  rlat1   = Math.toRadians(lat1);
158         double  rlon2   = Math.toRadians(lon2);
159         double  rlat2   = Math.toRadians(lat2);
160 
161         return(Polar.RADIUS*(Math.acos(Math.cos(rlat2)*Math.cos(rlat1)*Math.cos(rlon2-rlon1)+Math.sin(rlat2)*Math.sin(rlat1))));
162     }
163 
164 /**
165  * 功能：构造函数
166  * 参数：
167  *      pos     - 中心对应的屏幕位置
168  *      lon     - 中心对应的经度坐标
169  *      lat     - 中心对应的纬度坐标
170  * 返回值：
171  *      无
172  */
173     public Polar(Point pos, double lon, double lat) {
174         elevation       = 0.0;//无仰角
175         cosineElevation = 1.0;
176         setCenterPosition(pos);
177         setCenterCoordinate(lon, lat);
178     }
179 
180 /**
181  * 功能：构造函数
182  * 参数：
183  *      pos     - 中心对应的屏幕位置
184  *      lon     - 中心对应的经度坐标
185  *      lat     - 中心对应的纬度坐标
186  *      elv     = 仰角
187  * 返回值：
188  *      无
189  */
190     public Polar(Point pos, double lon, double lat, double elv) {
191         elevation       = Math.toRadians(Math.IEEEremainder(Math.abs(elv), 90.0));//在0-90度之间，但不能为90度
192         cosineElevation = Math.cos(elevation);//仰角的余弦值
193         setCenterPosition(pos);
194         setCenterCoordinate(lon, lat);
195     }
196 
197 /**
198  * 功能：获得极坐标中心位置
199  * 参数：
200  *      无
201  * 返回值：
202  *      极坐标中心对应的屏幕位置
203  */
204     public Point getCenterPosition() {
205         return(centerPosition);
206     }
207 
208 /**
209  * 功能：设置极坐标的中心位置(屏幕坐标)
210  * 参数：
211  *      pos     - 新的中心位置(屏幕坐标)
212  * 返回值：
213  *      无
214  */
215     public void setCenterPosition(Point pos) {
216         centerPosition  = pos;
217     }
218 
219 /**
220  * 功能：获得极坐标中心对应的经度坐标
221  * 参数：
222  *      无
223  * 返回值：
224  *      极坐标中心对应的经度坐标
225  */
226     public double getCenterLongitude() {
227         return(centerLongitude);
228     }
229 
230 /**
231  * 功能：获得极坐标中心对应的纬度坐标
232  * 参数：
233  *      无
234  * 返回值：
235  *      极坐标中心对应的纬度坐标
236  */
237     public double getCenterLatitude() {
238         return(centerLatitude);
239     }
240 
241 /**
242  * 功能：设置极坐标的中心位置(经纬度坐标)
243  * 参数：
244  *      lon     - 新的中心位置(经度值)
245  *      lat     - 新的中心位置(纬度值)
246  * 返回值：
247  *      无
248  */
249     public void setCenterCoordinate(double lon, double lat) {
250         centerLongitude = lon;
251         centerLatitude  = lat;
252         //中心经纬度或仰角发生改变，必须重新计算经向和纬向的1度对应的球面距离
253         kmPerDegreeX    = distanceOfSphere(centerLongitude, centerLatitude, centerLongitude+1.0, centerLatitude) / cosineElevation;
254         kmPerDegreeY    = distanceOfSphere(centerLongitude, centerLatitude, centerLongitude, centerLatitude+1.0) / cosineElevation;
255     }
256 
257 /**
258  * 功能：获得仰角
259  * 参数：
260  *      无
261  * 返回值：
262  *      仰角的度数
263  */
264     public double getElevation() {
265         return(Math.toDegrees(elevation));
266     }
267 
268 /**
269  * 功能：设置仰角
270  * 参数：
271  *      elv     - 新的仰角
272  * 返回值：
273  *      无
274  */
275     public void setElevation(double elv) {
276         elevation       = Math.toRadians(Math.IEEEremainder(Math.abs(elv), 90.0));//在0-90度之间，但不能为90度
277         cosineElevation = Math.cos(elevation);//仰角的余弦值
278         //中心经纬度或仰角发生改变，必须重新计算经向和纬向的1度对应的球面距离
279         kmPerDegreeX    = distanceOfSphere(centerLongitude, centerLatitude, centerLongitude+1.0, centerLatitude) / cosineElevation;
280         kmPerDegreeY    = distanceOfSphere(centerLongitude, centerLatitude, centerLongitude, centerLatitude+1.0) / cosineElevation;
281     }
282 
283 /**
284  * 功能：获得比例尺，即1公里对应的像素点数
285  * 参数：
286  *      无
287  * 返回值：
288  *      比例尺
289  */
290     public double getScale() {
291         return(perKilometer);
292     }
293 
294 /**
295  * 功能：设置比例尺，即1公里对应的像素点数
296  * 参数：
297  *      value   - 比例尺数值
298  * 返回值：
299  *      无
300  */
301     public void setScale(double value) {
302         perKilometer    = value;
303     }
304 
305 /**
306  * 功能：获得经纬度对应的屏幕像素坐标，与雷达仰角有关，主要用于体扫数据显示、底图叠加等。
307  * 参数：
308  *      lon - 经度
309  *      lat - 纬度
310  * 返回值：
311  *      对应的屏幕坐标
312  */
313     public Point getPosition(double lon, double lat) {
314         double  disX    = distanceOfSphere(lon, centerLatitude, centerLongitude, centerLatitude)/cosineElevation;
315         double  disY    = distanceOfSphere(centerLongitude, lat, centerLongitude, centerLatitude)/cosineElevation;
316         return(new Point(centerPosition.x+(lon>centerLongitude?1:-1)*(int)(disX*perKilometer),centerPosition.y-(lat>centerLatitude?1:-1)*(int)(disY*perKilometer)));
317     }
318 
319 /**
320  * 功能：获得极坐标对应的屏幕像素坐标，与雷达仰角无关，主要用于体扫数据显示、底图叠加等。
321  * 参数：
322  *      radius      - 极半径
323  *      angle       - 角度(以正北方向顺时针)
324  * 返回值：
325  *      对应的屏幕坐标
326  */
327 
328     public Point getXY(double radius, double angle) {
329         int x   = (new Double(radius * Math.sin(Math.toRadians(angle)))).intValue();
330         int y   = (new Double(radius * Math.cos(Math.toRadians(angle)))).intValue();
331         return(new Point(centerPosition.x+x, centerPosition.y-y));
332     }
333 
334 /**
335  * 功能：获得屏幕像素点位置的极坐标半径，由于是输入参数是扫描平面上的值，故与雷达仰角无关。
336  * 参数：
337  *      x   - 水平坐标
338  *      y   - 垂直坐标
339  * 返回值：
340  *      与极坐标中心的距离，即极半径
341  */
342     public double getRadius(int x, int y) {
343         return(Math.sqrt(1.0*(x-centerPosition.x)*(x-centerPosition.x)+1.0*(y-centerPosition.y)*(y-centerPosition.y)));
344     }
345 
346 /**
347  * 功能：获得经纬度位置的极坐标半径，与雷达仰角有关。
348  * 参数：
349  *      lon - 经度坐标
350  *      lat - 纬度坐标
351  * 返回值：
352  *      与极坐标中心的距离(象素点)，即极半径
353  */
354     public double getRadius(double lon, double lat) {
355         Point   pos = getPosition(lon, lat);//此函数已经考虑了仰角的影响
356         return(getRadius(pos.x, pos.y));
357     }
358 
359 /**
360  * 功能：获得屏幕像素点位置的极坐标角度(扫描平面与0度平面均相同)，与雷达仰角无关。
361  * 参数：
362  *      x   - 水平坐标
363  *      y   - 垂直坐标
364  * 返回值：
365  *      角度值，自正北方向顺时针
366  */
367     public double getAngle(int x, int y) {
368         double  agl = 0.0;
369         if( x == centerPosition.x && y == centerPosition.y ) {//重合
370             agl = 0.0;
371         }
372         else if( x == centerPosition.x ) {
373             agl = y > centerPosition.y ? 180.0 : 360.0;
374         }
375         else if( y == centerPosition.y ) {
376             agl = x > centerPosition.x ? 90.0 : 270.0;
377         }
378         else {
379             agl = Math.toDegrees(Math.atan(1.0*Math.abs(x-centerPosition.x)/Math.abs(y-centerPosition.y)));
380             agl =
381                     x > centerPosition.x && y < centerPosition.y ? agl :            //直角坐标的第一象限
382                     x < centerPosition.x && y < centerPosition.y ? 180.0 - agl :    //直角坐标的第二象限
383                     x < centerPosition.x && y > centerPosition.y ? 180.0 + agl :    //直角坐标的第三象限
384                     x > centerPosition.x && y > centerPosition.y ? 360.0 - agl :    //直角坐标的第四象限
385                     agl;
386         }
387         System.out.println(agl);
388         return(agl);
389     }
390 
391 /**
392  * 功能：获得经纬度位置的极坐标角度(扫描平面与0度平面均相同)，与雷达仰角无关。
393  * 参数：
394  *      lon - 水平坐标
395  *      lat - 垂直坐标
396  * 返回值：
397  *      角度值，自正北方向顺时针
398  */
399     public double getAngle(double lon, double lat) {
400 /*
401 //若通过获得屏幕坐标来计算角度，精度比较差，特别是在极坐标中心附近
402         Point   p   = getPosition(lon, lat);
403         return(getAngle(p.x, p.y);
404 */
405         double  agl = 0.0;
406         if( lon == centerLongitude && lat == centerLatitude ) {//重合
407             agl = 0.0;
408         }
409         else if( lon == centerLongitude ) {
410             agl = lat > centerLatitude ? 360.0 : 180.0;
411         }
412         else if( lat == centerLatitude ) {
413             agl = lon > centerLongitude ? 90.0 : 270.0;
414         }
415         else {
416             //注：由于经向和纬向的球面距离不等(华南，经向>纬向)，故点(1,1)与中心点(0,0)的极角不等45度，而应是略大于45度
417             agl = Math.toDegrees(Math.atan((Math.abs(lon-centerLongitude)*kmPerDegreeX)/(Math.abs(lat-centerLatitude)*kmPerDegreeY)));
418             agl =
419                     lon > centerLongitude && lat > centerLatitude ? agl :           //第一象限
420                     lon < centerLongitude && lat > centerLatitude ? 180.0 - agl :   //第二象限
421                     lon < centerLongitude && lat < centerLatitude ? 180.0 + agl :   //第三象限
422                     lon > centerLongitude && lat < centerLatitude ? 360.0 - agl :   //第四象限
423                     agl;
424         }
425         return(agl);
426     }
427 
428 /**
429  * 功能：获得屏幕坐标对应的经度值(根据目标点的经向球面距离来计算，雷达南面和北面的值略有差别)，与雷达仰角有关。
430  *       主要用于雷达产品的定位、底图叠加、转换为经纬度网格产品、拼图等。
431  * 参数：
432  *      x   - 横坐标，自西向东
433  *      y   - 纵坐标，自北向南
434  * 返回值：
435  *      对应的经度坐标
436  */
437     public double getLongitude(int x, int y) {
438         double  lat         = getLatitude(x, y);
439         double  disX0       = distanceOfSphere(centerLongitude, lat, centerLongitude+1.0, lat);//0度平面上1经度的球面距离
440         double  disX        = disX0 / cosineElevation;      //扫描平面上1经度的距离
441         double  perDegreeX  = disX * perKilometer;          //扫描平面上1经度的对应的像素点数
442         return(centerLongitude + (x - centerPosition.x) / perDegreeX);
443     }
444 
445 /**
446  * 功能：获得屏幕坐标对应的纬度值(根据极坐标中心点的纬向球面距离来计算)，与雷达仰角有关。
447  *       主要用于雷达产品的定位、底图叠加、转换为经纬度网格产品、拼图等。
448  * 参数：
449  *      x   - 横坐标，自西向东(未用到)
450  *      y   - 纵坐标，自北向南
451  * 返回值：
452  *      对应的纬度坐标
453  */
454     public double getLatitude(int x, int y) {
455         /*
456             目标点 A(X,Y) 弧度
457             中心点 B(A,B) 弧度
458             AB球面距离=R*{arccos[cos(B)*cos(Y)*cos(A-X)+sin(B)*sin(Y)]}, by Google
459             经度相同 => AB = R*{arccos[cos(B)*cos(Y)+sin(B)*sin(Y)]}
460             => AB = R*{arccos[cos(B-Y)]}
461             => AB = R * (B-Y)
462             => AB / R = B - Y
463             => Y = B - AB /R
464             => Y = B - (y-centerPosition.y)/perKilometer/R
465         */
466         return(Math.toDegrees(Math.toRadians(centerLatitude) + (centerPosition.y-y)*cosineElevation/perKilometer/Polar.RADIUS));
467     }
468 
469 /**
470  * 功能:
471  *      获得屏幕坐标对应的经纬度
472  * 参数:
473  *      x       - 屏幕水平坐标
474  *      y       - 屏幕垂直坐标
475  * 返回值:
476  *      对应的经纬度
477  */
478     public Point2D.Double getCoordinate(int x, int y) {
479         double  lat         = getLatitude(x, y);
480         double  disX0       = distanceOfSphere(centerLongitude, lat, centerLongitude+1.0, lat);//0度平面上1经度的球面距离
481         double  disX        = disX0 / cosineElevation;      //扫描平面上1经度的距离
482         double  perDegreeX  = disX * perKilometer;          //扫描平面上1经度的对应的像素点数
483         double  lon         = centerLongitude + (x - centerPosition.x) / perDegreeX;
484         return(new Point2D.Double(lon, lat));
485     }
486 
487 /**
488  * 功能：
489  *      画经线、纬线
490  * 参数：
491  *      g       - 图形设备
492  *      f       - 字体
493  *      c       - 画线颜色
494  *      inc_lon - 经线间隔
495  *      inc_lat - 纬线间隔
496  * 返回值：
497  *      无
498  */
499     public void drawGridLine(Graphics2D g, Font f, Color c, int inc_lon, int inc_lat) {
500         return;
501 /*
502         Color   saveColor   = g.getColor();
503         Font    saveFont    = g.getFont();
504         g.setFont(saveFont);
505         g.setColor(saveColor);
506 */
507     }
508 }