Polar坐标投影(C++)

//Polar.cpp

/**
*
*  Polar 投影(扫描方式，自正北方向顺时针)
*
* How to use Polar class:
*
* Polar polar = new Polar(Point(240, 240), 109.24, 24.35, 1.5);//构造函数
* polar->setScale(1.0);//设置比例尺，1公里对应1个像素点
* ...
*
**/

#include "Polar.h"

/**
 *
 *         扫描平面
 *            /
 *           /
 *          /
 *         /
 *        /
 *       /  仰角
 *      -------------------- 0度平面
 *
 * 如图所示：
 *          扫描平面=>0度平面需要乘以cos(仰角)
 *          0度平面=>扫描平面需要除以cos(仰角)
 *
 * 注意，日常显示的雷达图是扫描平面上的图。本类所说的屏幕指扫描平面。
 *
 */

/**
 * 雷达扫描示意图
 *
 *                    359 0
 *                        |     radius
 *                        |       /
 *                        |      /
 *                        |angle/
 *                        |    /
 *                        | ^ /
 *                        |  /
 *                        | /
 *                        |/
 * 270 -----------------中心----------------- 90
 *                        |
 *                        |
 *                        |
 *                        |
 *                        |
 *                        |
 *                        |
 *                        |
 *                        |
 *                       180
 */

/**
 * 功能：角度转换继弧度
 * 参数：
 *      degrees - 角度值
 * 返回值：
 *      弧度值
 */
double Polar::toRadians(double degrees) {
    return(PI/180.0*degrees);
}

/**
 * 功能：角度转换继弧度
 * 参数：
 *      degrees - 角度值
 * 返回值：
 *      弧度值
 */
double Polar::toDegrees(double radians) {
    return(radians*180.0/PI);
}

/**
 * 功能：计算球面上两点间的距离(单位：公在edu.gimdr.Atmos.Meteorology类中写有，为避免import过多的类，故重写一份
 * 参数：
 *  lon1,lat1   - 第1点的位置(经纬度)
 *  lon2,lat2   - 第2点的位置(经纬度)
 * 返回值：
 *      球面距离
 */
double Polar::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   = toRadians(lon1);
        double  rlat1   = toRadians(lat1);
        double  rlon2   = toRadians(lon2);
        double  rlat2   = toRadians(lat2);
        double  val     = acos(cos(rlat2)*cos(rlat1)*cos(rlon2-rlon1)+sin(rlat2)*sin(rlat1));
        val *= RADIUS;
        return(val);
    }

/**
 * 功能：构造函数
 * 参数：
 *      pos     - 中心对应的屏幕位置
 *      lon     - 中心对应的经度坐标
 *      lat     - 中心对应的纬度坐标
 * 返回值：
 *      无
 */
Polar::Polar(Point pos, double lon, double lat) {
        elevation       = 0.0;//无仰角
        cosineElevation = 1.0;
        setCenterPosition(pos);
        setCenterCoordinate(lon, lat);
    }

/**
 * 功能：构造函数
 * 参数：
 *      pos     - 中心对应的屏幕位置
 *      lon     - 中心对应的经度坐标
 *      lat     - 中心对应的纬度坐标
 *      elv     = 仰角
 * 返回值：
 *      无
 */
 Polar::Polar(Point pos, double lon, double lat, double elv) {
        elevation       = fabs(elv);
        while(elevation >= 90.0) {//在0-90度之间，但不能为90度
            elevation   -= 90.0;
        }
        cosineElevation = cos(elevation);//仰角的余弦值
        setCenterPosition(pos);
        setCenterCoordinate(lon, lat);
    }

/**
 * 功能：获得极坐标中心位置
 * 参数：
 *      无
 * 返回值：
 *      极坐标中心对应的屏幕位置
 */
Point Polar::getCenterPosition() {
        return(centerPosition);
    }

/**
 * 功能：设置极坐标的中心位置(屏幕坐标)
 * 参数：
 *      pos     - 新的中心位置(屏幕坐标)
 * 返回值：
 *      无
 */
void Polar::setCenterPosition(Point pos) {
        centerPosition  = pos;
    }

/**
 * 功能：获得极坐标中心对应的经度坐标
 * 参数：
 *      无
 * 返回值：
 *      极坐标中心对应的经度坐标
 */
