File: fixed_map_grid.pro

package info (click to toggle)
coyote 2019.01.29-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster
  • size: 6,316 kB
  • sloc: python: 184; makefile: 14; sh: 13
file content (788 lines) | stat: -rw-r--r-- 31,642 bytes parent folder | download | duplicates (3)
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
; $Id: //depot/idl/IDL_71/idldir/lib/map_grid.pro#1 $
;
; Copyright (c) 1996-2009, ITT Visual Information Solutions. All
;       rights reserved. Unauthorized reproduction is prohibited.
;+
; NAME:
;   FIXED_MAP_GRID
;
; PURPOSE:
;       The MAP_GRID procedure draws the graticule of parallels and meridians,
; according to the specifications established by MAP_SET. MAP_SET must be called
; before MAP_GRID to establish the projection type, the center of the
; projection, polar rotation and geographical limits.
;
; CATEGORY:
;   Mapping.
;
; CALLING SEQUENCE:
;       MAP_GRID
;
; INPUTS:
;   NONE
;
; OPTIONAL INPUTS:
;   NONE
;
; KEYWORD PARAMETERS:
;
;
; BOX_AXES: Surround the map window with a "box" style axes with
;         annotations, outside the box, where the parallels intersect the
;         sides, and the meridians intersect the bottom and top edges of the
;         box.  The border of the box is drawn in alternating foreground and
;         background colors, with color changes at each intersection with
;         a parallel or meridian.  This keyword determines the thickness of
;         the box's border, in millimeters.  If LABEL is not explicitly
;         specified, it defaults to 1 when this keyword is present.  If this
;         feature is selected, be sure to leave enough room around the map
;         window for the annotation, usually by specifying the XMARGIN and
;         YMARGIN keywords to MAP_SET.  See the example below.
;   CHARSIZE: The size of the characters used for the labels. The default is 1.
;      COLOR: The color index for the grid lines.
;FILL_HORIZON: Fills the current map_horizon.
;    HORIZON: Draws the current map horizon.
;  INCREMENT: Determines the spacing between graticle points.
; GLINESTYLE: If set, the line style used to draw the grid of parallels and
;             meridians. See the IDL User's Guide for options. The default is
;             dotted.
; GLINETHICK: The thickness of the grid lines. Default is 1.
;      LABEL: Set this keyword to label the parallels and meridians with their
;             corresponding latitudes and longitudes. Setting this keyword to
;             an integer will cause every LABEL gridline to be labeled (i.e,
;             if LABEL=3 then every third gridline will be labeled). The
;             starting point for determining which gridlines are labeled is the
;             minimum latitude or longitude (-180 to 180), unless the LATS or
;             LONS keyword is set to a single value. In this case, the starting
;             point is the value of LATS or LONS.
;   LATALIGN: This keyword controls the alignment of the text baseline for
;             latitude labels. A value of 0.0 left justifies the label,
;             1.0 right justifies it, and 0.5 centers it.
;     LATDEL: The spacing in degrees between parallels of latitude in the grid.
;             If not set, a suitable value is determined from the current map
;             projection.
;       LATS: The desired latitudes to be drawn (and optionally labeled). If
;             this keyword is omitted, appropriate latitudes will be generated
;             with consideration of the optional LATDEL keyword. If this
;             keyword is set to a single value, drawn (and optionally labeled)
;             latitudes will be automatically generated  as
;             [...,LATS-LATDEL,LATS,LATS+LATDEL,...] over the extent of the
;             map.  The single valued LATS is taken to be the starting point
;             for labelling (See the LABEL keyword).
;     LATLAB: The longitude at which to place latitude labels. The default is
;             the center longitude on the map.
;   LATNAMES: An array specifing the names to be used for the latitude labels.
;             By default, this array is automatically generated in units of
;             degrees. This array can be a string array or numeric, but should
;             not be of mixed type. When LATNAMES is specified, LATS must also
;             be specified, but the number of elments in the two arrays need not
;             be equal. Extra LATNAMES are ignored, while LATNAMES not supplied
;             are automatically reported in degrees. LATNAMES can be used when
;             LATS is set to a single value. It this case, the LATS used will
;             start at the specified latitude and progress northward, wrapping
;             around the globe if appropriate. Caution should be used when
;             using LATNAMES in conjunction with a single LATS, since the
;             number of visible latitude gridlines is dependent on many factors.
;   LONALIGN: This keyword controls the alignment of the text baseline for
;             longitude labels. A value of 0.0 left justifies the label,
;             1.0 right justifies it, and 0.5 centers it.
;     LONDEL: The spacing in degrees between meridians of longitude in the grid.
;             If not set, a suitable value is determined from the current map
;             projection.
;     LONLAB: The latitude at which to place longitude labels. The default is
;             the center latitude on the map.
;       LONS: The desired longitudes to be drawn (and optionally labeled). If
;             this keyword is omitted, appropriate longitudes will be generated
;             with consideration of the optional LONDEL keyword. If this
;             keyword is set to a single value, drawn (and optionally labeled)
;             longitudes will be automatically generated as
;             [...,LONS-LONDEL,LONS,LONS+LONDEL,...] over the extent of the map.
;             The single valued LONS is taken to be the starting point for
;             labeling (See the LABEL keyword).
;   LONNAMES: An array specifing the names to be used for the longitude labels.
;             By default, this array is automatically generated in units of
;             degrees. This array can be a string array or numeric, but should
;             not be of mixed type. When LONNAMES is specified, LATS must also
;             be specified, but the number of elments in the two arrays need not
;             be equal. Extra LONNAMES are ignored, while LATNAMES not supplied
;             are automatically reported in degrees. LONNAMES can be used when
;             LONS is set to a single value. It this case, the LONS used will
;             start at the specified longitude and progress eastward, wrapping
;             around the globe if appropriate. Caution should be used when
;             using LONNAMES in conjunction with a single LONS, since the number
;             of visible longitude gridlines is dependent on many factors.
;
; MAP_STRUCTURE: Set this keyword to a !MAP structure, as returned from
;       MAP_PROJ_INIT. If this keyword is set then the gridlines are passed
;       through MAP_PROJ_FORWARD and the resulting lines are drawn using
;       Cartesian coordinates. If this keyword is not set then it is assumed
;       that MAP_SET has been called to set up a map projection.
;
; NO_GRID: Set this keyword if you only want labels but not gridlines.
;
;ORIENTATION: Specifies the clockwise angle in degrees from horizontal
;             of the text baseline that the labels should be rotated. Note,
;             that the rotation used in MAP_SET is also clockwise, but XYOUTS
;             is normally counterclockwise.
;
; T3D: Set this keyword to indicate that the generalized transformation
;      matrix in !P.T is to be used. If not present, the user-supplied
;      coordinates are simply scaled to screen coordinates.
;
; ZVALUE: Sets the Z coordinate, in normalized coordinates in the
;         range of 0 to 1, at which to output the continents.
;
;      Note - This keyword has effect only if keyword T3D is set and the
;         transformation is stored in !P.T
;
;OUTPUTS:
;         Draws gridlines and labels on the current map.
;
; EXAMPLE:
;     map_set,10,20,23.5,/cont,/ortho,/horizon     ; establish a projection
;         lats=[-65,-52,-35,-20,-2.59,15,27.6,35,45,55,75,89]
;                     ; choose the parallels to draw
;         map_grid,lats=lats,londel=20,lons=-13,label=2,lonlab=-2.5,latlab=7
;                     ;draw the grid
;   Make a map with a grid surrounded by a box style axis:
;   map_set, /STEREO, 40, -90,scale=50e6,/CONTINENTS, XMARGIN=3, YMARGIN=3
;   map_grid, /BOX_AXES, COLOR=3, CHARSIZE=1.5  ;
;
; MODIFICATION HISTORY:
;   Written by: SVP, May, 23 1996.
;   DMS, Feb, 1996 Note that this version plots all gridlines
;           unless a 4 element LIMIT was specified to MAP_SET.
;       SVP, Nov 1996   Changed over !map1 (new IDL5 maps)
;   SVP, May 1996   Map_Grid used to live inside of map_set.pro, now
;                       it lives here.
;       SVP, May 1996   Changed LABEL to determine the skip between labelled
;                       gridlines.
;                       Also, added the LATS and LONS keywords to give complete
;                       control over where the labels go.
;       SVP, Sept 1996  Added LATNAMS,LONAMES, and ORIENTATION keywords
;   DMS, Nov, 1996  Rev 2 of maps.
;   DMS, Jul, 1997  Added Box Axes
;   CT, Jan 2004: Added MAP_STRUCTURE keyword.
;   Renamed the function FIXED_MAP_GRID with fixed at lines 555-556 and
;      664-665 as described here: http://www.dfanning.com/map_tips/missinggrid.html.
;      This is required to get proper grid output in some applications, but
;      it breaks other code.
;-
function map_grid_check_range, xy
;
; Determine if the xy label point is inside the plot coordinate axes.
; If not, return a single value, which is NOT drawn on the plot. This
; solves the problem of non-Hershey fonts not being clipped to the plot
; area.

    x = xy[0]
    y = xy[1]
    
    if x lt Min(!X.CRange) then return, xy[0]
    if x gt Max(!X.CRange) then return, xy[0]
    if y lt Min(!Y.CRange) then return, xy[0]
    if y gt Max(!Y.CRange) then return, xy[0]
    return, xy
