[Cmap-cvs] shp2mp tr.h,NONE,1.1 tr.cpp,NONE,1.1

[Denis P. <dy...@us...>]

Update of /cvsroot/cmap/shp2mp
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv1875

Added Files:
	tr.h tr.cpp 
Log Message:
Add projection conversion code

--- NEW FILE: tr.h ---
#ifndef __TR_H__
#define __TR_H__

const double c_PI = 3.14159265358979;

class dpoint_t {
public:
	dpoint_t () : x(0.0), y(0.0) {};
	dpoint_t (double _x, double _y) : x(_x), y(_y) {};

	double x;
	double y;
};

struct CEllipsoid {
	const char *strName;

	// The semi-major axis.
	const double a;
	// The flattening.
	const double f;

	// The eccentricity^2.
	const double e2;
	const double e1, ee2;

	CEllipsoid (const char *_strName, double _a, double _f);
};

class CTransformationBase {
  public:
	virtual void ApplyTo (dpoint_t & _point) const = 0;
};


class CMolodensky_To_WGS84 : public CTransformationBase {
  public:
	CMolodensky_To_WGS84 (const CEllipsoid & _SourceEllipsoid, double _fDX, double _fDY, double _fDZ, double _fElevation = 0);
	virtual void ApplyTo (dpoint_t & _point) const;

  protected:
	// The source ellipsoid.
	const CEllipsoid e;
	// The linear shifts, in meters.
	const double dx, dy, dz;

	// The average and difference of semin-major axes.
	const double a, da;
	// The average and difference of eccentricities.
	const double e2, de2;

	const double fElevation;
};

class CHelmert_To_WGS84 : public CMolodensky_To_WGS84 {
  public:
	CHelmert_To_WGS84 (
		const CEllipsoid & _SourceEllipsoid,
		double _fDX, double _fDY, double _fDZ,
		double _fWX, double _fWY, double _fWZ,
		double _fMS,
		double _fElevation = 0
	);
	virtual void ApplyTo (dpoint_t & _point) const;

  protected:
	// The angles between axis, in arc seconds.
	const double wx, wy, wz;
	// The scales difference.
	const double ms;
};

class CPulkovo42_To_WGS84 : public CHelmert_To_WGS84 {
  public:
	CPulkovo42_To_WGS84 (double _fElevation = 0);
};

class CTransverseMercator : public CTransformationBase {
  public:
	CTransverseMercator (const CEllipsoid & _e, double _k0, double _FE, double _FN, double _x0, double _y0);
	virtual void ApplyTo (dpoint_t & _point) const;

  private:
	const CEllipsoid m_Ellipsoid;
	// The scale factor.
	const double m_k0;
	// False Easting & Northing.
	const double m_FE, m_FN;
	// Longitude of the natural origin.
	const double m_x0;
	// Some parameter depending on the latitude of the natural origin.
	const double m_M0;
};

extern const CEllipsoid c_eKrasovsky1940;
extern const CEllipsoid c_eWGS84;

#endif

--- NEW FILE: tr.cpp ---
#include <math.h>

#include "tr.h"

	// Arc seconds per radian.
const double c_ro = 206264.8062;

const CEllipsoid c_eKrasovsky1940   ("Krasovsky 1940",          6378245.0,    1/298.3);
const CEllipsoid c_eWGS84           ("WGS 84",                  6378137.0,    1/298.257223563);

CEllipsoid::CEllipsoid (const char *_strName, double _a, double _f) :
		strName (_strName),
		a (_a),
		f (_f),
		e2 (2*_f - _f*_f),
		e1 ((1 - ::sqrt (1 - e2))/(1 + ::sqrt (1 - e2))),
		ee2 (e2/(1 - e2))
	{}

CMolodensky_To_WGS84::CMolodensky_To_WGS84 (
	const CEllipsoid & _SourceEllipsoid,
	double _fDX, double _fDY, double _fDZ,
	double _fElevation
) :
	e (_SourceEllipsoid),
	dx (_fDX), dy (_fDY), dz (_fDZ),
	fElevation (_fElevation),
	a  ((_SourceEllipsoid.a  + c_eWGS84.a)/2),
	e2 ((_SourceEllipsoid.e2 + c_eWGS84.e2)/2),
	da  (c_eWGS84.a  - _SourceEllipsoid.a),
	de2 (c_eWGS84.e2 - _SourceEllipsoid.e2)
{
}

