// Coords.h #ifndef __HEADER_COORDS_GEO__ #define __HEADER_COORDS_GEO__ namespace Geo ... {/**//*** 直角坐标点**/struct UPnt ...{ /**//** * x 坐标值 **/ double x; /**//** * y 坐标值 **/ double y;};/**//*** 大地坐标**/struct BLCoord...{ /**//** * 经度,单位为弧度 **/ double B; /**//** * 纬度,单位为弧度 **/ double L;};} #endif // __HEADER_COORDS_GEO__ // Ellipsoid.h #ifndef __HEADER_ELLIPSOID_GEO__ #define __HEADER_ELLIPSOID_GEO__ #include " Geo.h " namespace Geo ... {/**//*** 椭球参数**/struct Ellipsoid...{ /**//** * 长半坐标轴长 **/ double a; /**//** * 短半坐标轴长 **/ double b; /**//** * 返回扁率 **/ double Getf() ...{ return (a-b)/a; } /**//** * 返回第一偏心率 **/ double Gete1() const ...{ double rate = b/a; return sqrt(1.0-rate*rate); } /**//** * 返回第二偏心率 **/ double Gete2() const ...{ double rate = a/b; return sqrt(rate*rate-1.0); } /**//** * 计算卯酉圈曲率半径 * @param [in] B 纬度(弧度) **/ double GetN(double B) const ...{ double e1 = Gete1(); double e1sinB = e1*sin(B); return a/sqrt(1 - e1sinB*e1sinB); } /**//** * 计算子午圈曲率半径 * @param [in] B 纬度(弧度) **/ double GetR(double B) const ...{ double e1 = Gete1(); double e1sinB = e1*sin(B); double div = sqrt(1 - e1sinB*e1sinB); return a*(1-e1*e1)/(div*div*div); }};} #endif // __HEADER_ELLIPSOID_GEO__ // Geo.h #ifndef __HEADER_GEO__ #define __HEADER_GEO__ #ifdef GEO_EXPORTS #define GEO_API _declspec(dllexport) #else #define GEO_API _declspec(dllimport) #endif #include < math.h > #include < assert.h > #include < windows.h > #endif // __HEADER_GEO__ // Prj.h #ifndef __HEADER_PRJ_GEO__ #define __HEADER_PRJ_GEO__ #include " Ellipsoid.h " #include " Coords.h " namespace Geo ... {/**//*** 投影参数**/struct Prj...{ /**//** * 原点纬度 **/ double orgB; /**//** * 中央经度 **/ double centerL; /**//** * 东偏移 **/ double easteOffset; /**//** * 北偏移 **/ double northOffset; /**//** * 投影比例因子 **/ double k; /**//** * 保留参数 **/ double reserved0; double reserved1; double reserved2; };/**//*** 高斯投影正算* @param [in] ellipsoid 椭球参数* @param [in] prj 投影参数* @param [in] blCoords 经纬度坐标数组* @param [in] count 经纬度坐标数组的长度* @param [out] uPnts 传出计算的坐标点**/extern "C" GEO_API void __stdcall GaussProj(const Ellipsoid* ellipsoid, const Prj* prj, const BLCoord* blCoords, int count, UPnt* uPnts);/**//*** 高斯投影反算* @param [in] ellipsoid 椭球参数* @param [in] prj 投影参数* @param [in] blCoords 坐标点数组* @param [in] count 坐标点数组的长度* @param [out] uPnts 传出计算的经纬度坐标**/extern "C" GEO_API void __stdcall GaussProjInvert(const Ellipsoid* ellipsoid, const Prj* prj, const UPnt* uPnts, int count, BLCoord* blCoords);} #endif // __HEADER_PRJ_GEO__ // Publish.h #ifndef __HEADER_PUBLISH_GEO__ #define __HEADER_PUBLISH_GEO__ #include " Geo.h " #include " Ellipsoid.h " #include " Coords.h " #include " Prj.h " #include " Util.h " #endif // __HEADER_PUBLISH_GEO__ // Util.h #ifndef __HEADER_UTIL_GEO__ #define __HEADER_UTIL_GEO__ #include " Geo.h " namespace Geo ... {/**//*** 弧度转角度(如20.0940即20度,09分40秒)* @param [in] rad 弧度**/extern "C" GEO_API double __stdcall RadToDegree(double rad);/**//*** 角度转弧度* @param [in] degree 角度(格式同上) **/extern "C" GEO_API double __stdcall DegreeToRad(double degree);} #endif // __HEADER_UTIL_GEO__ // Geo.cpp #include " Geo.h " BOOL APIENTRY DllMain( HANDLE hModule, DWORD ul_reason_for_call, LPVOID lpReserved ) ... { return TRUE;} // Prj.cpp #include " Prj.h " namespace Geo ... {extern "C" GEO_API void __stdcall GaussProj(const Ellipsoid* ellipsoid, const Prj* prj, const BLCoord* blCoords, int count, UPnt* uPnts)...{ double e1 = ellipsoid->Gete1(); double e2 = ellipsoid->Gete2(); double ee2 = e1*e1; double ee4 = ee2*ee2; double ee6 = ee2*ee4; double k0 = (1 - ee2/4 - 3*ee4/64 - 5*ee6/256); double k1 = (-3*ee2/8 - 3*ee4/32 - 45*ee6/1024); double k2 = (15*ee4/256 + 45*ee6/1024); double k3 = (-35*ee6/3072); double B, L, T, C, A, M, N; for ( int cnt =0; cnt < count; cnt++ ) ...{ B = blCoords[cnt].B; L = blCoords[cnt].L; double tanB = tan(B); double cosB = cos(B); double e2cosB = e2*cosB; double sin2B = sin(2*B); double sin4B = sin(4*B); double sin6B = sin(6*B); double L0 = prj->centerL; T = tanB*tanB; C = e2cosB*e2cosB; A = (L-L0)*cosB; M = ellipsoid->a*(k0*B + k1*sin2B + k2*sin4B + k3*sin6B); N = ellipsoid->GetN(B); double AA2 = A*A; double AA4 = AA2*AA2; uPnts[cnt].x = prj->k*(M + N*tanB*(AA2/2 + (5 - T + 9*C + 4*C*C)*AA4/24) + (61 - 58*T + T*T + 270*C - 330*T*C)*AA2*AA4/720); uPnts[cnt].y = prj->easteOffset + prj->k*N*(A + (1 - T + C)*AA2*A/6 + (5 - 18*T + T*T + 14*C - 58*T*C)*AA4*A/120); }}extern "C" GEO_API void __stdcall GaussProjInvert(const Ellipsoid* ellipsoid, const Prj* prj, const UPnt* uPnts, int count, BLCoord* blCoords)...{ double e1 = ellipsoid->Gete1(); double e2 = ellipsoid->Gete2(); double ek = (1 - ellipsoid->b/ellipsoid->a)/(1 + ellipsoid->b/ellipsoid->a); double ee2 = e1*e1; double ee4 = ee2*ee2; double ee6 = ee2*ee4; double ek2 = ek*ek; double ek3 = ek2*ek; double ek4 = ek3*ek; double k0 = (3*ek/2 - 27*ek3/32); double k1 = (21*ek2/16 - 55*ek4/32); double k2 = (151*ek3/96); double XN, YE, Mf, Bf, Rf, Nf, Tf, Cf, D; for ( int cnt =0; cnt < count; cnt++ ) ...{ XN = uPnts->x; YE = uPnts->y; Mf = (XN - prj->northOffset)/prj->k; double rad = Mf/(ellipsoid->a*(1 - ee2/4 - 3*ee4/64 - 5*ee6/256)); double sin2rad = sin(2*rad); double sin4rad = sin(4*rad); double sin6rad = sin(6*rad); Bf = rad + k0*sin2rad + k1*sin4rad + k2*sin6rad; double div = sqrt(1 - ee2*sin(Bf)*sin(Bf)); Rf = ellipsoid->a*(1 - ee2)/(div*div*div); Nf = ellipsoid->a/div; Tf = tan(Bf); Tf *= Tf; Cf = e2*cos(Bf); Cf *= Cf; D = (YE - prj->easteOffset)/(prj->k*Nf); double D2 = D*D; double D4 = D2*D2; double D6 = D4*D2; blCoords[cnt].B = Bf - (Nf*tan(Bf)/Rf)*(D2/2 - (5 + 3*Tf + Cf - 9*Cf*Tf)*D4/24 + (61 + 90*Tf + 45*Tf*Tf)*D6/720); blCoords[cnt].L = prj->centerL + (D -(1 + 2*Tf + Cf)*D2*D/6 + (5 + 28*Tf + 6*Cf + 8*Tf*Cf + 24*Tf*Tf)*D4*D/120)/cos(Bf); }}} // Util.cpp #include " Util.h " namespace Geo ... {const double GEO_PI = 3.14159265358979323846;extern "C" GEO_API double __stdcall RadToDegree(double rad)...{ /**//*assert(rad);*/ double res = 0.0; double right = (rad*180.0)/GEO_PI; int left = int(right); res += left; right = right - left; right *= 60; left = int(right); res += left/100.0; right = right -left; right *= 60; left = int(right); res += left/10000.0; right = right -left; right *= 60; left = int(right); res += left/1000000.0; return res;}extern "C" GEO_API double __stdcall DegreeToRad(double degree)...{ /**//*assert(degree);*/ int left = int(degree); double right = degree - left; double res = (left/180.0)*GEO_PI; right *= 100.0; left = int(right); res += (left/10800.0)*GEO_PI; right = right - left; right *= 100.0; left = int(right); res += (left/648000.0)*GEO_PI; right = right - left; right *= 100.0; left = int(right); res += (left/38880000.0)*GEO_PI; return res;}}