43 #ifdef CHECK_MEMORY_LEAKS
45 #endif // CHECK_MEMORY_LEAKS
64 myInverseProjection(0),
68 myGeoScale(pow(10, (double) - shift)),
69 myProjectionMethod(NONE),
70 myUseInverseProjection(inverse),
72 myConvBoundary(conv) {
75 }
else if (proj ==
"-") {
77 }
else if (proj ==
"UTM") {
79 }
else if (proj ==
"DHDN") {
81 }
else if (proj ==
"DHDN_UTM") {
86 myProjection = pj_init_plus(proj.c_str());
87 if (myProjection == 0) {
98 if (myProjection != 0) {
99 pj_free(myProjection);
101 if (myInverseProjection != 0) {
102 pj_free(myInverseProjection);
104 if (myGeoProjection != 0) {
105 pj_free(myInverseProjection);
121 if (myProjection != 0) {
122 pj_free(myProjection);
125 if (myInverseProjection != 0) {
126 pj_free(myInverseProjection);
127 myInverseProjection = 0;
129 if (myGeoProjection != 0) {
130 pj_free(myGeoProjection);
133 if (orig.myProjection != 0) {
136 if (orig.myInverseProjection != 0) {
137 myInverseProjection = pj_init_plus(pj_get_def(orig.myInverseProjection, 0));
139 if (orig.myGeoProjection != 0) {
140 myGeoProjection = pj_init_plus(pj_get_def(orig.myGeoProjection, 0));
149 std::string proj =
"!";
150 int shift = oc.
getInt(
"proj.scale");
152 bool inverse = oc.
exists(
"proj.inverse") && oc.
getBool(
"proj.inverse");
154 if (oc.
getBool(
"simple-projection")) {
160 WRITE_ERROR(
"Inverse projection works only with explicit proj parameters.");
164 if (numProjections > 1) {
165 WRITE_ERROR(
"The projection method needs to be uniquely defined.");
171 }
else if (oc.
getBool(
"proj.dhdn")) {
173 }
else if (oc.
getBool(
"proj.dhdnutm")) {
201 oc.
addSynonyme(
"simple-projection",
"proj.simple",
true);
202 oc.
addDescription(
"simple-projection",
"Projection",
"Uses a simple method for projection");
206 oc.
addDescription(
"proj.scale",
"Projection",
"Number of places to shift decimal point to right in geo-coordinates");
210 oc.
addDescription(
"proj.utm",
"Projection",
"Determine the UTM zone (for a universal transversal mercator projection based on the WGS84 ellipsoid)");
213 oc.
addDescription(
"proj.dhdn",
"Projection",
"Determine the DHDN zone (for a transversal mercator projection based on the bessel ellipsoid, \"Gauss-Krueger\")");
216 oc.
addDescription(
"proj",
"Projection",
"Uses STR as proj.4 definition for projection");
219 oc.
addDescription(
"proj.inverse",
"Projection",
"Inverses projection");
222 oc.
addDescription(
"proj.dhdnutm",
"Projection",
"Convert from Gauss-Krueger to UTM");
249 p = pj_inv(p, myProjection);
260 if (includeInBoundary) {
265 if (myProjection == 0) {
269 int zone = (
int)((x - 500000.) / 1000000.);
270 if (zone < 1 || zone > 5) {
275 " +k=1 +x_0=" +
toString(zone * 1000000 + 500000) +
276 " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
277 myInverseProjection = pj_init_plus(
myProjString.c_str());
278 myGeoProjection = pj_init_plus(
"+proj=latlong +datum=WGS84");
280 x = ((x - 500000.) / 1000000.) * 3;
283 int zone = (
int)(x + 180) / 6 + 1;
285 " +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
291 int zone = (
int)(x / 3);
292 if (zone < 1 || zone > 5) {
297 " +k=1 +x_0=" +
toString(zone * 1000000 + 500000) +
298 " +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";
307 if (myInverseProjection != 0) {
310 if (pj_transform(myInverseProjection, myGeoProjection, 1, 1, &x, &y, NULL)) {
319 if (includeInBoundary) {
336 if (x > 180.1 || x < -180.1 || y > 90.1 || y < -90.1) {
340 if (myProjection != 0) {
342 p.u = x * DEG_TO_RAD;
343 p.v = y * DEG_TO_RAD;
344 p = pj_fwd(p, myProjection);
352 x *= 111320. * cos(
DEG2RAD(ys));