📄 projctns.cpp
字号:
//////////////////////////////////////////////////////
//
// NRDB Pro - Spatial database and mapping application
//
// Copyright (c) 1989-2004 Richard D. Alexander
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 2
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//
// NRDB Pro is part of the Natural Resources Database Project
//
// Homepage: http://www.nrdb.co.uk/
// Users' Forum: http://nrdb.mypalawan.info/
//
#include "stdafx.h"
#include "nrdb.h"
#include "projctns.h"
#include "bdimportexport.h"
#include "dlgprojection.h"
#ifdef _DEBUG
#undef THIS_FILE
static char BASED_CODE THIS_FILE[] = __FILE__;
#endif
////////////////////////////////////////////////////////////////////////////////
//
// Internal parameters
//
#undef PI
#define PI 3.14159265358979323844
#define EARTHRADIUS 6378137 // radius at equator in metres
////////////////////////////////////////////////////////////////////////////////
//
// Additional, non intrinsic/library maths functions
//
#define cotan(x) (1.0/tan(x))
inline double square(double n) {return n*n;};
inline double pow2(double n) {return n*n;};
inline double pow3(double n) {return pow2(n)*n;}
inline double pow4(double n) {return pow2(n)*pow2(n);}
inline double pow5(double n) {return pow2(n)*pow3(n);}
inline double pow6(double n) {return pow3(n)*pow3(n);}
inline double pow7(double n) {return pow3(n)*pow4(n);}
/////////////////////////////////////////////////////////////////////////////
CString g_sRegCode;
///////////////////////////////////////////////////////////////////////////////
ELLIPSOID CEllipsoid::m_aEllipsoids[] =
{
{"UNIVERSAL TRANSVERSE MERCATOR (UTM)",6378206.400, 294.9786982},
{"LATITUDE / LONGITUDE",LATLON_ID,LATLON_ID},
{"USER DEFINED",0,0},
{"Airy (1830)",6377563.396,299.3249647},
{"US/UK - Modified Airy",6377340.189,299.3249647},
{"Australian National (1966)",6378160.000,298.2500000},
{"APL 4.5 (1968)",6378144.000,298.2300000},
{"Average Terrestrial System 1977", 6378135.000, 298.2570000},
{"Airy (War Office)",6377542.178, 299.3250000},
{"Bessel (Modified)",6377492.018, 299.1528000},
{"Bessel 1841 (Namibia)", 6377483.865, 299.1528128},
{"US - Bessel 1841 (Ethiopia, Indonesia, Japan, Korea)",6377397.155, 299.1528128},
{"UK - Bessel (1841) Revised",6377397.155, 299.1528128},
{"Clarke 1858",6378235.600, 294.2606768},
{"Clarke 1858 (Modified)",6378293.645, 294.26000000},
{"Clarke 1866",6378206.400, 294.9786982},
{"US - Clarke 1880",6378249.145, 293.4650000},
{"Clarke 1880 (Cape)",6378249.145, 293.4663077},
{"Clarke 1880 (Palestine)",6378300.782, 293.4663077},
{"Clarke 1880 (IGN)",6378249.200, 293.4660208},
{"Clarke 1880 (Syria)",6378247.842, 293.4663517},
{"Clarke 1880 (Fiji)",6378301.000, 293.4650000},
{"Danish (1876) or Andrae",6377104.430, 300.0000000},
{"Delambre 1810",6376985.228, 308.6400000},
{"Delambre (Carte de France)",6376985.000, 308.6400000},
{"US - Everest (India 1830)",6377276.345, 300.8017000},
{"UK - Everest (1830)",6377276.345, 300.8017000},
{"US - Everest (Brunei and E. Malaysia (Sabah and Sarawak))",6377298.556, 300.8017000},
{"UK - Everest (Borneo)",6377298.556, 300.8017000},
{"US - Everest (India 1956)",6377301.243, 300.8017000},
{"UK - Everest (India)",6377301.243,300.8017255},
{"US - Everest (W. Malaysia 1969)",6377295.664, 300.8017000},
{"UK - Everest (Malaya RSO)",6377295.664, 300.8017000},
{"US - Everest (W. Malaysia and Singapore 1948)",6377304.063, 300.8017000},
{"UK - Everest (Malaya RKT)",6377304.063, 300.8017000},
{"Everest (Pakistan)",6377309.613, 300.8017000},
{"US - Modified Fischer 1960 (South Asia)",6378155.000, 298.3000000},
{"UK - Fischer 1960 (South Asia)",6378155.000, 298.3000000},
{"Fischer 1968",6378150.000, 298.3000000},
{"Fischer 1960 (Mercury)",6378166.000, 298.3000000},
{"Germaine (Djibouti)",6378284.000, 294.0000000},
{"Hayford 1909",6378388.000, 296.9592630},
{"Helmert 1906",6378200.000, 298.3000000},
{"Hough 1960",6378270.000, 297.0000000},
{"IAG Best Estimate 1975",6378140.000, 298.2570000},
{"Indonesian National (1974)",6378160.000, 298.2470000},
{"US - International 1924",6378388.000, 297.0000000},
{"UK - International",6378388.000, 297.0000000},
{"Krassovsky (1940)",6378245.000, 298.3000000},
{"Krayenhoff 1827",6376950.400, 309.6500000},
{"NWL-8E",6378145.000, 298.2500000},
{"Plessis Modified",6376523.000, 308.6400000},
{"Plessis Reconstituted",6376523.994, 308.6248070},
{"Geodetic Reference System 1967",6378160.000, 298.2471674},
{"Geodetic Reference System 1980",6378137.000, 298.2572221},
{"South American",6378160.000, 298.2500000},
{"Soviet Geodetic System 1985",6378136.000, 298.2570000},
{"Soviet Geodetic System 1990",6378136.000, 298.2578393},
{"Struve 1860",6378298.300, 294.7300000},
{"Svanberg",6376797.000, 304.2506000},
{"Walbeck 1819 (Planheft 1942)",6376895.000, 302.7821565},
{"Walbeck 1819 (AMS 1963)",6376896.000, 302.7800000},
{"World Geodetic System 1966",6378145.000, 298.2500000},
{"World Geodetic System 1972",6378135.000, 298.2600000},
{"World Geodetic System 1984",6378137.000, 298.2572236},
{"US - War Office 1924 (McCaw)",6378300.000, 296.0000000},
{"UK - War Office 1924",6378300.000, 296.0000000},
{"World Geodetic System 1960",6378165.000, 298.3000000}
};
////////////////////////////////////////////////////////////////////////////////
double CEllipsoid::GetSemiMinorAxis()
{
double F, ESQ, EPS;
double ER = m_dSemiMajorAxis;
double RF = m_dRF;
F=1.0/RF;
ESQ=(F+F-pow(F,2));
EPS=ESQ/(1.0-ESQ);
return (1.0-F)*ER;
}
///////////////////////////////////////////////////////////////////////////////
BOOL CProjection::IsUTMOld()
{
return fabs(m_dSemiMajorAxis - UTM_ID_OLD) < 0.001;
}
///////////////////////////////////////////////////////////////////////////////
BOOL CProjection::IsLatLon()
{
return fabs(m_dSemiMajorAxis - LATLON_ID) < 0.001;
}
///////////////////////////////////////////////////////////////////////////////
CProjection::CProjection()
{
m_lScaleCoords = 1;
}
///////////////////////////////////////////////////////////////////////////////
//
// Reloads the projection from the database
//
BOOL CProjection::InitialiseProjection(BDHANDLE hConnect)
{
BOOL bDefault = FALSE;
if (hConnect == NULL)
{
bDefault = TRUE;
hConnect = BDHandle();
};
do
{
// Determine the default projection
CBDProjection projection;
BOOL bFound = BDProjection(hConnect, &projection, BDGETINIT);
while (bFound)
{
if (projection.m_bDefault)
{
(CBDProjection&)*this = projection;
// If projection is UTM then update parameters (backwards compatibility)
if (((CProjection&)projection).IsUTMOld())
{
projection.m_dSemiMajorAxis = CEllipsoid::m_aEllipsoids[UTM_INDEX].m_dSemiMajorAxis;
projection.m_dSemiMinorAxis = ((CEllipsoid&)(CEllipsoid::m_aEllipsoids[UTM_INDEX])).GetSemiMinorAxis();
projection.m_dFalseEasting = 0;
projection.m_dFalseEasting = 500000;
projection.m_dScaleFactorAtOrigin = UTM_SCALEFACTOR;
// Update
BDEnd(hConnect);
BDProjection(hConnect, &projection, BDUPDATE);
}
break;
};
bFound = BDGetNext(hConnect);
}
BDEnd(hConnect);
// If no default projection is defined then display dialog
if (m_nID == 0)
{
if (bDefault)
{
if (AfxMessageBox(IDS_NOPROJECTION, MB_YESNO|MB_DEFBUTTON1) == IDYES)
{
CDlgProjection dlg;
dlg.DoModal();
// If the projection is latitude/longitude then set the default display
if (IsLatLon())
{
BDGetSettings().m_bCoordAsLatLon = TRUE;
}
// Prevent message re-appearing
bDefault = FALSE;
} else
{
return FALSE;
}
} else
{
return FALSE;
}
}
} while (m_nID == 0);
return TRUE;
};
///////////////////////////////////////////////////////////////////////////////
//
// Dummy header for supporting easting and northing as long integers
//
void CProjection::LatLonToTransMercator(double latitude, double longitude,
long *easting, long *northing)
{
double dEasting, dNorthing;
LatLonToTransMercator(latitude, longitude, &dEasting, &dNorthing);
*easting = (long)(dEasting+0.5);
*northing = (long)(dNorthing+0.5);
}
///////////////////////////////////////////////////////////////////////////////
//
// Converts from latitude longitude in decimal degrees to easting and northing
// using a tranverse mercator projection with the parameters defined in this
// class
void CProjection::LatLonToTransMercator(double latitude, double longitude,
double *easting, double *northing)
{
double a,b,e2,N,phi,v,rho,n2,phi1,M,P,IV,VI,I,II,III,IIIA;
// Initialise
a = m_dSemiMajorAxis * m_dScaleFactorAtOrigin;
b = m_dSemiMinorAxis * m_dScaleFactorAtOrigin;
e2 = ((a * a) - (b * b)) / (a * a);
N = (a - b) / (a + b);
phi = latitude * PI / 180.0;
v = a / sqrt(1.0 - e2 * (pow2(sin(phi))));
rho = (v * (1.0 - e2)) / (1.0 - e2 * pow2(sin(phi)));
n2 = (v / rho) - 1.0;
phi1 = m_dLatitudeOrigin * PI / 180.0;
M = b * (((1.0 + N + ((pow2(N)) * 5.0 / 4.0) + ((pow3(N)) * 5.0 / 4.0)) * (phi - phi1))
- ((3.0 * N + 3.0 * (pow2(N)) + (pow3(N)) * 21.0 / 8.0) * sin(phi - phi1) * cos(phi + phi1))
+ ((((pow2(N)) * 15.0 / 8.0) + ((pow3(N)) * 15.0 / 8.0)) * (sin(2.0 * (phi - phi1))) * (cos(2.0 * (phi + phi1))))
- (((pow3(N)) * 35.0 / 24.0) * (sin(3.0 * (phi - phi1))) * (cos(3.0 * (phi + phi1)))));
P = (longitude - m_dLongitudeOrigin) * PI / 180;
// Calculate northing
I = M + m_dFalseNorthing;
II = (v / 2.0) * sin(phi) * cos(phi);
III = (v / 24.0) * sin(phi) * (pow3(cos(phi))) * (5.0 - (pow2(tan(phi))) + 9.0 * n2);
IIIA = (v / 720.0) * sin(phi) * (pow5(cos(phi))) * (61.0 - 58.0 * (pow2(tan(phi))) + (pow4(tan(phi))));
*northing = I + ((pow2(P)) * II) + ((pow4(P)) * III) + ((pow6(P)) * IIIA);
// Calculating easting
IV = v * cos(phi);
v = (v / 6.0) * (pow3(cos(phi))) * ((v / rho) - (pow2(tan(phi))));
VI = (v / 120.0) * (pow5(cos(phi))) * (5.0 - 18.0 * (pow2(tan(phi))) + (pow4(tan(phi))) + 14.0 * n2 - 58.0 * (pow2(tan(phi))) * n2);
VI = (v / 120.0) * (pow5(cos(phi))) * (5.0 - 18.0 * (pow2(tan(phi))) + (pow4(tan(phi))) + 14 * n2 - 58.0 * (pow2(tan(phi))) * n2);
*easting = m_dFalseEasting + (P * IV) + ((pow3(P)) * v) + ((pow5(P)) * VI);
}
///////////////////////////////////////////////////////////////////////////////
//
// Converts from easting and northings to latitude longitude in decimal degrees
// using a tranverse mercator projection with the parameters defined in this
// class
void CProjection::TransMercatorToLatLon(double easting, double northing,
double *latitude, double *longitude)
{
double a, b, e2, N, phiDash, phi, phi1, M;
double v, rho, n2, VII, VIII, IX, Et;
double X, XI, XII, XIIA;
// Initialise parameters
a = m_dSemiMajorAxis * m_dScaleFactorAtOrigin;
b = m_dSemiMinorAxis * m_dScaleFactorAtOrigin;
e2 = ((a * a) - (b * b)) / (a * a);
N = (a - b) / (a + b);
// Calculate latitude
phiDash = ((northing - m_dFalseNorthing) / a) + (m_dLatitudeOrigin * PI / 180);
phi = phiDash;
phi1 = (m_dLatitudeOrigin * PI / 180.0);
M = b * (((1.0 + N + ((pow2(N)) * 5.0 / 4) + ((pow3(N)) * 5.0 / 4.0)) * (phi - phi1))
- ((3.0 * N + 3.0 * (pow2(N)) + (pow3(N)) * 21.0 / 8.0) * sin(phi - phi1) * cos(phi + phi1))
+ ((((pow2(N)) * 15.0 / 8.0) + ((pow3(N)) * 15.0 / 8.0)) * (sin(2 * (phi - phi1))) * (cos(2.0 * (phi + phi1))))
- (((pow3(N)) * 35.0 / 24.0) * (sin(3.0 * (phi - phi1))) * (cos(3.0 * (phi + phi1)))));
while ((northing - m_dFalseNorthing - M) > 0.001)
{
phi = phi + ((northing - m_dFalseNorthing - M) / a);
M = b * (((1 + N + ((pow2(N)) * 5.0 / 4.0) + ((pow3(N)) * 5.0 / 4.0)) * (phi - phi1))
- ((3.0 * N + 3.0 * (pow2(N)) + (pow3(N)) * 21.0 / 8.0) * sin(phi - phi1) * cos(phi + phi1))
+ ((((pow2(N)) * 15.0 / 8.0) + ((pow3(N)) * 15.0 / 8.0)) * (sin(2.0 * (phi - phi1))) * (cos(2 * (phi + phi1))))
- (((pow3(N)) * 35.0 / 24.0) * (sin(3.0 * (phi - phi1))) * (cos(3.0 * (phi + phi1)))));
};
v = a / sqrt(1.0 - e2 * (pow2((sin(phi)))));
rho = (v * (1.0 - e2)) / (1 - e2 * (pow2(sin(phi))));
n2 = (v / rho) - 1.0;
VII = (tan(phi)) / (2.0 * rho * v);
VIII = ((tan(phi)) / (24.0 * rho * (pow3(v)))) * (5.0 + (3.0 * (pow2(tan(phi)))) + n2 - (9.0 * n2 * (pow2(tan(phi)))));
IX = (tan(phi) / (720.0 * rho * (pow5(v)))) * (61.0 + (90.0 * (pow2(tan(phi)))) + (45.0 * (pow4(tan(phi)))));
Et = easting - m_dFalseEasting;
*latitude = (phi - ((pow2(Et)) * VII) + ((pow4(Et)) * VIII) - ((pow6(Et)) * IX)) * 180.0 / PI;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -