from:http://hi.baidu.com/sungaoyong/item/0c4584d25873f131e3108f05
///刘泽军java版本的极坐标投影c#版本的移植
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Drawing;
namespace MyPolar
{
/// <summary>
/// 自定义的Math类，支持角度到弧度和弧度到角度的计算
/// 对应Java Math中的Math.toDegrees和Math.toRadians
/// </summary>
class CMath
{
public static Double PI = Math.PI;
static public Double toDegrees(Double rad)
{
return (rad * (180.0f / PI));
}
static public Double toRadians(Double deg)
{
return (deg * (PI / 180.0f));
}
public static double cos(double d)
{
return Math.Cos(d);
}
public static double acos(double d)
{
return Math.Acos(d);
}
public static double sin(double d)
{
return Math.Sin(d);
}
public static double abs(double d)
{
return Math.Abs(d);
}
public static double IEEERemainder(double x, double y)
{
return Math.IEEERemainder(x, y);
}
public static double sqrt(double d)
{
return Math.Sqrt(d);
}
public static double atan(double d)
{
return Math.Atan(d);
}
}
/*
Polar 投影(扫描方式，自正北方向顺时针)
PACKAGE: cma.common.projection
FILENAME: Polar.java
LANGUAGE: Java2 v1.4
ORIGINAL: 无
DESCRIPTION: 极坐标投影(主要用于雷达图像处理)
RELATED: cma.common.projection.Lambert(兰勃特投影)
EDITOR: UltraEdit-32 v12.20a(Windows) NEdit(Linux)
CREATE: 2007-05-06 20:08:23
UPDATE: 2007-07-18 修改为抽象类Coordinate的扩展类
AUTHOR: 刘泽军 ()
广西气象减灾研究所
Guangxi Institude of Meteorology and Disaster-reducing Research(GIMDR)
Compile : javac Coordinate.java Polar.java
How to use Polar class:
Polar polar = new Polar(109.24, 24.35, 512, 384, 1.0, 0.0);//构造函数
...
孙高勇2011-02-10移植到DotNet版本。
*/
/**
*
* 扫描平面
* /
* /
* /
* /
* /
* / 仰角
* -------------------- 0度平面
*
* 如图所示：
* 扫描平面=>0度平面，需要乘以cos(仰角)
* 0度平面=>扫描平面，需要除以cos(仰角)
*
* 注意，日常显示的雷达图是扫描平面上的图。本类所说的屏幕指扫描平面。
*
*/
/**
* 雷达扫描示意图
*
* 359 0
* | radius
* | /
* | /
* |angle/
* | /
* | ^ /
* | /
* | /
* |/
* 270 -----------------中心----------------- 90
* |
* |
* |
* |
* |
* |
* |
* |
* |
* 180
*/
class Polar
{
//静态常量，地球半径，来源：《大气科学常用公式》，P601，附录
public static double RADIUS = 6371.004;//地球平均半径，单位：公里(Km)。
public static double RADIUS_POLAR = 6356.755;//地球两极半径，单位：公里(Km)。
public static double RADIUS_EQUATOR = 6373.140;//地球赤道半径，单位：公里(Km)。
//私有成员
private PointF center; //中心对应的屏幕位置
private Point place; //中心经纬度对应的屏幕坐标
private Point offset; //偏移量
//缩放系数
private double scale;
private double scaleOriginal;
//private double centerLongitude = 0.0; //中心经度
//private double centerLatitude = 0.0; //中心纬度
private double perKilometer = 1.0; //比例尺：一公里对应的像素点数(扫描平面)
private double elevation = 0.0; //仰角
private double cosineElevation = 1.0; //仰角的余弦值
private double kmPerDegreeX = 1.0; //1经度对应的距离(公里)，不同纬度数值不同
private double kmPerDegreeY = 1.0; //1纬度对应的距离(公里)，不同纬度数值不同
/**
* 功能：计算球面上两点间的距离(单位：公里)，原在edu.gimdr.Atmos.Meteorology类中写有，为避免import过多的类，故重写一份
* 参数：
* lon1,lat1 - 第1点的位置(经纬度)
* lon2,lat2 - 第2点的位置(经纬度)
* 返回值：
* 球面距离
*/
public static double distanceOfSphere(double lon1, double lat1, double lon2, double lat2)
{
/* 公式：
A(x,y) B(a,b)
AB点的球面距离=R*{arccos[cos(b)*cos(y)*cos(a-x)+sin(b)*sin(y)]}, by Google
*/
double rlon1 = CMath.toRadians(lon1);
double rlat1 = CMath.toRadians(lat1);
double rlon2 = CMath.toRadians(lon2);
double rlat2 = CMath.toRadians(lat2);
return (RADIUS * (CMath.acos(CMath.cos(rlat2) * CMath.cos(rlat1) * CMath.cos(rlon2 - rlon1) + CMath.sin(rlat2) * CMath.sin(rlat1))));
}
/**
* 功能:
* 重置参数
* 参数:
* lon,lat - 中心经纬度，
* px,py - 中心经纬度对应的屏幕坐标
* sc - 缩放系数
* agl - 仰角
* 返回值:
* 无
*/
public void reset(double lon, double lat, int px, int py, double sc, double agl)
{
//type = POLAR;
center = new PointF(
(float)(lon < 0.0 ? 0.0 : lon > 360.0 ? 360.0 : lon),
(float)(lat < -90.0 ? -90.0 : lat > 90.0 ? 90.0 : lat)
);
place = new Point(px, py);
elevation = CMath.toRadians(CMath.IEEERemainder(CMath.abs(agl), 90.0));//在0-90度之间，但不能为90度
cosineElevation = CMath.cos(elevation);//仰角的余弦值
scale = sc == 0.0 ? 1.0 : CMath.abs(sc);//缩放系数
scaleOriginal = scale;
offset = new Point(0, 0);
perKilometer = 1.0;//标准比例尺
////中心经纬度或仰角发生改变，必须重新计算经向和纬向的1度对应的球面距离
kmPerDegreeX = distanceOfSphere(center.X, center.Y, center.X + 1.0, center.Y) / cosineElevation;
kmPerDegreeY = distanceOfSphere(center.X, center.Y, center.X, center.Y + 1.0) / cosineElevation;
}
/**
* 功能：构造函数
* 参数：
* lon - 中心对应的经度坐标
* lat - 中心对应的纬度坐标
* x - 中心对应的屏幕位置x
* y - 中心对应的屏幕位置y
* sc - 缩放系数
* 返回值：
* 无
*/
public Polar(double lon, double lat, int x, int y, double sc)
{
reset(lon, lat, x, y, sc, 0.0);
}
/**
* 功能：构造函数
* 参数：
* lon - 中心对应的经度坐标
* lat - 中心对应的纬度坐标
* x - 中心对应的屏幕位置x
* y - 中心对应的屏幕位置y
* sc - 缩放系数
* agl = 仰角
* 返回值：
* 无
*/
public Polar(double lon, double lat, int x, int y, double sc, double agl)
{
reset(lon, lat, x, y, sc, agl);
}
/**
* 功能：获得仰角
* 参数：
* 无
* 返回值：
* 仰角的度数
*/
public double getElevation()
{
return (CMath.toDegrees(elevation));
}
/**
* 功能：获得经纬度对应的屏幕像素坐标，与雷达仰角有关，主要用于体扫数据显示、底图叠加等。
* 参数：
* lon - 经度
* lat - 纬度
* 返回值：
* 对应的屏幕坐标
*/
public Point getPosition(double lon, double lat)
{
double disX = distanceOfSphere(lon, center.Y, center.X, center.Y) / cosineElevation;
double disY = distanceOfSphere(center.X, lat, center.X, center.Y) / cosineElevation;
double x = (lon > center.X ? 1 : -1) * (disX * perKilometer * scale) + place.X + 0.5;
double y = -(lat > center.Y ? 1 : -1) * (disY * perKilometer * scale) + place.Y + 0.5;
return (new Point((int)x, (int)y));
}
/**
* 功能：获得极坐标对应的屏幕像素坐标，与雷达仰角无关，主要用于体扫数据显示、底图叠加等。
* 参数：
* radius - 极半径
* angle - 角度(以正北方向顺时针)
* 返回值：
* 对应的屏幕坐标
*/
public Point getXY(double radius, double angle)
{
int x = (int)(0.5 + radius * CMath.sin(CMath.toRadians(angle)));
int y = (int)(0.5 + radius * CMath.cos(CMath.toRadians(angle)));
return (new Point(place.X + x, place.Y - y));
}
/**
* 功能：获得屏幕像素点位置的极坐标半径，由于是输入参数是扫描平面上的值，故与雷达仰角无关。
* 参数：
* x - 水平坐标
* y - 垂直坐标
* 返回值：
* 与极坐标中心的距离，即极半径
*/
public double getRadius(int x, int y)
{
return (CMath.sqrt(1.0 * (x - place.X) * (x - place.X) + 1.0 * (y - place.Y) * (y - place.Y)));
}
/**
* 功能：获得经纬度位置的极坐标半径，与雷达仰角有关。
* 参数：
* lon - 经度坐标
* lat - 纬度坐标
* 返回值：
* 与极坐标中心的距离(象素点)，即极半径
*/
public double getRadius(double lon, double lat)
{
Point pos = getPosition(lon, lat);//此函数已经考虑了仰角的影响
return (getRadius(pos.X, pos.Y));
}
/**
* 功能：获得屏幕像素点位置的极坐标角度(扫描平面与0度平面均相同)，与雷达仰角无关。
* 参数：
* x - 水平坐标
* y - 垂直坐标
* 返回值：
* 角度值，自正北方向顺时针
*/
public double getAngle(int x, int y)
{
double agl = 0.0;
if (x == place.X && y == place.Y)
{
//重合
agl = 0.0;
}
else if (x == place.X)
{
agl = y > place.Y ? 180.0 : 360.0;
}
else if (y == place.Y)
{
agl = x > place.X ? 90.0 : 270.0;
}
else
{
agl = CMath.toDegrees(CMath.atan(1.0 * CMath.abs(x - place.X) / CMath.abs(y - place.Y)));
agl =
x > place.X && y < place.Y ? agl : //直角坐标的第一象限
x < place.X && y < place.Y ? 180.0 - agl : //直角坐标的第二象限
x < place.X && y > place.Y ? 180.0 + agl : //直角坐标的第三象限
x > place.X && y > place.Y ? 360.0 - agl : //直角坐标的第四象限
agl;
}
//System.out.println(agl);
return (agl);
}
/**
* 功能：获得经纬度位置的极坐标角度(扫描平面与0度平面均相同)，与雷达仰角无关。
* 参数：
* lon - 水平坐标
* lat - 垂直坐标
* 返回值：
* 角度值，自正北方向顺时针
*/
public double getAngle(double lon, double lat)
{
/*
//若通过获得屏幕坐标来计算角度，精度比较差，特别是在极坐标中心附近
Point p = getPosition(lon, lat);
return(getAngle(p.x, p.y);
*/
double agl = 0.0;
if( lon == center.X && lat == center.Y) //重合
{
agl = 0.0;
}
else if( lon == center.X )
{
agl = lat > center.Y ? 360.0 : 180.0;
}
else if( lat == center.Y )
{
agl = lon > center.X ? 90.0 : 270.0;
}
else
{
//注：由于经向和纬向的球面距离不等(华南，经向>纬向)，故点(1,1)与中心点(0,0)的极角不等45度，而应是略大于45度
agl = CMath.toDegrees(CMath.atan((CMath.abs(lon-center.X)*kmPerDegreeX)/(CMath.abs(lat-center.Y)*kmPerDegreeY)));
agl =
lon > center.X && lat > center.Y ? agl : //第一象限
lon < center.X && lat > center.Y ? 180.0 - agl : //第二象限
lon < center.X && lat < center.Y ? 180.0 + agl : //第三象限
lon > center.X && lat < center.Y ? 360.0 - agl : //第四象限
agl;
}
return(agl);
}
/**
* 功能:
* 获得屏幕坐标对应的经纬度
* 参数:
* x - 屏幕水平坐标
* y - 屏幕垂直坐标
* 返回值:
* 对应的经纬度
*/
public PointF getCoordinate(int x, int y)
{
/*
目标点 A(X,Y) 弧度
中心点 B(A,B) 弧度
AB球面距离=R*{arccos[cos(B)*cos(Y)*cos(A- X)+sin(B)*sin(Y)]}, by Google
经度相同 =& gt; AB = R*{arccos[cos(B)*cos(Y)+sin(B)*sin(Y)]}
=> AB = R*{arccos[cos(B-Y)]}
=> AB = R * (B-Y)
=> AB / R = B - Y
=> Y = B - AB /R
=> Y = B - (y-centerPosition.y)*cosineElevation/perKilometer/scale/R
*/
double lat = CMath.toDegrees(CMath.toRadians(center.Y) + (place.Y-y)*cosineElevation/perKilometer/scale/Polar.RADIUS);
double disX0 = distanceOfSphere(center.X, lat, center.X+1.0, lat);//0度平面上1经度的球面距离
double disX = disX0 / cosineElevation; //扫描平面上1经度的距离
double perDegreeX = disX * perKilometer * scale; //扫描平面上1经度的对应的像素点数
double lon = center.X + (x - place.X) / perDegreeX;
return (new PointF((float)lon, (float)lat));
}
/**
* 功能:
* 获得四角坐标对应的经纬度
* 参数:
* W - 图像高度
* H - 图像宽度
* 返回值:
* 对应的经纬度
* add by sungaoyong 2011-02-10
*/
public PointF[] getRecF(Double W, double H)
{
PointF PointLeftTop = getCoordinate((int)(0.5 + place.X - scale * W / 2), (int)(0.5 + place.Y - scale * H / 2));
PointF PointRightTop = getCoordinate((int)(0.5 + place.X + scale * W /2), (int)(0.5 + place.Y - scale * H / 2));
PointF PointLeftBottom = getCoordinate((int)(0.5 + place.X - scale * W / 2), (int)(0.5 + place.Y + scale * H / 2));
PointF PointRightBottom = getCoordinate((int)(0.5 + place.X + scale * W / 2), (int)(0.5 + place.Y + scale * H / 2));
return (new PointF[] { PointLeftTop, PointRightTop, PointLeftBottom, PointRightBottom });
}
/**
* 功能：
* 画经线、纬线
* 参数：
* g - 图形设备
* f - 字体
* c - 画线颜色
* inc_lon - 经线间隔//未使用
* inc_lat - 纬线间隔//未使用
* 返回值：
* 无
*/
public void drawGridLine(System.Drawing.Graphics g, Font f, Color c, int inc_lon, int inc_lat)
{
//Color saveColor = g.getColor();
//Font saveFont = g.getFont();
Pen blue = new Pen(Color.Blue);
//以下两行改进线条的锯齿
//RenderingHints renderHints = new RenderingHints(RenderingHints.KEY_ANTIALIASING,RenderingHints.VALUE_ANTIALIAS_ON);
//g.setRenderingHints(renderHints);
// g.setColor(Color.black);//背景色
// g.fillRect(c.x-(int)(z*240), c.y-(int)(z*240), (int)(z*240*2), (int)(z*240*2));
//g.setColor(c);//雷达图形区域的边框颜色
g.DrawArc(blue,new RectangleF((int)(0.5+place.X-scale*240), (int)(0.5+place.Y-scale*240), (int)(0.5+scale*240*2), (int)(0.5+scale*240*2)),0,0);
//画极径
Point pos1, pos2;
for(double i=0.0; i<180.0; i=i+30.0)
{
pos1 = getXY(scale*240.0, 0.0+i);
pos2 = getXY(scale*240.0, 180.0+i);
g.DrawLine(blue,pos1.X, pos1.Y, pos2.X, pos2.Y);
}
//画极圈
for(int i=50; i<=200; i=i+50) //每50公里画一个圈
{
g.DrawArc(blue,new RectangleF((int)(0.5+place.X-scale*i), (int)(0.5+place.Y-scale*i), (int)(0.5+scale*i*2), (int)(0.5+scale*i*2)), 0, 360);
}
g.DrawArc(blue,new RectangleF((int)(0.5 + place.X - scale * 240), (int)(0.5 + place.Y - scale * 240), (int)(0.5 + scale * 240 * 2), (int)(0.5 + scale * 240 * 2)), 0, 360);//外圈240公里
//g.setFont(saveFont);
//g.setColor(saveColor);
}
}
}