double Polar::getCenterLongitude() {
        return(centerLongitude);
    }

/**
 * 功能：获得极坐标中心对应的纬度坐标
 * 参数：
 *      无
 * 返回值：
 *      极坐标中心对应的纬度坐标
 */
double Polar::getCenterLatitude() {
        return(centerLatitude);
    }

/**
 * 功能：设置极坐标的中心位置(经纬度坐标)
 * 参数：
 *      lon     - 新的中心位置(经度值)
 *      lat     - 新的中心位置(纬度值)
 * 返回值：
 *      无
 */
void Polar::setCenterCoordinate(double lon, double lat) {
        centerLongitude = lon;
        centerLatitude  = lat;
        //中心经纬度或仰角发生改变，必须重新计算经向和纬向的1度对应的球面距离
        kmPerDegreeX    = distanceOfSphere(erLongitude, centerLatitude, centerLongitude+1.0, centerLatitude) / cosineElevation;
        kmPerDegreeY    = distanceOfSphere(erLongitude, centerLatitude, centerLongitude, centerLatitude+1.0) / cosineElevation;
    }

/**
 * 功能：获得仰角
 * 参数：
 *      无
 * 返回值：
 *      仰角的度数
 */
double Polar::getElevation() {
        return(toDegrees(elevation));
    }

/**
 * 功能：设置仰角
 * 参数：
 *      elv     - 新的仰角
 * 返回值：
 *      无
 */
void Polar::setElevation(double elv) {
        elevation       = fabs(elv);
        while(elevation >= 90.0) {//在0-90度之间，但不能为90度
            elevation   -= 90.0;
        }
        cosineElevation = cos(elevation);//仰角的余弦值
        //中心经纬度或仰角发生改变，必须重新计算经向和纬向的1度对应的球面距离
        kmPerDegreeX    = distanceOfSphere(erLongitude, centerLatitude, centerLongitude+1.0, centerLatitude) / cosineElevation;
        kmPerDegreeY    = distanceOfSphere(erLongitude, centerLatitude, centerLongitude, centerLatitude+1.0) / cosineElevation;
    }

/**
 * 功能：获得比例尺，即1公里对应的像素点数
 * 参数：
 *      无
 * 返回值：
 *      比例尺
 */
double Polar::getScale() {
        return(perKilometer);
    }

/**
 * 功能：设置比例尺，即1公里对应的像素点数
 * 参数：
 *      value   - 比例尺数值
 * 返回值：
 *      无
 */
    void Polar::setScale(double value) {
        perKilometer    = value;
    }

/**
 * 功能：获得经纬度对应的屏幕像素坐标，与雷达仰角有关，主要用于体扫数据显示、底图叠加等。
 * 参数：
 *      lon - 经度
 *      lat - 纬度
 * 返回值：
 *      对应的屏幕坐标
 */
    Point Polar::getPosition(double lon, double lat) {
        double  disX    = distanceOfSphere( centerLatitude, centerLongitude, centerLatitude)/cosineElevation;
        double  disY    = distanceOfSphere(erLongitude, lat, centerLongitude, centerLatitude)/cosineElevation;
        return(Point(centerPosition.x+(lon>erLongitude?1:-1)*(int)(disX*perKilometer),centerPosition.y-(lat>erLatitude?1:-1)*(int)(disY*perKilometer)));
    }

/**
 * 功能：获得极坐标对应的屏幕像素坐标，与雷达仰角无关，主要用于体扫数据显示、底图叠加等。
 * 参数：
 *      radius      - 极半径
 *      angle       - 角度(以正北方向顺时针)
 * 返回值：
 *      对应的屏幕坐标
 */

    Point Polar::getXY(double radius, double angle) {
        int x   = (int)(radius * sin(toRadians(angle)));
        int y   = (int)(radius * cos(toRadians(angle)));
        return(Point(centerPosition.x+x, centerPosition.y-y));
    }

/**
 * 功能：获得屏幕像素点位置的极坐标半径，由于是输入参数是扫描平面上的值，故与雷达仰角无关。
 * 参数：
 *      x   - 水平坐标
 *      y   - 垂直坐标
 * 返回值：
 *      与极坐标中心的距离，即极半径
 */
    double Polar::getRadius(int x, int y) {
        return((1.0*(x-centerPosition.x)*(x-centerPosition.x)+1.0*(y-centerPosition.y)*(y-centerPosition.y)));
    }

