; docformat = 'rst' ; ; NAME: ; cgMap_Grid ; ; PURPOSE: ; Provides a significantly modified MAP_GRID command that can be used together ; with other Coyote Graphics routines and in resizeable graphics windows. ; ;******************************************************************************************; ; ; ; Copyright (c) 2011, by Fanning Software Consulting, Inc. All rights reserved. ; ; ; ; Redistribution and use in source and binary forms, with or without ; ; modification, are permitted provided that the following conditions are met: ; ; ; ; * Redistributions of source code must retain the above copyright ; ; notice, this list of conditions and the following disclaimer. ; ; * Redistributions in binary form must reproduce the above copyright ; ; notice, this list of conditions and the following disclaimer in the ; ; documentation and/or other materials provided with the distribution. ; ; * Neither the name of Fanning Software Consulting, Inc. nor the names of its ; ; contributors may be used to endorse or promote products derived from this ; ; software without specific prior written permission. ; ; ; ; THIS SOFTWARE IS PROVIDED BY FANNING SOFTWARE CONSULTING, INC. ''AS IS'' AND ANY ; ; EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ; ; OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT ; ; SHALL FANNING SOFTWARE CONSULTING, INC. BE LIABLE FOR ANY DIRECT, INDIRECT, ; ; INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED ; ; TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; ; ; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ; ; ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ; ; (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ; ; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ; ;******************************************************************************************; ; ;+-------------------------------------------------------------------------- ; This program provides a significantly modified MAP_GRID command that ; can be used with other Coyote Graphics routines. ; ; Description of known MAP_GRID problems:: ; http://www.idlcoyote.com/map_tips/noclip.html ; http://www.idlcoyote.com/map_tips/missinggrid.html ; http://www.idlcoyote.com/map_tips/extralines.php ; :Categories: ; Graphics, Map Projections ; ; :Author: ; FANNING SOFTWARE CONSULTING:: ; David W. Fanning ; 1645 Sheely Drive ; Fort Collins, CO 80526 USA ; Phone: 970-221-0438 ; E-mail: david@idlcoyote.com ; Coyote's Guide to IDL Programming: http://www.idlcoyote.com ; ; :History: ; Change History:: ; Written by David W. Fanning, 7 November 2011. ; Significant modification of the MAP_GRID command in IDL to fix known ; problems, especially when used in conjunction with map projection spaces ; set up with MAP_PROJ_INIT. David W. Fanning, 7 November 2011. ; Added an ERASE=0 to the /NOGRAPHICS keyword on the Draw method call to cgMap. 27 Dec 2011. ; Changed the default line thickness to !P.Thick to better support PostScript files. 28 Dec 2011. DWF. ; Fixed a problem with grid labeling when values were passed with LATS or LONS keyword. 6 April 2012. DWF. ; Modified slightly to allow a three-element byte array to be used as the COLOR. 18 April 2012. DWF. ; If a Map object is available, I make sure to call DRAW method before drawing graphics. 12 Sept 2012. DWF. ; Added cgGRID keyword to allow the cgMap object to create latitude and longitude grid in its ; LatLonLabels method. Previously used by default, but it doesn't work well with global ; map projections. It works best with small map areas in UTM projection space. 3 Jan 2013. DWF. ; Removed some old code that was used to correct latitude and longitude values. No longer needed, ; I hope, with the new cgGRID keyword. 3 Jan 2013. ; Corrected bug in variable spelling that affect LONDELTA and LATDELTA keywords. 6 Jan 2013. DWF. ; Lost a piece of code that allows longitude box axes. Added back in. 23 Jan 2013. DWF. ; T3D keyword was not being applied. Fixed. 11 February 2013. DWF. ; Added NOCLIP keyword. 15 February 2013. DWF. ; Sometimes a longitude line is draw incorrectly due to the fact that the longitude vector does not ; have a point in the XRANGE of the projection. A fix to that problem has failed to work in all ; circumstances, so I have done more work on that algorithm to see if I can solve the problem is ; a better way. Now usine Value_Locate to test for point. 19 February 2013. DWF. ; Now implementing label formatting when using CGGRID keyword. 27 November 2013. DWF. ; ; :Copyright: ; Copyright (c) 2011-2013, Fanning Software Consulting, Inc. ;--------------------------------------------------------------------------- ; ; ;+------------------------------------------------------------------------ ; ; Find the grid increment find defining the latitude and longitude delta ; values, if they are not currently defined. ; ; :Params: ; span: in, required, type=float ; The data range. ; ;------------------------------------------------------------------------- function cgMap_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 = [ 0.1, 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 retvalue = span ge 0 ? t : -t return, retvalue 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. ; ; :Params: ; c0: in, required, type=integer ; Input coordinate? ; c1: in, required, type=integer ; Input coordinate? ; icoord: in, required, type=integer ; Index of input coordinate? ; gwant: in, required, type=integer ; Global wrapping? ; ; :Keywords: ; map_structure: in, optional, type=structure ; The map structure to covert XY projected meters to lat/lon space. ; ;------------------------------------------------------------------------- Function cgMap_Grid_Solve, c0, c1, icoord, gwant, MAP_STRUCTURE=mapStruct 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 ; ;+------------------------------------------------------------------------ ; This routine fixes a bug in MAP_GRID that causes map labels to be ; written outside the map boundary when using hardware or true-type ; fonts. It checks to be sure the label is inside the map boundary ; before it is written. Users can control how "exact" the boundary is ; when using GCTP map projections by setting the FUZZY keyword to ; a multiplication factor that is multiplied times the calculated ; data range of the map projection. ; ; If a point is determined to be outside the map boundary, a single ; data value is returned by the function. This is a signal that this ; data point should not be drawn. ; ; :Params: ; xy: in, required, type=float ; The input label point. Normally, a two element array. ; ; :Keywords: ; fuzzy: in, optional, type=float, default=0.0 ; This keyword applies only if the GCTP keyword is set. Set the ; keyword to a value that is a percentage of the current data range. ; This percentage of the range is added to or subtracted from the ; values used to determine if the label is "inside" the boundary. ; It allows you to be a little less exact when selecting inside ; points. There are occasional aesthetic reasons for allowing fuzzy ; boundaries. A reasonable value for fuzziness might be 0.0125. ; gctp: in, optional, type=boolean, default=0 ; Set this keyword to calculate boundaries in terms of the the ; current calculated data range variables !X.CRange and !Y.CRange. ; Otherwise, the boundaries are taken from !Map.LL_BOX. ; ;------------------------------------------------------------------------- function cgMap_Grid_Check_Range, xy, FUZZY=fuzzy, GCTP=gctp ; You need at least two points coming in here. IF N_Elements(xy) NE 2 THEN RETURN, xy x = xy[0] y = xy[1] IF Keyword_Set(gctp) THEN SetDefaultValue, fuzzy, 0.0 ; If this is a CGTP projection, then check the axis range values. Otherwise, ; use the LL_BOX form the map structure !MAP. The latter is more strict about ; being inside the box. IF Keyword_Set(gctp) THEN BEGIN IF x LT (Min(!X.CRange)- fuzzy*(!X.CRange[1]-!X.CRange[0])) THEN RETURN, xy[0] IF x GT (Max(!X.CRange)+ fuzzy*(!X.CRange[1]-!X.CRange[0])) THEN RETURN, xy[0] IF y LT (Min(!Y.CRange)- fuzzy*(!Y.CRange[1]-!Y.CRange[0])) THEN RETURN, xy[0] IF y GT (Max(!Y.CRange)+ fuzzy*(!Y.CRange[1]-!Y.CRange[0])) THEN RETURN, xy[0] ENDIF ELSE BEGIN ENDELSE RETURN, xy end ;+-------------------------------------------------------------------------- ; Provides a MAP_GRID command that can be used in conjunction with other ; Coyote Graphics programs and in resizeable graphics windows. Any keyword ; appropriate for the MAP_GRID command in IDL can be used. ; ; :Keywords: ; addcmd: in, optional, type=boolean, default=0 ; If this keyword is set, the object is added to the resizeable graphics ; window, cgWindow. Note that a map projection command must be added to the ; window before this command is added to be effective. ; bcolor: in, optional, type=string, default='opposite' ; The name of the color to draw box axes with. ; bthick: in, optional, type=float ; The thickness of the box axes in device units. ; blabel: in, optional, type=float, default=0 ; Set this keyword to 0 to label all four of the box axes. ; Set this keyword to 1 to label the left and bottom box axes only. ; Set this keyword to 2 to label the top and right box axes only. ; box_axes: in, optional, type=boolean, default=0 ; Set this keyword to draw box axes on the map projection. ; cggrid: in, optional, type=boolean, default=0 ; Set this keyword to allow the latitude and longitude values to be set by ; the LatLon_Labels method in the cgMap object. Previously this was used by ; default, but it caused a lot of problems with global or near global map projections. ; This really should be used ONLY if you are mapping a very small region of the Earth, ; and maybe if you are using a UTM map projection. Othersize, it is probably not ; needed, so I have made it an optional choice. ; charsize: in, optional, type=float ; The character size for the labels. Default is cgDefCharSize()*0.75. ; clip_text: in, optional, type=boolean, default=1 ; Set this keyword to a zero value to turn off clipping of text labels. ; By default, text labels are clipped. This keyword is ignored if the `Box_Axes` ; keyword is set. ; color: in, optional, type=string, default='opposite' ; The name of the drawing color for the program. ; fill_horizon: in, optional, type=boolean, default=0 ; format: in, optional, type=string ; Set this keyword to a format for the grid labels. ; fuzzy: in, optional, type=float, default=0.0 ; This keyword applies only if the MAP_STRUCTURE keyword is used. Set the ; keyword to a value that is a percentage of the current data range. ; This percentage of the range is added to or subtracted from the ; values used to determine if the label is "inside" the boundary. ; It allows you to be a little less exact when selecting inside ; points. There are occasional aesthetic reasons for allowing fuzzy ; boundaries. A reasonable value for fuzziness might be 0.0125. ; glinestyle: in, optional ; Not used. Use `Linestyle` instead. ; glinethick: in, optional ; Not used. Use `Thick` instead. ; horizon: in, optional, type=boolean, default=0 ; Set this keyword to draw the current map horizon. ; increment: in, optional, type=integer, default=4 ; Determines the smoothness of graticle lines by setting the spacing between ; them. A low number is a smooth number. ; label: in, optional, type=integer, default=0 ; An number that tells how to label the grid line. A 0 means no grid lines are ; labelled. A 1 means all grid lines are labelled. A 2 means label every 2nd ; grid line is labelled. A 3 means every 3rd grid line is labelled, and so on. ; latalign: in, optional, type=float, default=0.0 ; 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. This keyword is ignored if the `Box_Axes` keyword is set. ; latdelta: in, optional, type=float ; Set this keyword to the spacing in degrees between parallels of latitude in ; the grid. If this keyword is not set, a suitable value is determined from the ; current map projection. ; latlabel: in, optional, type=float ; The longitude at which to place latitude labels. The default is the center ; longitude on the map if using Map_Set, and a longitude line on the left of the ; plot if using Map_Proj_Init. This keyword is ignored if the `Box_Axes` keyword is set. ; latnames: in, optional, type=string ; Set this keyword equal to an array specifying the names to be used for the ; latitude labels. By default, this array is automatically generated in units ; of degrees. The LATNAMES array can be either type string or any single numeric ; type, but should not be of mixed type. When LATNAMES is specified, the LATS ; keyword must also be specified. The number of elements in the two arrays need ; not be equal. If there are more elements in the LATNAMES array than in the LATS ; array, the extra LATNAMES are ignored. If there are more elements in the LATS ; array than in the LATNAMES array, labels in degrees will be automatically ; provided for the missing latitude labels. The LATNAMES keyword can be also used ; when the LATS keyword is set to a single value. It this case, the first label ; supplied will be used at the specified latitude; subsequent names will be ; placed at the next latitude line to the north, wrapping around the globe if ; appropriate. Caution should be used when using LATNAMES in conjunction with a ; single LATS value, since the number of visible latitude gridlines is dependent ; on many factors. ; lats: in, optional, type=float ; Set this keyword equal to a one or more element vector of latitudes for which ; lines will be drawn (and optionally labeled). If LATS is omitted, appropriate ; latitudes will be generated based on the value of the (optional) LATDEL keyword. ; If LATS is set to a single value, that latitude and a series of automatically ; generated latitudes will be drawn (and optionally labeled). Automatically generated ; latitudes have values that extend over the map. If LATS is a single value, that ; value is taken to be the starting point for labelling (See the LABEL keyword). ; lcolor: in, optional, type=string ; Set this to the name of the label color to use in labeling grid lines. ; By default, the same as COLOR, or if BOX_AXIS is set, then same as BCOLOR. ; linestyle: in, optional, type=integer, default=1 ; This keyword is the same as the GLineStyle keyword, but is a more ; natural way of setting the value. ; lonalign: in, optional, type=float, default=0.0 ; 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. This keyword is ignored if the `Box_Axes` keyword is set. ; londelta: in, optional, type=float ; Set this keyword to the spacing in degrees between parallels of longitude in ; the grid. If this keyword is not set, a suitable value is determined from the ; current map projection. ; lonlabel: in, optional, type=float ; The latitude at which to place longitude labels. The default is the center ; longitude on the map if using Map_Set, and a longitude line on the left of the ; plot if using Map_Proj_Init. This keyword is ignored if the `Box_Axes` keyword is set. ; lonnames: in, optional, type=string ; Set this keyword equal to an array specifying the names to be used for the ; longitude labels. By default, this array is automatically generated in units ; of degrees. The LATNAMES array can be either type string or any single numeric ; type, but should not be of mixed type. When LATNAMES is specified, the LATS ; keyword must also be specified. The number of elements in the two arrays need ; not be equal. If there are more elements in the LATNAMES array than in the LATS ; array, the extra LATNAMES are ignored. If there are more elements in the LATS ; array than in the LATNAMES array, labels in degrees will be automatically ; provided for the missing longitude labels. The LATNAMES keyword can be also used ; when the LATS keyword is set to a single value. It this case, the first label ; supplied will be used at the specified longitude; subsequent names will be ; placed at the next longitude line to the north, wrapping around the globe if ; appropriate. Caution should be used when using LATNAMES in conjunction with a ; single LATS value, since the number of visible longitude gridlines is dependent ; on many factors. ; lons: in, optional, type=float ; Set this keyword equal to a one or more element vector of longitudes for which ; lines will be drawn (and optionally labeled). If LATS is omitted, appropriate ; longitudes will be generated based on the value of the (optional) LATDEL keyword. ; If LATS is set to a single value, that longitude and a series of automatically ; generated longitudes will be drawn (and optionally labeled). Automatically generated ; longitudes have values that extend over the map. If LATS is a single value, that ; value is taken to be the starting point for labelling (See the LABEL keyword). ; map_structure: in, optional ; This keyword is normally used to pass in a map structure of the type ; created by Map_Proj_Init. In this version of the program, it can also ; be used to pass in a cgMap object, from which the map structure and other ; pertinent information for creating map grid lines can be obtained. ; noclip: in, optional, type=boolean, default=0 ; Normally, output is clipped to the map projection boundaries. Set this keyword ; to be able to draw outside the map boundaries. ; no_grid: in, optional, type=boolean, default=0 ; Set this keyword if you only want labels but not grid lines. ; orientation: in, optional, type=float, default=0.0 ; Set this keyword equal to an angle in degrees from horizontal (in the clockwise ; direction) to rotate the labels. This keyword is ignored if the `Box_Axes` keyword is set. ; t3d: in, optional, type=boolean, default=0 ; Set this keyword if you are labeling in 3D space. ; thick: in, optional, type=integer, default=!P.Thick ; Sets the thickness of the grid lines. ; zvalue: in, optional, type=float ; Draw the grid lines at this Z-value. Implies the use of `T3D`. ; ;--------------------------------------------------------------------------- PRO cgMap_Grid, $ ADDCMD=addcmd, $ BCOLOR=bcolor, $ BLABEL=blabel, $ BTHICK=bthick, $ BOX_AXES=box_axes, $ CGGRID=cggrid, $ CHARSIZE=charsize, $ CLIP_TEXT=clip_text, $ COLOR=scolor, $ FILL_HORIZON=fill_horizon, $ FORMAT=format, $ FUZZY=fuzzy, $ GLINESTYLE=glinestyle, $ GLINETHICK=glinethick, $ HORIZON=horizon, $ INCREMENT=increment, $ LABEL=label, $ LATALIGN=latalign, $ LATDELTA = latdelta, $ LATLABEL=latlab, $ LATNAMES=latnames, $ LATS=lats, $ LCOLOR=lcolor, $ LINESTYLE=linestyle, $ LONALIGN=lonalign, $ LONDELTA=londelta, $ LONLABEL=lonlab, $ LONNAMES=lonnames, $ LONS=lons, $ MAP_STRUCTURE=mapStruct, $ NOCLIP=noclip, $ NO_GRID=no_grid, $ ORIENTATION=orientation, $ T3D=t3d, $ THICK=thick, $ ZVALUE=zvalue Compile_Opt strictarr Catch, theError IF theError NE 0 THEN BEGIN Catch, /CANCEL void = cgErrorMsg() IF N_Elements(thisState) NE 0 THEN cgSetColorState, thisState RETURN ENDIF ; Should this be added to a resizeable graphics window? IF (Keyword_Set(addcmd)) && ((!D.Flags && 256) NE 0) THEN BEGIN cgWindow, 'cgMap_Grid', $ BLABEL=blabel, $ BTHICK=bthick, $ BOX_AXES=box_axes, $ CHARSIZE=charsize, $ CGGRID=cgGrid, $ CLIP_TEXT=clip_text, $ COLOR=scolor, $ FILL_HORIZON=fill_horizon, $ FORMAT=format, $ FUZZY=fuzzy, $ GLINESTYLE=glinestyle, $ GLINETHICK=glinethick, $ HORIZON=horizon, $ INCREMENT=increment, $ LABEL=label, $ LATALIGN=latalign, $ LATDELTA = latdelta, $ LATLABEL=latlab, $ LATNAMES=latnames, $ LATS=lats, $ LINESTYLE=linestyle, $ LONALIGN=lonalign, $ LONDELTA=londelta, $ LONLABEL=lonlab, $ LONNAMES=lonnames, $ LONS=lons, $ MAP_STRUCTURE=mapStruct, $ NOCLIP=noclip, $ NO_GRID=no_grid, $ ORIENTATION=orientation, $ T3D=t3d, $ THICK=thick, $ WHOLE_MAP=obsolete_keyword, $ ZVALUE=zvalue, $ ADDCMD=1 RETURN ENDIF ; Set up PostScript device for working with colors. IF !D.Name EQ 'PS' THEN Device, COLOR=1, BITS_PER_PIXEL=8 SetDefaultValue, charsize, cgDefCharsize() * 0.75 SetDefaultValue, bcolor, "opposite" SetDefaultValue, blabel, 0 noclip = Keyword_Set(noclip) ; I want to use the more natural LINESTYLE and THICK keywords to this routine. IF (N_Elements(linestyle) EQ 0) && (N_Elements(glinestyle) NE 0) THEN BEGIN linestyle = glinestyle ENDIF ELSE IF (N_Elements(linestyle) EQ 0) THEN linestyle = 1 ; dotted IF (N_Elements(thick) EQ 0) && (N_Elements(gthick) NE 0) THEN BEGIN thick = gthick ENDIF ELSE IF (N_Elements(thick) EQ 0) THEN thick = !P.Thick ; Try to do this in decomposed color, if possible. cgSetColorState, 1, Current=thisState ; Need a color. IF N_Elements(scolor) NE 0 THEN BEGIN CASE Size(scolor, /TNAME) OF 'STRING': 'LONG': 'BYTE': BEGIN IF N_Elements(scolor) NE 3 THEN scolor = StrTrim(Fix(scolor), 2) END ELSE: scolor = StrTrim(scolor,2) ENDCASE ENDIF ELSE scolor = "opposite" IF N_Elements(scolor) EQ 0 THEN color = !P.Color ELSE color = sColor IF (Size(scolor, /TNAME) EQ 'BYTE') AND (N_Elements(scolor) EQ 3) THEN color = cgColor(scolor) IF Size(color, /TYPE) EQ 3 THEN IF cgGetColorState() EQ 0 THEN color = Byte(color) IF Size(color, /TYPE) LE 2 THEN color = StrTrim(Fix(color),2) ; Check for label color now. Depends on other colors and value of BOX_AXES. IF N_Elements(lcolor) EQ 0 THEN lcolor = (Keyword_Set(box_axes)) ? bcolor : color ; Determine if you have a map structure or a map object. IF N_Elements(mapStruct) EQ 0 THEN BEGIN hasMap = 0 hasMapObj = 0 ENDIF ELSE BEGIN IF Obj_Valid(mapStruct) THEN BEGIN hasMapObj = 1 mapObj = mapStruct mapObj -> Draw, /NoGraphics thisMapStruct = mapObj -> GetMapStruct() ; This routine is here because Map_Grid does not select good line for small areas. IF Keyword_Set(cgGrid) THEN BEGIN mapObj -> LatLonLabels, LATS=mlats, LATLAB=mlatlab, LATDELTA=latdelta, LATNAMES=mlatnames, $ LONS=mlons, LONLAB=mlonlab, LONDELTA=londelta, LONNAMES=mlonnames, $ FORMAT=format IF N_Elements(lats) EQ 0 THEN BEGIN lats = mlats IF (N_Elements(latnames) EQ 0) THEN latnames = mlatnames IF N_Elements(latlab) EQ 0 THEN latlab = mlatlab ENDIF IF N_Elements(lons) EQ 0 THEN BEGIN lons = mlons IF (N_Elements(lonnames) EQ 0) THEN lonnames = mlonnames IF N_Elements(lonlab) EQ 0 THEN lonlab = mlonlab ENDIF ENDIF ENDIF ELSE BEGIN hasMapObj = 0 thisMapStruct = mapStruct ENDELSE hasMap = 1 ENDELSE if ((!x.type NE 3) && ~hasMap) THEN $ message, 'cgMap_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. IF ~noclip THEN BEGIN noclip = (N_ELEMENTS(clip_text) gt 0) ? ~KEYWORD_SET(clip_text) : 0 ENDIF ; Append the graphics keywords: if n_elements(t3d) then map_struct_append, extra,'T3D',t3d ;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 latdelta or londelta 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,'cgMap_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,'cgMap_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 ? thisMapStruct : !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(latdelta) eq 0 THEN begin latdelta = cgMap_Grid_incr(latmax - latmin) latd = 1 endif else latd = latdelta IF n_elements(londelta) eq 0 THEN begin londelta = cgMap_Grid_incr(lonmax - lonmin) lond = 1 endif else lond = londelta ; 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? is_cyl = 0 if (hasMap) then begin IF Obj_Valid(mapObj) THEN is_cyl = mapObj -> Is_Cylindrical() endif else begin map_proj_info, iproj, CYLINDRICAL=is_cyl, /CURRENT endelse 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. IF lonmin EQ lonmax THEN lonmax = lonmin + 360. 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) / londelta) longitudes = lonmin - (lonmin mod londelta) + findgen(n_lons) * londelta endif else if n_elements(lons) eq 1 then begin i0 = ceil((lonmin - lons[0]) / float(londelta)) ;First tick i1 = floor((lonmax - lons[0]) / float(londelta)) ;Last tick n_lons = i1 - i0 + 1 > 1 longitudes = (findgen(n_lons) + i0) * londelta + 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(latdelta)) ;1st lat for grid n_lats = 1 + fix((latmax-lat0)/float(latdelta)) latitudes = lat0 + findgen(n_lats)*latdelta endif else if n_elements(lats) eq 1 then begin i0 = ceil((latmin - lats[0]) / float(latdelta)) ;First tick i1 = floor((latmax - lats[0]) / float(latdelta)) ;Last tick n_lats = i1 - i0 + 1 > 1 latitudes = (findgen(n_lats) + i0) * latdelta + 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) ne 0 then dc = dc * charsize IF hasMapObj THEN mapObj -> Draw, /NoGraphics, Erase=0 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 IF hasMapObj THEN BEGIN mapObj -> GetProperty, XRANGE=xrange, YRANGE=yrange xww = (xrange * !x.s[1] + !x.s[0]) * !d.x_size + del yww = (yrange * !y.s[1] + !y.s[0]) * !d.y_size + del ENDIF ELSE BEGIN 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 ENDELSE IF (N_Elements(bthick) NE 0) THEN bdel = bthick ELSE bdel = box_thick * !d.y_px_cm ;Thickness of box in device units 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 cgPlotS, xw[[0,1,1,0,0]], yw[[0,0,1,1,0]], /DEVICE, $ COLOR=bcolor, LINESTYLE=0, THICK=thick cgPlotS, 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, LINESTYLE=0, THICK=thick 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 fraught 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. hcolor = (Size(color, /TNAME) EQ 'STRING') ? cgColor(color) : color if keyword_set(horizon) then map_horizon, COLOR=hcolor, _EXTRA=e if keyword_set(fill_horizon) then map_horizon, COLOR=hcolor, _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=thisMapStruct, POLYLINES=polylines) ;;; DWF correction next two commands. uv = MAP_PROJ_FORWARD(REPLICATE(lon, N_ELEMENTS(lati)), lati, $ MAP_STRUCTURE=thisMapStruct) 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] cgPlotS, REFORM(uv[0,indices]), REFORM(uv[1,indices]), $ zvalue, $ NOCLIP=noclip, $ COLOR=color, LINESTYLE=linestyle, THICK=thick, $ _EXTRA=extra endif index += nline + 1 endwhile endif else begin cgPlotS, lon, lati, zvalue, NOCLIP=noclip, $ COLOR=color, LINESTYLE=linestyle, THICK=thick, _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=thisMapStruct) if (FINITE(uv[0]) && FINITE(uv[1])) then begin xy = uv[0:1] gctp = 1 endif endif else begin if (noclip eq 1) || map_point_valid(lon, LonLab) then $ xy = [lon, LonLab] endelse IF N_Elements(xy) NE 0 THEN xy = cgMap_Grid_check_range(xy, GCTP=gctp, FUZZY=fuzzy) IF N_Elements(xy) EQ 2 THEN BEGIN cgText, xy[0], xy[1], lonname, ALIGNMENT=lonalign, $ NOCLIP=1, Z=zvalue, COLOR=lcolor, CHARSIZE=charsize, _EXTRA=extra Undefine, xy ENDIF 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] = cgMap_Grid_solve( $ [xww[0], yw[j]+k*dy], $ [xww[1], yw[j]+k*dy], 0, lon, $ MAP_STRUCTURE=thisMapStruct) 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 begin latname=strtrim(string(lat, format=fmt),2) endelse 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=thisMapStruct) if (FINITE(uv[0]) && FINITE(uv[1])) then begin xy = uv[0:1] gctp = 1 endif endif else begin if (noclip eq 1) || map_point_valid(latlab, lat) then $ xy = [latlab, lat] endelse IF N_Elements(xy) NE 0 THEN xy = cgMap_Grid_check_range(xy, GCTP=gctp, FUZZY=fuzzy) IF N_Elements(xy) EQ 2 THEN BEGIN cgText, xy[0], xy[1], latname, CHARSIZE=charsize, $ alignment=latalign, NOCLIP=1, COLOR=lcolor, Z=zvalue, _EXTRA=extra Undefine, xy ENDIF endif ENDIF IF Abs(lat) EQ 90 THEN BEGIN oldlat = lat lat = -89.975 > lat < 89.975 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 IF (hasMapObj) THEN BEGIN ; Added this because if there is not a point in the longitude vector ; that is inside our actual data range, then the line can be drawn ; at the wrong place. I don't really understand why. But this is a ; check to be sure there is a point inside the data range. mapObj -> GetProperty, XRANGE=xrange, YRANGE=yrange index = Value_Locate(xrange, loni) ; If there is a zero in index vector than there is a point inside range ; and we continue normally. Otherwise, we try to construct a vector that ; has a point in the data range. void = Where(index EQ 0, count) IF count EQ 0 THEN BEGIN ; No point inside range, so construct one. thisMapStructure = mapObj -> GetMapStruct() ll = Map_Proj_Inverse(xrange, yrange, MAP_STRUCTURE=thisMapStructure) lonii = Reform(ll[0,*]) lonii = cgScaleVector(Findgen(N_Elements(loni)), lonii[0], lonii[1]) uv = MAP_PROJ_FORWARD(lonii, REPLICATE(lat, N_ELEMENTS(lonii)), $ MAP_STRUCTURE=thisMapStruct, POLYLINES=polylines) ; uv = MAP_PROJ_FORWARD(lonii, REPLICATE(lat, N_ELEMENTS(lonii)), $ ; MAP_STRUCTURE=thisMapStruct) ;;; Screws up near the poles, I think. ; polylines = [N_Elements(loni), Indgen(N_Elements(loni))] loni = lonii ENDIF ELSE BEGIN uv = MAP_PROJ_FORWARD(loni, REPLICATE(lat, N_ELEMENTS(loni)), $ MAP_STRUCTURE=thisMapStruct, POLYLINES=polylines) ; uv = MAP_PROJ_FORWARD(loni, REPLICATE(lat, N_ELEMENTS(loni)), $ ; MAP_STRUCTURE=thisMapStruct) ;;; Screws up near the poles, I think. ; polylines = [N_Elements(loni), Indgen(N_Elements(loni))] ENDELSE ENDIF ELSE BEGIN uv = MAP_PROJ_FORWARD(loni, REPLICATE(lat, N_ELEMENTS(loni)), $ MAP_STRUCTURE=thisMapStruct, POLYLINES=polylines) ; uv = MAP_PROJ_FORWARD(loni, REPLICATE(lat, N_ELEMENTS(lonii)), $ ; MAP_STRUCTURE=thisMapStruct) ;;; Screws up near the poles, I think. ; polylines = [N_Elements(loni), Indgen(N_Elements(loni))] ENDELSE 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] cgPLOTS, REFORM(uv[0,indices]), REFORM(uv[1,indices]), $ zvalue, $ NOCLIP=noclip, $ COLOR=color, LINESTYLE=linestyle, THICK=thick, $ _EXTRA=extra endif index += nline + 1 endwhile endif else $ cgPLOTS, loni, lat, zvalue, NOCLIP=noclip, $ COLOR=color, LINESTYLE=linestyle, THICK=thick, _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] = cgMap_Grid_Solve( $ [xw[j]+dx*k, yww[0]], $ [xw[j]+dx*k, yww[1]], 1, lat, $ MAP_STRUCTURE=thisMapStruct) 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 ; iaxis = 0 for lon axis, 1 for lat axis 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 $ cgColorFill, xp0, yp0, /DEVICE, COLOR=bcolor xychar[iaxis] = z if strlen(vtext[i]) gt 0 then begin CASE blabel OF 0: cgText, xychar[0], xychar[1], vtext[i], $ ORIENTATION=dy * (90-180*j), CHARSIZE=charsize, $ ALIGN=0.5, CLIP=0, Z=zvalue, COLOR=lcolor, /DEVICE, _EXTRA=extra 1: IF (j EQ 0) THEN cgText, xychar[0], xychar[1], vtext[i], $ ORIENTATION=dy * (90-180*j), CHARSIZE=charsize, $ ALIGN=0.5, CLIP=0, Z=zvalue, COLOR=lcolor, /DEVICE, _EXTRA=extra 2: IF (j EQ 1) THEN cgText, xychar[0], xychar[1], vtext[i], $ ORIENTATION=dy * (90-180*j), CHARSIZE=charsize, $ ALIGN=0.5, CLIP=0, Z=zvalue, COLOR=lcolor, /DEVICE, _EXTRA=extra ENDCASE 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]] cgColorFill, xp0, yp0, /DEVICE, COLOR=bcolor endif endfor endif ; box_thick ; Restore color mode cgSetColorState, thisState end