1 #include "GeomCoords.h"
2
3 #include "method.h"
4
5 namespace Upp {
6
7 #define LLOG(x) // RLOG(x)
8
GisCoordsLonLat(double prime_meridian)9 GisCoordsLonLat::GisCoordsLonLat(double prime_meridian)
10 : prime_meridian(prime_meridian)
11 {
12 }
13
LonLat(Pointf xy) const14 Pointf GisCoordsLonLat::LonLat(Pointf xy) const
15 {
16 return Pointf(xy.x + prime_meridian, xy.y);
17 }
18
LonLatExtent(const Rectf & xy_extent) const19 Rectf GisCoordsLonLat::LonLatExtent(const Rectf& xy_extent) const
20 {
21 return xy_extent.OffsetedHorz(prime_meridian);
22 }
23
Project(Pointf lonlat,int branch) const24 Pointf GisCoordsLonLat::Project(Pointf lonlat, int branch) const
25 {
26 return Pointf(lonlat.x - prime_meridian, lonlat.y);
27 }
28
ProjectExtent(const Rectf & lonlat_extent) const29 Rectf GisCoordsLonLat::ProjectExtent(const Rectf& lonlat_extent) const
30 {
31 return lonlat_extent.OffsetedHorz(-prime_meridian);
32 }
33
ProjectDeviation(Pointf lonlat1,Pointf lonlat2,int branch) const34 double GisCoordsLonLat::ProjectDeviation(Pointf lonlat1, Pointf lonlat2, int branch) const
35 {
36 return 0;
37 }
38
ProjectRatio(Pointf lonlat,int branch) const39 double GisCoordsLonLat::ProjectRatio(Pointf lonlat, int branch) const
40 {
41 return 1;
42 }
43
44 /*
45 String GisCoordsLonLat::ToString() const
46 {
47 return NFormat("%s, PM=%s", GetName(), FormatAngle(prime_meridian, 3));
48 }
49 */
50
EnumArgs()51 Array<GisCoords::Arg> GisCoordsLonLat::EnumArgs()
52 {
53 Array<GisCoords::Arg> args;
54 args.Add(GisCoords::Arg::Angle(prime_meridian, "PM", "Prime meridian", Null));
55 return args;
56 }
57
EnumConstArgs() const58 Array<GisCoords::ConstArg> GisCoordsLonLat::EnumConstArgs() const
59 {
60 Array<GisCoords::ConstArg> args;
61 return args;
62 }
63
SyncArgs()64 void GisCoordsLonLat::SyncArgs()
65 {
66 }
67
IsProjected() const68 bool GisCoordsLonLat::IsProjected() const
69 {
70 return false;
71 }
72
GisCoordsSpherical(double bp)73 GisCoordsSpherical::GisCoordsSpherical(double bp)
74 {
75 base_parallel = bp;
76 SyncArgs();
77 }
78
LonLat(Pointf xy) const79 Pointf GisCoordsSpherical::LonLat(Pointf xy) const
80 {
81 CheckInterpolator();
82 return Pointf(minmax(xy.x / alpha, -M_PI, +M_PI),
83 gauss_latitude(minmax(xy.y, -M_PI / 2, +M_PI / 2)));
84 }
85
LonLatExtent(const Rectf & xy_extent) const86 Rectf GisCoordsSpherical::LonLatExtent(const Rectf& xy_extent) const
87 {
88 CheckInterpolator();
89 return Rectf(xy_extent.left / alpha, gauss_latitude(xy_extent.top),
90 xy_extent.right / alpha, gauss_latitude(xy_extent.bottom));
91 }
92
Project(Pointf lonlat,int branch) const93 Pointf GisCoordsSpherical::Project(Pointf lonlat, int branch) const
94 {
95 CheckInterpolator();
96 return Pointf(lonlat.x * alpha, gauss_projected(lonlat.y));
97 }
98
ProjectExtent(const Rectf & lonlat_extent) const99 Rectf GisCoordsSpherical::ProjectExtent(const Rectf& lonlat_extent) const
100 {
101 CheckInterpolator();
102 return Rectf(lonlat_extent.left * alpha, gauss_projected(lonlat_extent.top),
103 lonlat_extent.right * alpha, gauss_projected(lonlat_extent.bottom));
104 }
105
ProjectDeviation(Pointf lonlat1,Pointf lonlat2,int branch) const106 double GisCoordsSpherical::ProjectDeviation(Pointf lonlat1, Pointf lonlat2, int branch) const
107 {
108 CheckInterpolator();
109 return 0;
110 }
111
ProjectRatio(Pointf lonlat,int branch) const112 double GisCoordsSpherical::ProjectRatio(Pointf lonlat, int branch) const
113 {
114 double U = gauss_projected(lonlat.y);
115 double coef = alpha * R * ellipsoid.N(lonlat.y);
116 double cu = cos(U), cp = cos(lonlat.y);
117 return (cp ? coef * cu / cp : 1);
118 }
119
120 /*
121 String GisCoordsSpherical::ToString() const
122 {
123 return NFormat("%s, base parallel = %s", GetName(), FormatAngle(base_parallel, 3));
124 }
125 */
126
EnumArgs()127 Array<GisCoords::Arg> GisCoordsSpherical::EnumArgs()
128 {
129 Array<GisCoords::Arg> out;
130 out.Add(GisCoords::Arg::Angle(base_parallel, "BP", "Base parallel", Null));
131 return out;
132 }
133
EnumConstArgs() const134 Array<GisCoords::ConstArg> GisCoordsSpherical::EnumConstArgs() const
135 {
136 Array<GisCoords::ConstArg> args;
137 args.Add(GisCoords::ConstArg("A", "alpha", alpha));
138 args.Add(GisCoords::ConstArg("K", "k", k));
139 args.Add(GisCoords::ConstArg("R", "R", R));
140 args.Add(GisCoords::ConstArg("e", "e", e));
141 args.Add(GisCoords::ConstArg("U0", "U0", FormatDegree(U0, 3)));
142 return args;
143 }
144
SyncArgs()145 void GisCoordsSpherical::SyncArgs()
146 {
147 e = sqrt(ellipsoid.e2);
148 double phi0 = base_parallel * DEGRAD;
149 alpha = sqrt(1 + (ellipsoid.e2 * sqr(sqr(cos(phi0)))) / (1 - ellipsoid.e2));
150 double sinphi = sin(phi0);
151 U0 = asin(sinphi / alpha);
152 k = exp(alpha * (log(tan(phi0 / 2 + M_PI / 4)) + e / 2 * log((1 - e * sinphi) / (1 + e * sinphi))))
153 / tan(U0 / 2 + M_PI / 4);
154 // k = pow(tan(base_parallel / 2 + M_PI / 4), alpha)
155 // * pow((1 - e * sinphi) / (1 + e * sinphi), alpha * e / 2)
156 // / tan(U0 / 2 + M_PI / 4);
157 R = ellipsoid.a * sqrt(1 - ellipsoid.e2) / (1 - ellipsoid.e2 * sqr(sinphi));
158 gauss_projected.Clear();
159 gauss_latitude.Clear();
160 }
161
SyncInterpolator() const162 void GisCoordsSpherical::SyncInterpolator() const
163 {
164 SphericalLatitudeFunction gslf(alpha, k, R, e, U0);
165 gslf.Dump(-1.58, +1.58, 1000);
166 gslf.Dump(-1.58, -1.56, 1000);
167 gslf.Dump(+1.56, +1.58, 1000);
168 gauss_projected.Create(-M_PI / 2, +M_PI / 2, gslf, 300, 5000, 4);
169 gauss_latitude.CreateInverse(-M_PI / 2, +M_PI / 2, gslf, 300, 5000, 4);
170 }
171
GisCoordsConical(double gauss_base_parallel_,double pole_meridian_,double pole_parallel_,double central_parallel_,double north_parallel_,double south_parallel_,double multiplier_,double xmeteroffset_,double ymeteroffset_)172 GisCoordsConical::GisCoordsConical(
173 // bool gauss_sphere_,
174 double gauss_base_parallel_,
175 double pole_meridian_,
176 double pole_parallel_,
177 double central_parallel_,
178 // double prime_meridian_,
179 double north_parallel_,
180 double south_parallel_,
181 double multiplier_,
182 double xmeteroffset_,
183 double ymeteroffset_)
184 {
185 // gauss_sphere = gauss_sphere_;
186 gauss_base_parallel = gauss_base_parallel_;
187 pole_meridian = pole_meridian_;
188 pole_parallel = pole_parallel_;
189 central_parallel = central_parallel_;
190 // prime_meridian = prime_meridian_;
191 north_parallel = north_parallel_;
192 south_parallel = south_parallel_;
193 multiplier = multiplier_;
194 xmeteroffset = xmeteroffset_;
195 ymeteroffset = ymeteroffset_;
196 SyncArgs();
197 }
198
LonLat(Pointf xy) const199 Pointf GisCoordsConical::LonLat(Pointf xy) const
200 {
201 CheckInterpolator();
202 xy.x -= xmeteroffset;
203 xy.y -= ymeteroffset;
204 double eps = (modulo(Bearing(xy) / DEGRAD - 90, 360) - 180) / n;
205 double rho = Length(xy);
206 double lat = tgsin_latitude(rho);
207 // double prcr = 2 / DEGRAD * atan(pow(rho_coef / rho, 1 / n)) - 90;
208 Pointf global = orientation.Global(eps, lat);
209 return Pointf(global.x / gauss_alpha, gauss_latitude(global.y));
210 }
211
LonLatExtent(const Rectf & xy_extent) const212 Rectf GisCoordsConical::LonLatExtent(const Rectf& xy_extent) const
213 {
214 // CheckInterpolator();
215 // xy_extent.Offset(-xmeteroffset, -ymeteroffset);
216
217 return GisCoords::Data::LonLatExtent(xy_extent);
218 }
219
GetBranchCount() const220 int GisCoordsConical::GetBranchCount() const
221 {
222 return 3;
223 }
224
LonLatBoundary(const Rectf & lonlat_extent,int branch) const225 Array<Pointf> GisCoordsConical::LonLatBoundary(const Rectf& lonlat_extent, int branch) const
226 {
227 Array<Pointf> out;
228 int em = DegreeMask(lonlat_extent.left /* - prime_meridian*/, lonlat_extent.right /*- prime_meridian*/);
229 if(em & AM_E2)
230 {
231 if(branch == 1)
232 {
233 out.SetCount(4);
234 double x0 = (em & 1 ? 0 : modulo(lonlat_extent.left + 180, 360) - 180);
235 out[0] = Pointf(x0, lonlat_extent.top);
236 out[1] = Pointf(180, lonlat_extent.top);
237 out[2] = Pointf(180, lonlat_extent.bottom);
238 out[3] = Pointf(x0, lonlat_extent.bottom);
239 }
240 if(branch == 2)
241 {
242 out.SetCount(4);
243 double x0 = (em & 1 ? 0 : modulo(lonlat_extent.right - 180, 360) + 180);
244 out[0] = Pointf(x0, lonlat_extent.top);
245 out[1] = Pointf(-180, lonlat_extent.top);
246 out[2] = Pointf(-180, lonlat_extent.bottom);
247 out[3] = Pointf(x0, lonlat_extent.bottom);
248 }
249 }
250 else if(!branch)
251 out = Data::LonLatBoundary(lonlat_extent, 0);
252 return out;
253 }
254
GetBranchTree(const Rectf & lonlat_extent) const255 GisBSPTree GisCoordsConical::GetBranchTree(const Rectf& lonlat_extent) const
256 {
257 double pm = modulo(pole_parallel - lonlat_extent.left, 360) + lonlat_extent.left;
258 double l = lonlat_extent.left - pm, r = lonlat_extent.right - pm;
259 if(l >= -90 + 1 / 60.0 && r <= +90 - 1 / 60.0)
260 return GisBSPTree::Node(0);
261 GisBSPTree::Node minus = (l >= -90 + 1 / 60.0
262 ? GisBSPTree::Node(2)
263 : GisBSPTree::Node(new GisBSPTree::Split(Pointf(1, 0), pm - 180,
264 GisBSPTree::Node(1),
265 GisBSPTree::Node(2))));
266 GisBSPTree::Node plus = (r <= 180 - 1 / 60.0
267 ? GisBSPTree::Node(1)
268 : GisBSPTree::Node(new GisBSPTree::Split(Pointf(1, 0), pm + 180,
269 GisBSPTree::Node(1),
270 GisBSPTree::Node(2))));
271 return GisBSPTree::Node(new GisBSPTree::Split(Pointf(1, 0), pm, pick(minus), pick(plus)));
272 }
273
Project(Pointf lonlat,int branch) const274 Pointf GisCoordsConical::Project(Pointf lonlat, int branch) const
275 {
276 CheckInterpolator();
277 lonlat = orientation.Local(lonlat.x * gauss_alpha, gauss_projected(lonlat.y));
278 double comp = (branch == 0 ? 180 : branch == 1 ? 90 : 270);
279 double l = modulo(lonlat.x /*- prime_meridian*/ + comp, 360) - comp;
280 double eps = n * l;
281 double rho = tgsin_projected(lonlat.y);
282 Pointf out = PolarPointf(rho, (eps - 90) * DEGRAD);
283 out.x += xmeteroffset;
284 out.y += ymeteroffset;
285 return out;
286 }
287
ProjectExtent(const Rectf & lonlat_extent) const288 Rectf GisCoordsConical::ProjectExtent(const Rectf& lonlat_extent) const
289 {
290 // Rectf local = orientation.LocalExtent(lonlat_extent);
291 // return PolarToExtent(local.left * n, tgsin_projected(local.top),
292 // local.right * n, tgsin_projected(local.bottom));
293 return GisCoords::Data::ProjectExtent(lonlat_extent);
294 }
295
ProjectDeviation(Pointf lonlat1,Pointf lonlat2,int branch) const296 double GisCoordsConical::ProjectDeviation(Pointf lonlat1, Pointf lonlat2, int branch) const
297 {
298 return ExtentDeviation(SortRectf(lonlat1, lonlat2));
299 }
300
ExtentDeviation(const Rectf & lonlat_extent) const301 double GisCoordsConical::ExtentDeviation(const Rectf& lonlat_extent) const
302 {
303 CheckInterpolator();
304 Rectf local = orientation.LocalExtent(lonlat_extent);
305 double rho = tgsin_projected(local.bottom);
306 return rho * (1 + sin(local.Width() * (DEGRAD / 2) - M_PI / 2));
307 }
308
ProjectRatio(Pointf lonlat,int branch) const309 double GisCoordsConical::ProjectRatio(Pointf lonlat, int branch) const
310 {
311 Pointf local = orientation.Local(lonlat);
312 double den = ellipsoid.N(local.y) * cos(local.y);
313 return n * tgsin_projected(local.y) / (den ? den : 1);
314 }
315
MinMaxRatio(const Rectf & lonlat_extent) const316 Sizef GisCoordsConical::MinMaxRatio(const Rectf& lonlat_extent) const
317 {
318 CheckInterpolator();
319 Rectf local = orientation.LocalExtent(lonlat_extent);
320 double a = n * tgsin_projected(local.top) * DEGRAD;
321 double b = n * tgsin_projected(local.bottom) * DEGRAD;
322 double c = a;
323 if(local.top < central_parallel && local.bottom > central_parallel)
324 c = n * tgsin_projected(central_parallel) * DEGRAD;
325 return Sizef(min(a, min(b, c)), max(a, max(b, c)));
326 }
327
328 /*
329 String GisCoordsConical::ToString() const
330 {
331 return NFormat("%s, central parallel = %s, "
332 "north parallel = %s, south parallel = %s, "
333 "x meter offset = %3n, y meter offset = %3n",
334 GetName(),
335 FormatAngle(central_parallel, 3), FormatAngle(prime_meridian, 3),
336 FormatAngle(north_parallel, 3), FormatAngle(south_parallel, 3),
337 xmeteroffset, ymeteroffset);
338 }
339 */
340
EnumArgs()341 Array<GisCoords::Arg> GisCoordsConical::EnumArgs()
342 {
343 Array<GisCoords::Arg> out;
344 // out.Add(GisCoords::Arg::Option(gauss_sphere, "GAUSS", "Gauss preprocess"));
345 out.Add(GisCoords::Arg::Angle (gauss_base_parallel, "GAUSSBP", "Gauss base parallel"));
346 out.Add(GisCoords::Arg::Angle (pole_meridian, "POLEM", "Cone pole meridian"));
347 out.Add(GisCoords::Arg::Angle (pole_parallel, "POLEP", "Cone pole parallel"));
348 out.Add(GisCoords::Arg::Angle (central_parallel, "CP", "Central parallel"));
349 // out.Add(GisCoords::Arg::Angle (prime_meridian, "PM", "Prime meridian"));
350 out.Add(GisCoords::Arg::Angle (north_parallel, "NP", "North parallel"));
351 out.Add(GisCoords::Arg::Angle (south_parallel, "SP", "South parallel"));
352 out.Add(GisCoords::Arg::Edit (multiplier, "MUL", "Multiplier"));
353 out.Add(GisCoords::Arg::Edit (xmeteroffset, "DX", "X meter offset"));
354 out.Add(GisCoords::Arg::Edit (ymeteroffset, "DY", "Y meter offset"));
355 return out;
356 }
357
EnumConstArgs() const358 Array<GisCoords::ConstArg> GisCoordsConical::EnumConstArgs() const
359 {
360 Array<GisCoords::ConstArg> args;
361 // if(gauss_sphere)
362 {
363 args.Add(GisCoords::ConstArg("GA", "gauss_alpha", gauss_alpha));
364 args.Add(GisCoords::ConstArg("GK", "gauss_k", gauss_k));
365 args.Add(GisCoords::ConstArg("GR", "gauss_R", gauss_R));
366 args.Add(GisCoords::ConstArg("GE", "gauss_e", gauss_e));
367 args.Add(GisCoords::ConstArg("GU0", "gauss_U0", FormatDegree(gauss_U0, 3)));
368 }
369
370 args.Add(GisCoords::ConstArg("R0", "rho0", rho0));
371 args.Add(GisCoords::ConstArg("N", "n", n));
372 // args.Add(GisCoords::ConstArg("e", e));
373 args.Add(GisCoords::ConstArg("RC", "rho_coef", rho_coef));
374 return args;
375 }
376
SyncArgs()377 void GisCoordsConical::SyncArgs()
378 {
379 // double NN, NS;
380 // if(gauss_sphere)
381 {
382 gauss_e = sqrt(ellipsoid.e2);
383 gauss_alpha = sqrt(1 + (ellipsoid.e2 * sqr(sqr(cos(gauss_base_parallel * DEGRAD)))) / (1 - ellipsoid.e2));
384 double sinphi = sin(gauss_base_parallel * DEGRAD);
385 gauss_U0 = asin(sinphi / gauss_alpha);
386 gauss_k = exp(gauss_alpha * (log(tan(gauss_base_parallel * (DEGRAD / 2) + M_PI / 4))
387 + gauss_e / 2 * log((1 - gauss_e * sinphi) / (1 + gauss_e * sinphi)))) / tan(gauss_U0 / 2 + M_PI / 4);
388 // k = pow(tan(base_parallel / 2 + M_PI / 4), alpha)
389 // * pow((1 - e * sinphi) / (1 + e * sinphi), alpha * e / 2)
390 // / tan(U0 / 2 + M_PI / 4);
391 gauss_R = ellipsoid.a * sqrt(1 - ellipsoid.e2) / (1 - ellipsoid.e2 * sqr(sinphi));
392 // e = 0;
393 // NN = NS = gauss_R;
394 }
395 // else
396 // {
397 // e = sqrt(ellipsoid.e2);
398 // NN = ellipsoid.N(north_parallel);
399 // NS = ellipsoid.N(south_parallel);
400 // }
401 /*
402 double esin_n = e * sin(north_parallel), esin_s = e * sin(south_parallel);
403 double esin_0 = e * sin(central_parallel);
404 double nc_n = NN * cos(north_parallel);
405 double nc_s = NS * cos(south_parallel);
406 double tan_n = tan(north_parallel / 2 + M_PI / 4);
407 double tan_s = tan(south_parallel / 2 + M_PI / 4);
408 double tan_0 = tan(central_parallel / 2 + M_PI / 4);
409 n = (log(nc_n) - log(nc_s)) / (log(tan_s) - log(tan_n)
410 + e / 2 * (log((1 + esin_n) * (1 - esin_s)) - log((1 - esin_n) * (1 + esin_s))));
411 rho0 = nc_n / n * pow(tan_n / tan_0, n)
412 * pow((1 + esin_0) * (1 - esin_n) / ((1 - esin_0) * (1 + esin_n)), n * e / 2);
413 rho_coef = rho0 * pow(tan_0, n) * pow((1 - esin_0) / (1 + esin_0), n * e / 2);
414 */
415 //SphericalLatitudeFunction gslf(gauss_alpha, gauss_k, gauss_R, gauss_e, gauss_U0);
416 double Us = south_parallel * DEGRAD, Un = north_parallel * DEGRAD, Uc = central_parallel * DEGRAD;
417 double Tn = tan(Un / 2 + M_PI / 4), Ts = tan(Us / 2 + M_PI / 4), Tc = tan(Uc / 2 + M_PI / 4);
418 if(fabs(Un - Us) <= SECRAD)
419 {
420 rho0 = 1.0 / tan(Us);
421 n = sin(Us);
422 }
423 else
424 {
425 n = log(cos(Us) / cos(Un)) / log(Tn / Ts);
426 rho0 = cos(Us) / n * pow(Ts / Tc, n);
427 }
428 rho_coef = gauss_R * rho0 * pow(Tc, n) * multiplier;
429
430 gauss_projected.Clear();
431 gauss_latitude.Clear();
432 tgsin_projected.Clear();
433 tgsin_latitude.Clear();
434
435 orientation = Pointf(pole_meridian, pole_parallel);
436 }
437
DumpProjectedDelta() const438 String GisCoordsConical::DumpProjectedDelta() const
439 {
440 ConicalRadiusFunction lcrf(n, rho_coef);
441 return GisInterpolatorDelta(-45, 90.1, lcrf, 100, 10000, 8, 80000);
442 }
443
SyncInterpolator() const444 void GisCoordsConical::SyncInterpolator() const
445 {
446 // if(gauss_sphere)
447 {
448 enum { BUCKETS = 50, SPLITS = 200 };
449 LLOG("gauss_projected");
450 SphericalLatitudeFunction gslf(gauss_alpha, gauss_k, gauss_R, gauss_e, gauss_U0);
451 gauss_projected.Create(lonlat_limits.top, lonlat_limits.bottom, gslf, BUCKETS, SPLITS, 4);
452 LLOG("gauss_latitude");
453 gauss_latitude.CreateInverse(lonlat_limits.top, lonlat_limits.bottom, gslf, BUCKETS, SPLITS, 4);
454 }
455 {
456 enum { BUCKETS = 50, SPLITS = 200 };
457 ConicalRadiusFunction lcrf(n /*, e*/, rho_coef);
458 // lcrf.Dump(+1.4, +1.6, 1000);
459 double max_phi = min(0.0, gauss_projected(lonlat_limits.bottom) - pole_parallel) + 90;
460 double min_phi = gauss_projected(lonlat_limits.top);
461 LLOG("tgsin_projected");
462 tgsin_projected.Create(min_phi, max_phi, lcrf, BUCKETS, SPLITS, 4);
463 LLOG("tgsin_inverse");
464 tgsin_latitude.CreateInverse(min_phi, max_phi, lcrf, BUCKETS, SPLITS, 4);
465 }
466 }
467
468 static const double LAMBERT_LIMDIFF = 1 * DEGRAD;
469
ConicalRadiusFunction(double n,double rho_coef)470 ConicalRadiusFunction::ConicalRadiusFunction(double n /*, double e*/, double rho_coef)
471 : n(n), rho_coef(rho_coef) //, e(e)
472 {
473 static const double DELTA = 1e-6;
474 static const double ANG = M_PI / 2 - LAMBERT_LIMDIFF, ANG1 = ANG + DELTA;
475 double out0 = pow(cos(ANG) / sin(ANG), n), out1 = pow(cos(ANG1) / sin(ANG1), n);
476 double slope = (out0 - out1) / DELTA;
477 k1 = (out0 - LAMBERT_LIMDIFF * slope) / (LAMBERT_LIMDIFF * LAMBERT_LIMDIFF);
478 k0 = slope + 2 * k1 * LAMBERT_LIMDIFF;
479 }
480
Get(double phi) const481 double ConicalRadiusFunction::Get(double phi) const
482 {
483 // RTIMING("GisCoordsConical::CalcTgSinProjected");
484 phi *= DEGRAD;
485 double ang = phi / 2 + M_PI / 4;
486 double diff = M_PI / 2 - ang;
487 double out;
488 // static double prev = 0;
489 // bool limit = false;
490 out = (diff <= LAMBERT_LIMDIFF ? diff * (k0 - k1 * diff) : pow(cos(ang) / sin(ang), n));
491 return out * rho_coef;
492 /*
493 {
494 // limit = true;
495 static const double tglim = tan(0.001);
496 out = pow(tglim, n) * (1e3 * diff);
497 }
498 else
499 if(phi >= -1.5)
500 */
501 // else
502 // out = slope * phi + sat;
503 // LLOG(NFormat("tg_sin(%5nf): limit = %[0:no;yes]s, out = %10nf, diff = %10nf",
504 // phi, limit, out, out - prev));
505 // prev = out;
506 // double esin = e * sin(phi);
507 // return out * pow((1 + esin) / (1 - esin), n * e / 2);
508 }
509
pow3(double x)510 static double pow3(double x) { return x * x * x; }
511
GisCoordsUTM(double central_meridian_,double multiplier_,double xmeteroffset_,double ymeteroffset_)512 GisCoordsUTM::GisCoordsUTM(
513 double central_meridian_,
514 double multiplier_,
515 double xmeteroffset_,
516 double ymeteroffset_)
517 {
518 central_meridian = central_meridian_;
519 multiplier = multiplier_;
520 xmeteroffset = xmeteroffset_;
521 ymeteroffset = ymeteroffset_;
522 SyncArgs();
523 }
524
LonLat(Pointf xy) const525 Pointf GisCoordsUTM::LonLat(Pointf xy) const
526 {
527 // int paspol = ffloor(xy.x / 1000000);
528 xy.x = xy.x - xmeteroffset;
529 double xyx2 = xy.x * xy.x;
530 double b1 = (xy.y - ymeteroffset) / 111134.861084 * DEGRAD, b = b1;
531 b = b1 - (-0.002518467884 * sin(2 * b) + 0.26428e-5 * sin(4 * b) - 3.681e-9 * sin(6 * b));
532 b = b1 - (-0.002518467884 * sin(2 * b) + 0.26428e-5 * sin(4 * b) - 3.681e-9 * sin(6 * b));
533 double sb = sin(b), cb = cos(b);
534 double T = sb / cb, T2 = T * T;
535 double ETA2 = E12 * sqr(cb);
536 double N = sqr(ellipsoid.a) / (ellipsoid.b * sqrt(ETA2 + 1)), N2 = N * N;
537 double L = central_meridian + (1 / DEGRAD) * xy.x * (
538 + 1 / (N * cb)
539 - xyx2 / (6 * N * N2 * cb) * (1 + 2 * T2 + ETA2
540 + xyx2 * (5 + T2 * (28 + 24 * T2 + 8 * ETA2) + 6 * ETA2) / (20 * N2)));
541 b -= (T * xyx2 / (2 * N2)) * (1 + ETA2
542 + xyx2 * (5 + 3 * ETA2 * (2 - ETA2) + T2 * (3 - 6 * ETA2 - 9 * sqr(ETA2))) / (12 * N2));
543 return Pointf(L, b / DEGRAD);
544 }
545
Project(Pointf lonlat,int branch) const546 Pointf GisCoordsUTM::Project(Pointf lonlat, int branch) const
547 {
548 double LL = (lonlat.x - central_meridian) * DEGRAD;
549 double LL2 = LL * LL;//, LL3 = LL2 * LL;
550 lonlat *= DEGRAD;
551 double sx = sin(lonlat.y), cx = cos(lonlat.y);
552 double cx2 = cx * cx;//, cx3 = cx2 * cx;
553 double B = 111134.861084 / DEGRAD * lonlat.y - 16036.480269 * sin(2 * lonlat.y)
554 + 16.828067 * sin(4 * lonlat.y) - 0.021975 * sin(6 * lonlat.y) + 0.000031 * sin(8 * lonlat.y);
555 double T = sx / cx, T2 = T * T;
556 double ETA2 = ellipsoid.e2 * cx2;
557 double N = ellipsoid.a / sqrt(1 - E21 * sqr(sx));
558
559 double xout = B + N * sx * cx * (LL2 / 2) * (1 + cx2 * (5 - T2 + (9 + 4 * ETA2) * ETA2) * (LL2 / 12))
560 + ymeteroffset;
561
562 double yout = N * cx * LL * (1
563 + cx2 * (LL2 / 6) * ((1 - T2 + ETA2) + cx2 * (5 - 18 * pow3(T2) + (14 - 58 * T2) * ETA2) * (LL2 / 20)))
564 + xmeteroffset;
565
566 return Pointf(yout, xout);
567 }
568
EnumArgs()569 Array<GisCoords::Arg> GisCoordsUTM::EnumArgs()
570 {
571 Array<GisCoords::Arg> out;
572 // out.Add(GisCoords::Arg::Option(gauss_sphere, "GAUSS", "Gauss preprocess"));
573 out.Add(GisCoords::Arg::Angle (central_meridian, "CM", "Central meridian"));
574 out.Add(GisCoords::Arg::Edit (multiplier, "MUL", "Multiplier"));
575 out.Add(GisCoords::Arg::Edit (xmeteroffset, "DX", "X meter offset"));
576 out.Add(GisCoords::Arg::Edit (ymeteroffset, "DY", "Y meter offset"));
577 return out;
578 }
579
EnumConstArgs() const580 Array<GisCoords::ConstArg> GisCoordsUTM::EnumConstArgs() const
581 {
582 Array<GisCoords::ConstArg> args;
583 args.Add(GisCoords::ConstArg("E21", "1 - b^2/a^2", E21));
584 args.Add(GisCoords::ConstArg("E12", "a^2/b^2 - 1", E12));
585 return args;
586 }
587
SyncArgs()588 void GisCoordsUTM::SyncArgs()
589 {
590 E21 = 1 - sqr(ellipsoid.b / ellipsoid.a);
591 E12 = sqr(ellipsoid.a / ellipsoid.b) - 1;
592 }
593
GisCoordsAzimuthal(const Pointf & p,const Sizef & sc,const Sizef & o)594 GisCoordsAzimuthal::GisCoordsAzimuthal(const Pointf& p, const Sizef& sc, const Sizef& o)
595 {
596 pole = p;
597 scale = sc;
598 offset = o;
599 }
600
DeepCopy() const601 GisCoords GisCoordsAzimuthal::DeepCopy() const
602 {
603 One<GisCoordsAzimuthal> out = new GisCoordsAzimuthal(pole, scale, offset);
604 out->ellipsoid = ellipsoid;
605 out->SyncArgs();
606 return GisCoords(out.Detach());
607 }
608
GetBranchCount() const609 int GisCoordsAzimuthal::GetBranchCount() const
610 {
611 return 1;
612 }
613
LonLatBoundary(const Rectf & lonlat_extent,int branch) const614 Array<Pointf> GisCoordsAzimuthal::LonLatBoundary(const Rectf& lonlat_extent, int branch) const
615 {
616 return Array<Pointf>();
617 }
618
GetBranchTree(const Rectf & lonlat_extent) const619 GisBSPTree GisCoordsAzimuthal::GetBranchTree(const Rectf& lonlat_extent) const
620 {
621 return 0;
622 }
623
LonLat(Pointf xy) const624 Pointf GisCoordsAzimuthal::LonLat(Pointf xy) const
625 {
626 xy -= offset;
627 xy /= scale;
628 double rho = Length(xy);
629 double eps = Bearing(xy) / DEGRAD;
630 double sine = min(rho / (2 * gauss.radius), 1.0);
631 double psi = 2 * asin(sine);
632 //double psi = rho / Rdeg;
633 //return Pointf(eps, M_PI / 2 - psi);
634 Pointf out = orientation.Global(eps, 90 - psi);
635 //out.y = gauss.Elliptical(out.y);
636 return out;
637 }
638
639 /*
640 Rectf GisCoordsAzimuthal::LonLatExtent(const Rectf& xy_extent) const
641 {
642 Rectf lonlat = ExtentToDegree(xy_extent);
643 return orientation.GlobalExtent(lonlat.left, 90 - lonlat.bottom / Rdeg, lonlat.right, 90 - lonlat.top / Rdeg);
644 }
645 */
646
Project(Pointf lonlat,int branch) const647 Pointf GisCoordsAzimuthal::Project(Pointf lonlat, int branch) const
648 {
649 lonlat.y = gauss.Spherical(lonlat.y);
650 lonlat = orientation.Local(lonlat);
651 double psi = DEGRAD * (90 - lonlat.y);
652 //double rho = R * E * sin(psi) / (c + R * cos(psi));
653 //double rho = Rdeg * psi;
654 double rho = 2 * gauss.radius * sin(psi / 2);
655 return scale * (Sizef)PolarPointf(rho, (lonlat.x - 90) * DEGRAD) + offset;
656 }
657
658 /*
659 Rectf GisCoordsAzimuthal::ProjectExtent(const Rectf& lonlat_extent) const
660 {
661 Rectf local = orientation.LocalExtent(lonlat_extent);
662 return DegreeToExtent(local.left, Rdeg * (90 - local.bottom), local.right, Rdeg * (90 - local.top));
663 }
664 */
665
666 /*
667 double GisCoordsAzimuthal::ProjectDeviation(Pointf lonlat1, Pointf lonlat2, int branch) const
668 {
669 return 0;
670 }
671 */
672
673 /*
674 double GisCoordsAzimuthal::ProjectRatio(Pointf lonlat, int branch) const
675 {
676 return 1;
677 }
678 */
679
EnumArgs()680 Array<GisCoords::Arg> GisCoordsAzimuthal::EnumArgs()
681 {
682 Array<GisCoords::Arg> out;
683 out.Add(GisCoords::Arg::Angle(pole.x, "POLE_M", "Pole meridian", "", -180, +180));
684 out.Add(GisCoords::Arg::Angle(pole.y, "POLE_P", "Pole parallel", "", -90, +90));
685 out.Add(GisCoords::Arg::Edit(offset.cx, "X_OFFSET", "False easting"));
686 out.Add(GisCoords::Arg::Edit(offset.cy, "Y_OFFSET", "False northing"));
687 return out;
688 }
689
EnumConstArgs() const690 Array<GisCoords::ConstArg> GisCoordsAzimuthal::EnumConstArgs() const
691 {
692 Array<GisCoords::ConstArg> out;
693 return out;
694 }
695
SyncArgs()696 void GisCoordsAzimuthal::SyncArgs()
697 {
698 gauss.Create(ellipsoid.a, ellipsoid.e2, pole.y);
699 gauss.radius = 6384000;
700 orientation = Pointf(pole.x, gauss.Spherical(pole.y));
701 }
702
SyncInterpolator() const703 void GisCoordsAzimuthal::SyncInterpolator() const
704 {
705 }
706
EnumEPSG()707 const Vector<int>& GisCoords::EnumEPSG()
708 {
709 static Vector<int> epsg;
710 if(epsg.IsEmpty())
711 epsg
712 << 2065 // JTSK
713 << 3035 // ETRS89 / ETRS-LAEA
714 << 4326 // WGS-84
715 << 4818 // JTSK-geo
716 << 28403 // S-42
717 << 102065 // JTSK
718 << 102067 // JTSK ?
719 ;
720 return epsg;
721 }
722
GetEPSG(int code)723 GisCoords GisCoords::GetEPSG(int code)
724 {
725 GeomRefPtr<GisCoords::Data> gc;
726 switch(code)
727 {
728 case 2065:
729 case 102065:
730 case 102067:
731 gc = new GisCoordsConical(
732 49.5,
733 ANGDEG(24, 50, 53.41635),
734 ANGDEG(59, 42, 42.6969),
735 78.5,
736 78.5,
737 78.5,
738 0.9999,
739 0, 0);
740 gc->lonlat_limits = Rectf(10, 45, 20, 55);
741 gc->ellipsoid = GisEllipsoid::GetEPSG(GisEllipsoid::BESSEL_1841);
742 // gc->ellipsoid = GisEllipsoid::GetEPSG(GisEllipsoid::WGS_1984);
743 gc->name = "S-JTSK (Ferro) / Krovak";
744 gc->coordsys =
745 "PROJCS[\"S-JTSK (Ferro) / Krovak\", "
746 "GEOGCS[\"S-JTSK (Ferro)\", "
747 "DATUM[\"D_S_JTSK_Ferro\", "
748 "SPHEROID[\"Bessel 1841\",6377397.155,299.1528128]], "
749 "PRIMEM[\"Ferro\",-17.66666666666667], "
750 "UNIT[\"degree\",0.0174532925199433]], "
751 "PROJECTION[\"Krovak\"], "
752 "PARAMETER[\"latitude_of_center\",49.5], "
753 "PARAMETER[\"longitude_of_center\",42.5], "
754 "PARAMETER[\"azimuth\",30.28813972222222], "
755 "PARAMETER[\"pseudo_standard_parallel_1\",78.5], "
756 "PARAMETER[\"scale_factor\",0.9999], "
757 "PARAMETER[\"false_easting\",0], "
758 "PARAMETER[\"false_northing\",0], "
759 "UNIT[\"metre\",1]]";
760 break;
761
762 case 3035: {
763 gc = new GisCoordsAzimuthal(
764 Pointf(10, 52),
765 Sizef(1, 1),
766 Sizef(4321000, 3210000));
767 gc->ellipsoid = GisEllipsoid::GetEPSG(GisEllipsoid::WGS_1984); //GRS_1980);
768 gc->name = "ETRS89 / ETRS-LAEA";
769 gc->description = "Single CRS for all Europe";
770 gc->coordsys =
771 "PROJCS[\"ETRS89 / ETRS-LAEA\","
772 "GEOGCS[\"ETRS89\","
773 "DATUM[\"European_Terrestrial_Reference_System_1989\","
774 "SPHEROID[\"GRS 1980\",6378137,298.257222101,"
775 "AUTHORITY[\"EPSG\",\"7019\"]],"
776 "AUTHORITY[\"EPSG\",\"6258\"]],"
777 "PRIMEM[\"Greenwich\",0,"
778 "AUTHORITY[\"EPSG\",\"8901\"]],"
779 "UNIT[\"degree\",0.01745329251994328,"
780 "AUTHORITY[\"EPSG\",\"9122\"]],"
781 "AUTHORITY[\"EPSG\",\"4258\"]],"
782 "UNIT[\"metre\",1,"
783 "AUTHORITY[\"EPSG\",\"9001\"]],"
784 "PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],"
785 "PARAMETER[\"latitude_of_center\",52],"
786 "PARAMETER[\"longitude_of_center\",10],"
787 "PARAMETER[\"false_easting\",4321000],"
788 "PARAMETER[\"false_northing\",3210000],"
789 "AUTHORITY[\"EPSG\",\"3035\"],"
790 "AXIS[\"X\",EAST],"
791 "AXIS[\"Y\",NORTH]]"
792 ;
793 break;
794 }
795
796 case 4326:
797 gc = new GisCoordsLonLat(0);
798 gc->ellipsoid = GisEllipsoid::GetEPSG(GisEllipsoid::WGS_1984);
799 gc->name = "WGS 84";
800 gc->description = "World Geodetic System 1984";
801 gc->coordsys =
802 "GEOGCS[\"WGS 84\", "
803 "DATUM[\"WGS_1984\", "
804 "SPHEROID[\"WGS 84\",6378137,298.257223563, "
805 "AUTHORITY[\"EPSG\",7030]], "
806 "TOWGS84[0,0,0,0,0,0,0], "
807 "AUTHORITY[\"EPSG\",6326]], "
808 "PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",8901]], "
809 "UNIT[\"DMSH\",0.0174532925199433,AUTHORITY[\"EPSG\",9108]], "
810 "AXIS[\"Lat\",NORTH], "
811 "AXIS[\"Long\",EAST], "
812 "AUTHORITY[\"EPSG\",4326]]";
813 break;
814
815 case 4818:
816 gc = new GisCoordsLonLat(-ANGDEG(17, 40, 0)); // Ferro
817 gc->ellipsoid = GisEllipsoid::GetEPSG(GisEllipsoid::BESSEL_1841);
818 gc->name = "S-JTSK (Ferro)";
819 gc->description = "Modification of Austrian MGI (Ferro) datum.";
820 break;
821
822 case 28403:
823 gc = new GisCoordsUTM(15, 1, 3500e3, 0);
824 gc->ellipsoid = GisEllipsoid::GetEPSG(GisEllipsoid::KRASSOWSKY_1940);
825 gc->name = "Pulkovo 1942 / Gauss-Kruger zone 3";
826 gc->description = "Military mapping. Czech Republic and Germany ( former DDR) - east of 12 deg East; "
827 "Poland and Slovakia - west of 18 deg East.";
828 gc->lonlat_limits = Rectf(12, 45, 18, 55);
829 gc->coordsys =
830 "PROJCS[\"Pulkovo 1942(83) / Gauss Kruger zone 3\", "
831 "GEOGCS[\"Pulkovo 1942(83)\", "
832 "DATUM[\"Pulkovo_1942_83\", "
833 "SPHEROID[\"Krassowsky 1940\",6378245,298.3, "
834 "AUTHORITY[\"EPSG\",\"7024\"]], "
835 "TOWGS84[24,-123,-94,0.02,-0.25,-0.13,1.1], "
836 "AUTHORITY[\"EPSG\",\"6178\"]], "
837 "PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]], "
838 "UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]], "
839 "AUTHORITY[\"EPSG\",\"4178\"]], "
840 "PROJECTION[\"Transverse_Mercator\"], "
841 "PARAMETER[\"latitude_of_origin\",0], "
842 "PARAMETER[\"central_meridian\",9], "
843 "PARAMETER[\"scale_factor\",1], "
844 "PARAMETER[\"false_easting\",3500000], "
845 "PARAMETER[\"false_northing\",0], "
846 "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]], "
847 "AUTHORITY[\"EPSG\",\"2166\"]]";
848 break;
849
850 default:
851 return GisCoords();
852 }
853 gc->ident = Format("EPSG:%d", code);
854 gc->SyncArgs();
855 return GisCoords(gc);
856 }
857
858 }
859