void CMolodensky_To_WGS84::ApplyTo (dpoint_t & _point) const {
	const double & Ld = _point.x;
	const double & Bd = _point.y;

	const double B = Bd*c_PI/180;
	const double L = Ld*c_PI/180;

	const double fSinB  = ::sin (B);
	const double fCosB  = ::cos (B);
	const double fSinL  = ::sin (L);
	const double fCosL  = ::cos (L);

	const double fNu = 1 - e2*fSinB*fSinB;
	const double fSqrtNu = ::sqrt (fNu);

	const double M = a*(1 - e2)/(fNu*fSqrtNu);
	const double N = a/fSqrtNu;

	//
	// Latitude.
	//
	const double dB = c_ro/(M + fElevation)*(
		N/a*e2*fSinB*fCosB*da
		+ (N*N/(a*a) + 1)*N*fSinB*fCosB*de2/2
		- (dx*fCosL + dy*fSinL)*fSinB
		+ dz*fCosB
	);
	const double fWGS84Lat = Bd + dB/3600;

	//
	// Longitude.
	//
	const double dL = c_ro/((N + fElevation)*fCosB)*(-dx*fSinL + dy*fCosL);
	const double fWGS84Long = Ld + dL/3600;

	_point.x = fWGS84Long;
	_point.y = fWGS84Lat;
}

///////////////////////////////////////////////////////

CHelmert_To_WGS84::CHelmert_To_WGS84 (
	const CEllipsoid & _SourceEllipsoid,
	double _fDX, double _fDY, double _fDZ,
	double _fWX, double _fWY, double _fWZ,
	double _fMS,
	double _fElevation
) :
	CMolodensky_To_WGS84 (_SourceEllipsoid, _fDX, _fDY, _fDZ, _fElevation),
	wx (_fWX), wy (_fWY), wz (_fWZ),
	ms (_fMS)
{}

void CHelmert_To_WGS84::ApplyTo (dpoint_t & _point) const {
	const double & Ld = _point.x;
	const double & Bd = _point.y;

	const double B = Bd*c_PI/180;
	const double L = Ld*c_PI/180;

	const double fSinB  = ::sin (B);
	const double fCosB  = ::cos (B);
	const double fTanB  = fSinB/fCosB;
	const double fCos2B = 2*fCosB*fCosB - 1;

	const double fSinL  = ::sin (L);
	const double fCosL  = ::cos (L);

	const double fNu = 1 - e2*fSinB*fSinB;
	const double fSqrtNu = ::sqrt (fNu);

	const double M = a*(1 - e2)/(fNu*fSqrtNu);
	const double N = a/fSqrtNu;

	//
	// Latitude.
	//
	const double dB = c_ro/(M + fElevation)*(
	    N/a*e2*fSinB*fCosB*da
		+ (N*N/(a*a) + 1)*N*fSinB*fCosB*de2/2
		- (dx*fCosL + dy*fSinL)*fSinB
		+ dz*fCosB
	)
		- wx*fSinL*(1 + e2*fCos2B)
		+ wy*fCosL*(1 + e2*fCos2B)
		- c_ro*ms*e2*fSinB*fCosB;

	const double fWGS84Lat = Bd + dB/3600;

	//
	// Longitude.
	//
	const double dL = c_ro/((N + fElevation)*fCosB)*(-dx*fSinL + dy*fCosL)
		+ fTanB*(1 - e2)*(wx*fCosL + wy*fSinL) - wz;

	const double fWGS84Long = Ld + dL/3600;
/*
	//
	// Elevation.
	//

	const double dH = -a/N*da + N*fSinB*fSinB*c_de2/2
		+ (dx*fCosL + dy*fSinL)*fCosB + dz*fSinB
		- N*e2*fSinB*fCosB*(wx/c_ro*fSinL - wy/c_ro*fCosL)
		+ (a*a/N + _fElevation)*ms;

	const double fWGS84Elevation = fElevation + dH;
*/
	_point.x = fWGS84Long;
	_point.y = fWGS84Lat;
}

