1 /*
2 * geoLnArrToMap.c --
3 *
4 * This file defines functions that convert the geopoints in geolinearrays
5 * to map-points.
6 *
7 * Copyright (c) 2004 Gordon D. Carrie. All rights reserved.
8 *
9 * Licensed under the Open Software License version 2.1
10 *
11 * Please address questions and feedback to user0@tkgeomap.org
12 *
13 * @(#) $Id: geoLnArrToMap.c,v 1.7 2007/06/29 15:59:53 tkgeomap Exp $
14 *
15 ********************************************
16 *
17 */
18
19 #include <stdio.h>
20 #include <math.h>
21
22 #include <tcl.h>
23 #include "geography.h"
24 #include "geoLines.h"
25 #include "mapLines.h"
26 #include "geoProj.h"
27 #include "geoLnArrToMap.h"
28
29 /*
30 * Declarations of local functions.
31 */
32
33 static GeoPt noAdjust _ANSI_ARGS_((GeoPt geoPt));
34 static GeoPt offPole _ANSI_ARGS_((GeoPt geoPt));
35 static GeoPt nToEq _ANSI_ARGS_((GeoPt geoPt));
36 static GeoPt sToEq _ANSI_ARGS_((GeoPt geoPt));
37 static MapLnArr geoLnArrToCylMap _ANSI_ARGS_((GeoLnArr geoLnArr,
38 GeoProj proj));
39 static MapLnArr geoLnArrToHemisphere _ANSI_ARGS_((GeoLnArr geoLnArr,
40 GeoProj proj));
41
42 /*
43 *----------------------------------------------------------------------
44 *
45 * GeoLnArrToMap --
46 *
47 * This procedure converts the geopoints in a geolinearray to map-points.
48 *
49 * Results:
50 * Return value is a new MapLnArr or NULL if something goes wrong.
51 *
52 * Side effects:
53 * Storage for the returned array is dynamically allocated and should
54 * eventually be freed with a call to MapLnArrDestroy;
55 *
56 *----------------------------------------------------------------------
57 */
58
59 MapLnArr
GeoLnArrToMap(geoLnArr,proj)60 GeoLnArrToMap(geoLnArr, proj)
61 GeoLnArr geoLnArr;
62 GeoProj proj;
63 {
64 GeoProjInfo projInfo; /* GeoProj information */
65 MapLnArr mapLnArr; /* Return value */
66
67 if (!proj) {
68 return NULL;
69 }
70 projInfo = GeoProjGetInfo(proj);
71 switch (projInfo.type) {
72 case CylEqDist:
73 case LambertConfConic:
74 case CylEqArea:
75 case Mercator:
76 mapLnArr = geoLnArrToCylMap(geoLnArr, proj);
77 break;
78 case Orthographic:
79 case Stereographic:
80 case LambertEqArea:
81 mapLnArr = geoLnArrToHemisphere(geoLnArr, proj);
82 break;
83 default:
84 return NULL;
85 }
86 if ( !mapLnArr ) {
87 return NULL;
88 }
89 MapLnArrSet(mapLnArr, geoLnArr, proj);
90 return mapLnArr;
91 }
92
93 /*
94 *----------------------------------------------------------------------
95 *
96 * noAdjust --
97 *
98 * This procedure just returns its input. It is a place holder.
99 *
100 * Results:
101 * The input value.
102 *
103 * Side effects:
104 * None.
105 *
106 *----------------------------------------------------------------------
107 */
108
109 static GeoPt
noAdjust(geoPt)110 noAdjust(geoPt)
111 GeoPt geoPt;
112 {
113 return geoPt;
114 }
115
116 /*
117 *----------------------------------------------------------------------
118 *
119 * offPole --
120 *
121 * This procedure moves a point off the North or South pole. This is
122 * necessary for cylindrical projections which are undefined at the poles.
123 *
124 * Results:
125 * If geoPt is not on a pole, return value is geoPt. If geoPt is on a
126 * pole, return value is a point near geoPt.
127 *
128 * Side effects:
129 * None.
130 *
131 *----------------------------------------------------------------------
132 */
133
134 static GeoPt
offPole(geoPt)135 offPole(geoPt)
136 GeoPt geoPt;
137 {
138 Angle nearAng = AngleFmDeg(89.0);
139 Angle d90 = AngleFmDeg(90.0);
140
141 geoPt.lat = (LatCmp(geoPt.lat, d90) == Eq) ? nearAng
142 : (LatCmp(geoPt.lat, -d90) == Eq) ? -nearAng : geoPt.lat;
143 return geoPt;
144 }
145
146 /*
147 *----------------------------------------------------------------------
148 *
149 * sToEq --
150 *
151 * This procedure moves a point in the Southern Hemisphere to the equator.
152 *
153 * Results:
154 * A point with the same longitude as the given geopoint, with a latitude
155 * not in the Southern Hemisphere.
156 *
157 * Side effects:
158 * None.
159 *
160 *----------------------------------------------------------------------
161 */
162
163 static GeoPt
sToEq(geoPt)164 sToEq(geoPt)
165 GeoPt geoPt;
166 {
167 geoPt.lat = (LatCmp(geoPt.lat, 0) == South) ? 0 : geoPt.lat;
168 return geoPt;
169 }
170
171 /*
172 *----------------------------------------------------------------------
173 *
174 * nToEq --
175 *
176 * This procedure moves a point in the Northern Hemisphere to the equator.
177 *
178 * Results:
179 * A point with the same longitude as the given geopoint, with a latitude
180 * not in the Northern Hemisphere.
181 *
182 * Side effects:
183 * None.
184 *
185 *----------------------------------------------------------------------
186 */
187
188 static GeoPt
nToEq(geoPt)189 nToEq(geoPt)
190 GeoPt geoPt;
191 {
192 geoPt.lat = (LatCmp(geoPt.lat, 0) == North) ? 0 : geoPt.lat;
193 return geoPt;
194 }
195
196 /*
197 *----------------------------------------------------------------------
198 *
199 * geoLnArrToCylMap --
200 *
201 * This procedure applies a cylindrical projection to the geopoints in a
202 * geolinearray. Then it rotates the projected points.
203 *
204 * In a cylindrical projection, the lat-lon's are projected onto a cylinder
205 * that wraps around the Earth with a longitude discontinuity 180 degrees
206 * from the reference longitude. Segments that cross the discontinuity are
207 * split.
208 *
209 * Results:
210 * Return value is a new mapLnArr or NULL if something goes wrong.
211 *
212 * Side effects:
213 * Storage for the returned array is dynamically allocated and should
214 * eventually be freed with a call to MapLnArrDestroy;
215 *
216 *----------------------------------------------------------------------
217 */
218
219 static MapLnArr
geoLnArrToCylMap(geoLnArr,proj)220 geoLnArrToCylMap(geoLnArr, proj)
221 GeoLnArr geoLnArr;
222 GeoProj proj;
223 {
224 MapLnArr mapLnArr = NULL; /* Return value */
225 unsigned nl, nls, np; /* Loop parameters */
226 GeoProjInfo projInfo; /* Projection information */
227 GeoPt (*adjust)(GeoPt) = noAdjust; /* Function to bring a point into
228 * projection domain, if necessary */
229 Angle refLon = 0.0; /* Reference longitude */
230 Angle iWest, iEast; /* West and east end of domain */
231 Angle lon0, lon1; /* Longitudes of last and current pt */
232 GeoPt geoPt0, geoPt1; /* Points in the current segment */
233 GeoPt cross; /* Point on domain boundary */
234 enum LonSgn side0, side1; /* Which side of rLon geoPt0 and geoPt1
235 * are on */
236 enum LonSgn side; /* If line has touched rLon, this stores
237 * which side it approached from */
238 GeoLn geoLn = NULL; /* Current line in latlon coordinates */
239 MapLn mapLn[2] = {NULL, NULL}; /* Current line in map coordinates */
240 MapLn currLn; /* Line we are adding points to */
241 int i; /* Counter to toggle mapLn index */
242 Angle d180 = AngleFmDeg(180.0);
243
244 nls = geoLnArr->nLines;
245 if ( !(mapLnArr = MapLnArrCreate(nls))
246 || !(mapLn[0] = MapLnCreate(1000))
247 || !(mapLn[1] = MapLnCreate(1000)) ) {
248 MapLnDestroy(mapLn[0]);
249 MapLnDestroy(mapLn[1]);
250 MapLnArrDestroy(mapLnArr);
251 return NULL;
252 }
253
254 projInfo = GeoProjGetInfo(proj);
255 switch (projInfo.type) {
256 case CylEqDist:
257 refLon = projInfo.prm.refPt.lon;
258 adjust = noAdjust;
259 break;
260 case LambertConfConic:
261 refLon = projInfo.prm.refPt.lon;
262 if (LatCmp(projInfo.prm.refPt.lat, 0) == North) {
263 adjust = sToEq;
264 } else {
265 adjust = nToEq;
266 }
267 break;
268 case CylEqArea:
269 refLon = projInfo.prm.refLon;
270 adjust = noAdjust;
271 break;
272 case Mercator:
273 refLon = projInfo.prm.refLon;
274 adjust = offPole;
275 break;
276 case Orthographic:
277 case Stereographic:
278 case LambertEqArea:
279 return NULL;
280 }
281 iWest = refLon - d180;
282 iEast = refLon + d180;
283 if (refLon > 0) {
284 /*
285 * If reference longitude is in the eastern hemisphere, split lines
286 * that cross the western boundary.
287 */
288
289 for (nl = 0; nl < nls; nl++) {
290 geoLn = GeoLnArrGetLine(geoLnArr, nl);
291 currLn = mapLn[(i = 0)];
292 side = PrMd;
293
294 /*
295 * Process the first segment.
296 */
297
298 geoPt0 = adjust(GeoLnGetPt(geoLn, 0));
299 geoPt1 = adjust(GeoLnGetPt(geoLn, 1));
300 lon0 = DomainLon(geoPt0.lon, iWest);
301 lon1 = DomainLon(geoPt1.lon, iWest);
302 side0 = LonCmp(geoPt0.lon, iWest);
303 side1 = LonCmp(geoPt1.lon, iWest);
304
305 if (side0 == PrMd && side1 == West) {
306 /*
307 * Line is departing from western boundary to the west.
308 * Draw the segment at the east end.
309 */
310
311 geoPt0.lon = iEast;
312 }
313
314 if (side0 != PrMd && side1 == PrMd) {
315 /*
316 * If line has just intersected west, make a note
317 * of which side it is approaching from.
318 */
319
320 side = side0;
321 if (side0 == West) {
322 /*
323 * If line has arrived at west side from beyond, draw
324 * second point at east end.
325 */
326 geoPt1.lon = iEast;
327 }
328 }
329
330 if (LonBtwn(iWest, lon0, lon1)) {
331 /*
332 * Segment is crossing western boundary.
333 */
334
335 cross.lat = geoPt0.lat + (double)(geoPt1.lat - geoPt0.lat)
336 / (lon1 - lon0) * (iWest - lon0);
337 if (LonCmp(lon0, iWest) == West) {
338 /*
339 * Crossing west to east.
340 */
341
342 cross.lon = iEast;
343 MapLnAddPt(LatLonToProj(cross, proj), currLn);
344 currLn = mapLn[++i % 2];
345 cross.lon = iWest;
346 MapLnAddPt(LatLonToProj(cross, proj), currLn);
347 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
348 } else {
349 /*
350 * Crossing east to west.
351 */
352
353 cross.lon = iWest;
354 MapLnAddPt(LatLonToProj(cross, proj), currLn);
355 currLn = mapLn[++i % 2];
356 cross.lon = iEast;
357 MapLnAddPt(LatLonToProj(cross, proj), currLn);
358 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
359 }
360 } else {
361 /*
362 * No crossing => no special treatment.
363 */
364
365 MapLnAddPt(LatLonToProj(geoPt0, proj), currLn);
366 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
367 }
368 geoPt0 = geoPt1;
369
370 for (np = 2; np < geoLn->nPts; np++) {
371 /*
372 * Process the rest of the segments.
373 */
374
375 geoPt1 = adjust(GeoLnGetPt(geoLn, np));
376 lon0 = DomainLon(geoPt0.lon, iWest);
377 lon1 = DomainLon(geoPt1.lon, iWest);
378 side0 = LonCmp(geoPt0.lon, iWest);
379 side1 = LonCmp(geoPt1.lon, iWest);
380
381 if (side0 != PrMd && side1 == PrMd) {
382 /*
383 * If line has just arrived at west, make a note
384 * of which side it is approaching from.
385 */
386
387 side = side0;
388 }
389
390 if (side != PrMd && side0 == PrMd
391 && side1 != PrMd && side1 != side) {
392 /*
393 * Line has intersected west and
394 * then left on the opposite side.
395 */
396
397 currLn = mapLn[++i % 2];
398 geoPt0.lon = (side == East) ? iEast : iWest;
399 MapLnAddPt(LatLonToProj(geoPt0, proj), currLn);
400 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
401 side = PrMd;
402 } else if (side != PrMd && side0 == PrMd
403 && side1 != PrMd && side1 == side) {
404 /*
405 * Line has traveled along rLon and
406 * then left on the same side.
407 */
408
409 if (side == West) {
410 geoPt1.lon = DomainLon(geoPt1.lon, iEast);
411 }
412 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
413 side = PrMd;
414 } else if (LonBtwn(iWest, lon0, lon1)) {
415 /*
416 * Segment is crossing western boundary.
417 */
418
419 cross.lat = geoPt0.lat
420 + (double)(geoPt1.lat - geoPt0.lat) / (lon1 - lon0)
421 * (iWest - lon0);
422 if (LonCmp(lon0, iWest) == West) {
423 /*
424 * Crossing west to east.
425 */
426
427 cross.lon = iEast;
428 MapLnAddPt(LatLonToProj(cross, proj), currLn);
429 currLn = mapLn[++i % 2];
430 cross.lon = iWest;
431 MapLnAddPt(LatLonToProj(cross, proj), currLn);
432 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
433 } else {
434 /*
435 * Crossing east to west.
436 */
437
438 cross.lon = iWest;
439 MapLnAddPt(LatLonToProj(cross, proj), currLn);
440 currLn = mapLn[++i % 2];
441 cross.lon = iEast;
442 MapLnAddPt(LatLonToProj(cross, proj), currLn);
443 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
444 }
445 } else {
446 /*
447 * Segment is not crossing a boundary. Just add the point.
448 */
449
450 if (side == West) {
451 geoPt1.lon = DomainLon(geoPt1.lon, iEast);
452 }
453 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
454 }
455 geoPt0 = geoPt1;
456 }
457
458 /*
459 * If projected rotated line has any points, add it to the
460 * output array.
461 */
462
463 if (mapLn[0]->nPts > 1) {
464 MapLnArrAddLine(mapLn[0], mapLnArr);
465 }
466 if (mapLn[1]->nPts > 1) {
467 MapLnArrAddLine(mapLn[1], mapLnArr);
468 }
469 MapLnClear(mapLn[0]);
470 MapLnClear(mapLn[1]);
471 }
472 } else {
473 /*
474 * Reference longitude <= 0, i.e. reference longitude is in the
475 * western hemisphere. Split lines that cross the eastern boundary.
476 */
477
478 for (nl = 0; nl < nls; nl++) {
479 geoLn = GeoLnArrGetLine(geoLnArr, nl);
480 currLn = mapLn[(i = 0)];
481 side = PrMd;
482
483 /*
484 * Process the first segment.
485 */
486
487 geoPt0 = adjust(GeoLnGetPt(geoLn, 0));
488 geoPt1 = adjust(GeoLnGetPt(geoLn, 1));
489 lon0 = DomainLon(geoPt0.lon, iEast);
490 lon1 = DomainLon(geoPt1.lon, iEast);
491 side0 = LonCmp(geoPt0.lon, iEast);
492 side1 = LonCmp(geoPt1.lon, iEast);
493
494 if (side0 == PrMd && side1 == East) {
495 /*
496 * Line is departing from eastern boundary to the east.
497 * Draw the segment at the west end.
498 */
499
500 geoPt0.lon = iWest;
501 }
502
503 if (side0 != PrMd && side1 == PrMd) {
504 /*
505 * If line has just intersected east, make a note
506 * of which side it is approaching from.
507 */
508
509 side = side0;
510 if (side0 == East) {
511 /*
512 * If line has arrived at east side from beyond, draw
513 * second point at west end.
514 */
515
516 geoPt1.lon = iWest;
517 }
518 }
519
520 if (LonBtwn(iEast, lon0, lon1)) {
521 /*
522 * Segment is crossing eastern boundary.
523 */
524
525 cross.lat = geoPt0.lat + (double)(geoPt1.lat - geoPt0.lat)
526 / (lon1 - lon0) * (iEast - lon0);
527 if (LonCmp(lon0, iEast) == West) {
528 /*
529 * Crossing west to east
530 */
531
532 cross.lon = iEast;
533 MapLnAddPt(LatLonToProj(cross, proj), currLn);
534 cross.lon = iWest;
535 currLn = mapLn[++i % 2];
536 MapLnAddPt(LatLonToProj(cross, proj), currLn);
537 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
538 } else {
539 /*
540 * Crossing east to west
541 */
542
543 cross.lon = iWest;
544 MapLnAddPt(LatLonToProj(cross, proj), currLn);
545 cross.lon = iEast;
546 currLn = mapLn[++i % 2];
547 MapLnAddPt(LatLonToProj(cross, proj), currLn);
548 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
549 }
550 } else {
551 /*
552 * No crossing => no special treatment.
553 */
554
555 MapLnAddPt(LatLonToProj(geoPt0, proj), currLn);
556 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
557 }
558 geoPt0 = geoPt1;
559
560 for (np = 2; np < geoLn->nPts; np++) {
561 /*
562 * Process the rest of the segments.
563 */
564
565 geoPt1 = adjust(GeoLnGetPt(geoLn, np));
566 lon0 = DomainLon(geoPt0.lon, iEast);
567 lon1 = DomainLon(geoPt1.lon, iEast);
568 side0 = LonCmp(geoPt0.lon, iEast);
569 side1 = LonCmp(geoPt1.lon, iEast);
570
571 if (side0 != PrMd && side1 == PrMd) {
572 /*
573 * If line has just intersected east, make a note
574 * of which side it is approaching from.
575 */
576
577 side = side0;
578 }
579
580 if (side != PrMd && side0 == PrMd
581 && side1 != PrMd && side1 != side) {
582 /*
583 * Line has intersected east and
584 * then left on the opposite side.
585 */
586
587 currLn = mapLn[++i % 2];
588 geoPt0.lon = (side == West) ? iWest : iEast;
589 MapLnAddPt(LatLonToProj(geoPt0, proj), currLn);
590 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
591 side = PrMd;
592 } else if (side != PrMd && side0 == PrMd
593 && side1 != PrMd && side1 == side) {
594 /*
595 * Line has intersected east and
596 * then left on the same side.
597 */
598
599 if (side == East) {
600 geoPt1.lon = DomainLon(geoPt1.lon, iWest);
601 }
602 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
603 side = PrMd;
604 } else if (LonBtwn(iEast, lon0, lon1)) {
605 /*
606 * Segment is crossing eastern boundary.
607 */
608
609 cross.lat = geoPt0.lat
610 + (double)(geoPt1.lat - geoPt0.lat)
611 / (lon1 - lon0) * (iEast - lon0);
612 if (LonCmp(lon0, iEast) == West) {
613 /*
614 * Crossing west to east
615 */
616
617 cross.lon = iEast;
618 MapLnAddPt(LatLonToProj(cross, proj), currLn);
619 cross.lon = iWest;
620 currLn = mapLn[++i % 2];
621 MapLnAddPt(LatLonToProj(cross, proj), currLn);
622 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
623 } else {
624 /*
625 * Crossing east to west
626 */
627
628 cross.lon = iWest;
629 MapLnAddPt(LatLonToProj(cross, proj), currLn);
630 cross.lon = iEast;
631 currLn = mapLn[++i % 2];
632 MapLnAddPt(LatLonToProj(cross, proj), currLn);
633 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
634 }
635 } else {
636 /*
637 * Segment is not crossing a boundary. Just add the point.
638 */
639
640 if (side == East) {
641 geoPt1.lon = DomainLon(geoPt1.lon, iWest);
642 }
643 MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
644 }
645 geoPt0 = geoPt1;
646 }
647
648 /*
649 * If projected rotated line has any points, add it to the
650 * output array.
651 */
652
653 if (mapLn[0]->nPts > 1) {
654 MapLnArrAddLine(mapLn[0], mapLnArr);
655 }
656 if (mapLn[1]->nPts > 1) {
657 MapLnArrAddLine(mapLn[1], mapLnArr);
658 }
659 MapLnClear(mapLn[0]);
660 MapLnClear(mapLn[1]);
661 }
662 }
663
664 MapLnDestroy(mapLn[0]);
665 MapLnDestroy(mapLn[1]);
666 MapLnArrSetAlloc(mapLnArr, mapLnArr->nLines);
667 return mapLnArr;
668 }
669
670 /*
671 *----------------------------------------------------------------------
672 *
673 * geoLnArrToHemisphere --
674 *
675 * This procedure applies a hemispheric projection to the geopoints in
676 * a geolinearray. Then it rotates the projected points.
677 *
678 * In a hemispheric projection, points are defined for the hemisphere
679 * centered on the projection's reference point. Segments that cross the
680 * great circle 90 degrees from the reference point are truncated.
681 *
682 * Results:
683 * Return value is a new mapLnArr or NULL if something goes wrong.
684 *
685 * Side effects:
686 * Storage for the returned array is dynamically allocated and should
687 * eventually be freed with a call to MapLnArrDestroy;
688 *
689 *----------------------------------------------------------------------
690 */
691
692 static MapLnArr
geoLnArrToHemisphere(geoLnArr,proj)693 geoLnArrToHemisphere(geoLnArr, proj)
694 GeoLnArr geoLnArr;
695 GeoProj proj;
696 {
697 MapLnArr mapLnArr = NULL; /* Return value */
698 unsigned nl, nls, np, nps; /* Loop parameters */
699 GeoProjInfo projInfo; /* GeoProj info */
700 GeoPt refPt; /* Reference point */
701 GeoPt ccl0, ccl1; /* Two points on domain boundary */
702 GeoPt geoPt0, geoPt1; /* Points in the current segment */
703 GeoPt bndPt, exitPt, entr0; /* Points on domain boundary */
704 int in0, in1 = 0; /* True if point in domain */
705 int in00; /* True if 1st point in line in domain */
706 Angle az0, az1, az; /* Arc angle and increment */
707 Angle dAz = AngleFmDeg(2.0);/* Angle increment in loop */
708 MapLn mapLn = NULL; /* Line being made right now */
709 MapPt mapPt; /* Point to add to mapLn */
710 Angle d90 = AngleFmDeg(90.0);
711
712 /*
713 * Get reference point from projection
714 */
715
716 projInfo = GeoProjGetInfo(proj);
717 refPt = projInfo.prm.refPt;
718 ccl0 = GeoStep(refPt, 0, d90);
719 ccl1 = GeoStep(refPt, d90, d90);
720
721 nls = geoLnArr->nLines;
722 if (!(mapLnArr = MapLnArrCreate(nls)) || !(mapLn = MapLnCreate(1000))) {
723 MapLnDestroy(mapLn);
724 MapLnArrDestroy(mapLnArr);
725 return NULL;
726 }
727
728 for (nl = 0; nl < nls; nl++) {
729 /*
730 * Convert a line from the array to a rotated projection coordinates
731 */
732
733 GeoLn geoPts = GeoLnArrGetLine(geoLnArr, nl);
734 nps = geoPts->nPts;
735 entr0 = exitPt = GeoPtNowhere();
736
737 /*
738 * Handle first point in line
739 */
740
741 geoPt0 = GeoLnGetPt(geoPts, 0);
742 in0 = in00 = (AngleCmp(GeoDistance(refPt, geoPt0), d90) == -1);
743 if (in0) {
744 mapPt = LatLonToProj(geoPt0, proj);
745 MapLnAddPt(mapPt, mapLn);
746 }
747
748 /*
749 * Loop through the rest of the points
750 */
751
752 for (np = 1; np < nps; np++) {
753 geoPt1 = GeoLnGetPt(geoPts, np);
754 in1 = (AngleCmp(GeoDistance(refPt, geoPt1), d90) == -1);
755 if (in0 && in1) {
756 /*
757 * Segment is in projection domain
758 */
759
760 mapPt = LatLonToProj(geoPt1, proj);
761 MapLnAddPt(mapPt, mapLn);
762 } else if (in0 && !in1) {
763 /*
764 * Segment is leaving projection domain.
765 */
766
767 bndPt = GCircleX(ccl0, ccl1, geoPt0, geoPt1);
768 az = Azimuth(bndPt, refPt);
769 bndPt = GeoStep(bndPt, az, 100);
770 exitPt = bndPt;
771 mapPt = LatLonToProj(bndPt, proj);
772 MapLnAddPt(mapPt, mapLn);
773 } else if (!in0 && in1) {
774 /*
775 * Segment is entering projection domain.
776 */
777
778 bndPt = GCircleX(ccl0, ccl1, geoPt0, geoPt1);
779 az = Azimuth(bndPt, refPt);
780 bndPt = GeoStep(bndPt, az, 100);
781
782 /*
783 * If entering for first time, note location so we can close
784 * the polygon when done
785 */
786
787 if ( !in00 && GeoPtIsNowhere(entr0) ) {
788 entr0 = bndPt;
789 }
790
791 /*
792 * If line has left domain and is coming back in,
793 * draw arc from exit point to entry point (bndPt)
794 */
795
796 if (GeoPtIsSomewhere(exitPt)) {
797 az0 = Azimuth(refPt, exitPt);
798 az1 = DomainLon(Azimuth(refPt, bndPt), az0);
799 if (az1 > az0) {
800 for (az = az0 + dAz; az < az1; az += dAz) {
801 GeoPt pt = GeoStep(refPt, az, d90 - 100);
802 mapPt = LatLonToProj(pt, proj);
803 MapLnAddPt(mapPt, mapLn);
804 }
805 } else {
806 for (az = az0 - dAz; az > az1; az -= dAz) {
807 GeoPt pt = GeoStep(refPt, az, d90 - 100);
808 mapPt = LatLonToProj(pt, proj);
809 MapLnAddPt(mapPt, mapLn);
810 }
811 }
812 }
813
814 /*
815 * Add entry point and current point
816 */
817
818 mapPt = LatLonToProj(bndPt, proj);
819 MapLnAddPt(mapPt, mapLn);
820 mapPt = LatLonToProj(geoPt1, proj);
821 MapLnAddPt(mapPt, mapLn);
822 }
823
824 geoPt0 = geoPt1;
825 in0 = in1;
826
827 }
828
829 /*
830 * If line started and finished outside domain
831 * draw arc from last exit to first entry to close the polygon
832 */
833
834 if ( !in00 && !in1 && GeoPtIsSomewhere(entr0) ) {
835 az0 = Azimuth(refPt, exitPt);
836 az1 = DomainLon(Azimuth(refPt, entr0), az0);
837 if (az1 > az0) {
838 for (az = az0 + dAz; az < az1; az += dAz) {
839 GeoPt pt = GeoStep(refPt, az, d90 - 100);
840 mapPt = LatLonToProj(pt, proj);
841 MapLnAddPt(mapPt, mapLn);
842 }
843 } else {
844 for (az = az0 - dAz; az > az1; az -= dAz) {
845 GeoPt pt = GeoStep(refPt, az, d90 - 100);
846 mapPt = LatLonToProj(pt, proj);
847 MapLnAddPt(mapPt, mapLn);
848 }
849 }
850 }
851
852 if (mapLn->nPts > 1) {
853 MapLnArrAddLine(mapLn, mapLnArr);
854 }
855 MapLnClear(mapLn);
856
857 }
858
859 MapLnDestroy(mapLn);
860 return mapLnArr;
861 }
862