////////////////////////////////////
// C++ による 地図投影法のクラス //
////////////////////////////////////
#include "Projection.h"
#include "stdlib.h"
#include "stdio.h"
#include "string.h"
#include "math.h"
/////////////////////////////////////////////////////////////////////////////
// Projectionクラスの構築子
Projection::Projection()
{
}
/////////////////////////////////////////////////////////////////////////////
// Projectionクラスの消滅子
Projection::~Projection()
{
}
/////////////////////////////////////////////////////////////////////////////
// 文字表現の緯度を実数値に変換して返す
// 89 --> 89.000000
// 89.999999 --> 89.999999
// N89 --> 89.000000
// N89.999999 --> 89.999999
// N89.59' --> 89.983333
// N89.59'99 --> 89.999833
// N89.59'59" --> 89.999722
// N89.59'59"99 --> 89.99999722
// -89 --> -89.000000
// -89.999999 --> -89.999999
// S89 --> -89.000000
// S89.999999 --> -89.999999
// S89.59' --> -89.983333
// S89.59'99 --> -89.999833
// S89.59'59" --> -89.999722
// S89.59'59"99 --> -89.99999722
double Projection::ConvLatitudeToDouble( char* Str )
{
char *p = Str;
int Sign = 1; // Northern=1, Southern=-1
char digits[64];
int cnt = 0;
char degrees[64]; // degree
char minutes[64]; // minute
char seconds[64]; // second
BOOL minutesFlg = FALSE;
BOOL secondsFlg = FALSE;
BOOL fructionMode = FALSE;
double value; // latitude
if( Str[0] == '\0' ){
return INT_MAX;
}
memset( digits, 0, sizeof(digits) );
memset( degrees, 0, sizeof(degrees) );
memset( minutes, 0, sizeof(minutes) );
memset( seconds, 0, sizeof(seconds) );
while( *p != '\0' ){
if( ( *p >= '0' )&&( *p <= '9' ) ){
digits[cnt++] = *p;
}else if( *p == CHAR_PUNCT ){
if( ( minutesFlg == TRUE )||( secondsFlg == TRUE ) ){
return INT_MAX; // illegal fruction point!
}
strcpy( degrees, digits );
memset( digits, 0, sizeof(digits) );
cnt = 0;
fructionMode = TRUE;
}else if( *p == '\'' ){
if( ( minutesFlg == TRUE )||( secondsFlg == TRUE ) ){
return INT_MAX; // illegal order!
}
minutesFlg = TRUE;
strcpy( minutes, digits );
memset( digits, 0, sizeof(digits) );
cnt = 0;
fructionMode = TRUE;
}else if( *p == '\"' ){
if( secondsFlg == TRUE ){
return INT_MAX; // illegal order!
}
secondsFlg = TRUE;
strcpy( seconds, digits );
memset( digits, 0, sizeof(digits) );
cnt = 0;
fructionMode = TRUE;
}else if( *p == 'N' ){ // Northern
if( p != Str ){
return INT_MAX; // expect as top char!
}
Sign = 1;
}else if( *p == 'S' ){ // Southern
if( p != Str ){
return INT_MAX; // expect as top char!
}
Sign = -1;
}else if( *p == '+' ){ // Northern
if( p != Str ){
return INT_MAX; // expect as top char!
}
Sign = 1;
}else if( *p == '-' ){ // Southern
if( p != Str ){
return INT_MAX; // expect as top char!
}
Sign = -1;
}else if( ( *p == ' ' )||( *p == '\t' ) ){ // 空白
if( cnt > 0 ){
return INT_MAX; // illegal space
}
}else{
return INT_MAX;
}
p++;
}
if( cnt > 0 ){
if( fructionMode == TRUE ){
if( secondsFlg == TRUE ){
seconds[strlen(seconds)] = CHAR_PUNCT;
strcat(seconds, digits );
}else if( minutesFlg == TRUE ){
minutes[strlen(minutes)] = CHAR_PUNCT;
strcat(minutes, digits );
}else if( ( minutesFlg == FALSE )&&( secondsFlg == FALSE ) ){
degrees[strlen(degrees)] = CHAR_PUNCT;
strcat(degrees, digits );
}
}else{
if( ( minutesFlg == FALSE )&&( secondsFlg == FALSE ) ){
strcpy( degrees, digits );
}
}
}
value = atof(degrees) + (atof(minutes) / 60 ) + ( atof(seconds) / 3600 );
value *= Sign;
return value;
}
/////////////////////////////////////////////////////////////////////////////
// 文字表現の経度を実数値に変換して返す
// 179 --> 179.000000
// 179.999999 --> 179.999999
// E179 --> 179.000000
// E179.999999 --> 179.999999
// E179.59' --> 179.983333
// E179.59'99 --> 179.999833
// E179.59'59" --> 179.999722
// E179.59'59"99 --> 179.99999722
// -179 --> -179.000000
// -179.999999 --> -179.999999
// W179 --> -179.000000
// W179.999999 --> -179.999999
// W179.59' --> -179.983333
// W179.59'99 --> -179.999833
// W179.59'59" --> -179.999722
// W179.59'59"99 --> -179.99999722
double Projection::ConvLongitudeToDouble( char* Str )
{
char *p = Str;
int Sign = 1; // Eastern=1, Western=-1
char digits[64];
int cnt = 0;
char degrees[64]; // degree
char minutes[64]; // minute
char seconds[64]; // second
BOOL minutesFlg = FALSE;
BOOL secondsFlg = FALSE;
BOOL fructionMode = FALSE;
double value; // latitude
if( Str[0] == '\0' ){
return INT_MAX;
}
memset( digits, 0, sizeof(digits) );
memset( degrees, 0, sizeof(degrees) );
memset( minutes, 0, sizeof(minutes) );
memset( seconds, 0, sizeof(seconds) );
while( *p != '\0' ){
if( ( *p >= '0' )&&( *p <= '9' ) ){
digits[cnt++] = *p;
}else if( *p == CHAR_PUNCT ){
if( ( minutesFlg == TRUE )||( secondsFlg == TRUE ) ){
return INT_MAX; // illegal fruction point!
}
strcpy( degrees, digits );
memset( digits, 0, sizeof(digits) );
cnt = 0;
fructionMode = TRUE;
}else if( *p == '\'' ){
if( ( minutesFlg == TRUE )||( secondsFlg == TRUE ) ){
return INT_MAX; // illegal order!
}
minutesFlg = TRUE;
strcpy( minutes, digits );
memset( digits, 0, sizeof(digits) );
cnt = 0;
fructionMode = TRUE;
}else if( *p == '\"' ){
if( secondsFlg == TRUE ){
return INT_MAX; // illegal order!
}
secondsFlg = TRUE;
strcpy( seconds, digits );
memset( digits, 0, sizeof(digits) );
cnt = 0;
fructionMode = TRUE;
}else if( *p == 'E' ){ // Northern
if( p != Str ){
return INT_MAX; // expect as top char!
}
Sign = 1;
}else if( *p == 'W' ){ // Southern
if( p != Str ){
return INT_MAX; // expect as top char!
}
Sign = -1;
}else if( *p == '+' ){ // Eastern
if( p != Str ){
return INT_MAX; // expect as top char!
}
Sign = 1;
}else if( *p == '-' ){ // Western
if( p != Str ){
return INT_MAX; // expect as top char!
}
Sign = -1;
}else if( ( *p == ' ' )||( *p == '\t' ) ){ // 空白
if( cnt > 0 ){
return INT_MAX; // illegal space
}
}else{
return INT_MAX;
}
p++;
}
if( cnt > 0 ){
if( fructionMode == TRUE ){
if( secondsFlg == TRUE ){
seconds[strlen(seconds)] = CHAR_PUNCT;
strcat(seconds, digits );
}else if( minutesFlg == TRUE ){
minutes[strlen(minutes)] = CHAR_PUNCT;
strcat(minutes, digits );
}else if( ( minutesFlg == FALSE )&&( secondsFlg == FALSE ) ){
degrees[strlen(degrees)] = CHAR_PUNCT;
strcat(degrees, digits );
}
}else{
if( ( minutesFlg == FALSE )&&( secondsFlg == FALSE ) ){
strcpy( degrees, digits );
}
}
}
value = atof(degrees) + (atof(minutes) / 60 ) + ( atof(seconds) / 3600 );
value *= Sign;
return value;
}
/////////////////////////////////////////////////////////////////////////////
// 経緯度をランベルト正積図法(Lambert Azimuthal equal-area projection)の座標系に変換
// IN/ lon 経度の実数値 (度単位)
// IN/ lat 緯度の実数値 (度単位)
// IN/ lon0 ランベルト正積図法の中心経度 (度単位)
// IN/ lat0 ランベルト正積図法の中心緯度 (度単位)
// /OUT *x 求められたＸ座標 -1.0 〜 +1.0 地球の半径を掛けると実距離になる )
// /OUT *y 求められたＹ座標 -1.0 〜 +1.0 地球の半径を掛けると実距離になる )
void Projection::GeoToLea( double lon, double lat, double lon0, double lat0, double *x, double *y )
{
double phi; // latitude radian
double phi1;
double lamda; // longitude radian
double lamda0;
double pi = 3.14159265358979;
double k;
phi = lat * pi / 180.0;
phi1 = lat0 * pi / 180.0;
lamda = lon * pi / 180.0;
lamda0 = lon0 * pi / 180.0;
k = 2 / ( 1 + sin(phi1) * sin(phi) + cos(phi1) * cos(phi) * cos(lamda-lamda0) );
if( k > 0 ){
k = sqrt(k);
}else{
k = 0;
}
*x = k * cos(phi) * sin(lamda - lamda0);
*y = k * ( cos(phi1) * sin(phi) - sin(phi1) * cos(phi) * cos(lamda-lamda0) );
}
/////////////////////////////////////////////////////////////////////////////
// LEA座標系をマップ画像の位置に変換
POINT Projection::ConvLeaToMapXy( double x, double y, double shift_x, double shift_y, double magnify )
{
POINT pnt;
pnt.x = (long)(( x + shift_x ) * magnify);
pnt.y = (long)(( y * (-1) + shift_y) * magnify);
return pnt;
}
/////////////////////////////////////////////////////////////////////////////
// メートル単位の距離をLEA座標系のピクセル数に換算
double Projection::lengthLEA( int length, double magnify )
{
return ( (double)length * magnify / 6365000.0 );
}
/////////////////////////////////////////////////////////////////////////////
// 経緯度を横メルカトル図法(Transverse Mercator's projection)の座標系に変換
// IN/ lon 経度の実数値 (度単位)
// IN/ lat 緯度の実数値 (度単位)
// IN/ lon0 横メルカトル図法の中心経度 (度単位)
// /OUT *x 求められたＸ座標 中心経線からの距離 メートル単位)
// /OUT *y 求められたＹ座標 赤道からの距離 メートル単位)
void Projection::GeoToTM( double lon, double lat, double lon0, double *x, double *y )
{
double phi; // latitude radian
double lamda; // longitude radian
double pi = 3.14159265358979;
double a = 6378137;
double f = 1 / 298.2572336;
double k0 = 0.9996;
double lamda0;
double e2;
double E2;
double n;
double M;
double T;
double C;
double A;
phi = lat * pi / 180.0;
lamda = lon * pi / 180.0;
//// ユニバーサル横メルカトル図法ではないので、ゾーンで分けない
//// int zone;
//// zone = ( (int)( lon + 180 ) ) / 6 + 1;
//// lamda0 = ( ( zone - 1 ) * 6 - 180 + 3 ) * pi / 180.0;
lamda0 = lon0 * pi / 180.0; // center longitude
e2 = 2 * ( f - f * f );
E2 = e2 / ( 1 - e2 );
//// b = a * sqrt( 1 - e2 );
n = a / sqrt( 1 - e2 * sin(phi) * sin(phi) );
M = a * (
( 1 - e2/4 - 3*e2*e2/64 - 5*e2*e2*e2/256 ) * phi -
( 3*e2/8 + 3*e2*e2/32 + 45*e2*e2*e2/1024 ) * sin(2*phi) +
( 15*e2*e2/256 + 45*e2*e2*e2/1024 ) * sin(4*phi) -
( 35*e2*e2*e2/3072 ) * sin(6*phi)
);
T = tan(phi) * tan(phi);
C = E2 * cos(phi) * cos(phi);
A = ( lamda - lamda0 ) * cos(phi);
*x = k0 * n * ( A +
(1-T+C)*A*A*A/6 +
(5 - 18*T + T*T + 72*C - 58*E2)*A*A*A*A*A*120
)
; ////+ 500000; // ユニバーサル横メルカトル図法ではないので、経線方向にシフトさせない。Ｘ座標は負値あり
*y = k0 * ( M +
n * tan(phi) * ( A*A/2 +
( 5 - T + 9*C + 4*C*C) * A*A*A*A*24 +
( 61 - 58 * T + T*T + 600*C - 330*E2) * A*A*A*A*A*A*720
)
);
}
/////////////////////////////////////////////////////////////////////////////
// TM座標系をマップ画像の位置に変換
POINT Projection::ConvTMToMapXy( double x, double y, double shift_x, double shift_y, double magnify )
{
POINT pnt;
pnt.x = (long)(( x + shift_x ) * magnify);
pnt.y = (long)(( y * (-1) + shift_y) * magnify);
return pnt;
}
/////////////////////////////////////////////////////////////////////////////
// メートル単位の距離をTM座標系のピクセル数に換算
double Projection::lengthTM( int length, double magnify )
{
return ( (double)length * magnify );
}