;+ ; NAME: ; CONTOUR_ON_MAP ; ; PURPOSE: ; ; The purpose of this procedure is to provide an ; easy way to plot color contour plots on a map ; projection with an attached colorbar. ; ; AUTHOR: ; ; Nathan Johnson ; University of Arizona ; Department of Atmospheric Sciences ; johnson AT atmo DOT arizona DOT edu ; ; CATEGORY: ; ; Graphics ; ; CALLING SEQUENCE (WITH REQUIRED INPUT PARAMETERS): ; ; CONTOUR_ON_MAP, x, y, image_in : x is an array of the longitudes ; y is an array of the latitudes ; image_in is a 2D array containing the data that ; will be plotted on the map. ; The number of elements in x must be the same as the ; number of elements in the first dimension of image_in. ; The number of elements in y must be the same as the ; number of elements in the second dimension of image_in. ; ; OPTIONAL INPUT PARAMETERS: ; ; p0lat, p0lon, rot : These are the same optional keywords that can be called with ; the MAP_SET procedure. See the online help for details. ; ; INPUT KEYWORD PARAMETERS: ; ; Too many to list. All of the input keywords for the map_set ; procedure are allowed. Others include: ; ; ct : The numerical value of the color table you wish to display on the map. ; Default is B-W linear (grey scale). ; min_value : The minimum value to contour. If not set, use the minimum value of image_in. ; max_value : The maximum value to contour. If not set, use the maximum value of image_in. ; print : Set this keyword to print to the postscript device. ; filename : If the print keyword is set, this is the filename to which it will be saved. ; wrap : Makes plots that span 360 degrees of longitude look nicer. ; latrange : A 2 element array with the minimum and maximum latitude to plot, e.g. [0,90] ; lonrange : A 2 element array with the minimum and maximum longitude to plot, e.g. [-180, 180] ; colorbar_min : The minimum value to display on the colorbar. If not set, use the value in min_value. ; colorbar_max : The minimum value to display on the colorbar. If not set, use the value in min_value. ; cbar_title : Title to use for the colorbar. ; colorstep : Use this keyword to control the step size used in the current color table. A larger step ; will have larger 'distance' between the colors, and likely to have more distinct colors. ; coloroffset : The offset to use for the color table. Should be at least 1 (default) since most color tables ; have black as the first element. ; ncolors : The number of colors you want in your colortable. Also the same as the number of contours ; in the plot. Default is 12. ; xsize : Size of horizontal scale in pixels. Only used when plotting to the screen. ; ysize : Size of vertical scale in pixels. Only used when plotting to the screen. ; format : Format to use for the colorbar. Default is two decimal places. ; cbar_title : Title to use for the colorbar. ; show_contours: Will show the contours on the map ; c_position : 4 Element array specifying the endpoints of the position of the colorbar. ; position : 4 Element array used to specify the endpoints of the map. ; ; OUTPUT KEYWORD PARAMETERS: ; ; None. ; ; COMMON BLOCKS: ; None. ; ; SIDE EFFECTS: ; None known. ; ; RESTRICTIONS: ; ; Required programs from the Coyote Library: ; (Written by David Fanning) ; ; http://www.dfanning.com/programs/error_message.pro ; http://www.dfanning.com/programs/pickcolorname.pro ; http://www.dfanning.com/programs/fsc_color.pro ; http://www.dfanning.com/programs/colorbar.pro ; ; I suggest downloading the entire library from: ; http://www.dfanning.com/documents/programs.html ; ; EXAMPLE: ; ; contour_on_map, x, y, image_in, 90, -90, /wrap, /orthographic, ct=3 ; ; This command will contour image_in on an orthographic projection centered ; at a latitude of 90 degrees (North Pole). The color table loaded is 3, ; which gives a 'red temperature' shading to the data. ; ; ; MODIFICATION HISTORY: ; ; Written by: Nathan Johnson, 14 November 2007. ; This is a complete rewrite of unreleased code that performed the ; same function, but was less user-friendly. This procedure is ; entirely backward compatible with the previous version. ;- PRO contour_on_map, x, y, image_in, $ ; these variables are needed for all map projections ; x is the longitude, y is latitude, and image_in is the data to be contoured p0lat, p0lon, rot, $ ; some map projections require center points and rotations. see the map_set ; procedure for details. These default to zero if not set. ; ; ********** some keywords to make the plot look nice ct=ct, $ ;color table min_value=min_value, max_value=max_value, $ ; max and min values to contour print=print, filename=filename, $ ; is this to be printed to a file? if so, to what file? wrap=wrap, $ ; makes full longitude plots look prettier latrange=latrange, lonrange=lonrange, $ ; latitude and longitude ranges e.g. [0,90] or [-180, 180] colorbar_min=colorbar_min, colorbar_max=colorbar_max, $ ; min and max ranges for the colorbar - if not set, use the min and max values colorstep=colorstep, $ ; the colorset, coloroffset, and ncolors keywords give the user coloroffset=coloroffset, $ ; almost complete control over what color table to display on the ncolors=ncolors, $ ; map. Future versions may allow users to specify their own color table. xsize=xsize, ysize=ysize, $ ; Controls the size of window when plotting to the screen. format=format, $ ; Format to use for the colorbar. cbar_title=cbar_title, $ ; Title to use for the colorbar. show_contours=show_contours, $ ; Will show the contours on the map c_position=c_position, $ ; 4 Element array specifying the endpoints of the position of the colorbar. ; The POSITION keyword is used to specify the endpoints of the map. color_array=color_array, $ ; TEST TEST TEST -- fill in when working cont_color=cont_color, $ ; color of the continents cbmin_zero=cbmin_zero, $ ; use this keyword if you want the minimum value in the colorbar to be zero ; ; ********** Projection keywords, these are also used by map_set PROJECTION = proj, NAME = name, STEREOGRAPHIC = stereographic, ORTHOGRAPHIC = orthographic, CONIC = conic, $ LAMBERT = lambert, GNOMIC = gnomic, AZIMUTHAL = azimuth, SATELLITE = satellite, CYLINDRICAL = cylindrical, $ MERCATOR = mercator, MILLER_CYLINDRICAL=miller, MOLLWEIDE = mollweide, SINUSOIDAL = sinusoidal, AITOFF = aitoff, $ HAMMER = hammer, ALBERS = albers, TRANSVERSE_MERCATOR = utm, ROBINSON = robinson, GOODESHOMOLOSINE = goodes, $ ; ; **** Projection specific keywords, see map_set documentation for details ELLIPSOID = ellips, CENTRAL_AZIMUTH=cent_azim, STANDARD_PARALLELS = std_p, SAT_P = Sat_p, $ ; ; ********** MAP_SET specific keywords: CLIP=clip, REVERSE=reverse, SCALE=scale, ISOTROPIC = iso, LIMIT = limit, $ ; ; ********** MAP_SET graphics keywords: NOERASE=noerase, TITLE=title, ADVANCE = advance, COLOR=color, POSITION = position, NOBORDER=noborder, T3D=t3d, $ ZVALUE=zvalue, CHARSIZE = charsize, XMARGIN=xmargin, YMARGIN=ymargin, $ ; ; ********** MAP_HORIZON keywords: HORIZON=horizon, E_HORIZON=ehorizon, $ ; ; ********** MAP_CONTINENTS keywords: CONTINENTS = continents, E_CONTINENTS=econt, USA=usa, HIRES = hires, $ MLINESTYLE=mlinestyle, MLINETHICK=mlinethick, CON_COLOR=con_color, $ ; ; ********** MAP_GRID keywords: GRID=grid, E_GRID=egrid, GLINESTYLE=glinestyle, GLINETHICK=glinethick, $ LABEL=label, LATALIGN=latalign, LATDEL=latdel, LATLAB=latlab, $ LONALIGN=lonalign, LONDEL=londel, LONLAB=lonlab, $ ; Ignored, but here for compatibility: WHOLE_MAP=whole_map ;- ; on_error, 2 if max(y) gt 90 or min(y) lt -90 then message, 'Latitude is out of allowed range.' ; Remember current device thisdevice = !d.name device, get_decomposed = get_decomposed ; Check for correct number of parameters, and fail here instead of somewhere else. if N_PARAMS() lt 3 then message, 'This procedure requires at least 3 parameters: x, y, image_in' if (size(image_in))(0) ne 2 then message, 'image_in must be a two dimensional array.' ; original procedure call to contour_on_map -> make sure it's backwards compatible. It is. ;PRO contour_on_map_new, x, y, image, ct=ct, min_value=min_value, max_value=max_value, filename=filename, wrap=wrap, latrange=latrange, lonrange=lonrange, print=print image = image_in ; don't mess up input variable ; Set default values if not specified. if not keyword_set(colorbar_min) then begin if keyword_set(min_value) then colorbar_min = double(min_value) else colorbar_min = double(min(image, /nan)) if keyword_set(cbmin_zero) then colorbar_min = 0 endif if not keyword_set(colorbar_max) then $ if keyword_set(max_value) then colorbar_max = double(max_value) else colorbar_max = double(max(image, /nan)) if not keyword_set(min_value) then min_value=-999.9 if not keyword_set(max_value) then max_value=32766 ; set values lower then the minimum to the minimum and higher than maximum to the maximum s = where(image le colorbar_min, count) if count ne 0 then image(s) = colorbar_min s = where(image ge colorbar_max, count) if count ne 0 then image(s) = colorbar_max if not keyword_set(ct) then ct = 0 if not keyword_set(filename) then filename = 'a.eps' if not keyword_set(latrange) then latrange=[-70,70] if not keyword_set(lonrange) then lonrange=[-180,180] if not keyword_set(limit) then LIMIT=[latrange[0],lonrange[0],latrange[1],lonrange[1]] if not keyword_set(position) then POSITION=[0.15, 0.15, 0.95, 0.75] if not keyword_set(c_position) then c_position=[0.15, 0.90, 0.95, 0.95] if not keyword_set(colorstep) then colorstep = 16. if not keyword_set(coloroffset) then coloroffset = 1. if coloroffset le 1 then coloroffset = 1. if not keyword_set(ncolors) then ncolors = 12. if not keyword_set(format) then format = '(F0.2)' if not keyword_set(cont_color) then cont_color='black' if (size(cont_color))(1) ne 7 then cont_color='black' ; ; The /WRAP keyword provides an easy mechanism to wrap the data around the globe. If this keyword is set for data that does ; not span a longitude range of 360 degrees it will seriously mess up the contour plot. However, it will not return an error. ; This little procedure was written before I secured the integrity of the input array. It should probably be reduced in complexity now. if keyword_set(wrap) then begin sze = size(image) newimage = fltarr(sze[1]+1, sze[2]) newimage[0:sze[1]-1,*] = image newimage[sze[1],*] = image[0,*] newx = fltarr(N_ELEMENTS(x) + 1) newx[0:sze[1]-1] = x newx[sze[1]] = x[0] endif else begin newimage = image newx = x endelse ; Set the device based upon the /PRINT keyword if keyword_set(print) then begin set_plot, 'ps' !p.font=0 device, filename=filename, /color endif else begin if thisdevice eq 'PS' then begin ; this will only work on linux (including new macs) and windows machines! if !version.os eq 'linux' then set_plot, 'X' else set_plot, 'win' endif else begin set_plot, thisdevice endelse !p.font=-1 device, decomposed=0 window, xsize=xsize, ysize=ysize endelse ; Load the color table. if not keyword_set(color_array) then loadct, ct, /silent else $ tvlct, color_array[0,*], color_array[1,*], color_array[2,*] ; The first 3 parameters are needed. If any of the last three are missing, default to zeros. if n_params() lt 6 then rot = 0.0d0 if n_params() lt 5 then p0lon = 0d0 if n_params() lt 4 then p0lat = 0d0 ; Make the background color white, instead of whatever IDL wants to use. Erase, COLOR=FSC_COLOR('white') ; Don't erase the newly created white background when calling map_set ; Note: this keyword, although accepted by the contour_on_map ; procedure, will not be used. noerase=1 ; Call the MAP_SET procedure... with all its keywords MAP_SET, p0lat, p0lon, rot, PROJECTION=proj, NAME=name, STEREOGRAPHIC = stereographic, ORTHOGRAPHIC = orthographic, CONIC = conic, $ LAMBERT = lambert, GNOMIC = gnomic, AZIMUTHAL = azimuth, SATELLITE = satellite, CYLINDRICAL = cylindrical, MERCATOR = mercator, $ MILLER_CYLINDRICAL=miller, MOLLWEIDE = mollweide, SINUSOIDAL = sinusoidal, AITOFF = aitoff, HAMMER = hammer, ALBERS = albers, $ TRANSVERSE_MERCATOR = utm, ROBINSON = robinson, GOODESHOMOLOSINE = goodes, ELLIPSOID = ellips, CENTRAL_AZIMUTH=cent_azim, $ STANDARD_PARALLELS = std_p, SAT_P = Sat_p, CLIP=clip, REVERSE=reverse, SCALE=scale, ISOTROPIC = iso, LIMIT = limit, NOERASE=noerase, $ TITLE=title, ADVANCE = advance, COLOR=color, POSITION = position, NOBORDER=noborder, T3D=t3d, ZVALUE=zvalue, CHARSIZE = charsize, $ XMARGIN=xmargin, YMARGIN=ymargin, HORIZON=horizon, E_HORIZON=ehorizon, CONTINENTS = continents, E_CONTINENTS=econt, USA=usa, $ HIRES = hires, MLINESTYLE=mlinestyle, MLINETHICK=mlinethick, CON_COLOR=con_color, GRID=grid, E_GRID=egrid, GLINESTYLE=glinestyle, $ GLINETHICK=glinethick, LABEL=label, LATALIGN=latalign, LATDEL=latdel, LATLAB=latlab, LONALIGN=lonalign, LONDEL=londel, LONLAB=lonlab, $ WHOLE_MAP=whole_map ; Manually figure out the steps and levels because the NLEVELS keyword is just a 'suggestion'. levels = double(ncolors) step = (colorbar_max - colorbar_min) / levels userLevels = dIndGen(levels) * step + colorbar_min ; c_colors = indgen(ncolors)*colorstep+coloroffset c_colors = (indgen(ncolors)+coloroffset)*colorstep ; Do the actual contour. Put some continents and the US states on the map. contour, newimage, newx, y, /overplot, min_value=colorbar_min, c_colors=c_colors, /cell_fill, levels=userlevels if keyword_set(show_contours) then $ contour, newimage, newx, y, /overplot, min_value=colorbar_min, levels=userlevels, c_labels=fltarr(ncolors)+1, c_colors=fsc_color('white') ; These two commands should probably be optional. Feel free to comment them out. MAP_GRID, LABEL=1, /BOX_AXES, Color=fsc_color('black') MAP_CONTINENTS, /COASTS, /CONTINENTS, /COUNTRIES, /USA, Color=fsc_color(cont_color), LIMIT=limit ; Get the correct color to put in the colorbar based on the device specified if not keyword_set(print) then begin TVLCT, R, G, B, /GET TVLCT, [[r[c_colors]], [g[c_colors]], [b[c_colors]]], 1 x_colorbar_colors = ncolors x_bottom = coloroffset endif else begin x_colorbar_colors = ncolors * colorstep x_bottom = coloroffset * colorstep endelse ;if not keyword_set(color_array) then loadct, ct else $ ; tvlct, color_array[0,*], color_array[1,*], color_array[2,*] ColorBar, minrange=colorbar_min, maxrange=colorbar_max, NCOLORS=x_colorbar_colors, bottom=x_bottom, FORMAT=format, $ POSITION=c_position, XMinor=1, XTicklen=1, AnnotateColor='black', Color=255, title=cbar_title, divisions = ncolors ; Clean up. Close the postscript device if necessary, and set the plot back to ; whatever it was set to before this procedure was called. if keyword_set(print) then device, /close set_plot, thisdevice device, decomposed=get_decomposed END