-
Notifications
You must be signed in to change notification settings - Fork 359
/
gmt_map.c
10249 lines (9040 loc) · 472 KB
/
gmt_map.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*--------------------------------------------------------------------
*
* Copyright (c) 1991-2024 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
* See LICENSE.TXT file for copying and redistribution conditions.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation; version 3 or any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* Contact info: www.generic-mapping-tools.org
*--------------------------------------------------------------------*/
/*
* G M T _ M A P . C
*
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
* GMT_map.c contains code related to generic coordinate transformation.
* For the actual projection functions, see gmt_proj.c
*
* Map Transformation Setup Routines
* These routines initializes the selected map transformation
* The names and main function are listed below
* NB! Note that the transformation function does not check that they are
* passed valid lon,lat numbers. I.e asking for log10 scaling using values
* <= 0 results in problems.
*
* The ellipsoid used is selectable by editing the gmt.conf in your
* home directory. If no such file, create one by running gmtdefaults.
*
* Usage: Initialize system by calling gmt_g (separate module), and
* then just use gmt_geo_to_xy() and gmt_xy_to_geo() functions.
*
* Author: Paul Wessel
* Date: 1-JAN-2010
* Version: 5
*/
/*!
* \file gmt_map.c
* \brief gmt_map.c contains code related to generic coordinate transformation.
*
* A) List of exported gmt_* functions available to modules and libraries via gmt_dev.h:
*
* gmt_ECEF_forward
* gmt_ECEF_init
* gmt_ECEF_inverse
* gmt_ECEF_inverse_dest_datum
* gmt_UTMzone_to_wesn
* gmt_auto_frame_interval
* gmt_az_backaz
* gmt_azim_to_angle
* gmt_cart_outside
* gmt_cart_to_xy_line
* gmt_circle_to_region
* gmt_clip_to_map
* gmt_compact_line
* gmt_conv_datum
* gmt_datum_init
* gmt_dist_array
* gmt_dist_array_2
* gmt_distance
* gmt_geo_to_xy
* gmt_geo_to_xy_line
* gmt_geo_to_xy_noshift
* gmt_geo_to_xy_noshiftscale
* gmt_geoz_to_xy
* gmt_get_az_dist_from_components
* gmt_get_geo_ellipse
* gmt_get_smallcircle
* gmt_graticule_path
* gmt_grd_project
* gmt_great_circle_dist_meter
* gmt_half_map_width
* gmt_img_project
* gmt_init_distaz
* gmt_lat_swap
* gmt_line_length
* gmt_map_clip_path
* gmt_map_outside
* gmt_map_perimeter_search
* gmt_map_setup
* gmt_map_truncate
* gmt_mindist_to_point
* gmt_near_a_line
* gmt_near_a_point
* gmt_near_lines
* gmt_proj_setup
* gmt_project_init
* gmt_segment_BB_outside_map_BB
* gmt_set_datum
* gmt_split_poly_at_dateline
* gmt_translate_point
* gmt_wesn_search
* gmt_x_to_xx
* gmt_xy_to_geo
* gmt_xy_to_geo_noshift
* gmt_xy_to_geo_noshiftscale
* gmt_xyz_to_xy
* gmt_y_to_yy
* gmt_z_to_zz
*
* B) List of exported gmtlib_* functions available to libraries via gmt_internals.h:
*
* gmtlib_adjust_we_if_central_lon_set
* gmtlib_cartesian_dist
* gmtlib_cartesian_dist_proj
* gmtlib_distance_type
* gmtlib_genper_reset
* gmtlib_get_point_from_r_az
* gmtlib_great_circle_dist_cos
* gmtlib_great_circle_dist_degree
* gmtlib_great_circle_intersection
* gmtlib_init_ellipsoid
* gmtlib_init_geodesic
* gmtlib_lat_swap_quick
* gmtlib_latpath
* gmtlib_left_boundary
* gmtlib_lonpath
* gmtlib_map_latcross
* gmtlib_map_loncross
* gmtlib_map_path
* gmtlib_right_boundary
* gmtlib_scale_eqrad
* gmtlib_set_oblique_pole_and_origin
* gmtlib_small_circle_intersection
* gmtlib_split_line
* gmtlib_translate_point
*/
#include "gmt_dev.h"
#include "gmt_internals.h"
/* We put all the declaration of external functions from gmt_proj.h here since
* (a) they are only called in gmt_map.c so no need to go in gmt_internals.h
* (b) this will change completely when proj4 is relacing gmt_proj in GMT 6.
*/
#include "gmt_proj.c"
/*! CCW order of side in some tests */
enum GMT_side {
GMT_BOTTOM = 0,
GMT_RIGHT = 1,
GMT_TOP = 2,
GMT_LEFT = 3};
/* Basic error reporting when things go badly wrong. This Error_and_return macro can be
* used in stead of regular return(code) to print out what the code is before
* it returns. We assume the GMT pointer is available in the function!
*/
#define Error_and_return(code,err) { GMT_Report(GMT->parent,GMT_MSG_ERROR,"%s\n",gmt_error_string[code]); return (err);}
/* Note by P. Wessel, 18-Oct-2012, updated 08-JAN-2014:
* In the olden days, GMT only did great circle distances. In GMT 4 we implemented geodesic
* distances by Rudoe's formula as given in Bomford [1971]. However, that geodesic is not
* exactly what we wanted as it is a normal section and does not strictly follow the geodesic.
* Other candidates are Vincenty [1975], which is widely used and Karney [2012], which is super-
* accurate. At this point their differences are in the micro-meter level. For GMT 5 we have
* now switched to the Vincenty algorithm as provided by Gerald Evenden, USGS [author of proj4],
* which is a modified translation of the NGS algorithm and not exactly what is in proj4's geod
* program (which Evenden thinks is inferior.) I ran a comparison between many algorithms that
* either were available via codes or had online calculators. I sought the geodesic distance
* from (0,0) to (10,10) on WGS-84; the results were (in meters):
*
* GMT4 (Rudoe): 1565109.099232116
* proj4: 1565109.095557918
* vdist(0,0,10,10) [0] 1565109.09921775
* Karney [1]: 1565109.09921789
* Vincenty [2]: 1565109.099218036
* NGS [3] 1565109.0992
* Andoyer [4] 1565092.276857755
*
* [0] via Joaquim Luis, supposedly Vincenty [2012]
* [1] via online calculator at max precision https://geographiclib.sourceforge.net/cgi-bin/Geod
* [2] downloading, compiling and running https://article.gmane.org/gmane.comp.gis.proj-4.devel/3478.
* This is not identical to Vincenty in proj4 but written by Evenden (proj.4 author)
* [3] via online calculator https://www.ngs.noaa.gov/cgi-bin/Inv_Fwd/inverse2.prl. Their notes says
* this is Vincenty; unfortunately I cannot control the output precision.
* [4] Andoyer approximate from Astronomical Algorithms, Jean Meeus, 2009, second edition, Willmann-Bell, Inc.
*
* Based on these comparisons we decided to implement the Vincenty [2] code as given. The older Rudoe
* code is also accessible, as is the approximation by Andoyer which is good to a few tens of m.
* The choice was based on the readily available C code versus having to reimplement Karney in C.
* When GMT 6 is started we expect to use proj4
*/
static char *GEOD_TEXT[3] = {"Vincenty", "Andoyer", "Rudoe"};
/*! . */
GMT_LOCAL double gmtmap_get_angle (struct GMT_CTRL *GMT, double lon1, double lat1, double lon2, double lat2) {
double x1, y1, x2, y2, angle, direction;
gmt_geo_to_xy (GMT, lon1, lat1, &x1, &y1);
gmt_geo_to_xy (GMT, lon2, lat2, &x2, &y2);
if (doubleAlmostEqualZero (y1, y2) && doubleAlmostEqualZero (x1, x2)) { /* Special case that only(?) occurs at N or S pole or r=0 for GMT_POLAR */
if (fabs (fmod (lon1 - GMT->common.R.wesn[XLO] + 360.0, 360.0)) > fabs (fmod (lon1 - GMT->common.R.wesn[XHI] + 360.0, 360.0))) { /* East */
gmt_geo_to_xy (GMT, GMT->common.R.wesn[XHI], GMT->common.R.wesn[YLO], &x1, &y1);
gmt_geo_to_xy (GMT, GMT->common.R.wesn[XHI], GMT->common.R.wesn[YHI], &x2, &y2);
GMT->current.map.corner = 1;
}
else {
gmt_geo_to_xy (GMT, GMT->common.R.wesn[XLO], GMT->common.R.wesn[YLO], &x1, &y1);
gmt_geo_to_xy (GMT, GMT->common.R.wesn[XLO], GMT->common.R.wesn[YHI], &x2, &y2);
GMT->current.map.corner = 3;
}
angle = d_atan2d (y2-y1, x2-x1) - 90.0;
if (GMT->current.proj.got_azimuths) angle += 180.0;
if (GMT->current.proj.flip) angle += 180.0;
}
else
angle = d_atan2d (y2 - y1, x2 - x1);
if (abs (GMT->current.map.prev_x_status) == 2 && abs (GMT->current.map.prev_y_status) == 2) /* Last point outside */
direction = angle + 180.0;
else if (GMT->current.map.prev_x_status == 0 && GMT->current.map.prev_y_status == 0) /* Last point inside */
direction = angle;
else {
if (abs (GMT->current.map.this_x_status) == 2 && abs (GMT->current.map.this_y_status) == 2) /* This point outside */
direction = angle;
else if (GMT->current.map.this_x_status == 0 && GMT->current.map.this_y_status == 0) /* This point inside */
direction = angle + 180.0;
else { /* Special case of corners and sides only */
if (GMT->current.map.prev_x_status == GMT->current.map.this_x_status)
direction = (GMT->current.map.prev_y_status == 0) ? angle : angle + 180.0;
else if (GMT->current.map.prev_y_status == GMT->current.map.this_y_status)
direction = (GMT->current.map.prev_x_status == 0) ? angle : angle + 180.0;
else
direction = angle;
}
}
if (direction < 0.0) direction += 360.0;
if (direction >= 360.0) direction -= 360.0;
return (direction);
}
int gmtlib_adjust_we_if_central_lon_set (struct GMT_CTRL *GMT, double *west, double *east) {
/* Try to arrange longitudes relative to central meridian if it has been set */
int way = 0; /* Default is no change */
if (gmt_M_is_cartesian (GMT, GMT_IN)) return way; /* Adjustment is only sensible for geographic limits */
if (gmt_M_is_dnan (GMT->current.proj.central_meridian)) return way; /* Not set yet so nothing to consider */
/* Here we may consider various cases - for now just a general case */
if (*west > GMT->current.proj.central_meridian) {
*west -= 360.0;
*east -= 360.0;
way = -1; /* We shifted westwards */
}
else if (*east < GMT->current.proj.central_meridian) {
*west += 360.0;
*east += 360.0;
way = +1; /* We shifted eastwards */
}
return (way);
}
/*! . */
double gmtlib_left_boundary (struct GMT_CTRL *GMT, double y) {
return ((*GMT->current.map.left_edge) (GMT, y));
}
/*! . */
double gmtlib_right_boundary (struct GMT_CTRL *GMT, double y) {
return ((*GMT->current.map.right_edge) (GMT, y));
}
/* Private functions internal to gmt_map.c */
/*! . */
GMT_LOCAL bool gmtmap_quickconic (struct GMT_CTRL *GMT) {
/* Returns true if area/scale are large/small enough
* so that we can use spherical equations with authalic
* or conformal latitudes instead of the full ellipsoidal
* equations.
*/
double s, dlon, width;
if (GMT->common.j.active && GMT->common.j.mode == GMT_GEODESIC) return (false); /* Used -je for full ellipsoidal terms */
if (GMT->common.j.active && GMT->common.j.mode == GMT_GREATCIRCLE) return (true); /* Used -jg for spherical mode */
if (GMT->common.R.active[RSET]) { /* Use longitudinal extent to decide best option */
if (GMT->current.proj.gave_map_width) { /* Gave width */
dlon = GMT->common.R.wesn[XHI] - GMT->common.R.wesn[XLO];
width = GMT->current.proj.pars[4] * GMT->session.u2u[GMT->current.setting.proj_length_unit][GMT_M]; /* Convert to meters */
s = (dlon * GMT->current.proj.M_PR_DEG) / width;
}
else if (GMT->current.proj.units_pr_degree) { /* Gave scale */
/* Convert to meters */
s = GMT->current.proj.M_PR_DEG / (GMT->current.proj.pars[4] * GMT->session.u2u[GMT->current.setting.proj_length_unit][GMT_M]);
}
else { /* Got 1:xxx that was changed */
s = (1.0 / GMT->current.proj.pars[4]) / GMT->current.proj.unit;
}
if (s > 1.0e7) { /* if s in 1:s exceeds 1e7 we do the quick thing */
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Using spherical projection with conformal latitudes\n");
return (true);
}
else /* Use full ellipsoidal terms */
return (false);
}
return (true); /* Do the fast and easy if neighter -je|g nor -R is given */
}
/*! . */
GMT_LOCAL bool gmtmap_quicktm (struct GMT_CTRL *GMT, double lon0, double limit) {
/* Returns true if the region chosen is too large for the
* ellipsoidal series to be valid; hence use spherical equations
* with conformal latitudes instead.
* We let +-limit degrees from central meridian be the cutoff.
*/
double d_left, d_right;
if (GMT->common.j.active && GMT->common.j.mode == GMT_GEODESIC) return (false); /* Used -je for full ellipsoidal terms */
if (GMT->common.j.active && GMT->common.j.mode == GMT_GREATCIRCLE) return (true); /* Used -jg for spherical mode */
if (GMT->common.R.active[RSET]) { /* Use longitudinal extent to decide best option */
d_left = lon0 - GMT->common.R.wesn[XLO] - 360.0;
d_right = lon0 - GMT->common.R.wesn[XHI] - 360.0;
while (d_left < -180.0) d_left += 360.0;
while (d_right < -180.0) d_right += 360.0;
if (fabs (d_left) > limit || fabs (d_right) > limit) {
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Using spherical projection with conformal latitudes since area too big for ellipsoidal series\n");
return (true);
}
else /* Use full ellipsoidal terms */
return (false);
}
return (true); /* If neither -je|g nor -R is set we just do spherical */
}
/*! . */
GMT_LOCAL void gmtmap_set_polar (struct GMT_CTRL *GMT) {
/* Determines if the projection pole is N or S pole */
if (doubleAlmostEqual (fabs (GMT->current.proj.pars[1]), 90.0)) {
GMT->current.proj.polar = true;
GMT->current.proj.north_pole = (GMT->current.proj.pars[1] > 0.0);
}
}
GMT_LOCAL bool gmtmap_central_meridian_not_set (struct GMT_CTRL *GMT) {
/* Just to make it clearer to understand the code. If NaN then we were never given a central meridian */
return (gmt_M_is_dnan (GMT->current.proj.pars[0]));
}
GMT_LOCAL void gmtmap_set_default_central_meridian (struct GMT_CTRL *GMT) {
/* Pick half-way between w and e, but watch for -R+r */
if (GMT->common.R.oblique && GMT->common.R.wesn[XHI] < GMT->common.R.wesn[XLO])
GMT->current.proj.pars[0] = 0.5 * (GMT->common.R.wesn[XLO] + GMT->common.R.wesn[XHI] + 360.0); /* Set to middle lon but watch for e < w */
else
GMT->current.proj.pars[0] = 0.5 * (GMT->common.R.wesn[XLO] + GMT->common.R.wesn[XHI]); /* Set to middle lon */
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Central meridian not given, default to %g\n", GMT->current.proj.pars[0]);
GMT->current.proj.central_meridian = GMT->current.proj.pars[0];
}
/*! . */
GMT_LOCAL int gmtmap_cyl_validate_clon (struct GMT_CTRL *GMT, unsigned int mode) {
/* Make sure that for global (360-range) cylindrical projections, the central meridian is neither west nor east.
* If so then we reset it to the middle value or we change -R:
* mode == 0: <clon> should be reset based on w/e mid-point
* mode == 1: -J<clon> is firm so w/e is centered on <c.lon>
* mode == 2: -J<clon> is firm and e-w < 360 so must give an error instead.
*/
int error = GMT_NOERROR;
if (gmtmap_central_meridian_not_set (GMT)) /* If not set then we pick halfway between w and e */
gmtmap_set_default_central_meridian (GMT);
if (GMT->current.map.is_world) { /* For full 360 range the central meridian must be in the middle */
double w = GMT->current.proj.pars[0] - 180.0, e = GMT->current.proj.pars[0] + 180.0;
if (!doubleAlmostEqualZero (GMT->common.R.wesn[XLO], w)) { /* Not yet aligned */
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Region for global cylindrical projection had to be reset from %g/%g to %g/%g\n",
GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI], w, e);
GMT->common.R.wesn[XLO] = w; GMT->common.R.wesn[XHI] = e;
}
}
else if (!GMT->common.R.oblique) { /* For regional (<360) areas we cannot have clon > 180 away from either boundary */
double de, dw, clon;
bool is_clon0360 = GMT->current.proj.pars[0] > 180; /* Is the central meridian in the [0 360] range */;
bool is_reg0360 = GMT->common.R.wesn[XLO] > 180 || GMT->common.R.wesn[XHI] > 180; /* Is the region in the [0 360] range */
clon = GMT->current.proj.pars[0];
if (is_clon0360 && !is_reg0360 && clon > 180) /* Must take care of not checking for clon < 180 from either boundary and mix -Rg & -Rd */
clon -= 360.0;
else if (!is_clon0360 && is_reg0360 && clon < 0)
clon += 360.0;
dw = fabs (clon - GMT->common.R.wesn[XLO]); if (dw >= 360) dw -= 360.0;
de = fabs (clon - GMT->common.R.wesn[XHI]); if (de >= 360) de -= 360.0;
if (dw > 180.0 || de > 180.0) {
if (mode == 2) { /* Yield an error if fixed central longitude, range < 360, and exceed 180 to the border from central longitude */
static char *border[2] = {"Western", "Eastern"};
unsigned int kase = (dw > 180.0) ? 0 : 1;
GMT_Report (GMT->parent, GMT_MSG_ERROR, "%s boundary is > 180 degrees from specified central meridian and thus your region is invalid\n", border[kase]);
error = GMT_MAP_EXCEEDS_360;
}
else { /* Else we just adapt to the situation */
double new_lon = 0.5 * (GMT->common.R.wesn[XLO] + GMT->common.R.wesn[XHI]);
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Central meridian for cylindrical projection had to be reset from %g to %g\n", GMT->current.proj.pars[0], new_lon);
GMT->current.proj.pars[0] = new_lon;
}
}
}
return (error);
}
/*! . */
GMT_LOCAL void gmtmap_lat_swap_init (struct GMT_CTRL *GMT) {
/* Initialize values in GMT->current.proj.lat_swap_vals based on GMT->current.proj.
First compute GMT->current.proj.lat_swap_vals.ra (and rm), the radii to use in
spherical formulae for area (respectively, N-S distance) when
using the authalic (respectively, meridional) latitude.
Then for each type of swap:
First load GMT->current.proj.lat_swap_vals.c[itype][k], k=0,1,2,3 with the
coefficient for sin(2 (k +1) phi), based on series from Adams.
Next reshuffle these coefficients so they form a nested
polynomial using equations (3-34) and (3-35) on page 19 of
Snyder.
References: J. P. Snyder, "Map projections - a working manual",
U. S. Geological Survey Professional Paper #1395, 1987.
O. S. Adams, "Latitude Developments Connected With Geodesy and
Cartography", U. S. Coast and Geodetic Survey Special Publication
number 67, 1949.
P. D. Thomas, "Conformal Projections in Geodesy and Cartography",
US CGS Special Pub #251, 1952.
See also other US CGS Special Pubs (#53, 57, 68, 193, and 251).
Latitudes are named as follows (this only partly conforms to
names in the literature, which are varied):
Geodetic = G, angle between ellipsoid normal and equator
geocentric = O, angle between radius from Origin and equator
Parametric = P, angle such that x=a*cos(phi), y=b*sin(phi) is on ellipse
Authalic = A, angle to use in equal area development of ellipsoid
Conformal = GMT, angle to use in conformal development of ellipsoid
Meridional = M, angle to use in N-S distance calculation of ellipsoid
(The parametric latitude is the one used in orthogonal curvilinear
coordinates and ellipsoidal harmonics. The term "authalic" was coined
by A. Tissot in "Memoire sur la Representations des Surfaces et les
projections des cartes geographiques", Gauthier Villars, Paris, 1881;
it comes from the Greek meaning equal area.)
The idea of latitude swaps is this: Conformal, equal-area, and other
developments of spherical surfaces usually lead to analytic formulae
for the forward and inverse projections which are stable over a wide
range of the values. It is handy to use the same formulae when
developing the surface of the ellipsoid. The authalic (respectively,
conformal) lat is such that when plugged into a spherical development
formula, it results in an authalic (meaning equal area) (respectively,
conformal) development of the ellipsoid. The meridional lat does the
same thing for measurement of N-S distances along a meridian.
Adams gives coefficients for series in sin(2 (k +1) phi), for k
up to 2 or 3. I have extended these to k=3 in all cases except
the authalic. I have sometimes multiplied his coefficients by
(-1) so that the sense here is always to give a correction to be
added to the input lat to get the output lat in gmt_lat_swap().
I have tested this code by checking that
fabs (geocentric) < fabs (parametric) < fabs (geodetic)
and also, for each pair of possible conversions, that the
forward followed by the inverse returns the original lat to
within a small tolerance. This tolerance is as follows:
geodetic <-> authalic: max error (degrees) = 1.253344e-08
geodetic <-> conformal: max error (degrees) = 2.321796e-07 should be better after 13nov07 fix
geodetic <-> meridional: max error (degrees) = 4.490630e-12
geodetic <-> geocentric: max error (degrees) = 1.350031e-13
geodetic <-> parametric: max error (degrees) = 1.421085e-14
geocentric <-> parametric: max error (degrees) = 1.421085e-14
Currently, (GMT v5) the only ones we anticipate using are
geodetic, authalic, and conformal. I have put others in here
for possible future convenience.
Also, I made this depend on GMT->current.setting.ref_ellipsoid[GMT->current.setting.proj_ellipsoid]
rather than on GMT->current.proj, so that it will be possible to
call gmt_lat_swap() without having to pass -R and -J to
gmt_map_setup(), so that in the future we will be able to use
lat conversions without plotting maps.
W H F Smith, 10--13 May 1999. */
/* PW notes: Projections only convert latitudes if GMT->current.proj.GMT_convert_latitudes is true.
* This is set by gmt_map_setup if the ellipsoid is not a sphere. Calling gmtmap_lat_swap_init by itself
* does not affect the mapping machinery. Since various situations call for the use
* of auxiliary latitudes we initialize gmtmap_lat_swap_init in gmt_begin. This means
* programs can use functions like gmt_lat_swap whenever needed.
*/
unsigned int i;
double x, xx[4], a, f, e2, e4, e6, e8;
f = GMT->current.setting.ref_ellipsoid[GMT->current.setting.proj_ellipsoid].flattening;
a = GMT->current.setting.ref_ellipsoid[GMT->current.setting.proj_ellipsoid].eq_radius;
if (gmt_M_is_zero (f)) {
gmt_M_memset (GMT->current.proj.lat_swap_vals.c, GMT_LATSWAP_N * 4, double);
GMT->current.proj.lat_swap_vals.ra = GMT->current.proj.lat_swap_vals.rm = a;
GMT->current.proj.lat_swap_vals.spherical = true;
return;
}
GMT->current.proj.lat_swap_vals.spherical = false;
/* Below are two sums for x to get the two radii. I have nested the
parentheses to add the terms in the order that would minimize roundoff
error. However, in double precision there may be no need to do this.
I have carried these to 4 terms (eccentricity to the 8th power) because
this is as far as Adams goes with anything, but it is not clear what
the truncation error is, since every term in the sum has the same sign. */
e2 = f * (2.0 - f);
e4 = e2 * e2;
e6 = e4 * e2;
e8 = e4 * e4;
/* This expression for the Authalic radius comes from Adams [1949] */
xx[0] = 2.0 / 3.0;
xx[1] = 3.0 / 5.0;
xx[2] = 4.0 / 7.0;
xx[3] = 5.0 / 9.0;
x = xx[0] * e2 + ( xx[1] * e4 + ( xx[2] * e6 + xx[3] * e8));
GMT->current.proj.lat_swap_vals.ra = a * sqrt( (1.0 + x) * (1.0 - e2));
/* This expression for the Meridional radius comes from Gradshteyn and Ryzhik, 8.114.1,
because Adams only gets the first two terms. This can be worked out by expressing the
meridian arc length in terms of an integral in parametric latitude, which reduces to
equatorial radius times Elliptic Integral of the Second Kind. Expanding this using
binomial theorem leads to Gradshteyn and Ryzhik's expression: */
xx[0] = 1.0 / 4.0;
xx[1] = xx[0] * 3.0 / 16.0;
xx[2] = xx[1] * 3.0 * 5.0 / 36.0;
xx[3] = xx[2] * 5.0 * 7.0 / 64.0;
x = xx[0] * e2 + ( xx[1] * e4 + ( xx[2] * e6 + xx[3] * e8));
GMT->current.proj.lat_swap_vals.rm = a * (1.0 - x);
/* Geodetic to authalic: */
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2A][0] = -(e2 / 3.0 + (31.0 * e4 / 180.0 + 59.0 * e6 / 560.0));
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2A][1] = 17.0 * e4 / 360.0 + 61.0 * e6 / 1260;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2A][2] = -383.0 * e6 / 45360.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2A][3] = 0.0;
/* Authalic to geodetic: */
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_A2G][0] = e2 / 3.0 + (31.0 * e4 / 180.0 + 517.0 * e6 / 5040.0);
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_A2G][1] = 23.0 * e4 / 360.0 + 251.0 * e6 / 3780;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_A2G][2] = 761.0 * e6 / 45360.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_A2G][3] = 0.0;
/* Geodetic to conformal: */
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2C][0] = -(e2 / 2.0 + (5.0 * e4 / 24.0 + (3.0 * e6 / 32.0 + 281.0 * e8 / 5760.0)));
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2C][1] = 5.0 * e4 / 48.0 + (7.0 * e6 / 80.0 + 697.0 * e8 / 11520.0);
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2C][2] = -(13.0 * e6 / 480.0 + 461.0 * e8 / 13440.0);
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2C][3] = 1237.0 * e8 / 161280.0;
/* Conformal to geodetic: */
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_C2G][0] = e2 / 2.0 + (5.0 * e4 / 24.0 + (e6 / 12.0 + 13.0 * e8 / 360.0)) ;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_C2G][1] = 7.0 * e4 / 48.0 + (29.0 * e6 / 240.0 + 811.0 * e8 / 11520.0);
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_C2G][2] = 7.0 * e6 / 120.0 + 81.0 * e8 / 1120.0; /* Bug fixed 13nov07 whfs */
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_C2G][3] = 4279.0 * e8 / 161280.0;
/* The meridional and parametric developments use this parameter: */
x = f/(2.0 - f); /* Adams calls this n. It is f/(2-f), or -betaJK in my notes. */
xx[0] = x; /* n */
xx[1] = x * x; /* n-squared */
xx[2] = xx[1] * x; /* n-cubed */
xx[3] = xx[2] * x; /* n to the 4th */
/* Geodetic to meridional: */
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2M][0] = -(3.0 * xx[0] / 2.0 - 9.0 * xx[2] / 16.0);
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2M][1] = 15.0 * xx[1] / 16.0 - 15.0 * xx[3] / 32.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2M][2] = -35.0 * xx[2] / 48.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2M][3] = 315.0 * xx[3] / 512.0;
/* Meridional to geodetic: */
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_M2G][0] = 3.0 * xx[0] / 2.0 - 27.0 * xx[2] / 32.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_M2G][1] = 21.0 * xx[1] / 16.0 - 55.0 * xx[3] / 32.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_M2G][2] = 151.0 * xx[2] / 96.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_M2G][3] = 1097.0 * xx[3] / 512.0;
/* Geodetic to parametric equals parametric to geocentric: */
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2P][0] = GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_P2O][0] = -xx[0];
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2P][1] = GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_P2O][1] = xx[1] / 2.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2P][2] = GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_P2O][2] = -xx[2] / 3.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2P][3] = GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_P2O][3] = xx[3] / 4.0;
/* Parametric to geodetic equals geocentric to parametric: */
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_P2G][0] = GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_O2P][0] = xx[0];
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_P2G][1] = GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_O2P][1] = xx[1] / 2.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_P2G][2] = GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_O2P][2] = xx[2] / 3.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_P2G][3] = GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_O2P][3] = xx[3] / 4.0;
/* The geodetic <->geocentric use this parameter: */
x = 1.0 - e2;
x = (1.0 - x)/(1.0 + x); /* Adams calls this m. It is e2/(2-e2), or -betaJK in my notes. */
xx[0] = x; /* m */
xx[1] = x * x; /* m-squared */
xx[2] = xx[1] * x; /* m-cubed */
xx[3] = xx[2] * x; /* m to the 4th */
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2O][0] = -xx[0];
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2O][1] = xx[1] / 2.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2O][2] = -xx[2] / 3.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_G2O][3] = xx[3] / 4.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_O2G][0] = xx[0];
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_O2G][1] = xx[1] / 2.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_O2G][2] = xx[2] / 3.0;
GMT->current.proj.lat_swap_vals.c[GMT_LATSWAP_O2G][3] = xx[3] / 4.0;
/* Now do the Snyder Shuffle: */
for (i = 0; i < GMT_LATSWAP_N; i++) {
GMT->current.proj.lat_swap_vals.c[i][0] = GMT->current.proj.lat_swap_vals.c[i][0] - GMT->current.proj.lat_swap_vals.c[i][2];
GMT->current.proj.lat_swap_vals.c[i][1] = 2.0 * GMT->current.proj.lat_swap_vals.c[i][1] - 4.0 * GMT->current.proj.lat_swap_vals.c[i][3];
GMT->current.proj.lat_swap_vals.c[i][2] *= 4.0;
GMT->current.proj.lat_swap_vals.c[i][3] *= 8.0;
}
}
/* The *_outside routines return the status of the current point.
* Status is the sum of x_status and y_status.
* x_status may be
* 0 w < lon < e
* -1 lon == w
* 1 lon == e
* -2 lon < w
* 2 lon > e
* y_status may be
* 0 s < lat < n
* -1 lat == s
* 1 lat == n
* -2 lat < s
* 2 lat > n
*/
/*! . */
GMT_LOCAL bool gmtmap_wesn_outside (struct GMT_CTRL *GMT, double lon, double lat) {
/* Determine if a point (lon,lat) is outside or on the rectangular lon/lat boundaries
* The check GMT->current.map.lon_wrap is include since we need to consider the 360
* degree periodicity of the longitude coordinate.
* When we are making basemaps and may want to ensure that a point is
* slightly outside the border without having it automatically flip by
* 360 degrees. In that case GMT->current.map.lon_wrap will be temporarily set to false.
*/
if (GMT->current.map.lon_wrap) {
while (lon < GMT->common.R.wesn[XLO] && (lon + GMT->current.map.lon_wrap_range) <= GMT->common.R.wesn[XHI]) lon += GMT->current.map.lon_wrap_range;
while (lon > GMT->common.R.wesn[XHI] && (lon - GMT->current.map.lon_wrap_range) >= GMT->common.R.wesn[XLO]) lon -= GMT->current.map.lon_wrap_range;
}
else if (GMT->current.map.lat_wrap) {
while (lat < GMT->common.R.wesn[YLO] && (lat + GMT->current.map.lon_wrap_range) <= GMT->common.R.wesn[YHI]) lat += GMT->current.map.lon_wrap_range;
while (lat > GMT->common.R.wesn[YHI] && (lat - GMT->current.map.lon_wrap_range) >= GMT->common.R.wesn[YLO]) lat -= GMT->current.map.lon_wrap_range;
}
/* Note PW: 8-20-2014: Was GMT_CONV4_LIMIT instead of GMT_CONV8_LIMIT. Trying the latter */
if (GMT->current.map.on_border_is_outside && fabs (lon - GMT->common.R.wesn[XLO]) < GMT_CONV8_LIMIT)
GMT->current.map.this_x_status = -1;
else if (GMT->current.map.on_border_is_outside && fabs (lon - GMT->common.R.wesn[XHI]) < GMT_CONV8_LIMIT)
GMT->current.map.this_x_status = 1;
else if (lon < GMT->common.R.wesn[XLO])
GMT->current.map.this_x_status = -2;
else if (lon > GMT->common.R.wesn[XHI])
GMT->current.map.this_x_status = 2;
else
GMT->current.map.this_x_status = 0;
if (GMT->current.map.on_border_is_outside && fabs (lat - GMT->common.R.wesn[YLO]) < GMT_CONV8_LIMIT)
GMT->current.map.this_y_status = -1;
else if (GMT->current.map.on_border_is_outside && fabs (lat - GMT->common.R.wesn[YHI]) < GMT_CONV8_LIMIT)
GMT->current.map.this_y_status = 1;
else if (lat < GMT->common.R.wesn[YLO])
GMT->current.map.this_y_status = -2;
else if (lat > GMT->common.R.wesn[YHI])
GMT->current.map.this_y_status = 2;
else
GMT->current.map.this_y_status = 0;
return (GMT->current.map.this_x_status != 0 || GMT->current.map.this_y_status != 0);
}
/*! . */
GMT_LOCAL bool gmtmap_polar_outside (struct GMT_CTRL *GMT, double lon, double lat) {
gmtmap_wesn_outside (GMT, lon, lat);
if (!GMT->current.proj.edge[1]) GMT->current.map.this_x_status = 0; /* 360 degrees, no edge */
if (GMT->current.map.this_y_status < 0 && !GMT->current.proj.edge[0]) GMT->current.map.this_y_status = 0; /* South pole enclosed */
if (GMT->current.map.this_y_status > 0 && !GMT->current.proj.edge[2]) GMT->current.map.this_y_status = 0; /* North pole enclosed */
return (GMT->current.map.this_x_status != 0 || GMT->current.map.this_y_status != 0);
}
/*! . */
GMT_LOCAL bool gmtmap_radial_outside (struct GMT_CTRL *GMT, double lon, double lat) {
double dist;
/* Test if point is more than horizon spherical degrees from origin. For global maps, let all borders be "south" */
/* Note PW: 8-20-2014: Was GMT_CONV4_LIMIT instead of GMT_CONV8_LIMIT. Trying the latter */
GMT->current.map.this_x_status = 0;
dist = gmtlib_great_circle_dist_degree (GMT, lon, lat, GMT->current.proj.central_meridian, GMT->current.proj.pole);
if (GMT->current.map.on_border_is_outside && fabs (dist - GMT->current.proj.f_horizon) < GMT_CONV8_LIMIT)
GMT->current.map.this_y_status = -1;
else if (dist > GMT->current.proj.f_horizon)
GMT->current.map.this_y_status = -2;
else
GMT->current.map.this_y_status = 0;
return (GMT->current.map.this_y_status != 0);
}
/*! . */
GMT_LOCAL bool gmtmap_rect_outside (struct GMT_CTRL *GMT, double lon, double lat) {
double x, y;
gmt_geo_to_xy (GMT, lon, lat, &x, &y);
return (gmt_cart_outside (GMT, x, y));
}
/*! . */
GMT_LOCAL bool gmtmap_rect_outside2 (struct GMT_CTRL *GMT, double lon, double lat) {
/* For Azimuthal proj with rect borders since gmtmap_rect_outside may fail for antipodal points */
if (gmtmap_radial_outside (GMT, lon, lat)) return (true); /* Point > 90 degrees away */
return (gmtmap_rect_outside (GMT, lon, lat)); /* Must check if inside box */
}
/*! . */
GMT_LOCAL void gmtmap_x_wesn_corner (struct GMT_CTRL *GMT, double *x) {
/* if (fabs (fmod (fabs (*x - GMT->common.R.wesn[XLO]), 360.0)) <= GMT_CONV4_LIMIT)
*x = GMT->common.R.wesn[XLO];
else if (fabs (fmod (fabs (*x - GMT->common.R.wesn[XHI]), 360.0)) <= GMT_CONV4_LIMIT)
*x = GMT->common.R.wesn[XHI]; */
/* Note PW: 8-20-2014: Was GMT_CONV4_LIMIT instead of GMT_CONV8_LIMIT. Trying the latter */
if (fabs (*x - GMT->common.R.wesn[XLO]) <= GMT_CONV8_LIMIT)
*x = GMT->common.R.wesn[XLO];
else if (fabs (*x - GMT->common.R.wesn[XHI]) <= GMT_CONV8_LIMIT)
*x = GMT->common.R.wesn[XHI];
}
/*! . */
GMT_LOCAL void gmtmap_y_wesn_corner (struct GMT_CTRL *GMT, double *y) {
/* Note PW: 8-20-2014: Was GMT_CONV4_LIMIT instead of GMT_CONV8_LIMIT. Trying the latter */
if (fabs (*y - GMT->common.R.wesn[YLO]) <= GMT_CONV8_LIMIT)
*y = GMT->common.R.wesn[YLO];
else if (fabs (*y - GMT->common.R.wesn[YHI]) <= GMT_CONV8_LIMIT)
*y = GMT->common.R.wesn[YHI];
}
/*! . */
GMT_LOCAL bool gmtmap_is_wesn_corner (struct GMT_CTRL *GMT, double x, double y) {
/* Checks if point is a corner */
GMT->current.map.corner = 0;
if (doubleAlmostEqualZero (fmod(fabs(x), 360.0), fmod(fabs(GMT->common.R.wesn[XLO]), 360.0))) {
if (doubleAlmostEqualZero (y, GMT->common.R.wesn[YLO]))
GMT->current.map.corner = 1;
else if (doubleAlmostEqualZero (y, GMT->common.R.wesn[YHI]))
GMT->current.map.corner = 4;
}
else if (doubleAlmostEqualZero (fmod(fabs(x), 360.0), fmod(fabs(GMT->common.R.wesn[XHI]), 360.0))) {
if (doubleAlmostEqualZero (y, GMT->common.R.wesn[YLO]))
GMT->current.map.corner = 2;
else if (doubleAlmostEqualZero (y, GMT->common.R.wesn[YHI]))
GMT->current.map.corner = 3;
}
return (GMT->current.map.corner > 0);
}
/*! . */
GMT_LOCAL bool gmtmap_lon_inside (struct GMT_CTRL *GMT, double lon, double w, double e) {
while (lon < GMT->common.R.wesn[XLO]) lon += 360.0;
while (lon > GMT->common.R.wesn[XHI]) lon -= 360.0;
if (lon < w) return (false);
if (lon > e) return (false);
return (true);
}
/*! . */
GMT_LOCAL unsigned int gmtmap_wesn_crossing (struct GMT_CTRL *GMT, double lon0, double lat0, double lon1, double lat1, double *clon, double *clat, double *xx, double *yy, unsigned int *sides) {
/* Compute all crossover points of a line segment with the rectangular lat/lon boundaries
* Since it may not be obvious which side the line may cross, and since in some cases the two points may be
* entirely outside the region but still cut through it, we first find all possible candidates and then decide
* which ones are valid crossings. We may find 0, 1, or 2 intersections */
unsigned int n = 0, i;
double d, x0, y0;
/* If wrapping is allowed: first bring both points between W and E boundaries,
* then move the western-most point east if it is further than 180 degrees away.
* This may cause the points to span the eastern boundary */
if (GMT->current.map.lon_wrap) {
while (lon0 < GMT->common.R.wesn[XLO]) lon0 += GMT->current.map.lon_wrap_range;
while (lon0 > GMT->common.R.wesn[XHI]) lon0 -= GMT->current.map.lon_wrap_range;
while (lon1 < GMT->common.R.wesn[XLO]) lon1 += GMT->current.map.lon_wrap_range;
while (lon1 > GMT->common.R.wesn[XHI]) lon1 -= GMT->current.map.lon_wrap_range;
if (fabs (lon0 - lon1) <= 0.5 * GMT->current.map.lon_wrap_range) { /* Nothing */ }
else if (lon0 < lon1)
lon0 += GMT->current.map.lon_wrap_range;
else
lon1 += GMT->current.map.lon_wrap_range;
}
else if (GMT->current.map.lat_wrap) {
while (lat0 < GMT->common.R.wesn[YLO]) lat0 += GMT->current.map.lat_wrap_range;
while (lat0 > GMT->common.R.wesn[YHI]) lat0 -= GMT->current.map.lat_wrap_range;
while (lat1 < GMT->common.R.wesn[YLO]) lat1 += GMT->current.map.lat_wrap_range;
while (lat1 > GMT->common.R.wesn[YHI]) lat1 -= GMT->current.map.lat_wrap_range;
if (fabs (lat0 - lat1) <= 0.5 * GMT->current.map.lat_wrap_range) { /* Nothing */ }
else if (lat0 < lat1)
lat0 += GMT->current.map.lat_wrap_range;
else
lat1 += GMT->current.map.lat_wrap_range;
}
/* Then set 'almost'-corners to corners */
gmtmap_x_wesn_corner (GMT, &lon0);
gmtmap_x_wesn_corner (GMT, &lon1);
gmtmap_y_wesn_corner (GMT, &lat0);
gmtmap_y_wesn_corner (GMT, &lat1);
/* Crossing South */
if ((lat0 >= GMT->common.R.wesn[YLO] && lat1 <= GMT->common.R.wesn[YLO]) || (lat1 >= GMT->common.R.wesn[YLO] && lat0 <= GMT->common.R.wesn[YLO])) {
sides[n] = 0;
clat[n] = GMT->common.R.wesn[YLO];
d = lat0 - lat1;
clon[n] = (doubleAlmostEqualZero (lat0, lat1)) ? lon1 : lon1 + (lon0 - lon1) * (clat[n] - lat1) / d;
gmtmap_x_wesn_corner (GMT, &clon[n]);
if (fabs (d) > 0.0 && gmtmap_lon_inside (GMT, clon[n], GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI])) n++;
}
/* Crossing East */
if ((lon0 >= GMT->common.R.wesn[XHI] && lon1 <= GMT->common.R.wesn[XHI]) || (lon1 >= GMT->common.R.wesn[XHI] && lon0 <= GMT->common.R.wesn[XHI])) {
sides[n] = 1;
clon[n] = GMT->common.R.wesn[XHI];
d = lon0 - lon1;
clat[n] = (doubleAlmostEqualZero (lon0, lon1)) ? lat1 : lat1 + (lat0 - lat1) * (clon[n] - lon1) / d;
gmtmap_y_wesn_corner (GMT, &clat[n]);
if (fabs (d) > 0.0 && clat[n] >= GMT->common.R.wesn[YLO] && clat[n] <= GMT->common.R.wesn[YHI]) n++;
}
/* Now adjust the longitudes so that they might span the western boundary */
if (GMT->current.map.lon_wrap && MAX(lon0, lon1) > GMT->common.R.wesn[XHI]) {
lon0 -= GMT->current.map.lon_wrap_range; lon1 -= GMT->current.map.lon_wrap_range;
}
else if (GMT->current.map.lat_wrap && MAX(lat0, lat1) > GMT->common.R.wesn[YHI]) {
lat0 -= GMT->current.map.lat_wrap_range; lat1 -= GMT->current.map.lat_wrap_range;
}
/* Crossing North */
if ((lat0 >= GMT->common.R.wesn[YHI] && lat1 <= GMT->common.R.wesn[YHI]) || (lat1 >= GMT->common.R.wesn[YHI] && lat0 <= GMT->common.R.wesn[YHI])) {
sides[n] = 2;
clat[n] = GMT->common.R.wesn[YHI];
d = lat0 - lat1;
clon[n] = (doubleAlmostEqualZero (lat0, lat1)) ? lon1 : lon1 + (lon0 - lon1) * (clat[n] - lat1) / d;
gmtmap_x_wesn_corner (GMT, &clon[n]);
if (fabs (d) > 0.0 && gmtmap_lon_inside (GMT, clon[n], GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI])) n++;
}
/* Crossing West */
if ((lon0 <= GMT->common.R.wesn[XLO] && lon1 >= GMT->common.R.wesn[XLO]) || (lon1 <= GMT->common.R.wesn[XLO] && lon0 >= GMT->common.R.wesn[XLO])) {
sides[n] = 3;
clon[n] = GMT->common.R.wesn[XLO];
d = lon0 - lon1;
clat[n] = (doubleAlmostEqualZero (lon0, lon1)) ? lat1 : lat1 + (lat0 - lat1) * (clon[n] - lon1) / d;
gmtmap_y_wesn_corner (GMT, &clat[n]);
if (fabs (d) > 0.0 && clat[n] >= GMT->common.R.wesn[YLO] && clat[n] <= GMT->common.R.wesn[YHI]) n++;
}
if (n == 0) return (0);
for (i = 0; i < n; i++) {
gmt_geo_to_xy (GMT, clon[i], clat[i], &xx[i], &yy[i]);
if (GMT->current.proj.projection_GMT == GMT_POLAR && sides[i]%2) sides[i] = 4 - sides[i]; /* toggle 1 <-> 3 */
}
if (n == 1) return (1);
/* Check for corner xover if n == 2 */
if (gmtmap_is_wesn_corner (GMT, clon[0], clat[0])) return (1);
if (gmtmap_is_wesn_corner (GMT, clon[1], clat[1])) {
clon[0] = clon[1];
clat[0] = clat[1];
xx[0] = xx[1];
yy[0] = yy[1];
sides[0] = sides[1];
return (1);
}
/* Sort the two intermediate points into the right order based on projected distances from the first point */
gmt_geo_to_xy (GMT, lon0, lat0, &x0, &y0);
if (hypot (x0 - xx[1], y0 - yy[1]) < hypot (x0 - xx[0], y0 - yy[0])) {
gmt_M_double_swap (clon[0], clon[1]);
gmt_M_double_swap (clat[0], clat[1]);
gmt_M_double_swap (xx[0], xx[1]);
gmt_M_double_swap (yy[0], yy[1]);
gmt_M_uint_swap (sides[0], sides[1]);
}
return (2);
}
#if 0
GMT_LOCAL void gmtmap_x_rect_corner (struct GMT_CTRL *GMT, double *x) {
if (fabs (*x) <= GMT_CONV4_LIMIT)
*x = 0.0;
else if (fabs (*x - GMT->current.proj.rect[XHI]) <= GMT_CONV4_LIMIT)
*x = GMT->current.proj.rect[XHI];
}
GMT_LOCAL void gmtmap_y_rect_corner (struct GMT_CTRL *GMT, double *y) {
if (fabs (*y) <= GMT_CONV4_LIMIT)
*y = 0.0;
else if (fabs (*y - GMT->current.proj.rect[YHI]) <= GMT_CONV4_LIMIT)
*y = GMT->current.proj.rect[YHI];
}
#endif
/*! . */
GMT_LOCAL void gmtmap_x_rect_corner (struct GMT_CTRL *GMT, double *x) {
if (fabs (*x) <= GMT_CONV8_LIMIT)
*x = 0.0;
else if (fabs (*x - GMT->current.proj.rect[XHI]) <= GMT_CONV8_LIMIT)
*x = GMT->current.proj.rect[XHI];
}
/*! . */
GMT_LOCAL void gmtmap_y_rect_corner (struct GMT_CTRL *GMT, double *y) {
if (fabs (*y) <= GMT_CONV8_LIMIT)
*y = 0.0;
else if (fabs (*y - GMT->current.proj.rect[YHI]) <= GMT_CONV8_LIMIT)
*y = GMT->current.proj.rect[YHI];
}
/*! . */
GMT_LOCAL bool gmtmap_is_rect_corner (struct GMT_CTRL *GMT, double x, double y) {
/* Checks if point is a corner */
GMT->current.map.corner = -1;
if (doubleAlmostEqualZero (x, GMT->current.proj.rect[XLO])) {
if (doubleAlmostEqualZero (y, GMT->current.proj.rect[YLO]))
GMT->current.map.corner = 1;
else if (doubleAlmostEqualZero (y, GMT->current.proj.rect[YHI]))
GMT->current.map.corner = 4;
}
else if (doubleAlmostEqualZero (x, GMT->current.proj.rect[XHI])) {
if (doubleAlmostEqualZero (y, GMT->current.proj.rect[YLO]))
GMT->current.map.corner = 2;
else if (doubleAlmostEqualZero (y, GMT->current.proj.rect[YHI]))
GMT->current.map.corner = 3;
}
return (GMT->current.map.corner > 0);
}
/*! . */
GMT_LOCAL unsigned int gmtmap_rect_crossing (struct GMT_CTRL *GMT, double lon0, double lat0, double lon1, double lat1, double *clon, double *clat, double *xx, double *yy, unsigned int *sides) {
/* Compute all crossover points of a line segment with the boundaries in a rectangular projection */
unsigned int i, j, n = 0;
double x0, x1, y0, y1, d;
/* Since it may not be obvious which side the line may cross, and since in some cases the two points may be
* entirely outside the region but still cut through it, we first find all possible candidates and then decide
* which ones are valid crossings. We may find 0, 1, or 2 intersections */
gmt_geo_to_xy (GMT, lon0, lat0, &x0, &y0);
gmt_geo_to_xy (GMT, lon1, lat1, &x1, &y1);
/* First set 'almost'-corners to corners */
gmtmap_x_rect_corner (GMT, &x0);
gmtmap_x_rect_corner (GMT, &x1);
gmtmap_y_rect_corner (GMT, &y0);
gmtmap_y_rect_corner (GMT, &y1);
if ((y0 >= GMT->current.proj.rect[YLO] && y1 <= GMT->current.proj.rect[YLO]) || (y1 >= GMT->current.proj.rect[YLO] && y0 <= GMT->current.proj.rect[YLO])) {
sides[n] = 0;
yy[n] = GMT->current.proj.rect[YLO];
d = y0 - y1;
xx[n] = (doubleAlmostEqualZero (y0, y1)) ? x0 : x1 + (x0 - x1) * (yy[n] - y1) / d;
gmtmap_x_rect_corner (GMT, &xx[n]);
if (fabs (d) > 0.0 && xx[n] >= GMT->current.proj.rect[XLO] && xx[n] <= GMT->current.proj.rect[XHI]) n++;
}