; $Id: fixed_map_grid.pro 31224 2022-11-07 00:11:46Z rlillis3 $
;
; 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