end


function fixed_map_grid_incr, span
;
; Determine LONDEL or LATDEL if not specified
;
COMPILE_OPT hidden, IDL2

IF span eq 0 THEN return, 45.
ipow = 0
t = abs(span) < 450.
WHILE t lt 5 DO BEGIN t = t * 10 & ipow = ipow +1 & ENDWHILE
increments = [ 1., 2., 4., 5., 10., 15., 30., 45.]
i = 0
WHILE t gt (increments[i] * 10) DO i = i + 1
t = increments[i] / 10^ipow
return, span ge 0 ? t : -t
end


;-------------------------------------------------------------------------
; Find the point on the line between points c0 and c1, expressed in
; DEVICE coordinates, where the longitude (Icoord = 0) or latitude
; (Icoord = 1) is equal to Gwant.  If the segment between c0 and c1
; doesn't intersect the given meridan or parallel, or either endpoint
; is not mappable, return NaN.
; Otherwise, return the device coordinate, X if Icoord = 0, or Y if
; Icoord = 1, of the intersection.
;
Function fixed_map_grid_SOLVE, c0, c1, Icoord, Gwant, $
    MAP_STRUCTURE=mapStruct, RECURSIVE=recursive

    COMPILE_OPT hidden, IDL2

    hasMap = N_TAGS(mapStruct) gt 0

    p0 = CONVERT_COORD(c0, /DEVICE, /TO_DATA)
    p1 = CONVERT_COORD(c1, /DEVICE, /TO_DATA)

    if (hasMap) then begin   ; Convert from UV to latlon
        p0 = MAP_PROJ_INVERSE(p0[0], p0[1], MAP_STRUCTURE=mapStruct)
        p1 = MAP_PROJ_INVERSE(p1[0], p1[1], MAP_STRUCTURE=mapStruct)
    endif

    p0 = p0[Icoord]
    p1 = p1[Icoord]

    ; Not mappable or zero interval.
    if ~finite(p0) || ~finite(p1) || (p1 eq p0) then $
        return, !values.f_nan

    if (Icoord eq 0) && (p0 gt p1) then begin ;Wrap if we cross dateline
        if gwant ge 0 then p1 += 360. $
        else p0 -= 360.
    endif

    t = (Gwant - p0) / (p1-p0)
    if (t lt 0.0) || (t gt 1.0) then $
        return, !values.f_nan

    low = 0.0
    high = 1.0
    tol = 1.0e-5
    del = c1 - c0
    while abs(high-low) gt tol do begin ;Binary chopping method
        t = (low + high) / 2.
        c = c0 + t * del
        p = CONVERT_COORD(c, /DEVICE, /TO_DATA)
        if (hasMap) then $   ; Convert from UV to latlon
            p = MAP_PROJ_INVERSE(p[0], p[1], MAP_STRUCTURE=mapStruct)
        p = p[Icoord]
        if (~FINITE(p)) then $
            return, !values.f_nan
        if (Icoord eq 0) then begin ;Check for dateline?
            if p lt p0 then p += 360. $ ;Wrap?  P should be in interval p0-p1.
            else if p gt p1 then p -= 360.
        endif
        if (Gwant-p0) * (Gwant - p) gt 0.0 then begin ;In same interval as p0 : low
            low = t
            p0 = p
        endif else high = t         ;Keep low, and fcn at low = p0.
    endwhile

    return, c[Icoord]