///////////////////////////////////////////////

	// Based on the GOST R 51794-2001 of Russian Federation
CPulkovo42_To_WGS84::CPulkovo42_To_WGS84 (double _fElevation) :
	CHelmert_To_WGS84 (
		c_eKrasovsky1940,
		23.92f, -141.27f, -80.9f,
		0, -0.35f, -0.82f,
		-0.12e-6f,
		_fElevation
	)
{}

//
// Translate from TM grid (meters) to latitude/longitude degrees.
// The formula is given from http://posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34h.html
//
static
void TM2LL (const CEllipsoid & _e, double _k0, double _FE, double _FN, double _x0, double _M0, double _x, double _y, dpoint_t & _p) {
	const double M1 = _M0 + (_y - _FN)/_k0;

	const double mu1 = M1/(_e.a*(1 - _e.e2*(1./4 + _e.e2*(3./64 + 5./256*_e.e2))));
	const double fSin2Mu1 = ::sin (2*mu1);
	const double fCos2Mu1 = ::cos (2*mu1);

	const double e1_2 = _e.e1*_e.e1;
	const double e1_3 = e1_2*_e.e1;
	const double fi1 = mu1 +
		_e.e1*fSin2Mu1*(
			( 3./2        - 27./32*e1_2) + 
			(21./16*_e.e1 - 55./32*e1_3)*2*fCos2Mu1 +
			(              151./96*e1_2)*(3 - 4*fSin2Mu1*fSin2Mu1) +
			(            1097./512*e1_3)*4*fCos2Mu1*(1 - 2*fSin2Mu1*fSin2Mu1)
		);
	const double fSinFi1 = ::sin (fi1);
	const double fCosFi1 = ::cos (fi1);
	const double fTanFi1 = fSinFi1/fCosFi1;

	const double fEFi      = 1 - _e.e2*fSinFi1*fSinFi1;
	const double fSqrtEFi2 = ::sqrt (fEFi);

	const double D  = (_x - _FE)*fSqrtEFi2/(_e.a*_k0);
	const double D2 = D*D;
	const double T1 = fTanFi1*fTanFi1;
	const double C1 = _e.ee2*fCosFi1*fCosFi1;

	const double fLat = fi1 - fEFi*fTanFi1/(1 - _e.e2)*D2*(
		1./2 + D2*(
			- (5 + 3*T1 + C1*(10 - 4*C1) - 9*_e.ee2)/24 +
			(61 + T1*(90 + 45*T1) + C1*(298 - 3*C1) - 252*_e.ee2)*D2/720
		)
	);

	const double fLon = D*(
		1 + D2*(
			- (1 + 2*T1 + C1)/6 +
			(5 - C1*(2 + 3*C1) + T1*(28 + 24*T1) + 8*_e.ee2)*D2/120
		)
	)/fCosFi1;

	_p.y = fLat*(180/c_PI);
	_p.x = fLon*(180/c_PI) + _x0;
}

static
double GetM0 (const CEllipsoid & _e, double _y0) {
	const double fi = _y0*c_PI/180;
	return _e.a*(
		(1 -   _e.e2/4 -  3*_e.e2*_e.e2/64  -  5*_e.e2*_e.e2*_e.e2/256)*fi -
		(    3*_e.e2/8 +  3*_e.e2*_e.e2/32  + 45*_e.e2*_e.e2*_e.e2/1024)*::sin (2*fi) +
		(                15*_e.e2*_e.e2/256 + 45*_e.e2*_e.e2*_e.e2/1024)*::sin (4*fi) -
		(                                     35*_e.e2*_e.e2*_e.e2/3072)*::sin (6*fi)
	);
}

CTransverseMercator::CTransverseMercator (const CEllipsoid & _e, double _k0, double _FE, double _FN, double _x0, double _y0) :
	m_Ellipsoid (_e),
	m_k0 (_k0),
	m_FE (_FE),
	m_FN (_FN),
	m_x0 (_x0),
	m_M0 (GetM0 (_e, _y0))
{}

void CTransverseMercator::ApplyTo (dpoint_t & _point) const {
	TM2LL (m_Ellipsoid, m_k0, m_FE, m_FN, m_x0, m_M0, _point.x, _point.y, _point);
}