/**
 * 功能：获得经纬度位置的极坐标半径，与雷达仰角有关。
 * 参数：
 *      lon - 经度坐标
 *      lat - 纬度坐标
 * 返回值：
 *      与极坐标中心的距离(象素点)，即极半径
 */
    double Polar::getRadius(double lon, double lat) {
        Point   pos = getPosition(lon, lat);//此函数已经考虑了仰角的影响
        return(getRadius(pos.x, pos.y));
    }

/**
 * 功能：获得屏幕像素点位置的极坐标角度(扫描平面与0度平面均相同)，与雷达仰角无关。
 * 参数：
 *      x   - 水平坐标
 *      y   - 垂直坐标
 * 返回值：
 *      角度值，自正北方向顺时针
 */
    double Polar::getAngle(int x, int y) {
        double  agl = 0.0;
        if( x == centerPosition.x && y == centerPosition.y ) {//重合
            agl = 0.0;
        }
        else if( x == centerPosition.x ) {
            agl = y > centerPosition.y ? 180.0 : 360.0;
        }
        else if( y == centerPosition.y ) {
            agl = x > centerPosition.x ? 90.0 : 270.0;
        }
        else {
            agl = toDegrees(atan(1.0*fabs(x-centerPosition.x)/fabs(y-centerPosition.y)));
            agl =
                    x > centerPosition.x && y < centerPosition.y ? agl :            //直角坐标的第一象限
                    x < centerPosition.x && y < centerPosition.y ? 180.0 - agl :    //直角坐标的第二象限
                    x < centerPosition.x && y > centerPosition.y ? 180.0 + agl :    //直角坐标的第三象限
                    x > centerPosition.x && y > centerPosition.y ? 360.0 - agl :    //直角坐标的第四象限
                    agl;
        }
        return(agl);
    }

/**
 * 功能：获得经纬度位置的极坐标角度(扫描平面与0度平面均相同)，与雷达仰角无关。
 * 参数：
 *      lon - 水平坐标
 *      lat - 垂直坐标
 * 返回值：
 *      角度值，自正北方向顺时针
 */
    double Polar::getAngle(double lon, double lat) {
/*
//若通过获得屏幕坐标来计算角度，精度比较差，特别是在极坐标中心附近
        Point   p   = getPosition(lon, lat);
        return(getAngle(p.x, p.y);
*/
        double  agl = 0.0;
        if( lon == centerLongitude && lat == centerLatitude ) {//重合
            agl = 0.0;
        }
        else if( lon == centerLongitude ) {
            agl = lat > centerLatitude ? 360.0 : 180.0;
        }
        else if( lat == centerLatitude ) {
            agl = lon > centerLongitude ? 90.0 : 270.0;
        }
        else {
            //注：由于经向和纬向的球面距离不等(华南，经向>纬向)，故点(1,1)与中心点(0,角不等45度，而应是略大于45度
            agl = toDegrees(atan((fabs(lon-centerLongitude)*kmPerDegreeX)/(fabs(centerLatitude)*kmPerDegreeY)));
            agl =
                    lon > centerLongitude && lat > centerLatitude ? agl :           //第一象限
                    lon < centerLongitude && lat > centerLatitude ? 180.0 - agl :   //第二象限
                    lon < centerLongitude && lat < centerLatitude ? 180.0 + agl :   //第三象限
                    lon > centerLongitude && lat < centerLatitude ? 360.0 - agl :   //第四象限
                    agl;
        }
        //System.out.println(agl);
        return(agl);
    }

/**
 * 功能：获得屏幕坐标对应的经度值(根据目标点的经向球面距离来计算，雷达南面和北面的值略有差别)，与雷达仰角有关。
 *       主要用于雷达产品的定位、底图叠加、转换为经纬度网格产品、拼图等。
 * 参数：
 *      x   - 横坐标，自西向东
 *      y   - 纵坐标，自北向南
 * 返回值：
 *      对应的经度坐标
 */
    double Polar::getLongitude(int x, int y) {
        double  lat         = getLatitude(x, y);
        double  disX0       = distanceOfSphere(erLongitude, lat, centerLongitude+1.0, lat);//0度平面上1经度的球面距离
        double  disX        = disX0 / cosineElevation;      //扫描平面上1经度的距离
        double  perDegreeX  = disX * perKilometer;          //扫描平面上1经度的对应的像素点数
        return(centerLongitude + (x - centerPosition.x) / perDegreeX);
    }

/**
 * 功能：获得屏幕坐标对应的纬度值(根据极坐标中心点的纬向球面距离来计算)，与雷达仰角有关。
 *       主要用于雷达产品的定位、底图叠加、转换为经纬度网格产品、拼图等。
 * 参数：
 *      x   - 横坐标，自西向东(未用到)
 *      y   - 纵坐标，自北向南
 * 返回值：
 *      对应的纬度坐标
 */
    double Polar::getLatitude(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
 * 经度相同 => 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)/perKilometer/R
 */
        double  val = toRadians(centerLatitude) + (centerPosition.y-y)*cosineElevation/perKilometer;
        double  fRadius = RADIUS;
        val = val / fRadius;
        return(toDegrees(val));
    }

//Polar.h

#ifndef PolarH
#define PolarH

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

//静态常量，地球半径，来源：《大气科学常用公式》，P601，附录
#define RADIUS          6371.004;//地球平均半径，单位：公里(Km)。
#define RADIUS_POLAR    6356.755;//地球两极半径，单位：公里(Km)。
#define RADIUS_EQUATOR  6373.140;//地球赤道半径，单位：公里(Km)。

#define PI              3.14159265358979323846264338327950288

typedef struct _tagPOINT {
    int x;
    int y;
    _tagPOINT() {}
    _tagPOINT(int _x, int _y) : x(_x), y(_y) {}
}   Point;

class Polar {

private:    //私有成员
    Point   centerPosition;     //中心对应的屏幕位置
    double  centerLongitude;    //中心经度
    double  centerLatitude;     //中心纬度
    double  perKilometer;       //比例尺：一公里对应的像素点数(扫描平面)
    double  elevation;          //仰角
    double  cosineElevation;    //仰角的余弦值
    double  kmPerDegreeX;       //1经度对应的距离(公里)，不同纬度数值不同
    double  kmPerDegreeY;       //1纬度对应的距离(公里)，不同纬度数值不同

public: //公共成员
    //1、角度=>板度
    double toRadians(double degrees);
    //2、弧度=>角度
    double toDegrees(double radians);
    //3、计算球面距离
    double distanceOfSphere(double lon1, double lat1, double lon2, double lat2);
    //4、构造函数，不考虑仰角
    Polar(Point pos, double lon, double lat);
    //5、构造函数，考虑仰角
    Polar(Point pos, double lon, double lat, double elv);
    //6、获得极坐标中心点的屏幕位置
    Point   getCenterPosition();
    //7、设置极坐标中心点的屏幕位置
    void    setCenterPosition(Point pos);
    //8、获得极坐标中心点的经度
    double  getCenterLongitude();
    //9、获得极坐标中心点的纬度
    double  getCenterLatitude();
    //10、设置极坐标中心点的经纬度
    void    setCenterCoordinate(double lon, double lat);
    //11、获得仰角
    double  getElevation();
    //12、设置仰角
    void    setElevation(double elv);
    //13、获得缩放比例尺
    double  getScale();
    //14、设置缩放比例尺
    void    setScale(double value);
    //15、获得经纬度对应的屏幕坐标(扫描平面)，主要用于体扫数据显示、底图叠加等。
    Point   getPosition(double lon, double lat);
    //16、获得极坐标对应的屏幕坐标(扫描平面)，主要用于体扫数据显示、底图叠加等。
    Point   getXY(double radius, double angle);
    //17、根据屏幕坐标获得极半径
    double  getRadius(int x, int y);
    //18、根据经纬度坐标获得
    double  getRadius(double lon, double lat);
    //19、根据屏幕坐标获得极角
    double  getAngle(int x, int y);
    //20、根据经纬度坐标获得极角
    double  getAngle(double lon, double lat);
    //21、根据屏幕坐标获得对应的经度值，主要用于雷达产品的定位、底图叠加、转换为经纬度网格产品、拼图等。
    double  getLongitude(int x, int y);
    //22、根据屏幕坐标获得对应的纬度值，主要用于雷达产品的定位、底图叠加、转换为经纬度网格产品、拼图等。
    double  getLatitude(int x, int y);

};

#endif