end


;-------------------------------------------------------------------------
PRO fixed_map_grid, LABEL=label, LATDEL = latdel, NO_GRID=no_grid, $
       LONDEL=londel, GLINESTYLE=glinestyle, GLINETHICK=glinethick, $
       LONLAB=lonlab, LATLAB=latlab, LONALIGN=lonalign, $
       LATALIGN=latalign, LONNAMES=lonnames, LATNAMES=latnames, $
       LATS=lats, LONS=lons, ZVALUE=zvalue, $
       COLOR=color, T3D=t3d, CHARSIZE=charsize, ORIENTATION=orientation, $
       HORIZON=horizon, FILL_HORIZON=fill_horizon, _EXTRA=extra, $
       INCREMENT=increment, CLIP_TEXT=clip_text, BOX_AXES=box_axes, $
       MAP_STRUCTURE=mapStruct, FORMAT=format, $
       WHOLE_MAP=obsolete_keyword


compile_opt idl2

ON_ERROR, 2

hasMap = N_TAGS(mapStruct) gt 0

if ((!x.type NE 3) && ~hasMap) THEN $
   message, 'fixed_map_grid---Current ploting device must have mapping coordinates'

; Put a grid on a previously established map projection.
;
; no grid? - in case someone wants just to put labels
no_grid = keyword_set(no_grid)

