Polar.cpp
001 /**
002 *
003 * Polar 投影(扫描方式，自正北方向顺时针)
004 *
005 * PACKAGE:
006 * FILENAME: Polar.cpp
007 * LANGUAGE: C++
008 * ORIGINAL: Java2 v1.4
009 * DESCRIPTION: 极坐标投影(主要用于雷达图像处理)
010 * RELATED: Lambert.cpp Lambert.h
011 * EDITOR: UltraEdit-32 v12.20a(Windows) NEdit(Linux)
012 * CREATE: 2007-05-06 20:08:23
013 * UPDATE: 2007-06-06 改写为 C++ 版
014 * AUTHOR: 刘泽军 (BJ0773@126.com)
015 * 广西气象减灾研究所
016 * Guangxi Institude of Meteorology and Disaster-reducing Research(GIMDR)
017 *
018 * Compile : g++ -c Polar.cpp
019 *
020 * How to use Polar class:
021 *
022 * Polar polar = new Polar(Point(240, 240), 109.24, 24.35, 1.5);//构造函数
023 * polar->setScale(1.0);//设置比例尺，1公里对应1个像素点
024 * ...
025 *
026 **/
027
028 #include "Polar.h"
029
030 /**
031 *
032 * 扫描平面
033 * /
034 * /
035 * /
036 * /
037 * /
038 * / 仰角
039 * -------------------- 0度平面
040 *
041 * 如图所示：
042 * 扫描平面=>0度平面需要乘以cos(仰角)
043 * 0度平面=>扫描平面需要除以cos(仰角)
044 *
045 * 注意，日常显示的雷达图是扫描平面上的图。本类所说的屏幕指扫描平面。
046 *
047 */
048
049 /**
050 * 雷达扫描示意图
051 *
052 * 359 0
053 * | radius
054 * | /
055 * | /
056 * |angle/
057 * | /
058 * | ^ /
059 * | /
060 * | /
061 * |/
062 * 270 -----------------中心----------------- 90
063 * |
064 * |
065 * |
066 * |
067 * |
068 * |
069 * |
070 * |
071 * |
072 * 180
073 */
074
075 /**
076 * 功能：角度转换继弧度
077 * 参数：
078 * degrees - 角度值
079 * 返回值：
080 * 弧度值
081 */
082 double Polar::toRadians ( double degrees ) {
083 return ( PI/ 180.0 *degrees ) ;
084 }
085
086 /**
087 * 功能：角度转换继弧度
088 * 参数：
089 * degrees - 角度值
090 * 返回值：
091 * 弧度值
092 */
093 double Polar::toDegrees ( double radians ) {
094 return ( radians* 180.0 /PI ) ;
095 }
096
097 /**
098 * 功能：计算球面上两点间的距离(单位：公里)，原在edu.gimdr.Atmos.Meteorology类中写有，为避免import过多的类，故重写一份
099 * 参数：
100 * lon1,lat1 - 第1点的位置(经纬度)
101 * lon2,lat2 - 第2点的位置(经纬度)
102 * 返回值：
103 * 球面距离
104 */
105 double Polar::distanceOfSphere ( double lon1, double lat1, double lon2, double lat2 ) {
106 //A(x,y)
107 //B(a,b)
108 //AB球面距离=R*{arccos[cos(b)*cos(y)*cos(a-x)+sin(b)*sin(y)]}, by Google
109
110 double rlon1 = toRadians ( lon1 ) ;
111 double rlat1 = toRadians ( lat1 ) ;
112 double rlon2 = toRadians ( lon2 ) ;
113 double rlat2 = toRadians ( lat2 ) ;
114 double val = acos ( cos ( rlat2 ) *cos ( rlat1 ) *cos ( rlon2-rlon1 ) +sin ( rlat2 ) *sin ( rlat1 )) ;
115 val *= RADIUS;
116 return ( val ) ;
117 }
118
119 /**
120 * 功能：构造函数
121 * 参数：
122 * pos - 中心对应的屏幕位置
123 * lon - 中心对应的经度坐标
124 * lat - 中心对应的纬度坐标
125 * 返回值：
126 * 无
127 */
128 Polar::Polar ( Point pos, double lon, double lat ) {
129 elevation = 0.0 ; //无仰角
130 cosineElevation = 1.0 ;
131 setCenterPosition ( pos ) ;
132 setCenterCoordinate ( lon, lat ) ;
133 }
134
135 /**
136 * 功能：构造函数
137 * 参数：
138 * pos - 中心对应的屏幕位置
139 * lon - 中心对应的经度坐标
140 * lat - 中心对应的纬度坐标
141 * elv = 仰角
142 * 返回值：
143 * 无
144 */
145 Polar::Polar ( Point pos, double lon, double lat, double elv ) {
146 elevation = fabs ( elv ) ;
147 while ( elevation >= 90.0 ) { //在0-90度之间，但不能为90度
148 elevation -= 90.0 ;
149 }
150 cosineElevation = cos ( elevation ) ; //仰角的余弦值
151 setCenterPosition ( pos ) ;
152 setCenterCoordinate ( lon, lat ) ;
153 }
154
155 /**
156 * 功能：获得极坐标中心位置
157 * 参数：
158 * 无
159 * 返回值：
160 * 极坐标中心对应的屏幕位置
161 */
162 Point Polar::getCenterPosition () {
163 return ( centerPosition ) ;
164 }
165
166 /**
167 * 功能：设置极坐标的中心位置(屏幕坐标)
168 * 参数：
169 * pos - 新的中心位置(屏幕坐标)
170 * 返回值：
171 * 无
172 */
173 void Polar::setCenterPosition ( Point pos ) {
174 centerPosition = pos;
175 }
176
177 /**
178 * 功能：获得极坐标中心对应的经度坐标
179 * 参数：
180 * 无
181 * 返回值：
182 * 极坐标中心对应的经度坐标
183 */
184 double Polar::getCenterLongitude () {
185 return ( centerLongitude ) ;
186 }
187
188 /**
189 * 功能：获得极坐标中心对应的纬度坐标
190 * 参数：
191 * 无
192 * 返回值：
193 * 极坐标中心对应的纬度坐标
194 */
195 double Polar::getCenterLatitude () {
196 return ( centerLatitude ) ;
197 }
198
199 /**
200 * 功能：设置极坐标的中心位置(经纬度坐标)
201 * 参数：
202 * lon - 新的中心位置(经度值)
203 * lat - 新的中心位置(纬度值)
204 * 返回值：
205 * 无
206 */
207 void Polar::setCenterCoordinate ( double lon, double lat ) {
208 centerLongitude = lon;
209 centerLatitude = lat;
210 //中心经纬度或仰角发生改变，必须重新计算经向和纬向的1度对应的球面距离
211 kmPerDegreeX = distanceOfSphere ( centerLongitude, centerLatitude, centerLongitude+ 1.0 , centerLatitude ) / cosineElevation;
212 kmPerDegreeY = distanceOfSphere ( centerLongitude, centerLatitude, centerLongitude, centerLatitude+ 1.0 ) / cosineElevation;
213 }
214
215 /**
216 * 功能：获得仰角
217 * 参数：
218 * 无
219 * 返回值：
220 * 仰角的度数
221 */
222 double Polar::getElevation () {
223 return ( toDegrees ( elevation )) ;
224 }
225
226 /**
227 * 功能：设置仰角
228 * 参数：
229 * elv - 新的仰角
230 * 返回值：
231 * 无
232 */
233 void Polar::setElevation ( double elv ) {
234 elevation = fabs ( elv ) ;
235 while ( elevation >= 90.0 ) { //在0-90度之间，但不能为90度
236 elevation -= 90.0 ;
237 }
238 cosineElevation = cos ( elevation ) ; //仰角的余弦值
239 //中心经纬度或仰角发生改变，必须重新计算经向和纬向的1度对应的球面距离
240 kmPerDegreeX = distanceOfSphere ( centerLongitude, centerLatitude, centerLongitude+ 1.0 , centerLatitude ) / cosineElevation;
241 kmPerDegreeY = distanceOfSphere ( centerLongitude, centerLatitude, centerLongitude, centerLatitude+ 1.0 ) / cosineElevation;
242 }
243
244 /**
245 * 功能：获得比例尺，即1公里对应的像素点数
246 * 参数：
247 * 无
248 * 返回值：
249 * 比例尺
250 */
251 double Polar::getScale () {
252 return ( perKilometer ) ;
253 }
254
255 /**
256 * 功能：设置比例尺，即1公里对应的像素点数
257 * 参数：
258 * value - 比例尺数值
259 * 返回值：
260 * 无
261 */
262 void Polar::setScale ( double value ) {
263 perKilometer = value;
264 }
265
266 /**
267 * 功能：获得经纬度对应的屏幕像素坐标，与雷达仰角有关，主要用于体扫数据显示、底图叠加等。
268 * 参数：
269 * lon - 经度
270 * lat - 纬度
271 * 返回值：
272 * 对应的屏幕坐标
273 */
274 Point Polar::getPosition ( double lon, double lat ) {
275 double disX = distanceOfSphere ( lon, centerLatitude, centerLongitude, centerLatitude ) /cosineElevation;
276 double disY = distanceOfSphere ( centerLongitude, lat, centerLongitude, centerLatitude ) /cosineElevation;
277 return ( Point ( centerPosition.x+ ( lon>centerLongitude? 1 :- 1 ) * ( int )( disX*perKilometer ) ,centerPosition.y- ( lat>centerLatitude? 1 :- 1 ) * ( int )( disY*perKilometer ))) ;
278 }
279
280 /**
281 * 功能：获得极坐标对应的屏幕像素坐标，与雷达仰角无关，主要用于体扫数据显示、底图叠加等。
282 * 参数：
283 * radius - 极半径
284 * angle - 角度(以正北方向顺时针)
285 * 返回值：
286 * 对应的屏幕坐标
287 */
288
289 Point Polar::getXY ( double radius, double angle ) {
290 int x = ( int )( radius * sin ( toRadians ( angle ))) ;
291 int y = ( int )( radius * cos ( toRadians ( angle ))) ;
292 return ( Point ( centerPosition.x+x, centerPosition.y-y )) ;
293 }
294
295 /**
296 * 功能：获得屏幕像素点位置的极坐标半径，由于是输入参数是扫描平面上的值，故与雷达仰角无关。
297 * 参数：
298 * x - 水平坐标
299 * y - 垂直坐标
300 * 返回值：
301 * 与极坐标中心的距离，即极半径
302 */
303 double Polar::getRadius ( int x, int y ) {
304 return ( sqrt ( 1.0 * ( x-centerPosition.x ) * ( x-centerPosition.x ) + 1.0 * ( y-centerPosition.y ) * ( y-centerPosition.y ))) ;
305 }
306
307 /**
308 * 功能：获得经纬度位置的极坐标半径，与雷达仰角有关。
309 * 参数：
310 * lon - 经度坐标
311 * lat - 纬度坐标
312 * 返回值：
313 * 与极坐标中心的距离(象素点)，即极半径
314 */
315 double Polar::getRadius ( double lon, double lat ) {
316 Point pos = getPosition ( lon, lat ) ; //此函数已经考虑了仰角的影响
317 return ( getRadius ( pos.x, pos.y )) ;
318 }
319
320 /**
321 * 功能：获得屏幕像素点位置的极坐标角度(扫描平面与0度平面均相同)，与雷达仰角无关。
322 * 参数：
323 * x - 水平坐标
324 * y - 垂直坐标
325 * 返回值：
326 * 角度值，自正北方向顺时针
327 */
328 double Polar::getAngle ( int x, int y ) {
329 double agl = 0.0 ;
330 if ( x == centerPosition.x && y == centerPosition.y ) { //重合
331 agl = 0.0 ;
332 }
333 else if ( x == centerPosition.x ) {
334 agl = y > centerPosition.y ? 180.0 : 360.0 ;
335 }
336 else if ( y == centerPosition.y ) {
337 agl = x > centerPosition.x ? 90.0 : 270.0 ;
338 }
339 else {
340 agl = toDegrees ( atan ( 1.0 *fabs ( x-centerPosition.x ) /fabs ( y-centerPosition.y ))) ;
341 agl =
342 x > centerPosition.x && y < centerPosition.y ? agl : //直角坐标的第一象限
343 x < centerPosition.x && y < centerPosition.y ? 180.0 - agl : //直角坐标的第二象限
344 x < centerPosition.x && y > centerPosition.y ? 180.0 + agl : //直角坐标的第三象限
345 x > centerPosition.x && y > centerPosition.y ? 360.0 - agl : //直角坐标的第四象限
346 agl;
347 }
348 return ( agl ) ;
349 }
350
351 /**
352 * 功能：获得经纬度位置的极坐标角度(扫描平面与0度平面均相同)，与雷达仰角无关。
353 * 参数：
354 * lon - 水平坐标
355 * lat - 垂直坐标
356 * 返回值：
357 * 角度值，自正北方向顺时针
358 */
359 double Polar::getAngle ( double lon, double lat ) {
360 /*
361 //若通过获得屏幕坐标来计算角度，精度比较差，特别是在极坐标中心附近
362 Point p = getPosition(lon, lat);
363 return(getAngle(p.x, p.y);
364 */
365 double agl = 0.0 ;
366 if ( lon == centerLongitude && lat == centerLatitude ) { //重合
367 agl = 0.0 ;
368 }
369 else if ( lon == centerLongitude ) {
370 agl = lat > centerLatitude ? 360.0 : 180.0 ;
371 }
372 else if ( lat == centerLatitude ) {
373 agl = lon > centerLongitude ? 90.0 : 270.0 ;
374 }
375 else {
376 //注：由于经向和纬向的球面距离不等(华南，经向>纬向)，故点(1,1)与中心点(0,0)的极角不等45度，而应是略大于45度
377 agl = toDegrees ( atan (( fabs ( lon-centerLongitude ) *kmPerDegreeX ) / ( fabs ( lat-centerLatitude ) *kmPerDegreeY ))) ;
378 agl =
379 lon > centerLongitude && lat > centerLatitude ? agl : //第一象限
380 lon < centerLongitude && lat > centerLatitude ? 180.0 - agl : //第二象限
381 lon < centerLongitude && lat < centerLatitude ? 180.0 + agl : //第三象限
382 lon > centerLongitude && lat < centerLatitude ? 360.0 - agl : //第四象限
383 agl;
384 }
385 //System.out.println(agl);
386 return ( agl ) ;
387 }
388
389 /**
390 * 功能：获得屏幕坐标对应的经度值(根据目标点的经向球面距离来计算，雷达南面和北面的值略有差别)，与雷达仰角有关。
391 * 主要用于雷达产品的定位、底图叠加、转换为经纬度网格产品、拼图等。
392 * 参数：
393 * x - 横坐标，自西向东
394 * y - 纵坐标，自北向南
395 * 返回值：
396 * 对应的经度坐标
397 */
398 double Polar::getLongitude ( int x, int y ) {
399 double lat = getLatitude ( x, y ) ;
400 double disX0 = distanceOfSphere ( centerLongitude, lat, centerLongitude+ 1.0 , lat ) ; //0度平面上1经度的球面距离
401 double disX = disX0 / cosineElevation; //扫描平面上1经度的距离
402 double perDegreeX = disX * perKilometer; //扫描平面上1经度的对应的像素点数
403 return ( centerLongitude + ( x - centerPosition.x ) / perDegreeX ) ;
404 }
405
406 /**
407 * 功能：获得屏幕坐标对应的纬度值(根据极坐标中心点的纬向球面距离来计算)，与雷达仰角有关。
408 * 主要用于雷达产品的定位、底图叠加、转换为经纬度网格产品、拼图等。
409 * 参数：
410 * x - 横坐标，自西向东(未用到)
411 * y - 纵坐标，自北向南
412 * 返回值：
413 * 对应的纬度坐标
414 */
415 double Polar::getLatitude ( int x, int y ) {
416 /**
417 * 目标点 A(X,Y) 弧度
418 * 中心点 B(A,B) 弧度
419 * AB球面距离=R*{arccos[cos(B)*cos(Y)*cos(A-X)+sin(B)*sin(Y)]}, by Google
420 * 经度相同 => AB = R*{arccos[cos(B)*cos(Y)+sin(B)*sin(Y)]}
421 * => AB = R*{arccos[cos(B-Y)]}
422 * => AB = R * (B-Y)
423 * => AB / R = B - Y
424 * => Y = B - AB /R
425 * => Y = B - (y-centerPosition.y)/perKilometer/R
426 */
427 double val = toRadians ( centerLatitude ) + ( centerPosition.y-y ) *cosineElevation/perKilometer;
428 double fRadius = RADIUS;
429 val = val / fRadius;
430 return ( toDegrees ( val )) ;
431 }
Polar.h
01 #ifndef PolarH
02 #define PolarH
03
04 #include <stdio.h>
05 #include <stdlib.h>
06 #include <math.h>
07
08 //静态常量，地球半径，来源：《大气科学常用公式》，P601，附录
09 #define RADIUS 6371.004 ; //地球平均半径，单位：公里(Km)。
10 #define RADIUS_POLAR 6356.755 ; //地球两极半径，单位：公里(Km)。
11 #define RADIUS_EQUATOR 6373.140 ; //地球赤道半径，单位：公里(Km)。
12
13 #define PI 3.14159265358979323846264338327950288
14
15 typedef struct _tagPOINT {
16 int x;
17 int y;
18 _tagPOINT () {}
19 _tagPOINT ( int _x, int _y ) : x ( _x ) , y ( _y ) {}
20 } Point;
21
22 class Polar {
23
24 private : //私有成员
25 Point centerPosition; //中心对应的屏幕位置
26 double centerLongitude; //中心经度
27 double centerLatitude; //中心纬度
28 double perKilometer; //比例尺：一公里对应的像素点数(扫描平面)
29 double elevation; //仰角
30 double cosineElevation; //仰角的余弦值
31 double kmPerDegreeX; //1经度对应的距离(公里)，不同纬度数值不同
32 double kmPerDegreeY; //1纬度对应的距离(公里)，不同纬度数值不同
33
34 public : //公共成员
35 //1、角度=>板度
36 double toRadians ( double degrees ) ;
37 //2、弧度=>角度
38 double toDegrees ( double radians ) ;
39 //3、计算球面距离
40 double distanceOfSphere ( double lon1, double lat1, double lon2, double lat2 ) ;
41 //4、构造函数，不考虑仰角
42 Polar ( Point pos, double lon, double lat ) ;
43 //5、构造函数，考虑仰角
44 Polar ( Point pos, double lon, double lat, double elv ) ;
45 //6、获得极坐标中心点的屏幕位置
46 Point getCenterPosition () ;
47 //7、设置极坐标中心点的屏幕位置
48 void setCenterPosition ( Point pos ) ;
49 //8、获得极坐标中心点的经度
50 double getCenterLongitude () ;
51 //9、获得极坐标中心点的纬度
52 double getCenterLatitude () ;
53 //10、设置极坐标中心点的经纬度
54 void setCenterCoordinate ( double lon, double lat ) ;
55 //11、获得仰角
56 double getElevation () ;
57 //12、设置仰角
58 void setElevation ( double elv ) ;
59 //13、获得缩放比例尺
60 double getScale () ;
61 //14、设置缩放比例尺
62 void setScale ( double value ) ;
63 //15、获得经纬度对应的屏幕坐标(扫描平面)，主要用于体扫数据显示、底图叠加等。
64 Point getPosition ( double lon, double lat ) ;
65 //16、获得极坐标对应的屏幕坐标(扫描平面)，主要用于体扫数据显示、底图叠加等。
66 Point getXY ( double radius, double angle ) ;
67 //17、根据屏幕坐标获得极半径
68 double getRadius ( int x, int y ) ;
69 //18、根据经纬度坐标获得
70 double getRadius ( double lon, double lat ) ;
71 //19、根据屏幕坐标获得极角
72 double getAngle ( int x, int y ) ;
73 //20、根据经纬度坐标获得极角
74 double getAngle ( double lon, double lat ) ;
75 //21、根据屏幕坐标获得对应的经度值，主要用于雷达产品的定位、底图叠加、转换为经纬度网格产品、拼图等。
76 double getLongitude ( int x, int y ) ;
77 //22、根据屏幕坐标获得对应的纬度值，主要用于雷达产品的定位、底图叠加、转换为经纬度网格产品、拼图等。
78 double getLatitude ( int x, int y ) ;
79
80 } ;
81
82 #endif
5709
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
