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: 刘泽军 (BJ0773@gmail.com)
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  }