; if Label = n, then Labels are added every n gridlines
;   If box_axes is set, and LABEL isn't explicitly specified, set label.
;
nlabel = (N_ELEMENTS(label) gt 0) ? $
    FIX(ABS(label[0])) : KEYWORD_SET(box_axes)

have_lons =  n_elements(lons) gt 0
have_lats =  n_elements(lats) gt 0

if n_elements(zvalue) eq 0 THEN zvalue = 0

; CLIP_TEXT (default value = 1) = 1 to clip text within the map area,
; 0 to not clip text.
noclip = (N_ELEMENTS(clip_text) gt 0) ? ~KEYWORD_SET(clip_text) : 0

;   Append the graphics keywords:
if n_elements(t3d) then map_struct_append, extra,'T3D',t3d
if n_elements(color) then map_struct_append, extra, 'COLOR',color

if n_elements(glinestyle) eq 0 then glinestyle = 1 ;Default = dotted
map_struct_append, extra, 'LINESTYLE', glinestyle ;Append it

if n_elements(glinethick) then map_struct_append, extra,'THICK',glinethick
if n_elements(charsize) then map_struct_append, extra,'CHARSIZE', charsize

                                ;Orientation is reversed & conflicts w/box_axes
if n_elements(orientation) and (keyword_set(box_axes) eq 0) then $
  map_struct_append, extra,'ORIENTATION', -1 * orientation


;
; The gridlines can be specified by
;
;  1) an array of lats and/or lons
;  2) a single lats or lons which is taken to be the center
;     or 'for sure' lat or lon with gridlines every latdel or londel from it
;  3) automatically calculated if lats or lons are not specified.
;
;
; Require that LATS be specified when LATNAMES is ALSO SPECIFIED
;
if (n_elements(latnames) gt 0) and n_elements(lats) le 1 then $
  message,'fixed_map_grid---The LATNAMES keyword MUST be used in conjuction '+$
  'with the LATS keyword.'
if n_elements(lonnames) gt 0 and have_lons eq 0 then $
  message,'fixed_map_grid---The LONNAMES keyword MUST be used in conjuction '+$
  'with the LONS keyword.'

; Get lat/lon ranges from !MAP. Did MAP_SET specify 4 element limit?
map = hasMap ? mapStruct : !MAP
if n_elements(lats) gt 1 then latmin = min(lats, max=latmax) $
else if map.ll_box[0] ne map.ll_box[2] then begin
    latmin = map.ll_box[0]
    latmax = map.ll_box[2]
endif else begin
    latmin = -90
    latmax = 90
endelse

if have_lons then begin         ;Lons directly specified?
    lonmin = lons[0]
    lonmax = lons[n_elements(lons)-1]
endif else if (map.ll_box[1] ne map.ll_box[3]) and $ ;Lon limit specified?
  (latmax lt 90.) and (latmin gt -90.) then begin ; and poles not included
    lonmin = map.ll_box[1]
    lonmax = map.ll_box[3] ;Copy limits
endif else begin                ;If not, use entire globe
    lonmin = -180
    lonmax = 180
endelse

IF lonmax le lonmin THEN lonmax = lonmax + 360.

            ;Default grid spacings...
IF n_elements(latdel) eq 0 THEN begin
    latdel = fixed_map_grid_incr(latmax - latmin)
    latd = 1
endif else latd = latdel

IF n_elements(londel) eq 0 THEN begin
    londel = fixed_map_grid_incr(lonmax - lonmin)
    lond = 1
