1 /*
2 * geoProj.c --
3 *
4 * This file defines functions that convert back and forth between
5 * geopoints and map-points. See the user documentation for more
6 * information.
7 *
8 * Copyright (c) 2004 Gordon D. Carrie. All rights reserved.
9 *
10 * Licensed under the Open Software License version 2.1
11 *
12 * Please address questions and feedback to user0@tkgeomap.org
13 *
14 * @(#) $Id: geoProj.c,v 1.9 2007/06/27 18:38:56 tkgeomap Exp $
15 *
16 ********************************************
17 *
18 */
19
20 /*
21 * Formulas and notation for most projections are from:
22 * Snyder, John P.
23 * Map Projections used by the U.S. Geological Survey.
24 * (Geological Survey bulletin ; 1532)
25 * United States Government Printing Office, Washington: 1982.
26 */
27
28 #include <math.h>
29 #include <limits.h>
30 #include <stdlib.h>
31 #include <string.h>
32 #include <stdio.h>
33 #include <assert.h>
34 #include <tcl.h>
35 #include "geography.h"
36 #include "geoProj.h"
37
38 /*
39 * Local functions.
40 */
41
42 static MapPt latLonToCylEqDist _ANSI_ARGS_((GeoPt geoPt,
43 struct GeoProj *projPtr));
44 static GeoPt cylEqDistToLatLon _ANSI_ARGS_((MapPt mapPt,
45 struct GeoProj *projPtr));
46 static GeoProjInfo cylEqDistInfo _ANSI_ARGS_((struct GeoProj *projPtr));
47 static MapPt latLonToCylEqArea _ANSI_ARGS_((GeoPt geoPt,
48 struct GeoProj *projPtr));
49 static GeoPt cylEqAreaToLatLon _ANSI_ARGS_((MapPt mapPt,
50 struct GeoProj *projPtr));
51 static MapPt latLonToMercator _ANSI_ARGS_((GeoPt geoPt,
52 struct GeoProj *projPtr));
53 static GeoPt mercatorToLatLon _ANSI_ARGS_((MapPt mapPt,
54 struct GeoProj *projPtr));
55 static MapPt latLonToLambertConfConic _ANSI_ARGS_((GeoPt geoPt,
56 struct GeoProj *projPtr));
57 static GeoPt lambrtConfConicToLatLon _ANSI_ARGS_((MapPt mapPt,
58 struct GeoProj *projPtr));
59 static GeoProjInfo lambrtConfConicInfo _ANSI_ARGS_((
60 struct GeoProj *projPtr));
61 static MapPt latLonToLambrtEqArea _ANSI_ARGS_((GeoPt geoPt,
62 struct GeoProj *projPtr));
63 static GeoPt lambrtEqAreaToLatLon _ANSI_ARGS_((MapPt mapPt,
64 struct GeoProj *projPtr));
65 static MapPt latLonToStereographic _ANSI_ARGS_((GeoPt geoPt,
66 struct GeoProj *projPtr));
67 static GeoPt stereographicToLatLon _ANSI_ARGS_((MapPt mapPt,
68 struct GeoProj *projPtr));
69 static MapPt latLonToOrthographic _ANSI_ARGS_((GeoPt geoPt,
70 struct GeoProj *projPtr));
71 static GeoPt orthographicToLatLon _ANSI_ARGS_((MapPt mapPt,
72 struct GeoProj *projPtr));
73 static GeoProjInfo refLonProjInfo _ANSI_ARGS_((struct GeoProj *projPtr));
74 static GeoProjInfo refPtProjInfo _ANSI_ARGS_((struct GeoProj *projPtr));
75
76 /*
77 *----------------------------------------------------------------------
78 *
79 * GeoProjAlloc --
80 *
81 * This procedure allocates a new GeoProj structure.
82 *
83 * Results:
84 * See the user documentation.
85 *
86 * Side effects:
87 * A GeoProj structure is allocated and initialized.
88 * The memory should eventually be freed with a call to GeoProjDestroy.
89 *
90 *----------------------------------------------------------------------
91 */
92
93 GeoProj
GeoProjAlloc()94 GeoProjAlloc()
95 {
96 struct GeoProj *projPtr;
97
98 projPtr = (GeoProj )CKALLOC(sizeof(*projPtr));
99 GeoProjInit(projPtr);
100 return projPtr;
101 }
102
103 /*
104 *----------------------------------------------------------------------
105 *
106 * GeoProjInit --
107 *
108 * This procedure initializes a GeoProj structure.
109 *
110 * Results:
111 * None.
112 *
113 * Side effects:
114 * Modifies the member of a GeoProj structure. Previous contents,
115 * assumed to be garbage, are overwritten.
116 *
117 *----------------------------------------------------------------------
118 */
119
120 void
GeoProjInit(projPtr)121 GeoProjInit(projPtr)
122 struct GeoProj *projPtr;
123 {
124 projPtr->params = NULL;
125 SetCylEqDist(projPtr, 0, 0);
126 GeoProjSetRotation(projPtr, AngleFmDeg(0.0));
127 }
128
129 /*
130 *----------------------------------------------------------------------
131 *
132 * GeoProjFree --
133 *
134 * This procedure frees the internal allocations in a GeoProj structure.
135 *
136 * Results:
137 * None.
138 *
139 * Side effects:
140 * See the user documentation.
141 *
142 *----------------------------------------------------------------------
143 */
144
145 void
GeoProjFree(projPtr)146 GeoProjFree(projPtr)
147 struct GeoProj *projPtr;
148 {
149 if( projPtr ) {
150 CKFREE((char *)projPtr->params);
151 }
152 }
153
154 /*
155 *----------------------------------------------------------------------
156 *
157 * GeoProjDestroy --
158 *
159 * This procedure frees internal storage of a GeoProj structure and
160 * frees the structure.
161 *
162 * Results:
163 * None.
164 *
165 * Side effects:
166 * See the user documentation.
167 *
168 *----------------------------------------------------------------------
169 */
170
171 void
GeoProjDestroy(projPtr)172 GeoProjDestroy(projPtr)
173 struct GeoProj *projPtr;
174 {
175 GeoProjFree(projPtr);
176 CKFREE((char *)projPtr);
177 }
178
179 /*
180 *----------------------------------------------------------------------
181 *
182 * ProjToLatLon --
183 *
184 * This procedure converts a map-point to a geopoint.
185 *
186 * Results:
187 * See the user documentation.
188 *
189 * Side effects:
190 * See the user documentation.
191 *
192 *----------------------------------------------------------------------
193 */
194
195 GeoPt
ProjToLatLon(mapPt,projPtr)196 ProjToLatLon(mapPt, projPtr)
197 MapPt mapPt;
198 struct GeoProj *projPtr;
199 {
200 if ( MapPtIsNowhere(mapPt) ) {
201 return GeoPtNowhere();
202 }
203 if (projPtr->rotation != 0) {
204 double abs, ord;
205 abs = mapPt.abs * projPtr->cosr - mapPt.ord * projPtr->sinr;
206 ord = mapPt.abs * projPtr->sinr + mapPt.ord * projPtr->cosr;
207 mapPt.abs = abs;
208 mapPt.ord = ord;
209 }
210 return projPtr->projToLatLonProc(mapPt, projPtr);
211 }
212
213 /*
214 *----------------------------------------------------------------------
215 *
216 * LatLonToProj --
217 *
218 * This procedure converts a geopoint to a map-point.
219 *
220 * Results:
221 * See the user documentation.
222 *
223 * Side effects:
224 * See the user documentation.
225 *
226 *----------------------------------------------------------------------
227 */
228
229 MapPt
LatLonToProj(geoPt,projPtr)230 LatLonToProj(geoPt, projPtr)
231 GeoPt geoPt;
232 struct GeoProj *projPtr;
233 {
234 MapPt mapPt;
235
236 if ( GeoPtIsNowhere(geoPt) ) {
237 return MapPtNowhere();
238 }
239 mapPt = projPtr->latLonToProjProc(geoPt, projPtr);
240 if (projPtr->rotation != 0) {
241 double abs, ord;
242 abs = mapPt.abs * projPtr->cosr + mapPt.ord * projPtr->sinr;
243 ord = mapPt.ord * projPtr->cosr - mapPt.abs * projPtr->sinr;
244 mapPt.abs = abs;
245 mapPt.ord = ord;
246 }
247 return mapPt;
248 }
249
250 /*
251 *----------------------------------------------------------------------
252 *
253 * GeoProjGetInfo --
254 *
255 * This is a member access function.
256 *
257 * Results:
258 * See the user documentation.
259 *
260 * Side effects:
261 * See the user documentation.
262 *
263 *----------------------------------------------------------------------
264 */
265
266 GeoProjInfo
GeoProjGetInfo(projPtr)267 GeoProjGetInfo(projPtr)
268 struct GeoProj *projPtr;
269 {
270 return projPtr->infoProc(projPtr);
271 }
272
273 /*
274 *----------------------------------------------------------------------
275 *
276 * GeoProjGetType --
277 *
278 * This is a member access function.
279 *
280 * Results:
281 * See the user documentation.
282 *
283 * Side effects:
284 * See the user documentation.
285 *
286 *----------------------------------------------------------------------
287 */
288
289 enum ProjType
GeoProjGetType(projPtr)290 GeoProjGetType(projPtr)
291 struct GeoProj *projPtr;
292 {
293 return projPtr->type;
294 }
295
296 /*
297 *----------------------------------------------------------------------
298 *
299 * GeoProjDescriptor --
300 *
301 * This is a member access function.
302 *
303 * Results:
304 * See the user documentation.
305 *
306 * Side effects:
307 * See the user documentation.
308 *
309 *----------------------------------------------------------------------
310 */
311
312 char *
GeoProjDescriptor(projPtr)313 GeoProjDescriptor(projPtr)
314 struct GeoProj *projPtr;
315 {
316 return projPtr->descriptor;
317 }
318
319 /*
320 *----------------------------------------------------------------------
321 *
322 * GeoProjSetRotation --
323 *
324 * This function determines how much a map is rotated.
325 *
326 * Results:
327 * None.
328 *
329 * Side effects:
330 * The rotation angle an associated members are set in the structure.
331 *
332 *----------------------------------------------------------------------
333 */
334
335 void
GeoProjSetRotation(projPtr,angle)336 GeoProjSetRotation(projPtr, angle)
337 struct GeoProj *projPtr; /* Projection to modify */
338 Angle angle; /* Projection angle */
339 {
340 projPtr->rotation = angle;
341 projPtr->cosr = ICos(angle);
342 projPtr->sinr = ISin(angle);
343 }
344
345 /*
346 *----------------------------------------------------------------------
347 *
348 * SetCylEqDist --
349 *
350 * This procedure sets a projection to Cylindrical Equidistant.
351 *
352 * Results:
353 * None.
354 *
355 * Side effects:
356 * The members of a GeoProj structure are modified.
357 *
358 *----------------------------------------------------------------------
359 */
360
361 void
SetCylEqDist(projPtr,refLat,refLon)362 SetCylEqDist(projPtr, refLat, refLon)
363 struct GeoProj *projPtr; /* Projection to modify */
364 Angle refLat; /* Reference latitude */
365 Angle refLon; /* Reference longitude */
366 {
367 GeoProjCylEqDistParams *params;
368
369 projPtr->type = CylEqDist;
370 params = (GeoProjCylEqDistParams *)CKALLOC(sizeof(GeoProjCylEqDistParams));
371 params->refLat = DomainLat(refLat);
372 params->cosRLat = ICos(refLat);
373 params->refLon = GwchLon(refLon);
374 if (projPtr->params) {
375 CKFREE(projPtr->params);
376 }
377 projPtr->params = params;
378 projPtr->latLonToProjProc = latLonToCylEqDist;
379 projPtr->projToLatLonProc = cylEqDistToLatLon;
380 sprintf(projPtr->descriptor, "CylEqDist {%9.3f %-9.3f}",
381 AngleToDeg(params->refLat), AngleToDeg(params->refLon));
382 projPtr->infoProc = cylEqDistInfo;
383 }
384
385 /*
386 *----------------------------------------------------------------------
387 *
388 * cylEqDistInfo --
389 *
390 * This is a value for an infoProc member in a GeoProj structure.
391 * It is usually called by the GeoProjGetInfo procedure
392 *
393 * Results:
394 * Projection information is copied into a GeoProjInfo structure.
395 *
396 * Side effects:
397 * None.
398 *
399 *----------------------------------------------------------------------
400 */
401
402 static GeoProjInfo
cylEqDistInfo(projPtr)403 cylEqDistInfo(projPtr)
404 struct GeoProj *projPtr;
405 {
406 GeoProjInfo geoProjInfo;
407 GeoProjCylEqDistParams *params = projPtr->params;
408
409 geoProjInfo.type = projPtr->type;
410 geoProjInfo.prm.refPt.lat = params->refLat;
411 geoProjInfo.prm.refPt.lon = params->refLon;
412 geoProjInfo.rotation = projPtr->rotation;
413
414 return geoProjInfo;
415 }
416
417 /*
418 *----------------------------------------------------------------------
419 *
420 * latLonToCylEqDist --
421 *
422 * This is a value for a latLonToProjProc member of a GeoProj structure.
423 * It is usually called by the LatLonToProj procedure.
424 *
425 * Results:
426 * A map-point corresponding to a geopoint.
427 *
428 * Side effects:
429 * None.
430 *
431 *----------------------------------------------------------------------
432 */
433
434 static MapPt
latLonToCylEqDist(geoPt,projPtr)435 latLonToCylEqDist(geoPt, projPtr)
436 GeoPt geoPt;
437 struct GeoProj *projPtr;
438 {
439 MapPt mapPt;
440 double lat, lon;
441 GeoProjCylEqDistParams *params = projPtr->params;
442
443 GeoPtGetDeg(DomainLonPt(geoPt, params->refLon), &lat, &lon);
444 mapPt.ord = MPERDEG * lat;
445 mapPt.abs = MPERDEG * params->cosRLat * lon;
446 return mapPt;
447 }
448
449 /*
450 *----------------------------------------------------------------------
451 *
452 * cylEqDistToLatLon --
453 *
454 * This is a value for a projToLatLonProc member of a GeoProj structure.
455 * It is usually called by the ProjToLatLon procedure.
456 *
457 * Results:
458 * A geopoint corresponding to a map-point.
459 *
460 * Side effects:
461 * None.
462 *
463 *----------------------------------------------------------------------
464 */
465
466 static GeoPt
cylEqDistToLatLon(mapPt,projPtr)467 cylEqDistToLatLon(mapPt, projPtr)
468 MapPt mapPt;
469 struct GeoProj *projPtr;
470 {
471 GeoPt geoPt;
472 GeoProjCylEqDistParams *params = projPtr->params;
473
474 geoPt = GeoPtFmDeg(DEGPERM * mapPt.ord,
475 DEGPERM * mapPt.abs / params->cosRLat);
476 return DomainLonPt(geoPt, params->refLon);
477 }
478
479 /*
480 *----------------------------------------------------------------------
481 *
482 * SetMercator --
483 *
484 * This procedure sets a projection to Mercator.
485 *
486 * Results:
487 * None.
488 *
489 * Side effects:
490 * The members of a GeoProj structure are modified.
491 *
492 *----------------------------------------------------------------------
493 */
494
495 void
SetMercator(projPtr,refLon)496 SetMercator(projPtr, refLon)
497 struct GeoProj *projPtr; /* Projection to modify */
498 Angle refLon; /* Reference longitude */
499 {
500 Angle *refLonPtr;
501
502 refLonPtr = (Angle *)CKALLOC(sizeof(Angle));
503 projPtr->type = Mercator;
504 *refLonPtr = GwchLon(refLon);
505 if (projPtr->params) {
506 CKFREE(projPtr->params);
507 }
508 projPtr->params = refLonPtr;
509 projPtr->latLonToProjProc = latLonToMercator;
510 projPtr->projToLatLonProc = mercatorToLatLon;
511 sprintf(projPtr->descriptor, "Mercator %-9.3f", AngleToDeg(*refLonPtr));
512 projPtr->infoProc = refLonProjInfo;
513 }
514
515 /*
516 *----------------------------------------------------------------------
517 *
518 * latLonToMercator --
519 *
520 * This is a value for a latLonToProjProc member of a GeoProj structure.
521 * It is usually called by the LatLonToProj procedure.
522 *
523 * Results:
524 * A map-point corresponding to a geopoint.
525 *
526 * Side effects:
527 * None.
528 *
529 *----------------------------------------------------------------------
530 */
531
532 static MapPt
latLonToMercator(geoPt,projPtr)533 latLonToMercator(geoPt, projPtr)
534 GeoPt geoPt;
535 struct GeoProj *projPtr;
536 {
537 Angle *refLon = projPtr->params;
538 double lat, lon;
539 MapPt mapPt;
540 Angle d90 = AngleFmDeg(90.0);
541
542 /*
543 * Undefined at pole
544 */
545
546 if (LatCmp(geoPt.lat, d90) == Eq || LatCmp(geoPt.lat, -d90) == Eq) {
547 return MapPtNowhere();
548 }
549
550 GeoPtGetRad(DomainLonPt(geoPt, *refLon), &lat, &lon);
551 mapPt.ord = REarth() * log(tan(M_PI_4 + 0.5 * lat));
552 mapPt.abs = REarth() * lon;
553 return mapPt;
554 }
555
556 /*
557 *----------------------------------------------------------------------
558 *
559 * mercatorToLatLon --
560 *
561 * This is a value for a projToLatLonProc member of a GeoProj structure.
562 * It is usually called by the ProjToLatLon procedure.
563 *
564 * Results:
565 * A geopoint corresponding to a map-point.
566 *
567 * Side effects:
568 * None.
569 *
570 *----------------------------------------------------------------------
571 */
572
573 static GeoPt
mercatorToLatLon(mapPt,projPtr)574 mercatorToLatLon(mapPt, projPtr)
575 MapPt mapPt;
576 struct GeoProj *projPtr;
577 {
578 double lat, lon; /* Latitude and longitude in radians */
579 Angle *refLon = projPtr->params;
580
581 lat = 2.0 * (atan(exp(mapPt.ord / REarth())) - M_PI_4);
582 lon = mapPt.abs / REarth();
583 return DomainLonPt(GeoPtFmRad(lat, lon), *refLon);
584 }
585
586 /*
587 *----------------------------------------------------------------------
588 *
589 * SetCylEqArea --
590 *
591 * This procedure sets a projection to Cylindrical Equal Area.
592 *
593 * Results:
594 * None.
595 *
596 * Side effects:
597 * The members of a GeoProj structure are modified.
598 *
599 *----------------------------------------------------------------------
600 */
601
602 void
SetCylEqArea(projPtr,refLon)603 SetCylEqArea(projPtr, refLon)
604 struct GeoProj *projPtr; /* Projection to modify */
605 Angle refLon; /* Reference longitude */
606 {
607 Angle *refLonPtr;
608
609 refLonPtr = (Angle *)CKALLOC(sizeof(Angle));
610 projPtr->type = CylEqArea;
611 *refLonPtr = GwchLon(refLon);
612 projPtr->latLonToProjProc = latLonToCylEqArea;
613 projPtr->projToLatLonProc = cylEqAreaToLatLon;
614 sprintf(projPtr->descriptor, "CylEqArea %-9.3f", AngleToDeg(*refLonPtr));
615 projPtr->infoProc = refLonProjInfo;
616 }
617
618 /*
619 *----------------------------------------------------------------------
620 *
621 * latLonToCylEqArea --
622 *
623 * This is a value for a latLonToProjProc member of a GeoProj structure.
624 * It is usually called by the LatLonToProj procedure.
625 *
626 * Results:
627 * A map-point corresponding to a geopoint.
628 *
629 * Side effects:
630 * None.
631 *
632 *----------------------------------------------------------------------
633 */
634
635 static MapPt
latLonToCylEqArea(geoPt,projPtr)636 latLonToCylEqArea(geoPt, projPtr)
637 GeoPt geoPt;
638 struct GeoProj *projPtr;
639 {
640 MapPt mapPt;
641 Angle *refLon = projPtr->params;
642 double lat, lon;
643
644 GeoPtGetRad(DomainLonPt(geoPt, *refLon), &lat, &lon);
645 mapPt.ord = REarth() * sin(lat);
646 mapPt.abs = REarth() * lon;
647 return mapPt;
648 }
649
650 /*
651 *----------------------------------------------------------------------
652 *
653 * cylEqAreaToLatLon --
654 *
655 * This is a value for a projToLatLonProc member of a GeoProj structure.
656 * It is usually called by the ProjToLatLon procedure.
657 *
658 * Results:
659 * A geopoint corresponding to a map-point.
660 *
661 * Side effects:
662 * None.
663 *
664 *----------------------------------------------------------------------
665 */
666
667 static GeoPt
cylEqAreaToLatLon(mapPt,projPtr)668 cylEqAreaToLatLon(mapPt, projPtr)
669 MapPt mapPt;
670 struct GeoProj *projPtr;
671 {
672 Angle *refLon = projPtr->params;
673 double r, lat, lon;
674
675 r = mapPt.ord / REarth();
676 if (fabs(r) > 1.0) {
677 return GeoPtNowhere();
678 }
679 lat = asin(r);
680 lon = mapPt.abs / REarth();
681 return DomainLonPt(GeoPtFmRad(lat, lon), *refLon);
682 }
683
684 /*
685 *----------------------------------------------------------------------
686 *
687 * SetLambertConfConic --
688 *
689 * This procedure sets a projection to Lambert Conformal Conic.
690 *
691 * Results:
692 * None.
693 *
694 * Side effects:
695 * The members of a GeoProj structure are modified.
696 *
697 *----------------------------------------------------------------------
698 */
699
700 void
SetLambertConfConic(projPtr,refLat,refLon)701 SetLambertConfConic(projPtr, refLat, refLon)
702 struct GeoProj *projPtr; /* Projection to modify */
703 Angle refLat; /* Reference latitude */
704 Angle refLon; /* Reference longitude */
705 {
706 /*
707 * Notation from Snyder, p. 105
708 */
709
710 double phi0, n, tan_n, RF;
711 GeoProjLambertConfConicParams *params;
712
713 params = (GeoProjLambertConfConicParams *)
714 CKALLOC(sizeof(GeoProjLambertConfConicParams));
715 refLat = DomainLat(refLat);
716 if (AngleCmp(refLat, 0) == 0) {
717 SetMercator(projPtr, refLat);
718 return;
719 }
720 phi0 = AngleToRad(refLat);
721 n = sin(phi0);
722 tan_n = pow(tan(M_PI_4 + 0.5 * phi0), n);
723 RF = REarth() * cos(phi0) * tan_n / n;
724
725 projPtr->type = LambertConfConic;
726 params->refLat = refLat;
727 params->refLon = GwchLon(refLon);
728 params->n = n;
729 params->RF = RF;
730 params->rho0 = REarth() / tan(phi0);
731 if (projPtr->params) {
732 CKFREE(projPtr->params);
733 }
734 projPtr->params = params;
735 projPtr->projToLatLonProc = lambrtConfConicToLatLon;
736 projPtr->latLonToProjProc = latLonToLambertConfConic;
737 sprintf(projPtr->descriptor, "LambertConfConic {%9.3f %-9.3f}",
738 AngleToDeg(params->refLat), AngleToDeg(params->refLon));
739 projPtr->infoProc = lambrtConfConicInfo;
740 }
741
742 /*
743 *----------------------------------------------------------------------
744 *
745 * lambrtConfConicInfo --
746 *
747 * This is a value for an infoProc member in a GeoProj structure.
748 * It is usually called by the GeoProjGetInfo procedure
749 *
750 * Results:
751 * Projection information is copied into a GeoProjInfo structure.
752 *
753 * Side effects:
754 * None.
755 *
756 *----------------------------------------------------------------------
757 */
758
759 static GeoProjInfo
lambrtConfConicInfo(projPtr)760 lambrtConfConicInfo(projPtr)
761 struct GeoProj *projPtr;
762 {
763 GeoProjInfo geoProjInfo;
764 GeoProjLambertConfConicParams *params = projPtr->params;
765
766 geoProjInfo.type = projPtr->type;
767 geoProjInfo.prm.refPt.lat = params->refLat;
768 geoProjInfo.prm.refPt.lon = params->refLon;
769 geoProjInfo.rotation = projPtr->rotation;
770 return geoProjInfo;
771 }
772
773 /*
774 *----------------------------------------------------------------------
775 *
776 * latLonToLambertConfConic --
777 *
778 * This is a value for a latLonToProjProc member of a GeoProj structure.
779 * It is usually called by the LatLonToProj procedure.
780 *
781 * Results:
782 * A map-point corresponding to a geopoint.
783 *
784 * Side effects:
785 * None.
786 *
787 *----------------------------------------------------------------------
788 */
789
790 static MapPt
latLonToLambertConfConic(geoPt,projPtr)791 latLonToLambertConfConic(geoPt, projPtr)
792 GeoPt geoPt;
793 struct GeoProj *projPtr;
794 {
795 /*
796 * Formulas and notation from Snyder, p. 105
797 */
798
799 MapPt mapPt;
800 double theta, rho;
801 GeoProjLambertConfConicParams *params = projPtr->params;
802 double lat, lon;
803 Angle d90 = AngleFmDeg(90.0);
804
805 GeoPtGetRad(DomainLonPt(geoPt, params->refLon), &lat, &lon);
806
807 if (AngleCmp(geoPt.lat, d90) == 0) {
808 /*
809 * Point at North Pole
810 */
811
812 if (AngleCmp(params->refLat, 0) == South) {
813 return MapPtNowhere();
814 }
815 rho = 0.0;
816 } else if (AngleCmp(geoPt.lat, -d90) == 0) {
817 /*
818 * Point at South Pole
819 */
820
821 if (AngleCmp(params->refLat, 0) == North) {
822 return MapPtNowhere();
823 }
824 rho = 0.0;
825 } else {
826 /*
827 * Point between poles
828 */
829
830 rho = params->RF / pow(tan(M_PI_4 + 0.5 * lat), params->n);
831 }
832
833 theta = params->n * (lon - AngleToRad(params->refLon));
834 mapPt.abs = rho * sin(theta);
835 mapPt.ord = params->rho0 - rho * cos(theta);
836 return mapPt;
837 }
838
839 /*
840 *----------------------------------------------------------------------
841 *
842 * lambrtConfConicToLatLon --
843 *
844 * This is a value for a projToLatLonProc member of a GeoProj structure.
845 * It is usually called by the ProjToLatLon procedure.
846 *
847 * Results:
848 * A geopoint corresponding to a map-point.
849 *
850 * Side effects:
851 * None.
852 *
853 *----------------------------------------------------------------------
854 */
855
856 static GeoPt
lambrtConfConicToLatLon(mapPt,projPtr)857 lambrtConfConicToLatLon(mapPt, projPtr)
858 MapPt mapPt;
859 struct GeoProj *projPtr;
860 {
861 /*
862 * Formulas and notation from Snyder, p. 105
863 */
864
865 double lat, lon; /* Latitude and longitude in radians */
866 GeoProjLambertConfConicParams *params = projPtr->params;
867 double RF = params->RF;
868 double n = params->n;
869 double x = n > 0.0 ? mapPt.abs : -mapPt.abs;
870 double dy = n > 0.0 ? params->rho0 - mapPt.ord : mapPt.ord - params->rho0;
871 double rho = n > 0.0
872 ? sqrt(mapPt.abs * mapPt.abs + dy * dy)
873 : -sqrt(mapPt.abs * mapPt.abs + dy * dy);
874 double theta = atan2(x, dy);
875
876 lat = 2.0 * atan(pow(RF/rho, 1.0/n)) - M_PI_2;
877 lon = AngleToRad(params->refLon) + theta / n;
878 return DomainLonPt(GeoPtFmRad(lat, lon), params->refLon);
879 }
880
881 /*
882 *----------------------------------------------------------------------
883 *
884 * SetLambertEqArea --
885 *
886 * This procedure sets a projection to Lambert Equal Area.
887 *
888 * Results:
889 * None.
890 *
891 * Side effects:
892 * The members of a GeoProj structure are modified.
893 *
894 *----------------------------------------------------------------------
895 */
896
897 void
SetLambertEqArea(projPtr,refPt)898 SetLambertEqArea(projPtr, refPt)
899 struct GeoProj *projPtr; /* Projection to modify */
900 GeoPt refPt; /* Reference point */
901 {
902 GeoProjRefPtParams *params;
903 double lat, lon;
904 Angle d90 = AngleFmDeg(90.0);
905
906 params = (GeoProjRefPtParams *)CKALLOC(sizeof(GeoProjRefPtParams));
907 projPtr->type = LambertEqArea;
908 params->refPt = GwchLonPt(refPt);
909 GeoPtGetRad(params->refPt, &lat, &lon);
910 params->cosRLat = cos(lat);
911 params->sinRLat = sin(lat);
912 if (projPtr->params) {
913 CKFREE(projPtr->params);
914 }
915 projPtr->params = params;
916 projPtr->projToLatLonProc = lambrtEqAreaToLatLon;
917 projPtr->latLonToProjProc = latLonToLambrtEqArea;
918 if (AngleCmp(refPt.lat, d90) == 0) {
919 sprintf(projPtr->descriptor, "LambertEqArea {90.0 0.0}");
920 } else if (AngleCmp(refPt.lat, -d90) == 0) {
921 sprintf(projPtr->descriptor, "LambertEqArea {-90.0 0.0}");
922 } else {
923 sprintf(projPtr->descriptor, "LambertEqArea {%9.3f %-9.3f}",
924 AngleToDeg(params->refPt.lat), AngleToDeg(params->refPt.lon));
925 }
926 projPtr->infoProc = refPtProjInfo;
927 }
928
929 /*
930 *----------------------------------------------------------------------
931 *
932 * latLonToLambrtEqArea --
933 *
934 * This is a value for a latLonToProjProc member of a GeoProj structure.
935 * It is usually called by the LatLonToProj procedure.
936 *
937 * Results:
938 * A map-point corresponding to a geopoint.
939 *
940 * Side effects:
941 * None.
942 *
943 *----------------------------------------------------------------------
944 */
945
946 static MapPt
latLonToLambrtEqArea(geoPt,projPtr)947 latLonToLambrtEqArea(geoPt, projPtr)
948 GeoPt geoPt;
949 struct GeoProj *projPtr;
950 {
951 MapPt mapPt;
952 GeoProjRefPtParams *params = projPtr->params;
953 double lat, lon;
954 double cosRLat = params->cosRLat, sinRLat = params->sinRLat;
955 double dlon;
956 double k;
957 Angle dist;
958 double cosLat, sinLat, cosDLon;
959 Angle d90 = AngleFmDeg(90.0);
960
961 GeoPtGetRad(geoPt, &lat, &lon);
962 cosLat = cos(lat);
963 sinLat = sin(lat);
964 dlon = AngleToRad(geoPt.lon - params->refPt.lon);
965 cosDLon = cos(dlon);
966 dist = GeoDistance(params->refPt, geoPt);
967 if (AngleCmp(dist, d90) == 1) {
968 /*
969 * Point is further than 90 deg from reference point.
970 */
971
972 return MapPtNowhere();
973 }
974 k = 1.0 + sinRLat * sinLat + cosRLat * cosLat * cosDLon;
975 k = sqrt(2.0 / k);
976 mapPt.abs = REarth() * k * cosLat * sin(dlon);
977 mapPt.ord = REarth() * k * (cosRLat * sinLat - sinRLat * cosLat * cosDLon);
978 return mapPt;
979 }
980
981 /*
982 *----------------------------------------------------------------------
983 *
984 * lambrtEqAreaToLatLon --
985 *
986 * This is a value for a projToLatLonProc member of a GeoProj structure.
987 * It is usually called by the ProjToLatLon procedure.
988 *
989 * Results:
990 * A geopoint corresponding to a map-point.
991 *
992 * Side effects:
993 * None.
994 *
995 *----------------------------------------------------------------------
996 */
997
998 static GeoPt
lambrtEqAreaToLatLon(mapPt,projPtr)999 lambrtEqAreaToLatLon(mapPt, projPtr)
1000 MapPt mapPt;
1001 struct GeoProj *projPtr;
1002 {
1003 GeoProjRefPtParams *params = projPtr->params;
1004 double rho, c, cosC, sinC, ord;
1005 double cosRLat = params->cosRLat, sinRLat = params->sinRLat;
1006 double lat, lon;
1007
1008 rho = sqrt(mapPt.abs * mapPt.abs + mapPt.ord * mapPt.ord);
1009 if (rho == 0.0) {
1010 return params->refPt;
1011 }
1012 if (rho > 2.0 * REarth()) {
1013 return GeoPtNowhere();
1014 }
1015 c = 2.0 * asin(rho / (2.0 * REarth()));
1016 cosC = cos(c);
1017 sinC = sin(c);
1018 ord = cosC * sinRLat + (mapPt.ord * sinC * cosRLat / rho);
1019 if (ord > 1.0) {
1020 return GeoPtNowhere();
1021 }
1022 lat = asin(ord);
1023 lon = AngleToRad(params->refPt.lon)
1024 + atan2(mapPt.abs * sinC,
1025 (rho * cosRLat * cosC - mapPt.ord * sinRLat * sinC));
1026 return DomainLonPt(GeoPtFmRad(lat, lon), params->refPt.lon);
1027 }
1028
1029 /*
1030 *----------------------------------------------------------------------
1031 *
1032 * SetStereographic --
1033 *
1034 * This procedure sets a projection to Stereographic.
1035 *
1036 * Results:
1037 * None.
1038 *
1039 * Side effects:
1040 * The members of a GeoProj structure are modified.
1041 *
1042 *----------------------------------------------------------------------
1043 */
1044
1045 void
SetStereographic(projPtr,refPt)1046 SetStereographic(projPtr, refPt)
1047 struct GeoProj *projPtr; /* Projection to modify */
1048 GeoPt refPt; /* Reference point */
1049 {
1050 GeoProjRefPtParams *params;
1051 double lat, lon;
1052 Angle d90 = AngleFmDeg(90.0);
1053
1054 params = (GeoProjRefPtParams *)CKALLOC(sizeof(GeoProjRefPtParams));
1055 projPtr->type = Stereographic;
1056 params->refPt = GwchLonPt(refPt);
1057 GeoPtGetRad(params->refPt, &lat, &lon);
1058 params->cosRLat = cos(lat);
1059 params->sinRLat = sin(lat);
1060 if (projPtr->params) {
1061 CKFREE(projPtr->params);
1062 }
1063 projPtr->params = params;
1064 projPtr->latLonToProjProc = latLonToStereographic;
1065 projPtr->projToLatLonProc = stereographicToLatLon;
1066 if (LatCmp(refPt.lat, d90) == Eq) {
1067 sprintf(projPtr->descriptor, "Stereographic {90.0 0.0}");
1068 } else if (LatCmp(refPt.lat, -d90) == Eq) {
1069 sprintf(projPtr->descriptor, "Stereographic {-90.0 0.0}");
1070 } else {
1071 sprintf(projPtr->descriptor, "Stereographic {%9.3f %-9.3f}",
1072 AngleToDeg(params->refPt.lat), AngleToDeg(params->refPt.lon));
1073 }
1074 projPtr->infoProc = refPtProjInfo;
1075 }
1076
1077 /*
1078 *----------------------------------------------------------------------
1079 *
1080 * latLonToStereographic --
1081 *
1082 * This is a value for a latLonToProjProc member of a GeoProj structure.
1083 * It is usually called by the LatLonToProj procedure.
1084 *
1085 * Results:
1086 * A map-point corresponding to a geopoint.
1087 *
1088 * Side effects:
1089 * None.
1090 *
1091 *----------------------------------------------------------------------
1092 */
1093
1094 static MapPt
latLonToStereographic(geoPt,projPtr)1095 latLonToStereographic(geoPt, projPtr)
1096 GeoPt geoPt;
1097 struct GeoProj *projPtr;
1098 {
1099 MapPt mapPt;
1100 GeoProjRefPtParams *params = projPtr->params;
1101 GeoPt refPt = params->refPt;
1102 double cosRLat = params->cosRLat, sinRLat = params->sinRLat;
1103 double lat, lon;
1104 double dlon, cosDLon;
1105 double k;
1106 Angle dist;
1107 double cosLat, sinLat;
1108 Angle d90 = AngleFmDeg(90.0);
1109
1110 /*
1111 * Go to 90 degrees from reference point (i.e. follow convention and
1112 * treat as hemisphere projection)
1113 */
1114
1115 dist = GeoDistance(refPt, geoPt);
1116 if (AngleCmp(dist, d90) == 1) {
1117 return MapPtNowhere();
1118 }
1119 GeoPtGetRad(geoPt, &lat, &lon);
1120 cosLat = cos(lat);
1121 sinLat = sin(lat);
1122 dlon = AngleToRad(geoPt.lon - params->refPt.lon);
1123 cosDLon = cos(dlon);
1124 k = 2.0 / (1.0 + sinRLat * sinLat + cosRLat * cosLat * cosDLon);
1125 mapPt.abs = REarth() * k * cosLat * sin(dlon);
1126 mapPt.ord = REarth() * k * (cosRLat * sinLat - sinRLat * cosLat * cosDLon);
1127 return mapPt;
1128 }
1129
1130 /*
1131 *----------------------------------------------------------------------
1132 *
1133 * stereographicToLatLon --
1134 *
1135 * This is a value for a projToLatLonProc member of a GeoProj structure.
1136 * It is usually called by the ProjToLatLon procedure.
1137 *
1138 * Results:
1139 * A geopoint corresponding to a map-point.
1140 *
1141 * Side effects:
1142 * None.
1143 *
1144 *----------------------------------------------------------------------
1145 */
1146
1147 static GeoPt
stereographicToLatLon(mapPt,projPtr)1148 stereographicToLatLon(mapPt, projPtr)
1149 MapPt mapPt;
1150 struct GeoProj *projPtr;
1151 {
1152 GeoProjRefPtParams *params = projPtr->params;
1153 double rho, c, cosC, sinC, ord;
1154 double cosRLat = params->cosRLat, sinRLat = params->sinRLat;
1155 double lat, lon;
1156
1157 rho = sqrt(mapPt.abs * mapPt.abs + mapPt.ord * mapPt.ord);
1158 if (rho == 0.0) {
1159 return params->refPt;
1160 }
1161 c = 2.0 * atan2(rho, 2.0 * REarth());
1162 cosC = cos(c);
1163 sinC = sin(c);
1164 ord = cosC * sinRLat + (mapPt.ord * sinC * cosRLat / rho);
1165 if (ord > 1.0) {
1166 return GeoPtNowhere();
1167 }
1168 lat = asin(ord);
1169 lon = AngleToRad(params->refPt.lon) + atan2(mapPt.abs * sinC,
1170 (rho * cosRLat * cosC - mapPt.ord * sinRLat * sinC));
1171 return DomainLonPt(GeoPtFmRad(lat, lon), params->refPt.lon);
1172 }
1173
1174 /*
1175 *----------------------------------------------------------------------
1176 *
1177 * SetOrthographic --
1178 *
1179 * This procedure sets a projection to Orthographic.
1180 *
1181 * Results:
1182 * None.
1183 *
1184 * Side effects:
1185 * The members of a GeoProj structure are modified.
1186 *
1187 *----------------------------------------------------------------------
1188 */
1189
1190 void
SetOrthographic(projPtr,refPt)1191 SetOrthographic(projPtr, refPt)
1192 struct GeoProj *projPtr; /* Projection to modify */
1193 GeoPt refPt; /* Reference point */
1194 {
1195 double lat, lon;
1196 GeoProjRefPtParams *params;
1197
1198 params = (GeoProjRefPtParams *)CKALLOC(sizeof(GeoProjRefPtParams));
1199 projPtr->type = Orthographic;
1200 params->refPt = GwchLonPt(refPt);
1201 GeoPtGetRad(params->refPt, &lat, &lon);
1202 params->cosRLat = cos(lat);
1203 params->sinRLat = sin(lat);
1204 if (projPtr->params) {
1205 CKFREE(projPtr->params);
1206 }
1207 projPtr->params = params;
1208 projPtr->latLonToProjProc = latLonToOrthographic;
1209 projPtr->projToLatLonProc = orthographicToLatLon;
1210 sprintf(projPtr->descriptor, "Orthographic {%9.3f %-9.3f}",
1211 AngleToDeg(params->refPt.lat), AngleToDeg(params->refPt.lon));
1212 projPtr->infoProc = refPtProjInfo;
1213 }
1214
1215 /*
1216 *----------------------------------------------------------------------
1217 *
1218 * latLonToOrthographic --
1219 *
1220 * This is a value for a latLonToProjProc member of a GeoProj structure.
1221 * It is usually called by the LatLonToProj procedure.
1222 *
1223 * Results:
1224 * A map-point corresponding to a geopoint.
1225 *
1226 * Side effects:
1227 * None.
1228 *
1229 *----------------------------------------------------------------------
1230 */
1231
1232 static MapPt
latLonToOrthographic(geoPt,projPtr)1233 latLonToOrthographic(geoPt, projPtr)
1234 GeoPt geoPt;
1235 struct GeoProj *projPtr;
1236 {
1237 MapPt mapPt;
1238 GeoProjRefPtParams *params = projPtr->params;
1239 GeoPt refPt = params->refPt;
1240 Angle dist;
1241 double lat, lon;
1242 double coslat;
1243 double cosRLat = params->cosRLat, sinRLat = params->sinRLat;
1244 double dlon;
1245 Angle d90 = AngleFmDeg(90.0);
1246
1247 dist = GeoDistance(refPt, geoPt);
1248 if (AngleCmp(dist, d90) == 1) {
1249 /*
1250 * Point further than 90 deg from refPt
1251 */
1252
1253 return MapPtNowhere();
1254 }
1255 GeoPtGetRad(geoPt, &lat, &lon);
1256 coslat = cos(lat);
1257 dlon = AngleToRad(geoPt.lon - params->refPt.lon);
1258 mapPt.abs = REarth() * coslat * sin(dlon);
1259 mapPt.ord = REarth() * (cosRLat * sin(lat) - sinRLat * coslat * cos(dlon));
1260 return mapPt;
1261 }
1262
1263 /*
1264 *----------------------------------------------------------------------
1265 *
1266 * orthographicToLatLon --
1267 *
1268 * This is a value for a projToLatLonProc member of a GeoProj structure.
1269 * It is usually called by the ProjToLatLon procedure.
1270 *
1271 * Results:
1272 * A geopoint corresponding to a map-point.
1273 *
1274 * Side effects:
1275 * None.
1276 *
1277 *----------------------------------------------------------------------
1278 */
1279
1280 static GeoPt
orthographicToLatLon(mapPt,projPtr)1281 orthographicToLatLon(mapPt, projPtr)
1282 MapPt mapPt;
1283 struct GeoProj *projPtr;
1284 {
1285 GeoProjRefPtParams *params = projPtr->params;
1286 double r; /* Polar radius */
1287 double c; /* Angular distance from geoPt to ref point */
1288 double cosC, sinC, ord; /* Computational constants */
1289 double cosRLat = params->cosRLat, sinRLat = params->sinRLat;
1290 double lat, lon;
1291
1292 r = sqrt(mapPt.abs * mapPt.abs + mapPt.ord * mapPt.ord);
1293 if (r == 0.0) {
1294 return params->refPt;
1295 }
1296 if (r / REarth() > 1.0) {
1297 return GeoPtNowhere();
1298 }
1299 c = asin(r / REarth());
1300 cosC = cos(c);
1301 sinC = sin(c);
1302 ord = cosC * sinRLat + (mapPt.ord * sinC * cosRLat / r);
1303 if (ord > 1.0) {
1304 return GeoPtNowhere();
1305 }
1306 lat = asin(ord);
1307 lon = AngleToRad(params->refPt.lon) + atan2(mapPt.abs * sinC,
1308 (r * cosRLat * cosC - mapPt.ord * sinRLat * sinC));
1309 return DomainLonPt(GeoPtFmRad(lat, lon), params->refPt.lon);
1310 }
1311
1312 /*
1313 *----------------------------------------------------------------------
1314 *
1315 * refLonProjInfo --
1316 *
1317 * This is a value for an infoProc member in a GeoProj structure.
1318 * It is usually called by the GeoProjGetInfo procedure
1319 *
1320 * Results:
1321 * Projection information is copied into a GeoProjInfo structure.
1322 *
1323 * Side effects:
1324 * None.
1325 *
1326 *----------------------------------------------------------------------
1327 */
1328
1329 static GeoProjInfo
refLonProjInfo(projPtr)1330 refLonProjInfo(projPtr)
1331 struct GeoProj *projPtr;
1332 {
1333 GeoProjInfo geoProjInfo;
1334 Angle *refLonPtr = projPtr->params;
1335
1336 geoProjInfo.type = projPtr->type;
1337 geoProjInfo.prm.refLon = *refLonPtr;
1338 geoProjInfo.rotation = projPtr->rotation;
1339 return geoProjInfo;
1340 }
1341
1342 /*
1343 *----------------------------------------------------------------------
1344 *
1345 * refPtProjInfo --
1346 *
1347 * This is a value for an infoProc member in a GeoProj structure.
1348 * It is usually called by the GeoProjGetInfo procedure
1349 *
1350 * Results:
1351 * Projection information is copied into a GeoProjInfo structure.
1352 *
1353 * Side effects:
1354 * None.
1355 *
1356 *----------------------------------------------------------------------
1357 */
1358
1359 static GeoProjInfo
refPtProjInfo(projPtr)1360 refPtProjInfo(projPtr)
1361 struct GeoProj *projPtr;
1362 {
1363 GeoProjRefPtParams *params = projPtr->params;
1364 GeoProjInfo geoProjInfo;
1365
1366 geoProjInfo.type = projPtr->type;
1367 geoProjInfo.prm.refPt = params->refPt;
1368 geoProjInfo.rotation = projPtr->rotation;
1369
1370 return geoProjInfo;
1371 }
1372