Polar坐标投影(C++)
来源：互联网 发布：数据备份管理办法 编辑：程序博客网 时间：2024/06/12 01:44
Polar.cpp
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
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
- Polar坐标投影(C++)
- Polar极坐标投影(Java)
- Polar(极坐标)投影--主要用于天气雷达图
- 坐标＆投影
- 坐标投影
- 图像的Log-Polar极坐标变换
- 坐标、投影及坐标转换
- 地理坐标与投影坐标
- Linear线性坐标投影
- GIS坐标系统、投影
- 地理坐标与投影
- OpenGL 坐标系统 投影
- 利用proj.net进行投影坐标变换（c#，北京54坐标）
- MATLAB描绘极坐标图像——polar
- Log-Polar——关于对数极坐标
- OpenCV学习之Log-Polar极坐标变换
- 由投影坐标计算地理坐标
- 由投影坐标计算地理坐标
- 求职面试自我介绍一分钟
- HTTP Code简表
- 责任
- 关于循环标号的问题
- 2007年7月10日 星期二 晴
- Polar坐标投影(C++)
- 关于GBA的道歉声明
- 简单与争用simpleness & contention
- 数据库下载漏洞攻击技术
- 从原理上解决Tomcat中文问题
- 如何利用JSP的9种基本内置组件
- Java Class Loader
- TOMCAT用https替换http的方法
- Linear线性坐标投影