endif else lond = londel


; IF the deltas are < 1,
; do not convert the limits into integers
IF abs(latmax - latmin) gt 5. and latd ge 1 THEN BEGIN ;Make range integers
    latmin = float(floor(latmin))
    latmax = ceil(latmax)
ENDIF

IF abs(lonmax - lonmin) gt 5 and lond ge 1 THEN BEGIN ;Integerize long spans
    lonmin = float(floor(lonmin))
    lonmax = ceil(lonmax)
ENDIF

; Where we label things...
IF N_Elements(Latlab) eq 0 THEN Latlab = (lonmin + lonmax)/2
IF N_ELements(LonLab) eq 0 THEN LonLab = (latmin +latmax)/2

IF n_elements(latalign) eq 0 THEN latalign = .5 ;Text alignment of lat labels
IF n_elements(lonalign) eq 0 THEN lonalign = .5 ;Text alignment of lon labels

; Is this a cylindrical proj?
if (hasMap) then $
    is_cyl = 0 $
else $
    map_proj_info, iproj, CYLINDRICAL=is_cyl, /CURRENT

if keyword_set(increment) then step = increment $
else step = 4 < (latmax - latmin)/10.

len = long(float((latmax-latmin)) / float(step) + 1.0)

; Clip to avoid roundoff errors which can cause the latitude to exceed
; 90 degs by a very small amount.
lati = (float(latmax-latmin) / (len-1.)) * findgen(len) + latmin > (-90) < 90

; This fudge avoids curved meridians at the poles because of the split planes
if is_cyl and map.p0lat eq 0 then begin
    del = 2.0e-2
    if lati[0] eq -90 then lati[0] = del - 90.
    if lati[len-1] eq 90 then lati[len-1] = 90. - del
endif


;Compute longit distance between points for latitude lines.
step = 4 < (lonmax - lonmin)/10. ;At most 4 degrees
len = (lonmax-lonmin)/step + 1
loni = findgen(len) * step + lonmin
IF (loni[len-1] NE lonmax) THEN loni = [loni, lonmax]


;
; Determine the number of lons and the lon array
;
if n_elements(lons) eq 0 then begin
    n_lons = 1+fix((lonmax-lonmin) / londel)
    longitudes = lonmin - (lonmin mod londel) + findgen(n_lons) * londel
endif else if n_elements(lons) eq 1 then begin
    i0 = ceil((lonmin - lons[0]) / float(londel)) ;First tick
    i1 = floor((lonmax - lons[0]) / float(londel)) ;Last tick
    n_lons = i1 - i0 + 1 > 1
    longitudes = (findgen(n_lons) + i0) * londel + lons[0]
endif else begin
    n_lons=n_elements(lons)
    longitudes=lons
endelse

;
; Determine the number of lats and the lat array
;
if n_elements(lats) eq 0 then begin
    lat0 = latmin - (latmin mod float(latdel)) ;1st lat for grid
    n_lats = 1 + fix((latmax-lat0)/float(latdel))
    latitudes = lat0 + findgen(n_lats)*latdel
endif else if n_elements(lats) eq 1 then begin
    i0 = ceil((latmin - lats[0]) / float(latdel)) ;First tick
    i1 = floor((latmax - lats[0]) / float(latdel)) ;Last tick
    n_lats = i1 - i0 + 1 > 1
    latitudes = (findgen(n_lats) + i0) * latdel + lats[0]
endif else begin
    n_lats=n_elements(lats)
    latitudes=lats
endelse

;
; Build the Latitude/Longitude Label Flags
;
lon_label = bytarr(n_lons)
lat_label = bytarr(n_lats)
if nlabel ne 0 then begin
    if n_elements(lons) eq 1 then begin ; Ensure center is set and then go out
        index=where(longitudes eq lons[0])
        for i=(index[0] > 0) mod nlabel, n_lons-1, nlabel do lon_label[i] = 1
    endif else begin
        for i=0, n_lons-1, nlabel do lon_label[i] = 1
    endelse

    if n_elements(lats) eq 1 then begin ; Make sure the center one is set
                                ; and go out from there
        index=where(latitudes eq lats[0], count)
        for i=(index[0] > 0) mod nlabel, n_lats-1, nlabel do lat_label[i] = 1
    endif else begin            ; Start with latmin and label each nlabel point
        for i=0, n_lats-1, nlabel do lat_label[i] = 1
    endelse

