/* =============================================================================== FILE: geoprojectionconverter.cpp CONTENTS: see corresponding header file PROGRAMMERS: martin isenburg@cs.unc.edu chuck.gantz@globalstar.com gpotts@imagelinks.com COPYRIGHT: copyright (C) 2007 martin isenburg@cs.unc.edu This software 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. This 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. CHANGE HISTORY: see corresponding header file =============================================================================== */ #include "geoprojectionconverter.h" #include #include #include #include static const double PI = 3.141592653589793238462643383279502884197169399375105820974944592; static const double TWO_PI = PI * 2; static const double PI_OVER_2 = PI / 2; static const double PI_OVER_4 = PI / 4; static const double deg2rad = PI / 180.0; static const double rad2deg = 180.0 / PI; static const double feet2meter = 0.3048; static const double surveyfeet2meter = 0.3048006096012; class ReferenceEllipsoid { public: ReferenceEllipsoid(int id, char* name, double equatorialRadius, double eccentricitySquared, double inverseFlattening) { this->id = id; this->name = name; this->equatorialRadius = equatorialRadius; this->eccentricitySquared = eccentricitySquared; this->inverseFlattening = inverseFlattening; } int id; char* name; double equatorialRadius; double eccentricitySquared; double inverseFlattening; }; static const ReferenceEllipsoid ellipsoid_list[] = { // d, Ellipsoid name, Equatorial Radius, square of eccentricity ReferenceEllipsoid( -1, "Placeholder", 0, 0, 0), //placeholder to allow array indices to match id numbers ReferenceEllipsoid( 1, "Airy", 6377563.396, 0.00667054, 299.3249646), ReferenceEllipsoid( 2, "Australian National", 6378160.0, 0.006694542, 298.25), ReferenceEllipsoid( 3, "Bessel 1841", 6377397.155, 0.006674372, 299.1528128), ReferenceEllipsoid( 4, "Bessel 1841 (Nambia) ", 6377483.865, 0.006674372, 299.1528128), ReferenceEllipsoid( 5, "Clarke 1866 (NAD-27)", 6378206.4, 0.006768658, 294.9786982), ReferenceEllipsoid( 6, "Clarke 1880", 6378249.145, 0.006803511, 293.465), ReferenceEllipsoid( 7, "Everest 1830", 6377276.345, 0.006637847, 300.8017), ReferenceEllipsoid( 8, "Fischer 1960 (Mercury) ", 6378166, 0.006693422, 298.3), ReferenceEllipsoid( 9, "Fischer 1968", 6378150, 0.006693422, 298.3), ReferenceEllipsoid( 10, "GRS 1967", 6378160, 0.006694605, 298.247167427), ReferenceEllipsoid( 11, "GRS 1980 (NAD-83)", 6378137, 0.00669438002290, 298.257222101), ReferenceEllipsoid( 12, "Helmert 1906", 6378200, 0.006693422, 298.3), ReferenceEllipsoid( 13, "Hough", 6378270, 0.00672267, 297.0), ReferenceEllipsoid( 14, "International", 6378388, 0.00672267, 297.0), ReferenceEllipsoid( 15, "Krassovsky", 6378245, 0.006693422, 298.3), ReferenceEllipsoid( 16, "Modified Airy", 6377340.189, 0.00667054, 299.3249646), ReferenceEllipsoid( 17, "Modified Everest", 6377304.063, 0.006637847, 300.8017), ReferenceEllipsoid( 18, "Modified Fischer 1960", 6378155, 0.006693422, 298.3), ReferenceEllipsoid( 19, "South American 1969", 6378160, 0.006694542, 298.25), ReferenceEllipsoid( 20, "WGS 60", 6378165, 0.006693422, 298.3), ReferenceEllipsoid( 21, "WGS 66", 6378145, 0.006694542, 298.25), ReferenceEllipsoid( 22, "WGS-72", 6378135, 0.006694318, 298.26), ReferenceEllipsoid( 23, "WGS-84", 6378137, 0.00669437999013, 298.257223563) }; class StatePlaneLCC { public: StatePlaneLCC(char* zone, double falseEastingMeter, double falseNorthingMeter, double latOriginDegree, double longMeridianDegree, double firstStdParallelDegree, double secondStdParallelDegree) { this->zone = zone; this->falseEastingMeter = falseEastingMeter; this->falseNorthingMeter = falseNorthingMeter; this->latOriginDegree = latOriginDegree; this->longMeridianDegree = longMeridianDegree; this->firstStdParallelDegree = firstStdParallelDegree; this->secondStdParallelDegree = secondStdParallelDegree; } char* zone; double falseEastingMeter; double falseNorthingMeter; double latOriginDegree; double longMeridianDegree; double firstStdParallelDegree; double secondStdParallelDegree; }; static const StatePlaneLCC state_plane_lcc_nad27_list[] = { // zone, false east [m], false north [m], ProjOrig(Lat), CentMerid(Long), 1st std para, 2nd std para StatePlaneLCC("AK_10",914401.8288,0,51,-176,51.83333333,53.83333333), StatePlaneLCC("AR_N",609601.2192,0,34.33333333,-92,34.93333333,36.23333333), StatePlaneLCC("AR_S",609601.2192,0,32.66666667,-92,33.3,34.76666667), StatePlaneLCC("CA_I",609601.2192,0,39.33333333,-122,40,41.66666667), StatePlaneLCC("CA_II",609601.2192,0,37.66666667,-122,38.33333333,39.83333333), StatePlaneLCC("CA_III",609601.2192,0,36.5,-120.5,37.06666667,38.43333333), StatePlaneLCC("CA_IV",609601.2192,0,35.33333333,-119,36,37.25), StatePlaneLCC("CA_V",609601.2192,0,33.5,-118,34.03333333,35.46666667), StatePlaneLCC("CA_VI",609601.2192,0,32.16666667,-116.25,32.78333333,33.88333333), StatePlaneLCC("CA_VII",1276106.451,1268253.007,34.13333333,-118.3333333,33.86666667,34.41666667), StatePlaneLCC("CO_N",609601.2192,0,39.33333333,-105.5,39.71666667,40.78333333), StatePlaneLCC("CO_C",609601.2192,0,37.83333333,-105.5,38.45,39.75), StatePlaneLCC("CO_S",609601.2192,0,36.66666667,-105.5,37.23333333,38.43333333), StatePlaneLCC("CT",182880.3658,0,40.83333333,-72.75,41.2,41.86666667), StatePlaneLCC("FL_N",609601.2192,0,29,-84.5,29.58333333,30.75), StatePlaneLCC("IA_N",609601.2192,0,41.5,-93.5,42.06666667,43.26666667), StatePlaneLCC("IA_S",609601.2192,0,40,-93.5,40.61666667,41.78333333), StatePlaneLCC("KS_N",609601.2192,0,38.33333333,-98,38.71666667,39.78333333), StatePlaneLCC("KS_S",609601.2192,0,36.66666667,-98.5,37.26666667,38.56666667), StatePlaneLCC("KY_N",609601.2192,0,37.5,-84.25,37.96666667,38.96666667), StatePlaneLCC("KY_S",609601.2192,0,36.33333333,-85.75,36.73333333,37.93333333), StatePlaneLCC("LA_N",609601.2192,0,30.66666667,-92.5,31.16666667,32.66666667), StatePlaneLCC("LA_S",609601.2192,0,28.66666667,-91.33333333,29.3,30.7), StatePlaneLCC("LA_O",609601.2192,0,25.66666667,-91.33333333,26.16666667,27.83333333), StatePlaneLCC("MD",243840.4877,0,37.83333333,-77,38.3,39.45), StatePlaneLCC("MA_M",182880.3658,0,41,-71.5,41.71666667,42.68333333), StatePlaneLCC("MA_I",60960.12192,0,41,-70.5,41.28333333,41.48333333), StatePlaneLCC("MI_N",609601.2192,0,44.78333333,-87,45.48333333,47.08333333), StatePlaneLCC("MI_C",609601.2192,0,43.31666667,-84.33333333,44.18333333,45.7), StatePlaneLCC("MI_S",609601.2192,0,41.5,-84.33333333,42.1,43.66666667), StatePlaneLCC("MN_N",609601.2192,0,46.5,-93.1,47.03333333,48.63333333), StatePlaneLCC("MN_C",609601.2192,0,45,-94.25,45.61666667,47.05), StatePlaneLCC("MN_S",609601.2192,0,43,-94,43.78333333,45.21666667), StatePlaneLCC("MT_N",609601.2192,0,47,-109.5,47.85,48.71666667), StatePlaneLCC("MT_C",609601.2192,0,45.83333333,-109.5,46.45,47.88333333), StatePlaneLCC("MT_S",609601.2192,0,44,-109.5,44.86666667,46.4), StatePlaneLCC("NE_N",609601.2192,0,41.33333333,-100,41.85,42.81666667), StatePlaneLCC("NE_S",609601.2192,0,39.66666667,-99.5,40.28333333,41.71666667), StatePlaneLCC("NY_LI",609601.2192,30480.06096,40.5,-74,40.66666667,41.03333333), StatePlaneLCC("NC",609601.2192,0,33.75,-79,34.33333333,36.16666667), StatePlaneLCC("ND_N",609601.2192,0,47,-100.5,47.43333333,48.73333333), StatePlaneLCC("ND_S",609601.2192,0,45.66666667,-100.5,46.18333333,47.48333333), StatePlaneLCC("OH_N",609601.2192,0,39.66666667,-82.5,40.43333333,41.7), StatePlaneLCC("OH_S",609601.2192,0,38,-82.5,38.73333333,40.03333333), StatePlaneLCC("OK_N",609601.2192,0,35,-98,35.56666667,36.76666667), StatePlaneLCC("OK_S",609601.2192,0,33.33333333,-98,33.93333333,35.23333333), StatePlaneLCC("OR_N",609601.2192,0,43.66666667,-120.5,44.33333333,46), StatePlaneLCC("OR_S",609601.2192,0,41.66666667,-120.5,42.33333333,44), StatePlaneLCC("PA_N",609601.2192,0,40.16666667,-77.75,40.88333333,41.95), StatePlaneLCC("PA_S",609601.2192,0,39.33333333,-77.75,39.93333333,40.96666667), StatePlaneLCC("PR",152400.3048,0,17.83333333,-66.43333333,18.03333333,18.43333333), StatePlaneLCC("St.Croix",152400.3048,30480.06096,17.83333333,-66.43333333,18.03333333,18.43333333), StatePlaneLCC("SC_N",609601.2192,0,33,-81,33.76666667,34.96666667), StatePlaneLCC("SC_S",609601.2192,0,31.83333333,-81,32.33333333,33.66666667), StatePlaneLCC("SD_N",609601.2192,0,43.83333333,-100,44.41666667,45.68333333), StatePlaneLCC("SD_S",609601.2192,0,42.33333333,-100.3333333,42.83333333,44.4), StatePlaneLCC("TN",609601.2192,30480.06096,34.66666667,-86,35.25,36.41666667), StatePlaneLCC("TX_N",609601.2192,0,34,-101.5,34.65,36.18333333), StatePlaneLCC("TX_NC",609601.2192,0,31.66666667,-97.5,32.13333333,33.96666667), StatePlaneLCC("TX_C",609601.2192,0,29.66666667,-100.3333333,30.11666667,31.88333333), StatePlaneLCC("TX_SC",609601.2192,0,27.83333333,-99,28.38333333,30.28333333), StatePlaneLCC("TX_S",609601.2192,0,25.66666667,-98.5,26.16666667,27.83333333), StatePlaneLCC("UT_N",609601.2192,0,40.33333333,-111.5,40.71666667,41.78333333), StatePlaneLCC("UT_C",609601.2192,0,38.33333333,-111.5,39.01666667,40.65), StatePlaneLCC("UT_S",609601.2192,0,36.66666667,-111.5,37.21666667,38.35), StatePlaneLCC("VA_N",609601.2192,0,37.66666667,-78.5,38.03333333,39.2), StatePlaneLCC("VA_S",609601.2192,0,36.33333333,-78.5,36.76666667,37.96666667), StatePlaneLCC("WA_N",609601.2192,0,47,-120.8333333,47.5,48.73333333), StatePlaneLCC("WA_S",609601.2192,0,45.33333333,-120.5,45.83333333,47.33333333), StatePlaneLCC("WV_N",609601.2192,0,38.5,-79.5,39,40.25), StatePlaneLCC("WV_S",609601.2192,0,37,-81,37.48333333,38.88333333), StatePlaneLCC("WI_N",609601.2192,0,45.16666667,-90,45.56666667,46.76666667), StatePlaneLCC("WI_C",609601.2192,0,43.83333333,-90,44.25,45.5), StatePlaneLCC("WI_S",609601.2192,0,42,-90,42.73333333,44.06666667), StatePlaneLCC(0,-1,-1,-1,-1,-1,-1) }; static const StatePlaneLCC state_plane_lcc_nad83_list[] = { // zone, false east [m], false north [m], ProjOrig(Lat), CentMerid(Long), 1st std para, 2nd std para StatePlaneLCC("AK_10",1000000,0,51.000000,-176.000000,51.833333,53.833333), StatePlaneLCC("AR_N",400000,0,34.333333,-92.000000,34.933333,36.233333), StatePlaneLCC("AR_S",400000,400000,32.666667,-92.000000,33.300000,34.766667), StatePlaneLCC("CA_I",2000000,500000,39.333333,-122.000000,40.000000,41.666667), StatePlaneLCC("CA_II",2000000,500000,37.666667,-122.000000,38.333333,39.833333), StatePlaneLCC("CA_III",2000000,500000,36.500000,-120.500000,37.066667,38.433333), StatePlaneLCC("CA_IV",2000000,500000,35.333333,-119.000000,36.000000,37.250000), StatePlaneLCC("CA_V",2000000,500000,33.500000,-118.000000,34.033333,35.466667), StatePlaneLCC("CA_VI",2000000,500000,32.166667,-116.250000,32.783333,33.883333), StatePlaneLCC("CO_N",914401.8289,304800.6096,39.333333,-105.500000,39.716667,40.783333), StatePlaneLCC("CO_C",914401.8289,304800.6096,37.833333,-105.500000,38.450000,39.750000), StatePlaneLCC("CO_S",914401.8289,304800.6096,36.666667,-105.500000,37.233333,38.433333), StatePlaneLCC("CT",304800.6096,152400.3048,40.833333,-72.750000,41.200000,41.866667), StatePlaneLCC("FL_N",600000,0,29.000000,-84.500000,29.583333,30.750000), StatePlaneLCC("IA_N",1500000,1000000,41.500000,-93.500000,42.066667,43.266667), StatePlaneLCC("IA_S",500000,0,40.000000,-93.500000,40.616667,41.783333), StatePlaneLCC("KS_N",400000,0,38.333333,-98.000000,38.716667,39.783333), StatePlaneLCC("KS_S",400000,400000,36.666667,-98.500000,37.266667,38.566667), StatePlaneLCC("KY_N",500000,0,37.500000,-84.250000,37.966667,38.966667), StatePlaneLCC("KY_S",500000,500000,36.333333,-85.750000,36.733333,37.933333), StatePlaneLCC("LA_N",1000000,0,30.500000,-92.500000,31.166667,32.666667), StatePlaneLCC("LA_S",1000000,0,28.500000,-91.333333,29.300000,30.700000), StatePlaneLCC("LA_O",1000000,0,25.500000,-91.333333,26.166667,27.833333), StatePlaneLCC("MD",400000,0,37.666667,-77.000000,38.300000,39.450000), StatePlaneLCC("MA_M",200000,750000,41.000000,-71.500000,41.716667,42.683333), StatePlaneLCC("MA_I",500000,0,41.000000,-70.500000,41.283333,41.483333), StatePlaneLCC("MI_N",8000000,0,44.783333,-87.000000,45.483333,47.083333), StatePlaneLCC("MI_C",6000000,0,43.316667,-84.366667,44.183333,45.700000), StatePlaneLCC("MI_S",4000000,0,41.500000,-84.366667,42.100000,43.666667), StatePlaneLCC("MN_N",800000,100000,46.500000,-93.100000,47.033333,48.633333), StatePlaneLCC("MN_C",800000,100000,45.000000,-94.250000,45.616667,47.050000), StatePlaneLCC("MN_S",800000,100000,43.000000,-94.000000,43.783333,45.216667), StatePlaneLCC("MT",600000,0,44.250000,-109.500000,45.000000,49.000000), StatePlaneLCC("NE",500000,0,39.833333,-100.000000,40.000000,43.000000), StatePlaneLCC("NY_LI",300000,0,40.166667,-74.000000,40.666667,41.033333), StatePlaneLCC("NC",609601.22,0,33.750000,-79.000000,34.333333,36.166667), StatePlaneLCC("ND_N",600000,0,47.000000,-100.500000,47.433333,48.733333), StatePlaneLCC("ND_S",600000,0,45.666667,-100.500000,46.183333,47.483333), StatePlaneLCC("OH_N",600000,0,39.666667,-82.500000,40.433333,41.700000), StatePlaneLCC("OH_S",600000,0,38.000000,-82.500000,38.733333,40.033333), StatePlaneLCC("OK_N",600000,0,35.000000,-98.000000,35.566667,36.766667), StatePlaneLCC("OK_S",600000,0,33.333333,-98.000000,33.933333,35.233333), StatePlaneLCC("OR_N",2500000,0,43.666667,-120.500000,44.333333,46.000000), StatePlaneLCC("OR_S",1500000,0,41.666667,-120.500000,42.333333,44.000000), StatePlaneLCC("PA_N",600000,0,40.166667,-77.750000,40.883333,41.950000), StatePlaneLCC("PA_S",600000,0,39.333333,-77.750000,39.933333,40.966667), StatePlaneLCC("PR",200000,200000,17.833333,-66.433333,18.033333,18.433333), StatePlaneLCC("SC",609600,0,31.833333,-81.000000,32.500000,34.833333), StatePlaneLCC("SD_N",600000,0,43.833333,-100.000000,44.416667,45.683333), StatePlaneLCC("SD_S",600000,0,42.333333,-100.333333,42.833333,44.400000), StatePlaneLCC("TN",600000,0,34.333333,-86.000000,35.250000,36.416667), StatePlaneLCC("TX_N",200000,1000000,34.000000,-101.500000,34.650000,36.183333), StatePlaneLCC("TX_NC",600000,2000000,31.666667,-98.500000,32.133333,33.966667), StatePlaneLCC("TX_C",700000,3000000,29.666667,-100.333333,30.116667,31.883333), StatePlaneLCC("TX_SC",600000,4000000,27.833333,-99.000000,28.383333,30.283333), StatePlaneLCC("TX_S",300000,5000000,25.666667,-98.500000,26.166667,27.833333), StatePlaneLCC("UT_N",500000,1000000,40.333333,-111.500000,40.716667,41.783333), StatePlaneLCC("UT_C",500000,2000000,38.333333,-111.500000,39.016667,40.650000), StatePlaneLCC("UT_S",500000,3000000,36.666667,-111.500000,37.216667,38.350000), StatePlaneLCC("VA_N",3500000,2000000,37.666667,-78.500000,38.033333,39.200000), StatePlaneLCC("VA_S",3500000,1000000,36.333333,-78.500000,36.766667,37.966667), StatePlaneLCC("WA_N",500000,0,47.000000,-120.833333,47.500000,48.733333), StatePlaneLCC("WA_S",500000,0,45.333333,-120.500000,45.833333,47.333333), StatePlaneLCC("WV_N",600000,0,38.500000,-79.500000,39.000000,40.250000), StatePlaneLCC("WV_S",600000,0,37.000000,-81.000000,37.483333,38.883333), StatePlaneLCC("WI_N",600000,0,45.166667,-90.000000,45.566667,46.766667), StatePlaneLCC("WI_C",600000,0,43.833333,-90.000000,44.250000,45.500000), StatePlaneLCC("WI_S",600000,0,42.000000,-90.000000,42.733333,44.066667), StatePlaneLCC(0,-1,-1,-1,-1,-1,-1) }; class StatePlaneTM { public: StatePlaneTM(char* zone, double falseEastingMeter, double falseNorthingMeter, double latOriginDegree, double longMeridianDegree, double scaleFactor) { this->zone = zone; this->falseEastingMeter = falseEastingMeter; this->falseNorthingMeter = falseNorthingMeter; this->latOriginDegree = latOriginDegree; this->longMeridianDegree = longMeridianDegree; this->scaleFactor = scaleFactor; } char* zone; double falseEastingMeter; double falseNorthingMeter; double latOriginDegree; double longMeridianDegree; double scaleFactor; }; static const StatePlaneTM state_plane_tm_nad27_list[] = { // zone, false east [m], false north [m], ProjOrig(Lat), CentMerid(Long), scale factor StatePlaneTM("AL_E",152400.3048,0,30.5,-85.83333333,0.99996), StatePlaneTM("AL_W",152400.3048,0,30,-87.5,0.999933333), StatePlaneTM("AK_2",152400.3048,0,54,-142,0.9999), StatePlaneTM("AK_3",152400.3048,0,54,-146,0.9999), StatePlaneTM("AK_4",152400.3048,0,54,-150,0.9999), StatePlaneTM("AK_5",152400.3048,0,54,-154,0.9999), StatePlaneTM("AK_6",152400.3048,0,54,-158,0.9999), StatePlaneTM("AK_7",213360.4267,0,54,-162,0.9999), StatePlaneTM("AK_8",152400.3048,0,54,-166,0.9999), StatePlaneTM("AK_9",182880.3658,0,54,-170,0.9999), StatePlaneTM("AZ_E",152400.3048,0,31,-110.1666667,0.9999), StatePlaneTM("AZ_C",152400.3048,0,31,-111.9166667,0.9999), StatePlaneTM("AZ_W",152400.3048,0,31,-113.75,0.999933333), StatePlaneTM("DE",152400.3048,0,38,-75.41666667,0.999995), StatePlaneTM("FL_E",152400.3048,0,24.33333333,-81,0.999941177), StatePlaneTM("FL_W",152400.3048,0,24.33333333,-82,0.999941177), StatePlaneTM("GA_E",152400.3048,0,30,-82.16666667,0.9999), StatePlaneTM("GA_W",152400.3048,0,30,-84.16666667,0.9999), StatePlaneTM("HI_1",152400.3048,0,18.83333333,-155.5,0.999966667), StatePlaneTM("HI_2",152400.3048,0,20.33333333,-156.6666667,0.999966667), StatePlaneTM("HI_3",152400.3048,0,21.16666667,-158,0.99999), StatePlaneTM("HI_4",152400.3048,0,21.83333333,-159.5,0.99999), StatePlaneTM("HI_5",152400.3048,0,21.66666667,-160.1666667,1), StatePlaneTM("ID_E",152400.3048,0,41.66666667,-112.1666667,0.999947368), StatePlaneTM("ID_C",152400.3048,0,41.66666667,-114,0.999947368), StatePlaneTM("ID_W",152400.3048,0,41.66666667,-115.75,0.999933333), StatePlaneTM("IL_E",152400.3048,0,36.66666667,-88.33333333,0.999975), StatePlaneTM("IL_W",152400.3048,0,36.66666667,-90.16666667,0.999941177), StatePlaneTM("IN_E",152400.3048,0,37.5,-85.66666667,0.999966667), StatePlaneTM("IN_W",152400.3048,0,37.5,-87.08333333,0.999966667), StatePlaneTM("ME_E",152400.3048,0,43.83333333,-68.5,0.9999), StatePlaneTM("ME_W",152400.3048,0,42.83333333,-70.16666667,0.999966667), StatePlaneTM("MI_E",152400.3048,0,41.5,-83.66666667,0.999942857), StatePlaneTM("MI_C",152400.3048,0,41.5,-85.75,0.999909091), StatePlaneTM("MI_W",152400.3048,0,41.5,-88.75,0.999909091), StatePlaneTM("MS_E",152400.3048,0,29.66666667,-88.83333333,0.99996), StatePlaneTM("MS_W",152400.3048,0,30.5,-90.33333333,0.999941177), StatePlaneTM("MO_E",152400.3048,0,35.83333333,-90.5,0.999933333), StatePlaneTM("MO_C",152400.3048,0,35.83333333,-92.5,0.999933333), StatePlaneTM("MO_W",152400.3048,0,36.16666667,-94.5,0.999941177), StatePlaneTM("NV_E",152400.3048,0,34.75,-115.5833333,0.9999), StatePlaneTM("NV_C",152400.3048,0,34.75,-116.6666667,0.9999), StatePlaneTM("NV_W",152400.3048,0,34.75,-118.5833333,0.9999), StatePlaneTM("NH",152400.3048,0,42.5,-71.66666667,0.999966667), StatePlaneTM("NJ",609601.2192,0,38.83333333,-74.66666667,0.999975), StatePlaneTM("NM_E",152400.3048,0,31,-104.3333333,0.999909091), StatePlaneTM("NM_C",152400.3048,0,31,-106.25,0.9999), StatePlaneTM("NM_W",152400.3048,0,31,-107.8333333,0.999916667), StatePlaneTM("NY_E",152400.3048,0,40,-74.33333333,0.999966667), StatePlaneTM("NY_C",152400.3048,0,40,-76.58333333,0.9999375), StatePlaneTM("NY_W",152400.3048,0,40,-78.58333333,0.9999375), StatePlaneTM("RI",152400.3048,0,41.08333333,-71.5,0.99999375), StatePlaneTM("VT",152400.3048,0,42.5,-72.5,0.999964286), StatePlaneTM("WY_E",152400.3048,0,40.66666667,-105.1666667,0.999941177), StatePlaneTM("WY_EC",152400.3048,0,40.66666667,-107.3333333,0.999941177), StatePlaneTM("WY_WC",152400.3048,0,40.66666667,-108.75,0.999941177), StatePlaneTM("WY_W",152400.3048,0,40.66666667,-110.0833333,0.999941177), StatePlaneTM(0,-1,-1,-1,-1,-1) }; static const StatePlaneTM state_plane_tm_nad83_list[] = { // zone, false east [m], false north [m], ProjOrig(Lat), CentMerid(Long), scale factor StatePlaneTM("AL_E",200000,0,30.5,-85.83333333,0.99996), StatePlaneTM("AL_W",600000,0,30,-87.5,0.999933333), StatePlaneTM("AK_2",500000,0,54,-142,0.9999), StatePlaneTM("AK_3",500000,0,54,-146,0.9999), StatePlaneTM("AK_4",500000,0,54,-150,0.9999), StatePlaneTM("AK_5",500000,0,54,-154,0.9999), StatePlaneTM("AK_6",500000,0,54,-158,0.9999), StatePlaneTM("AK_7",500000,0,54,-162,0.9999), StatePlaneTM("AK_8",500000,0,54,-166,0.9999), StatePlaneTM("AK_9",500000,0,54,-170,0.9999), StatePlaneTM("AZ_E",213360,0,31,-110.1666667,0.9999), StatePlaneTM("AZ_C",213360,0,31,-111.9166667,0.9999), StatePlaneTM("AZ_W",213360,0,31,-113.75,0.999933333), StatePlaneTM("DE",200000,0,38,-75.41666667,0.999995), StatePlaneTM("FL_E",200000,0,24.33333333,-81,0.999941177), StatePlaneTM("FL_W",200000,0,24.33333333,-82,0.999941177), StatePlaneTM("GA_E",200000,0,30,-82.16666667,0.9999), StatePlaneTM("GA_W",700000,0,30,-84.16666667,0.9999), StatePlaneTM("HI_1",500000,0,18.83333333,-155.5,0.999966667), StatePlaneTM("HI_2",500000,0,20.33333333,-156.6666667,0.999966667), StatePlaneTM("HI_3",500000,0,21.16666667,-158,0.99999), StatePlaneTM("HI_4",500000,0,21.83333333,-159.5,0.99999), StatePlaneTM("HI_5",500000,0,21.66666667,-160.1666667,1), StatePlaneTM("ID_E",200000,0,41.66666667,-112.1666667,0.999947368), StatePlaneTM("ID_C",500000,0,41.66666667,-114,0.999947368), StatePlaneTM("ID_W",800000,0,41.66666667,-115.75,0.999933333), StatePlaneTM("IL_E",300000,0,36.66666667,-88.33333333,0.999975), StatePlaneTM("IL_W",700000,0,36.66666667,-90.16666667,0.999941177), StatePlaneTM("IN_E",100000,250000,37.5,-85.66666667,0.999966667), StatePlaneTM("IN_W",900000,250000,37.5,-87.08333333,0.999966667), StatePlaneTM("ME_E",300000,0,43.66666667,-68.5,0.9999), StatePlaneTM("ME_W",900000,0,42.83333333,-70.16666667,0.999966667), StatePlaneTM("MI_E",500000,0,41.5,-83.66666667,0.999942857), StatePlaneTM("MI_C",500000,0,41.5,-85.75,0.999909091), StatePlaneTM("MI_W",500000,0,41.5,-88.75,0.999909091), StatePlaneTM("MS_E",300000,0,29.5,-88.83333333,0.99995), StatePlaneTM("MS_W",700000,0,29.5,-90.33333333,0.99995), StatePlaneTM("MO_E",250000,0,35.83333333,-90.5,0.999933333), StatePlaneTM("MO_C",500000,0,35.83333333,-92.5,0.999933333), StatePlaneTM("MO_W",850000,0,36.16666667,-94.5,0.999941177), StatePlaneTM("NV_E",200000,8000000,34.75,-115.5833333,0.9999), StatePlaneTM("NV_C",500000,6000000,34.75,-116.6666667,0.9999), StatePlaneTM("NV_W",800000,4000000,34.75,-118.5833333,0.9999), StatePlaneTM("NH",300000,0,42.5,-71.66666667,0.999966667), StatePlaneTM("NJ",150000,0,38.83333333,-74.5,0.9999), StatePlaneTM("NM_E",165000,0,31,-104.3333333,0.999909091), StatePlaneTM("NM_C",500000,0,31,-106.25,0.9999), StatePlaneTM("NM_W",830000,0,31,-107.8333333,0.999916667), StatePlaneTM("NY_E",150000,0,38.83333333,-74.5,0.9999), StatePlaneTM("NY_C",250000,0,40,-76.58333333,0.9999375), StatePlaneTM("NY_W",350000,0,40,-78.58333333,0.9999375), StatePlaneTM("RI",100000,0,41.08333333,-71.5,0.99999375), StatePlaneTM("VT",500000,0,42.5,-72.5,0.999964286), StatePlaneTM("WY_E",200000,0,40.5,-105.1666667,0.9999375), StatePlaneTM("WY_EC",400000,100000,40.5,-107.3333333,0.9999375), StatePlaneTM("WY_WC",600000,0,40.5,-108.75,0.9999375), StatePlaneTM("WY_W",800000,100000,40.5,-110.0833333,0.9999375), StatePlaneTM(0,-1,-1,-1,-1,-1) }; bool GeoProjectionConverter::set_utm_projection(char* zone, char* description) { char* letter; utm_zone_number = strtoul(zone, &letter, 10); utm_zone_letter = *letter; if (utm_zone_letter < 'C' || utm_zone_letter > 'X') { return false; } sprintf(projection_name,"UTM zone %d%c",utm_zone_number,utm_zone_letter); if((*letter - 'N') >= 0) { utm_northern_hemisphere = true; // point is in northern hemisphere } else { utm_northern_hemisphere = false; //point is in southern hemisphere } utm_long_origin = (utm_zone_number - 1) * 6 - 180 + 3; // + 3 puts origin in middle of zone if (description) { sprintf(description, "%d%c - %s", utm_zone_number, utm_zone_letter, (utm_northern_hemisphere ? "northern hemisphere" : "southern hemisphere")); } return true; } bool GeoProjectionConverter::set_reference_ellipsoid(int id, char* description) { if (id <= 0 || id >= 24) { return false; } ellipsoid_id = id; ellipsoid_name = ellipsoid_list[id].name; equatorial_radius = ellipsoid_list[id].equatorialRadius; eccentricity_squared = ellipsoid_list[id].eccentricitySquared; eccentricity_prime_squared = (eccentricity_squared)/(1-eccentricity_squared); polar_radius = equatorial_radius*sqrt(1-eccentricity_squared); eccentricity = sqrt(eccentricity_squared); eccentricity_e1 = (1-sqrt(1-eccentricity_squared))/(1+sqrt(1-eccentricity_squared)); // recompute LCC parameters compute_lcc_parameters(); // recompute TM parameters compute_tm_parameters(); if (description) { sprintf(description, "%2d - %s (%g %g)", ellipsoid_id, ellipsoid_name, equatorial_radius, eccentricity_squared); } return true; } const char* GeoProjectionConverter::get_ellipsoid_name() { return ellipsoid_name; } // Configure a Lambert Conic Conformal Projection // // The function Set_Lambert_Parameters receives the ellipsoid parameters and // Lambert Conformal Conic projection parameters as inputs, and sets the // corresponding state variables. // // falseEastingMeter & falseNorthingMeter are just an offset in meters added // to the final coordinate calculated. // // latOriginDegree & longMeridianDegree are the "center" latitiude and // longitude in decimal degrees of the area being projected. All coordinates // will be calculated in meters relative to this point on the earth. // // firstStdParallelDegree & secondStdParallelDegree are the two lines of // longitude in decimal degrees (that is they run east-west) that define // where the "cone" intersects the earth. They bracket the area being projected. bool GeoProjectionConverter::set_lambert_conformal_conic_projection(double falseEastingMeter, double falseNorthingMeter, double latOriginDegree, double longMeridianDegree, double firstStdParallelDegree, double secondStdParallelDegree, char* description) { lcc_false_easting_meter = falseEastingMeter; lcc_false_northing_meter = falseNorthingMeter; lcc_lat_origin_degree = latOriginDegree; lcc_long_meridian_degree = longMeridianDegree; lcc_first_std_parallel_degree = firstStdParallelDegree; lcc_second_std_parallel_degree = secondStdParallelDegree; lcc_lat_origin_radian = deg2rad*lcc_lat_origin_degree; lcc_long_meridian_radian = deg2rad*lcc_long_meridian_degree; lcc_first_std_parallel_radian = deg2rad*lcc_first_std_parallel_degree; lcc_second_std_parallel_radian = deg2rad*lcc_second_std_parallel_degree; compute_lcc_parameters(); sprintf(projection_name,"Lambert Conformal Conic"); if (description) { sprintf(description, "false east/north: %g/%g [m], origin lat/ meridian long: %g/%g, parallel 1st/2nd: %g/%g", lcc_false_easting_meter, lcc_false_northing_meter, lcc_lat_origin_degree, lcc_long_meridian_degree, lcc_first_std_parallel_degree, lcc_second_std_parallel_degree); } return true; } /* * The function set_transverse_mercator_projection receives the Tranverse * Mercator projection parameters as input and sets the corresponding state * variables. * falseEastingMeter : Easting/X in meters at the center of the projection * falseNorthingMeter : Northing/Y in meters at the center of the projection * latOriginDegree : Latitude in decimal degree at the origin of the projection * longMeridianDegree : Longitude n decimal degree at the center of the projection * scaleFactor : Projection scale factor */ bool GeoProjectionConverter::set_transverse_mercator_projection(double falseEastingMeter, double falseNorthingMeter, double latOriginDegree, double longMeridianDegree, double scaleFactor, char* description) { tm_false_easting_meter = falseEastingMeter; tm_false_northing_meter = falseNorthingMeter; tm_lat_origin_degree = latOriginDegree; tm_long_meridian_degree = longMeridianDegree; tm_scale_factor = scaleFactor; tm_lat_origin_radian = deg2rad*tm_lat_origin_degree; tm_long_meridian_radian = deg2rad*tm_long_meridian_degree; compute_tm_parameters(); sprintf(projection_name,"Transverse Mercator"); if (description) { sprintf(description, "false east/north: %g/%g [m], origin lat/ meridian long: %g/%g, scale factor: %g", tm_false_easting_meter, tm_false_northing_meter, tm_lat_origin_degree, tm_long_meridian_degree, tm_scale_factor); } return true; } bool GeoProjectionConverter::set_state_plane_nad27_lcc(const char* zone, char* description) { int i = 0; while (state_plane_lcc_nad27_list[i].zone) { if (strcmp(zone, state_plane_lcc_nad27_list[i].zone) == 0) { set_reference_ellipsoid(5); return set_lambert_conformal_conic_projection(state_plane_lcc_nad27_list[i].falseEastingMeter, state_plane_lcc_nad27_list[i].falseNorthingMeter,state_plane_lcc_nad27_list[i].latOriginDegree,state_plane_lcc_nad27_list[i].longMeridianDegree,state_plane_lcc_nad27_list[i].firstStdParallelDegree,state_plane_lcc_nad27_list[i].secondStdParallelDegree,description); } i++; } return false; } void GeoProjectionConverter::print_all_state_plane_nad27_lcc() { int i = 0; while (state_plane_lcc_nad27_list[i].zone) { fprintf(stderr, "%s - false east/north: %g/%g [m], origin lat/meridian long: %g/%g, parallel 1st/2nd: %g/%g\n", state_plane_lcc_nad27_list[i].zone, state_plane_lcc_nad27_list[i].falseEastingMeter, state_plane_lcc_nad27_list[i].falseNorthingMeter,state_plane_lcc_nad27_list[i].latOriginDegree,state_plane_lcc_nad27_list[i].longMeridianDegree,state_plane_lcc_nad27_list[i].firstStdParallelDegree,state_plane_lcc_nad27_list[i].secondStdParallelDegree); i++; } } bool GeoProjectionConverter::set_state_plane_nad83_lcc(const char* zone, char* description) { int i = 0; while (state_plane_lcc_nad83_list[i].zone) { if (strcmp(zone, state_plane_lcc_nad83_list[i].zone) == 0) { set_reference_ellipsoid(11); return set_lambert_conformal_conic_projection(state_plane_lcc_nad83_list[i].falseEastingMeter, state_plane_lcc_nad83_list[i].falseNorthingMeter,state_plane_lcc_nad83_list[i].latOriginDegree,state_plane_lcc_nad83_list[i].longMeridianDegree,state_plane_lcc_nad83_list[i].firstStdParallelDegree,state_plane_lcc_nad83_list[i].secondStdParallelDegree,description); } i++; } return false; } void GeoProjectionConverter::print_all_state_plane_nad83_lcc() { int i = 0; while (state_plane_lcc_nad83_list[i].zone) { fprintf(stderr, "%s - false east/north: %g/%g [m], origin lat/meridian long: %g/%g, parallel 1st/2nd: %g/%g\n", state_plane_lcc_nad83_list[i].zone, state_plane_lcc_nad83_list[i].falseEastingMeter, state_plane_lcc_nad83_list[i].falseNorthingMeter,state_plane_lcc_nad83_list[i].latOriginDegree,state_plane_lcc_nad83_list[i].longMeridianDegree,state_plane_lcc_nad83_list[i].firstStdParallelDegree,state_plane_lcc_nad83_list[i].secondStdParallelDegree); i++; } } bool GeoProjectionConverter::set_state_plane_nad27_tm(const char* zone, char* description) { int i = 0; while (state_plane_tm_nad27_list[i].zone) { if (strcmp(zone, state_plane_tm_nad27_list[i].zone) == 0) { set_reference_ellipsoid(5); return set_transverse_mercator_projection(state_plane_tm_nad27_list[i].falseEastingMeter, state_plane_tm_nad27_list[i].falseNorthingMeter,state_plane_tm_nad27_list[i].latOriginDegree,state_plane_tm_nad27_list[i].longMeridianDegree,state_plane_tm_nad27_list[i].scaleFactor,description); } i++; } return false; } void GeoProjectionConverter::print_all_state_plane_nad27_tm() { int i = 0; while (state_plane_tm_nad27_list[i].zone) { fprintf(stderr, "%s - false east/north: %g/%g [m], origin lat/meridian long: %g/%g, scale factor: %g\n", state_plane_tm_nad27_list[i].zone, state_plane_tm_nad27_list[i].falseEastingMeter, state_plane_tm_nad27_list[i].falseNorthingMeter,state_plane_tm_nad27_list[i].latOriginDegree,state_plane_tm_nad27_list[i].longMeridianDegree,state_plane_tm_nad27_list[i].scaleFactor); i++; } } bool GeoProjectionConverter::set_state_plane_nad83_tm(const char* zone, char* description) { int i = 0; while (state_plane_tm_nad83_list[i].zone) { if (strcmp(zone, state_plane_tm_nad83_list[i].zone) == 0) { set_reference_ellipsoid(11); return set_transverse_mercator_projection(state_plane_tm_nad83_list[i].falseEastingMeter, state_plane_tm_nad83_list[i].falseNorthingMeter,state_plane_tm_nad83_list[i].latOriginDegree,state_plane_tm_nad83_list[i].longMeridianDegree,state_plane_tm_nad83_list[i].scaleFactor,description); } i++; } return false; } void GeoProjectionConverter::print_all_state_plane_nad83_tm() { int i = 0; while (state_plane_tm_nad83_list[i].zone) { fprintf(stderr, "%s - false east/north: %g/%g [m], origin lat/meridian long: %g/%g, scale factor: %g\n", state_plane_tm_nad83_list[i].zone, state_plane_tm_nad83_list[i].falseEastingMeter, state_plane_tm_nad83_list[i].falseNorthingMeter,state_plane_tm_nad83_list[i].latOriginDegree,state_plane_tm_nad83_list[i].longMeridianDegree,state_plane_tm_nad83_list[i].scaleFactor); i++; } } bool GeoProjectionConverter::has_projection() { return (projection_name[0] != '\0'); } const char* GeoProjectionConverter::get_projection_name() { return projection_name; } void GeoProjectionConverter::compute_lcc_parameters() { double es_sin; double t0,t1,t2; double m1,m2; es_sin = eccentricity * sin(lcc_lat_origin_radian); t0 = tan(PI_OVER_4 - lcc_lat_origin_radian / 2) / pow((1.0 - es_sin) / (1.0 + es_sin), (eccentricity / 2.0)); es_sin = eccentricity * sin(lcc_first_std_parallel_radian); t1 = tan(PI_OVER_4 - lcc_first_std_parallel_radian / 2) / pow((1.0 - es_sin) / (1.0 + es_sin), (eccentricity / 2.0)); m1 = cos(lcc_first_std_parallel_radian) / sqrt(1.0 - es_sin * es_sin); // compute lcc_n if (fabs(lcc_first_std_parallel_radian - lcc_second_std_parallel_radian) > 1.0e-10) { es_sin = eccentricity * sin(lcc_second_std_parallel_radian); t2 = tan(PI_OVER_4 - lcc_second_std_parallel_radian / 2) / pow((1.0 - es_sin) / (1.0 + es_sin), (eccentricity / 2.0)); m2 = cos(lcc_second_std_parallel_radian) / sqrt(1.0 - es_sin * es_sin); lcc_n = log(m1 / m2) / log(t1 / t2); } else { lcc_n = sin(lcc_first_std_parallel_radian); } // compute lcc_aF lcc_aF = equatorial_radius * m1 / (lcc_n * pow(t1, lcc_n)); // compute lcc_rho0; if ((t0 == 0) && (lcc_n < 0)) lcc_rho0 = 0.0; else lcc_rho0 = lcc_aF * pow(t0, lcc_n); } void GeoProjectionConverter::compute_tm_parameters() { /*True meridianal constants */ double tn = (equatorial_radius - polar_radius) / (equatorial_radius + polar_radius); double tn_other = eccentricity_e1; double tn2 = tn * tn; double tn3 = tn2 * tn; double tn4 = tn3 * tn; double tn5 = tn4 * tn; tm_ap = equatorial_radius * (1.e0 - tn + 5.e0 * (tn2 - tn3)/4.e0 + 81.e0 * (tn4 - tn5)/64.e0 ); tm_bp = 3.e0 * equatorial_radius * (tn - tn2 + 7.e0 * (tn3 - tn4) /8.e0 + 55.e0 * tn5/64.e0 )/2.e0; tm_cp = 15.e0 * equatorial_radius * (tn2 - tn3 + 3.e0 * (tn4 - tn5 )/4.e0) /16.0; tm_dp = 35.e0 * equatorial_radius * (tn3 - tn4 + 11.e0 * tn5 / 16.e0) / 48.e0; tm_ep = 315.e0 * equatorial_radius * (tn4 - tn5) / 512.e0; } // converts lat/long to UTM coords. Equations from USGS Bulletin 1532 // East Longitudes are positive, West longitudes are negative. // North latitudes are positive, South latitudes are negative // LatDegree and LongDegree are in decimal degrees // adapted from code written by Chuck Gantz- chuck.gantz@globalstar.com bool GeoProjectionConverter::LLtoUTM(const double LatDegree, const double LongDegree, double &UTMEastingMeter, double &UTMNorthingMeter, char* UTMZone) { const double k0 = 0.9996; // Make sure the longitude is between -180.00 .. 179.9 double LongTemp = (LongDegree+180)-int((LongDegree+180)/360)*360-180; // -180.00 .. 179.9; double LatRad = LatDegree*deg2rad; double LongRad = LongTemp*deg2rad; int ZoneNumber = int((LongTemp + 180)/6) + 1; if( LatDegree >= 56.0 && LatDegree < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0 ) ZoneNumber = 32; // Special zones for Svalbard if( LatDegree >= 72.0 && LatDegree < 84.0 ) { if( LongTemp >= 0.0 && LongTemp < 9.0 ) ZoneNumber = 31; else if( LongTemp >= 9.0 && LongTemp < 21.0 ) ZoneNumber = 33; else if( LongTemp >= 21.0 && LongTemp < 33.0 ) ZoneNumber = 35; else if( LongTemp >= 33.0 && LongTemp < 42.0 ) ZoneNumber = 37; } double LongOriginRad = ((ZoneNumber - 1)*6 - 180 + 3) * deg2rad; // + 3 puts origin in middle of zone // compute the UTM Zone from the latitude and longitude sprintf(UTMZone, "%d%c", ZoneNumber, UTMLetterDesignator(LatDegree)); double N = equatorial_radius/sqrt(1-eccentricity_squared*sin(LatRad)*sin(LatRad)); double T = tan(LatRad)*tan(LatRad); double C = eccentricity_prime_squared*cos(LatRad)*cos(LatRad); double A = cos(LatRad)*(LongRad-LongOriginRad); double M = equatorial_radius*((1 - eccentricity_squared/4 - 3*eccentricity_squared*eccentricity_squared/64 - 5*eccentricity_squared*eccentricity_squared*eccentricity_squared/256)*LatRad - (3*eccentricity_squared/8 + 3*eccentricity_squared*eccentricity_squared/32 + 45*eccentricity_squared*eccentricity_squared*eccentricity_squared/1024)*sin(2*LatRad) + (15*eccentricity_squared*eccentricity_squared/256 + 45*eccentricity_squared*eccentricity_squared*eccentricity_squared/1024)*sin(4*LatRad) - (35*eccentricity_squared*eccentricity_squared*eccentricity_squared/3072)*sin(6*LatRad)); UTMEastingMeter = (double)(k0*N*(A+(1-T+C)*A*A*A/6 + (5-18*T+T*T+72*C-58*eccentricity_prime_squared)*A*A*A*A*A/120) + 500000.0); UTMNorthingMeter = (double)(k0*(M+N*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24 + (61-58*T+T*T+600*C-330*eccentricity_prime_squared)*A*A*A*A*A*A/720))); if(LatDegree < 0) { UTMNorthingMeter += 10000000.0; //10000000 meter offset for southern hemisphere } return true; } // This routine determines the correct UTM letter designator for the given latitude // returns 'Z' if latitude is outside the UTM limits of 84N to 80S // Written by Chuck Gantz- chuck.gantz@globalstar.com char GeoProjectionConverter::UTMLetterDesignator(double LatDegree) { char LetterDesignator; if((84 >= LatDegree) && (LatDegree >= 72)) LetterDesignator = 'X'; else if((72 > LatDegree) && (LatDegree >= 64)) LetterDesignator = 'W'; else if((64 > LatDegree) && (LatDegree >= 56)) LetterDesignator = 'V'; else if((56 > LatDegree) && (LatDegree >= 48)) LetterDesignator = 'U'; else if((48 > LatDegree) && (LatDegree >= 40)) LetterDesignator = 'T'; else if((40 > LatDegree) && (LatDegree >= 32)) LetterDesignator = 'S'; else if((32 > LatDegree) && (LatDegree >= 24)) LetterDesignator = 'R'; else if((24 > LatDegree) && (LatDegree >= 16)) LetterDesignator = 'Q'; else if((16 > LatDegree) && (LatDegree >= 8)) LetterDesignator = 'P'; else if(( 8 > LatDegree) && (LatDegree >= 0)) LetterDesignator = 'N'; else if(( 0 > LatDegree) && (LatDegree >= -8)) LetterDesignator = 'M'; else if((-8> LatDegree) && (LatDegree >= -16)) LetterDesignator = 'L'; else if((-16 > LatDegree) && (LatDegree >= -24)) LetterDesignator = 'K'; else if((-24 > LatDegree) && (LatDegree >= -32)) LetterDesignator = 'J'; else if((-32 > LatDegree) && (LatDegree >= -40)) LetterDesignator = 'H'; else if((-40 > LatDegree) && (LatDegree >= -48)) LetterDesignator = 'G'; else if((-48 > LatDegree) && (LatDegree >= -56)) LetterDesignator = 'F'; else if((-56 > LatDegree) && (LatDegree >= -64)) LetterDesignator = 'E'; else if((-64 > LatDegree) && (LatDegree >= -72)) LetterDesignator = 'D'; else if((-72 > LatDegree) && (LatDegree >= -80)) LetterDesignator = 'C'; else LetterDesignator = 'Z'; //This is here as an error flag to show that the Latitude is outside the UTM limits return LetterDesignator; } // converts UTM coords to lat/long. Equations from USGS Bulletin 1532 // East Longitudes are positive, West longitudes are negative. // North latitudes are positive, South latitudes are negative // Lat and LongDegree are in decimal degrees. // adapted from code written by Chuck Gantz- chuck.gantz@globalstar.com bool GeoProjectionConverter::UTMtoLL(const double UTMEastingMeter, const double UTMNorthingMeter, double& LatDegree, double& LongDegree) { const double k0 = 0.9996; double x = UTMEastingMeter - 500000.0; // remove 500,000 meter offset for longitude double y = UTMNorthingMeter; if (!utm_northern_hemisphere) { y -= 10000000.0; //remove 10,000,000 meter offset used for southern hemisphere } double M = y / k0; double mu = M/(equatorial_radius*(1-eccentricity_squared/4-3*eccentricity_squared*eccentricity_squared/64-5*eccentricity_squared*eccentricity_squared*eccentricity_squared/256)); double phi1Rad = mu + (3*eccentricity_e1/2-27*eccentricity_e1*eccentricity_e1*eccentricity_e1/32)*sin(2*mu) + (21*eccentricity_e1*eccentricity_e1/16-55*eccentricity_e1*eccentricity_e1*eccentricity_e1*eccentricity_e1/32)*sin(4*mu) + (151*eccentricity_e1*eccentricity_e1*eccentricity_e1/96)*sin(6*mu); double phi1 = phi1Rad*rad2deg; double N1 = equatorial_radius/sqrt(1-eccentricity_squared*sin(phi1Rad)*sin(phi1Rad)); double T1 = tan(phi1Rad)*tan(phi1Rad); double C1 = eccentricity_prime_squared*cos(phi1Rad)*cos(phi1Rad); double R1 = equatorial_radius*(1-eccentricity_squared)/pow(1-eccentricity_squared*sin(phi1Rad)*sin(phi1Rad), 1.5); double D = x/(N1*k0); LatDegree = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccentricity_prime_squared)*D*D*D*D/24 + (61+90*T1+298*C1+45*T1*T1-252*eccentricity_prime_squared-3*C1*C1)*D*D*D*D*D*D/720); LatDegree = LatDegree * rad2deg; LongDegree = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccentricity_prime_squared+24*T1*T1)*D*D*D*D*D/120)/cos(phi1Rad); LongDegree = LongDegree * rad2deg + utm_long_origin; return true; } /* An alternate way to convert Lambert Conic Conformal Northing/Easting coordinates into Latitude & Longitude coordinates. The code adapted from Brenor Brophy (brenor dot brophy at gmail dot com) Homepage: www.brenorbrophy.com */ void lcc2ll( double e2, // Square of eccentricity double a, // Equatorial Radius double firstStdParallel, double secondStdParallel, double latOfOrigin, double longOfOrigin, double lccEasting, // Lambert coordinates of the point double lccNorthing, double falseEasting, double falseNorthing, double &LatDegree, // Latitude & Longitude of the point double &LongDegree) { double e = sqrt(e2); double phi1 = deg2rad*firstStdParallel; // Latitude of 1st std parallel double phi2 = deg2rad*secondStdParallel; // Latitude of 2nd std parallel double phio = deg2rad*latOfOrigin; // Latitude of Origin double lamdao = deg2rad*longOfOrigin; // Longitude of Origin double E = lccEasting; double N = lccNorthing; double Ef = falseEasting; double Nf = falseNorthing; double m1 = cos(phi1) / sqrt((1 - e2*sin(phi1)*sin(phi1))); double m2 = cos(phi2) / sqrt(( 1 - e2*sin(phi2)*sin(phi2))); double t1 = tan(PI_OVER_4-(phi1/2)) / pow(( ( 1 - e*sin(phi1) ) / ( 1 + e*sin(phi1) )),e/2); double t2 = tan(PI_OVER_4-(phi2/2)) / pow(( ( 1 - e*sin(phi2) ) / ( 1 + e*sin(phi2) )),e/2); double to = tan(PI_OVER_4-(phio/2)) / pow(( ( 1 - e*sin(phio) ) / ( 1 + e*sin(phio) )),e/2); double n = (log(m1)-log(m2)) / (log(t1)-log(t2)); double F = m1/(n*pow(t1,n)); double rf = a*F*pow(to,n); double r_ = sqrt( pow((E-Ef),2) + pow((rf-(N-Nf)),2) ); double t_ = pow(r_/(a*F),(1/n)); double theta_ = atan((E-Ef)/(rf-(N-Nf))); double lamda = theta_/n + lamdao; double phi0 = PI_OVER_2 - 2*atan(t_); phi1 = PI_OVER_2 - 2*atan(t_*pow(((1-e*sin(phi0))/(1+e*sin(phi0))),e/2)); phi2 = PI_OVER_2 - 2*atan(t_*pow(((1-e*sin(phi1))/(1+e*sin(phi1))),e/2)); double phi = PI_OVER_2 - 2*atan(t_*pow(((1-e*sin(phi2))/(1+e*sin(phi2))),e/2)); LatDegree = rad2deg*phi; LongDegree = rad2deg*lamda; } /* An alternate way to convert a Latitude/Longitude coordinate to an Northing/ Easting coordinate on a Lambert Conic Projection. The Northing/Easting parameters are calculated are in meters (because the datum used is in meters) and are relative to the falseNorthing/falseEasting coordinate. Which in turn is relative to the Lat/Long of origin. The formula were obtained from URL: http://www.ihsenergy.com/epsg/guid7_2.html. The code adapted from Brenor Brophy (brenor dot brophy at gmail dot com) Homepage: www.brenorbrophy.com */ void ll2lcc( double e2, // Square of eccentricity double a, // Equatorial Radius double firstStdParallel, double secondStdParallel, double latOfOrigin, double longOfOrigin, double &lccEasting, // Lambert coordinates of the point double &lccNorthing, double falseEasting, double falseNorthing, double LatDegree, // Latitude & Longitude of the point double LongDegree) { double e = sqrt(e2); double phi = deg2rad*LatDegree; // Latitude to convert double phi1 = deg2rad*firstStdParallel; // Latitude of 1st std parallel double phi2 = deg2rad*secondStdParallel; // Latitude of 2nd std parallel double lamda = deg2rad*LongDegree; // Lonitude to convert double phio = deg2rad*latOfOrigin; // Latitude of Origin double lamdao = deg2rad*longOfOrigin; // Longitude of Origin double m1 = cos(phi1) / sqrt((1 - e2*sin(phi1)*sin(phi1))); double m2 = cos(phi2) / sqrt((1 - e2*sin(phi2)*sin(phi2))); double t1 = tan(PI_OVER_4-(phi1/2)) / pow(( ( 1 - e*sin(phi1) ) / ( 1 + e*sin(phi1) )),e/2); double t2 = tan(PI_OVER_4-(phi2/2)) / pow(( ( 1 - e*sin(phi2) ) / ( 1 + e*sin(phi2) )),e/2); double to = tan(PI_OVER_4-(phio/2)) / pow(( ( 1 - e*sin(phio) ) / ( 1 + e*sin(phio) )),e/2); double t = tan(PI_OVER_4-(phi /2)) / pow(( ( 1 - e*sin(phi) ) / ( 1 + e*sin(phi) )),e/2); double n = (log(m1)-log(m2)) / (log(t1)-log(t2)); double F = m1/(n*pow(t1,n)); double rf = a*F*pow(to,n); double r = a*F*pow(t,n); double theta = n*(lamda - lamdao); lccEasting = falseEasting + r*sin(theta); lccNorthing = falseNorthing + rf - r*cos(theta); } /* * The function LCCtoLL() converts Lambert Conformal Conic projection * (easting and northing) coordinates to Geodetic (latitude and longitude) * coordinates, according to the current ellipsoid and Lambert Conformal * Conic projection parameters. * * LCCEastingMeter : input Easting/X in meters * LLCNorthingMeter : input Northing/Y in meters * LatDegree : output Latitude in decimal degrees * LongDegree : output Longitude in decimal degrees * * adapted from code by Garrett Potts ((C) 2000 ImageLinks Inc.) */ bool GeoProjectionConverter::LCCtoLL(const double LCCEastingMeter, const double LCCNorthingMeter, double& LatDegree, double& LongDegree) const { /* >>> alternate way to compute (but seems less precise) <<< lcc2ll(eccentricity_squared, equatorial_radius, lcc_first_std_parallel_degree, lcc_second_std_parallel_degree, lcc_lat_origin_degree, lcc_long_meridian_degree, LCCEastingMeter, LCCNorthingMeter, lcc_false_easting_meter, lcc_false_northing_meter, LatDegree, LongDegree); return true; */ double dx = LCCEastingMeter - lcc_false_easting_meter; double dy = LCCNorthingMeter - lcc_false_northing_meter; double rho0_MINUS_dy = lcc_rho0 - dy; double rho = sqrt(dx * dx + (rho0_MINUS_dy) * (rho0_MINUS_dy)); if (lcc_n < 0.0) { rho *= -1.0; dy *= -1.0; dx *= -1.0; rho0_MINUS_dy *= -1.0; } if (rho != 0.0) { double theta = atan2(dx, rho0_MINUS_dy); double t = pow(rho / (lcc_aF), 1.0 / lcc_n); double PHI = PI_OVER_2 - 2.0 * atan(t); double tempPHI = 0.0; while (fabs(PHI - tempPHI) > 4.85e-10) { double es_sin= eccentricity * sin(PHI); tempPHI = PHI; PHI = PI_OVER_2 - 2.0 * atan(t * pow((1.0 - es_sin) / (1.0 + es_sin), eccentricity / 2.0)); } LatDegree = PHI; LongDegree = theta / lcc_n + lcc_long_meridian_radian; if (fabs(LatDegree) < 2.0e-7) /* force tiny lat to 0 */ LatDegree = 0.0; else if (LatDegree > PI_OVER_2) /* force distorted lat to 90, -90 degrees */ LatDegree = 90.0; else if (LatDegree < -PI_OVER_2) LatDegree = -90.0; else LatDegree = rad2deg*LatDegree; if (fabs(LongDegree) < 2.0e-7) /* force tiny long to 0 */ LongDegree = 0.0; else if (LongDegree > PI) /* force distorted long to 180, -180 degrees */ LongDegree = 180.0; else if (LongDegree < -PI) LongDegree = -180.0; else LongDegree = rad2deg*LongDegree; } else { if (lcc_n > 0.0) LatDegree = 90.0; else LatDegree = -90.0; LongDegree = lcc_long_meridian_degree; } return true; } /* * The function LLtoLCC() converts Geodetic (latitude and longitude) * coordinates to Lambert Conformal Conic projection (easting and * northing) coordinates, according to the current ellipsoid and * Lambert Conformal Conic projection parameters. * * LatDegree : input Latitude in decimal degrees * LongDegree : input Longitude in decimal degrees * LCCEastingMeter : output Easting/X in meters * LLCNorthingMeter : output Northing/Y in meters * * adapted from code by Garrett Potts ((C) 2000 ImageLinks Inc.) */ bool GeoProjectionConverter::LLtoLCC(const double LatDegree, const double LongDegree, double& LCCEastingMeter, double& LCCNorthingMeter) const { /* >>> alternate way to compute (but seems less precise) <<< ll2lcc(eccentricity_squared, equatorial_radius, lcc_first_std_parallel_degree, lcc_second_std_parallel_degree, lcc_lat_origin_degree, lcc_long_meridian_degree, LCCEastingMeter, LCCNorthingMeter, lcc_false_easting_meter, lcc_false_northing_meter, LatDegreeAlt, LongDegreeAlt); return true; */ double rho = 0.0; double Latitude = LatDegree*deg2rad; double Longitude = LongDegree*deg2rad; if (fabs(fabs(Latitude) - PI_OVER_2) > 1.0e-10) { double slat = sin(Latitude); double es_sin = eccentricity*slat; double t = tan(PI_OVER_4 - Latitude / 2) / pow((1.0 - es_sin) / (1.0 + es_sin), eccentricity/2); rho = lcc_aF * pow(t, lcc_n); } else { if ((Latitude * lcc_n) <= 0) { // Point can not be projected return false; } } double dlam = Longitude - lcc_long_meridian_radian; double theta = lcc_n * dlam; LCCEastingMeter = rho * sin(theta) + lcc_false_easting_meter; LCCNorthingMeter = lcc_rho0 - rho * cos(theta) + lcc_false_northing_meter; return true; } #define SPHSN(Latitude) ((double) (equatorial_radius / sqrt( 1.e0 - eccentricity_squared * pow(sin(Latitude), 2)))) #define SPHTMD(Latitude) ((double) (tm_ap * Latitude \ - tm_bp * sin(2.e0 * Latitude) + tm_cp * sin(4.e0 * Latitude) \ - tm_dp * sin(6.e0 * Latitude) + tm_ep * sin(8.e0 * Latitude) ) ) #define DENOM(Latitude) ((double) (sqrt(1.e0 - eccentricity_squared * pow(sin(Latitude),2)))) #define SPHSR(Latitude) ((double) (equatorial_radius * (1.e0 - eccentricity_squared) / pow(DENOM(Latitude), 3))) /* * The function LLtoTM() converts geodetic (latitude and longitude) * coordinates to Transverse Mercator projection (easting and northing) * coordinates, according to the current ellipsoid and Transverse Mercator * projection parameters. * * LatDegree : input Latitude in decimal degrees * LongDegree : input Longitude in decimal degrees * TMEastingMeter : output Easting/X in meters * TMNorthingMeter : output Northing/Y in meters * * adapted from code by Garrett Potts ((C) 2000 ImageLinks Inc.) */ bool GeoProjectionConverter::LLtoTM(const double LatDegree, const double LongDegree, double& TMEastingMeter, double& TMNorthingMeter) const { double Latitude = LatDegree*deg2rad; double Longitude = LongDegree*deg2rad; double c; /* Cosine of latitude */ double c2; double c3; double c5; double c7; double dlam; /* Delta longitude - Difference in Longitude */ double eta; /* constant - eccentricity_prime_squared *c *c */ double eta2; double eta3; double eta4; double s; /* Sine of latitude */ double sn; /* Radius of curvature in the prime vertical */ double t; /* Tangent of latitude */ double tan2; double tan3; double tan4; double tan5; double tan6; double t1; /* Term in coordinate conversion formula - GP to Y */ double t2; /* Term in coordinate conversion formula - GP to Y */ double t3; /* Term in coordinate conversion formula - GP to Y */ double t4; /* Term in coordinate conversion formula - GP to Y */ double t5; /* Term in coordinate conversion formula - GP to Y */ double t6; /* Term in coordinate conversion formula - GP to Y */ double t7; /* Term in coordinate conversion formula - GP to Y */ double t8; /* Term in coordinate conversion formula - GP to Y */ double t9; /* Term in coordinate conversion formula - GP to Y */ double tmd; /* True Meridional distance */ double tmdo; /* True Meridional distance for latitude of origin */ if (Longitude > PI) Longitude -= TWO_PI; dlam = Longitude - tm_long_meridian_radian; if (dlam > PI) dlam -= TWO_PI; if (dlam < -PI) dlam += TWO_PI; if (fabs(dlam) < 2.e-10) dlam = 0.0; s = sin(Latitude); c = cos(Latitude); c2 = c * c; c3 = c2 * c; c5 = c3 * c2; c7 = c5 * c2; t = tan (Latitude); tan2 = t * t; tan3 = tan2 * t; tan4 = tan3 * t; tan5 = tan4 * t; tan6 = tan5 * t; eta = eccentricity_prime_squared * c2; eta2 = eta * eta; eta3 = eta2 * eta; eta4 = eta3 * eta; /* radius of curvature in prime vertical */ sn = SPHSN(Latitude); /* True Meridianal Distances */ tmd = SPHTMD(Latitude); /* Origin */ tmdo = SPHTMD (tm_lat_origin_radian); /* northing */ t1 = (tmd - tmdo) * tm_scale_factor; t2 = sn * s * c * tm_scale_factor/ 2.e0; t3 = sn * s * c3 * tm_scale_factor * (5.e0 - tan2 + 9.e0 * eta + 4.e0 * eta2) /24.e0; t4 = sn * s * c5 * tm_scale_factor * (61.e0 - 58.e0 * tan2 + tan4 + 270.e0 * eta - 330.e0 * tan2 * eta + 445.e0 * eta2 + 324.e0 * eta3 -680.e0 * tan2 * eta2 + 88.e0 * eta4 -600.e0 * tan2 * eta3 - 192.e0 * tan2 * eta4) / 720.e0; t5 = sn * s * c7 * tm_scale_factor * (1385.e0 - 3111.e0 * tan2 + 543.e0 * tan4 - tan6) / 40320.e0; TMNorthingMeter = tm_false_northing_meter + t1 + pow(dlam,2.e0) * t2 + pow(dlam,4.e0) * t3 + pow(dlam,6.e0) * t4 + pow(dlam,8.e0) * t5; /* Easting */ t6 = sn * c * tm_scale_factor; t7 = sn * c3 * tm_scale_factor * (1.e0 - tan2 + eta ) /6.e0; t8 = sn * c5 * tm_scale_factor * (5.e0 - 18.e0 * tan2 + tan4 + 14.e0 * eta - 58.e0 * tan2 * eta + 13.e0 * eta2 + 4.e0 * eta3 - 64.e0 * tan2 * eta2 - 24.e0 * tan2 * eta3 )/ 120.e0; t9 = sn * c7 * tm_scale_factor * ( 61.e0 - 479.e0 * tan2 + 179.e0 * tan4 - tan6 ) /5040.e0; TMEastingMeter = tm_false_easting_meter + dlam * t6 + pow(dlam,3.e0) * t7 + pow(dlam,5.e0) * t8 + pow(dlam,7.e0) * t9; return true; } /* * The function TMtoLL() converts Transverse Mercator projection (easting and * northing) coordinates to geodetic (latitude and longitude) coordinates, * according to the current ellipsoid and Transverse Mercator projection * parameters. * * TMEastingMeter : input Easting/X in meters * TMNorthingMeter : input Northing/Y in meters * LatDegree : output Latitude in decimal degrees * LongDegree : output Longitude in decimal degrees * * adapted from code by Garrett Potts ((C) 2000 ImageLinks Inc.) */ bool GeoProjectionConverter::TMtoLL(const double TMEastingMeter, const double TMNorthingMeter, double& LatDegree, double& LongDegree) const { double c; /* Cosine of latitude */ double de; /* Delta easting - Difference in Easting (Easting-Fe) */ double dlam; /* Delta longitude - Difference in Longitude */ double eta; /* constant - eccentricity_prime_squared *c *c */ double eta2; double eta3; double eta4; double ftphi; /* Footpoint latitude */ int i; /* Loop iterator */ double s; /* Sine of latitude */ double sn; /* Radius of curvature in the prime vertical */ double sr; /* Radius of curvature in the meridian */ double t; /* Tangent of latitude */ double tan2; double tan4; double t10; /* Term in coordinate conversion formula - GP to Y */ double t11; /* Term in coordinate conversion formula - GP to Y */ double t12; /* Term in coordinate conversion formula - GP to Y */ double t13; /* Term in coordinate conversion formula - GP to Y */ double t14; /* Term in coordinate conversion formula - GP to Y */ double t15; /* Term in coordinate conversion formula - GP to Y */ double t16; /* Term in coordinate conversion formula - GP to Y */ double t17; /* Term in coordinate conversion formula - GP to Y */ double tmd; /* True Meridional distance */ double tmdo; /* True Meridional distance for latitude of origin */ /* True Meridional Distances for latitude of origin */ tmdo = SPHTMD(tm_lat_origin_radian); /* Origin */ tmd = tmdo + (TMNorthingMeter - tm_false_northing_meter) / tm_scale_factor; /* First Estimate */ sr = SPHSR(0.e0); ftphi = tmd/sr; for (i = 0; i < 5 ; i++) { t10 = SPHTMD (ftphi); sr = SPHSR(ftphi); ftphi = ftphi + (tmd - t10) / sr; } /* Radius of Curvature in the meridian */ sr = SPHSR(ftphi); /* Radius of Curvature in the meridian */ sn = SPHSN(ftphi); /* Sine Cosine terms */ s = sin(ftphi); c = cos(ftphi); /* Tangent Value */ t = tan(ftphi); tan2 = t * t; tan4 = tan2 * tan2; eta = eccentricity_prime_squared * pow(c,2); eta2 = eta * eta; eta3 = eta2 * eta; eta4 = eta3 * eta; de = TMEastingMeter - tm_false_easting_meter; if (fabs(de) < 0.0001) de = 0.0; /* Latitude */ double Latitude; t10 = t / (2.e0 * sr * sn * pow(tm_scale_factor, 2)); t11 = t * (5.e0 + 3.e0 * tan2 + eta - 4.e0 * pow(eta,2) - 9.e0 * tan2 * eta) / (24.e0 * sr * pow(sn,3) * pow(tm_scale_factor,4)); t12 = t * (61.e0 + 90.e0 * tan2 + 46.e0 * eta + 45.E0 * tan4 - 252.e0 * tan2 * eta - 3.e0 * eta2 + 100.e0 * eta3 - 66.e0 * tan2 * eta2 - 90.e0 * tan4 * eta + 88.e0 * eta4 + 225.e0 * tan4 * eta2 + 84.e0 * tan2* eta3 - 192.e0 * tan2 * eta4) / ( 720.e0 * sr * pow(sn,5) * pow(tm_scale_factor, 6) ); t13 = t * ( 1385.e0 + 3633.e0 * tan2 + 4095.e0 * tan4 + 1575.e0 * pow(t,6))/ (40320.e0 * sr * pow(sn,7) * pow(tm_scale_factor,8)); Latitude = ftphi - pow(de,2) * t10 + pow(de,4) * t11 - pow(de,6) * t12 + pow(de,8) * t13; t14 = 1.e0 / (sn * c * tm_scale_factor); t15 = (1.e0 + 2.e0 * tan2 + eta) / (6.e0 * pow(sn,3) * c * pow(tm_scale_factor,3)); t16 = (5.e0 + 6.e0 * eta + 28.e0 * tan2 - 3.e0 * eta2 + 8.e0 * tan2 * eta + 24.e0 * tan4 - 4.e0 * eta3 + 4.e0 * tan2 * eta2 + 24.e0 * tan2 * eta3) / (120.e0 * pow(sn,5) * c * pow(tm_scale_factor,5)); t17 = (61.e0 + 662.e0 * tan2 + 1320.e0 * tan4 + 720.e0 * pow(t,6)) / (5040.e0 * pow(sn,7) * c * pow(tm_scale_factor,7)); /* Difference in Longitude */ dlam = de * t14 - pow(de,3) * t15 + pow(de,5) * t16 - pow(de,7) * t17; /* Longitude */ double Longitude = tm_long_meridian_radian + dlam; while (Latitude > PI_OVER_2) { Latitude = PI - Latitude; Longitude += PI; if (Longitude > PI) Longitude -= TWO_PI; } while (Latitude < -PI_OVER_2) { Latitude = - (Latitude + PI); Longitude += PI; if (Longitude > PI) Longitude -= TWO_PI; } if (Longitude > TWO_PI) Longitude -= TWO_PI; if (Longitude < -PI) Longitude += TWO_PI; LatDegree = rad2deg*Latitude; LongDegree = rad2deg*Longitude; return true; } GeoProjectionConverter::GeoProjectionConverter() { ellipsoid_id = 0; ellipsoid_name = 0; equatorial_radius = 0; eccentricity_squared = 0; eccentricity = 0; projection_name[0] = '\0'; utm_zone_number = 0; utm_zone_letter = 0; utm_northern_hemisphere = false; lcc_false_easting_meter = 0; lcc_false_northing_meter = 0; lcc_lat_origin_degree = 0; lcc_long_meridian_degree = 0; lcc_first_std_parallel_degree = 0; lcc_second_std_parallel_degree = 0; lcc_lat_origin_radian = 0; lcc_long_meridian_radian = 0; lcc_first_std_parallel_radian = 0; lcc_second_std_parallel_radian = 0; lcc_n = 0; lcc_aF = 0; lcc_rho0 = 0; // parameters for convenience coordinates2meter = 1.0; elevation2meter = 1.0; } GeoProjectionConverter::~GeoProjectionConverter() { } void GeoProjectionConverter::set_coordinates_in_survey_feet() { coordinates2meter = 0.3048006096012; } void GeoProjectionConverter::set_coordinates_in_feet() { coordinates2meter = 0.3048; } void GeoProjectionConverter::set_coordinates_in_meter() { coordinates2meter = 1.0; } const char* GeoProjectionConverter::get_coordinate_unit_description_string(bool abrev) { return (coordinates2meter == 1.0 ? (abrev ? "m" : "meter") : (coordinates2meter == 0.3048 ? (abrev ? "ft" : "feet") : (abrev ? "sft" : "surveyfeet"))); } void GeoProjectionConverter::set_elevation_in_feet() { elevation2meter = 0.3048; } void GeoProjectionConverter::set_elevation_in_meter() { elevation2meter = 1.0; } const char* GeoProjectionConverter::get_elevation_unit_description_string(bool abrev) { return (elevation2meter == 1.0 ? (abrev ? "m" : "meter") : (abrev ? "ft" : "feet")); } void GeoProjectionConverter::to_kml_style_lat_long_elevation(const double* point, double& latDeg, double& longDeg, float& elevation) { if (projection_name[0] == 'U') UTMtoLL(point[0], point[1], latDeg, longDeg); else if (projection_name[0] == 'L') LCCtoLL(coordinates2meter*point[0], coordinates2meter*point[1], latDeg, longDeg); else TMtoLL(coordinates2meter*point[0], coordinates2meter*point[1], latDeg, longDeg); elevation = (float)(elevation2meter*point[2]); } void GeoProjectionConverter::to_kml_style_lat_long_elevation(const float* point, double& latDeg, double& longDeg, float& elevation) { if (projection_name[0] == 'U') UTMtoLL(point[0], point[1], latDeg, longDeg); else if (projection_name[0] == 'L') LCCtoLL(coordinates2meter*point[0], coordinates2meter*point[1], latDeg, longDeg); else TMtoLL(coordinates2meter*point[0], coordinates2meter*point[1], latDeg, longDeg); elevation = (float)(elevation2meter*point[2]); } void GeoProjectionConverter::to_kml_style_lat_long_elevation(const int* point, double& latDeg, double& longDeg, float& elevation) { if (projection_name[0] == 'U') UTMtoLL(point[0], point[1], latDeg, longDeg); else if (projection_name[0] == 'L') LCCtoLL(coordinates2meter*point[0], coordinates2meter*point[1], latDeg, longDeg); else TMtoLL(coordinates2meter*point[0], coordinates2meter*point[1], latDeg, longDeg); elevation = (float)(elevation2meter*point[2]); }