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);
}
}
}