endif
;
;   Dont repeat 180 labelling if the projection is cylindrical or
;   polar and both 180 and -180 are present. This can be defeated by using
;   LONS=-180
;
if is_cyl or (abs(map.p0lat) eq 90) then begin
    id_180 = where(longitudes eq 180,count)
    id_m180 = where(longitudes eq -180,mcount)
    if count gt 0 and mcount gt 0 then begin
        if n_elements(lons) eq 1 then begin
            if lons[0] eq -180 then lon_label[id_180]=0
        endif else lon_label[id_m180]=0
    endif
endif

n = n_lons > n_lats             ;
latlontxt = strarr(n, 2)

if keyword_set(box_axes) then begin ;Draw a Box legend?
    box_thick = box_axes * 0.1  ;From mm to Thickness in cm
    dc = !d.y_ch_size           ;Char height to draw
    if n_elements(charsize) then dc = dc * charsize
    xw = !x.window * !d.x_size  ;Window coords in x & y
    yw = !y.window * !d.y_size
; xww and yww = corners of the uv_range that is mappable.  If NOBORDER
; was set for MAP_SET, this is the same as the window coords (xw,yw),
; otherwise, this rectangle is smaller than the window rectangle.
; Fudge factor for window to ensure that the edges are mappable.
    del = [1,-1]* 0.01
    xww = (map.uv_box[[0,2]] * !x.s[1] + !x.s[0]) * !d.x_size + del
    yww = (map.uv_box[[1,3]] * !y.s[1] + !y.s[0]) * !d.y_size + del

    bdel = box_thick * !d.y_px_cm ;Thickness of box in device units

    if n_elements(color) eq 0 then bcolor = !p.color $  ;Box color
    else bcolor = color

    xp = xw[0] - [0,bdel, bdel,0] ;X  & Y polygon coords for outer box
    yp = yw[0] - [0,0,bdel,bdel]
                                ;Draw the outline of the box
    plots, xw[[0,1,1,0,0]], yw[[0,0,1,1,0]], /DEVICE, COLOR=bcolor
    plots, xw[[0,1,1,0,0]]+[-bdel, bdel, bdel, -bdel, -bdel], $
      yw[[0,0,1,1,0]]+[-bdel, -bdel, bdel, bdel, -bdel], /DEVICE, COLOR=bcolor

    ychar = [yw[0]-bdel-dc, yw[1]+bdel+dc/4.]
    xchar = [xw[0] - bdel - dc/4., xw[1]+bdel+dc/4.]
    boxpos = replicate(!values.f_nan, n, 2,2)

;Device coordinates for box annotations. Go in to avoid edges of map &
;border which are frought with singularities.  Also to avoid effects
;of MAP_SET,/NOBORDER.  For box axes to be annotated, all the edges of the
;map rectangle must be mappable.
;
endif else box_thick = 0

                                ;  Do the horizon if specified.
if keyword_set(horizon) then map_horizon, _EXTRA=e
if keyword_set(fill_horizon) then map_horizon, _EXTRA=e, /FILL

;
;   ****************** Draw/Label the meridians ******************
;
FOR i=0,n_lons-1 DO BEGIN
    lon=longitudes[i]
    lon2 = (lon lt -180) ? (lon + 360) : $
        ((lon gt 180) ? (lon - 360) : lon)

; This block of code draws longitude lines that are at the center + or
; - 180 degrees, at center + or - (180-eps) to ensure that the grid
;   appears on the correct side.  Its really not necessary if people
;   would use the /HORIZON keyword, but they don't.

    if is_cyl then begin
        del = lon - map.p0lon
        while del gt 180 do del -= 360.
        while del lt -180 do del += 360.
        if abs(del) eq 180 then $
            lon -= 1.0e-5 * ((del ge 0) ? 1 : -1) ;fudge it (sign(1.0e-5, del))
    endif

    IF (~no_grid) THEN begin
        if (hasMap) then begin
            ; Make sure to clear out variable in case all points are clipped.
            polylines = -1
;            uv = MAP_PROJ_FORWARD(REPLICATE(lon, N_ELEMENTS(lati)), lati, $
;                MAP_STRUCTURE=mapStruct, POLYLINES=polylines)
            uv = MAP_PROJ_FORWARD(REPLICATE(lon, N_ELEMENTS(lati)), lati, $
                MAP_STRUCTURE=mapStruct)
            polylines = [N_Elements(lati), Indgen(N_Elements(lati))]
            index = 0L
            npoly = N_ELEMENTS(polylines)
            ; Loop thru our polylines connectivity array.
            while (index lt npoly) do begin
                nline = polylines[index]
                if (nline eq -1) then $
                    break
                if (nline gt 0) then begin
                    indices = polylines[index + 1 : index + nline]
                    PLOTS, REFORM(uv[0,indices]), REFORM(uv[1,indices]), $
                        zvalue, $
                        NOCLIP=1, $
                        _EXTRA=extra
                endif
                index += nline + 1
            endwhile
        endif else begin
            PLOTS, lon, lati, zvalue, NOCLIP=0, _EXTRA=extra
        endelse
    endif

    IF N_Elements(format) EQ 0 THEN BEGIN
        fmt = (lon2 ne long(lon2)) ? '(f7.2)' : '(i4)' 
    ENDIF ELSE BEGIN 
        IF format EQ "" THEN fmt = (lon2 ne long(lon2)) ? '(f7.2)' : '(i4)' ELSE fmt = format
    ENDELSE

    IF lon_label[i] THEN BEGIN

        IF i lt n_elements(lonnames) then begin ;User specified label?
            IF (reverse(size(lonnames[i])))[1] eq 7 then $ ;String?
              lonname=lonnames[i] else $
              lonname=strtrim(string(lonnames[i], FORMAT=fmt),2)
        endif else $
            lonname=strtrim(string(lon2, format=fmt),2)

        latlontxt[i,0] = lonname
        if (box_thick eq 0) then begin
            xy = 0
            if (hasMap) then begin
                ; We need to convert from latlon to UV ourself.
                uv = MAP_PROJ_FORWARD(lon, LonLab, MAP_STRUCTURE=mapStruct)
                if (FINITE(uv[0]) && FINITE(uv[1])) then begin
                    xy = uv[0:1]
                    xy = map_grid_check_range(xy)
                endif
            endif else begin
                if (noclip eq 1) || map_point_valid(lon, LonLab) then $
                    xy = [lon, LonLab]
            endelse
            if (N_ELEMENTS(xy) eq 2) then $
                XYOUTS, xy[0], xy[1], lonname, ALIGNMENT=lonalign, $
                    NOCLIP=noclip, Z=zvalue, _EXTRA=extra
        endif

    ENDIF

    if box_thick ne 0 then begin
        dy = (yw[1] - yw[0]) * 0.01 ;1% of the height
        for j=0,1 do begin      ;Save longitude crossings, try for edge...
            k = 0
            ; Try to find longitude crossings.  If it doesn't cross exactly
            ; at the edge, try going in until it crosses and is valid.
            while ~finite(boxpos[i,j,0]) && abs(k) lt 3 do begin
                boxpos[i, j, 0] = fixed_map_grid_solve( $
                    [xww[0], yw[j]+k*dy], $
                    [xww[1], yw[j]+k*dy], 0, lon, $
                    MAP_STRUCTURE=mapStruct)
                k = k + (j ? -1 : 1)
            endwhile
        endfor
    endif
ENDFOR

;
; Draw/Label the parallels of latitude  ******************
;

FOR i=0,n_lats-1 DO BEGIN

    lat=latitudes[i]

    IF N_Elements(format) EQ 0 THEN BEGIN
        fmt = (lat ne long(lat)) ? '(f7.2)' : '(i4)' 
    ENDIF ELSE BEGIN 
        IF format EQ "" THEN fmt = (lat ne long(lat)) ? '(f7.2)' : '(i4)'  ELSE fmt = format
    ENDELSE

    IF lat_label[i] THEN BEGIN
        IF i lt n_elements(latnames) then begin ;User specified latname?
            IF (reverse(size(latnames[i])))[1] eq 7 then $
              latname=latnames[i] else $
              latname=strtrim(string(latnames[i],format=fmt),2)
        endif else latname=strtrim(string(lat, format=fmt),2)
        latlontxt[i, 1] = latname
        if (box_thick eq 0) then begin
            xy = 0
            if (hasMap) then begin
                ; We need to convert from latlon to UV ourself.
                uv = MAP_PROJ_FORWARD(latlab, lat, MAP_STRUCTURE=mapStruct)
                if (FINITE(uv[0]) && FINITE(uv[1])) then begin
                    xy = uv[0:1]
                    xy = map_grid_check_range(xy)
                endif
            endif else begin
                if (noclip eq 1) || map_point_valid(latlab, lat) then $
                    xy = [latlab, lat]
            endelse
            if (N_ELEMENTS(xy) eq 2) then $
                XYOUTS, xy[0], xy[1], latname, $
                    alignment=latalign, NOCLIP=noclip, Z=zvalue, _EXTRA=extra
        endif
    ENDIF

    if (~no_grid && (ABS(lat) ne 90)) then begin
        if (hasMap) then begin
            ; Make sure to clear out variable in case all points are clipped.
            polylines = -1
;            uv = MAP_PROJ_FORWARD(loni, REPLICATE(lat, N_ELEMENTS(loni)), $
;                MAP_STRUCTURE=mapStruct, POLYLINES=polylines)
            uv = MAP_PROJ_FORWARD(loni, REPLICATE(lat, N_ELEMENTS(loni)), $
                MAP_STRUCTURE=mapStruct)
            polylines = [N_Elements(loni), Indgen(N_Elements(loni))]
            index = 0L
            npoly = N_ELEMENTS(polylines)
            ; Loop thru our polylines connectivity array.
            while (index lt npoly) do begin
                nline = polylines[index]
                if (nline eq -1) then $
                    break
                if (nline gt 0) then begin
                    indices = polylines[index + 1 : index + nline]
                    PLOTS, REFORM(uv[0,indices]), REFORM(uv[1,indices]), $
                        zvalue, $
                        NOCLIP=0, $
                        _EXTRA=extra
                endif
                index += nline + 1
            endwhile
        endif else $
            PLOTS, loni, lat, zvalue, NOCLIP=0, _EXTRA=extra
    endif

    if box_thick ne 0 then begin
        for j=0,1 do begin ;Save latitude crossings
            ; Start at edge and try for an intersection.
            ; If that doesn't work, go in some.
            k = 0
            dx = (xw[1] - xw[0]) * 0.01
            while ~finite(boxpos[i,j,1]) && abs(k) lt 3 do begin
                boxpos[i, j, 1] = fixed_map_grid_solve( $
                    [xw[j]+dx*k, yww[0]], $
                    [xw[j]+dx*k, yww[1]], 1, lat, $
                    MAP_STRUCTURE=mapStruct)
                k = k + (j ? -1 : 1)
            endwhile
        endfor
    endif

endfor


; ******************************** Do the box axes **************************
if box_thick ne 0 then begin
    for iaxis=0,1 do for j=0,1 do begin
        v = boxpos[*,j,iaxis]       ;Values along axes
        good = where(finite(v), count) ;Ignore bad values
        if (count eq 0) then $  ;Anything there?
            continue
        dy = iaxis eq 1
        v = v[good]             ;Remove unmappable elements
        subs = sort(v)          ;Sort the axis crossings
        v = v[subs]             ;Sorted Y values
        vtext = (latlontxt[good,iaxis])[subs]
        v0 = ([xw[0], yw[0]])[iaxis] ;Starting value on axis
        xp0 = xp + j * (xw[1]-xw[0] + bdel) ;Polygon X coords
        yp0 = yp + j * (yw[1]-yw[0] + bdel) ;Y coords
        xychar = [xchar[j], ychar[j]] ;Char position

        for i=0, count-1 do begin ;Draw each item
            z = v[i]            ;Axis crossing value
            if iaxis eq 0 then begin
                xp0 = (i eq (count-1) && (i and 1) && ~vtext[i]) ? $
                    [v0, xw[1], xw[1], v0] : [v0, z, z, v0]
            endif else yp0 = [v0, v0, z, z]
            if (i and 1) then $
                polyfill, xp0, yp0, /DEVICE, COLOR=bcolor
            xychar[iaxis] = z
            if strlen(vtext[i]) gt 0 then begin
                xyouts, xychar[0], xychar[1], vtext[i], $
                    ORIENTATION=dy * (90-180*j), $
                    ALIGN=0.5, CLIP=0, Z=zvalue, /DEVICE, _EXTRA=extra
            endif
            v0 = z
        endfor
                                ;Fill to the end of the axis
        if i and 1 then begin
            if iaxis eq 0 then xp0 = [v0, xw[1], xw[1], v0] $
            else yp0 = [v0, v0, yw[1], yw[1]]
            polyfill, xp0, yp0, /DEVICE, COLOR=bcolor
        endif
    endfor
endif   ; box_thick


end