; $Id: map_set.pro,v 1.93 2001/03/14 15:43:57 dave Exp $ ; ; Copyright (c) 1993-2001, Research Systems, Inc. All rights reserved. ; Unauthorized reproduction prohibited. ; MAP_SET ; ; KEYWORDS: (partial list) ; ; NAME: Set this keyword to a string indicating the projection ; that you wish to use. A list of available projections can ; be found using MAP_PROJ_INFO, PROJ_NAMES=names. ; This keyword will override any of the individual projection ; keywords. ; ; PROJECTION: Set this keyword to a projection number indicating the ; projection that you wish to use. This keyword will override NAME ; or any of the individual projection keywords. ; ; Note - This keyword is intended for internal use only, and is subject ; to change in the future. ; ; REVERSE: Set this keyword to reverse the X and/or Y axes: ; 0 = no reversal (the default) ; 1 = reverse X ; 2 = reverse Y ; 3 = reverse both. PRO map_struct_append, In_struct, Name, Value, SUPERCEDE = super ; ; Append/Replace the tagname Name, with a given value to the In_Struct. ; If SUPERCEDE is set, and a tag with the given Name already exists, ; replace its value. ; If SUPERCEDE is NOT set, and the given tag exists, do nothing. ; COMPILE_OPT hidden ntags=N_TAGS(In_struct) if ntags eq 0 then begin ;If In_struct is undef, just return the tag In_struct = CREATE_STRUCT(name, value) endif else begin tnames=tag_names(IN_STRUCT) index=where(tnames eq NAME, count) if count eq 0 then $ ;No match, append In_Struct = CREATE_STRUCT(in_struct, name, value) $ else if keyword_set(super) then $ in_struct.(index[0]) = value ;Overwrite value? ;Otherwise, tag is already there, don't add endelse end PRO MAP_STRUCT_MERGE, orig, add, SUPERCEDE = super ; Add the keywords in the struct add to the original struct orig. ; ADD must be either undefined or defined as a structure. ; If a given tag name is present in both structures: ; If SUPERCEDE is set, then the value of add.tag replaces the ; original value. If SUPERCEDE is not set, the original value ; remains unchanged. COMPILE_OPT hidden if n_elements(orig) eq 0 then begin if n_elements(add) ne 0 then orig = add return endif otag = tag_names(orig) atag = tag_names(add) for i=0, n_elements(atag)-1 do begin ;Add the new ones match = where(otag eq atag[i], count) if count eq 0 then $ orig = create_struct(orig, atag[i], add.(i)) $ ;Disjoint tag name, add it else if keyword_set(super) then $ ;Got a match, if super is set, overwrite. orig.(match[0]) = add.(i) ;Common, add if set endfor end FUNCTION map_point_valid, lon, lat, u, v ; Return 1 if the point at geo-coordinates (lon, lat) is mappable, ; 0 if not. ; Mappable means that the point is transformable, and is not clipped ; by the mapping pipeline. ; Input: lon, lat = input coordinates in degrees. ; Output: u, v = uv coordinates if point is mappable. The contents of ; u and v are undefined if the function returns 0, indicating that the ; point is not mappable. ; COMPILE_OPT hidden u = 0.0 & v = 0.0 p = convert_coord(lon, lat, /DATA, /TO_NORM) ;Cvt to normalized coords if finite(p[0]) eq 0 then return, 0 u = (p[0] - !x.s[0]) / !x.s[1] ;To UV coords v = (p[1] - !y.s[0]) / !y.s[1] clat = cos(lat * !dtor) ;Cvt lon lat to cartesian 3D coords xyz = [ clat * cos(lon * !dtor), clat * sin(lon * !dtor), sin(lat * !dtor)] p = !map.pipeline s = size(p) for i=0, s[2]-1 do begin ;# of stages case p[0,i] of 0 : return, 1 ;0 = END, means we're successful. 1 : dummy = 0 ;Ignore splits 2 : if total(xyz * p[1:3,i]) + p[4,i] lt 0 then $ ;Outside clip plane? return, 0 3 : dummy = 0 ;Ignore transforms 4 : if u lt p[1,i] or v lt p[2,i] or $ ;Within UV bounding box? u gt p[3,i] or v gt p[4,i] then return, 0 endcase endfor return, 1 ;Should never hit here. end ;map_point_valid PRO map_set_ll_box COMPILE_OPT hidden ;Try to determine a lon/lat range from the existing parameters, and ;save in !MAP.LL_BOX = [latmin, lonmin, latmax, lonmax]. ; If we can't reliably determine the lat/lon limits then the ; coorresponding min and max are both set to zero. ; If either pole is mappable, then all longitudes are mappable north = map_point_valid(0, 90) south = map_point_valid(0, -90) dateline = (north + south eq 0) and map_point_valid(-180., !map.p0lat) and $ (map_point_valid(0., !map.p0lat) eq 0) ; Dateline visible if north and south then begin ;Both poles visible sum = 0 for i=0,270,90 do sum = sum + map_point_valid(!map.p0lon+i, 0) if sum eq 4 then begin ;Entire globe is visible lonmin = !map.p0lon - 180 ;Center window about projection center while lonmin lt (-180.) do lonmin = lonmin + 360. !map.LL_BOX = [-90, lonmin, 90, lonmin + 360.] endif return ;The whole globe [might be] visible endif ; Edges are in order: bottom, right, top, left. sidex = [0, 2, 2, 0 ,0] ;X index of UV corner sidey = [1, 1, 3, 3, 1] ;Y index ; convert to float to truncate extraneous digits ; (otherwise Mercator will have incorrect LL_BOX for non-zero p0lon) box = float(!map.uv_box) uvcen = (box[[0,1]] + box[[2,3]])/2. del = 8 ;# of intervals again: found = bytarr(4) if north then begin lonmax = 180. & lonmin = -180. ;All longitudes are visible latmax = 90. & latmin = 90. endif else if south then begin lonmax = 180. & lonmin = -180. latmax = -90. & latmin = -90. endif else begin if dateline then begin lonmax = 0. & lonmin = 360. endif else begin lonmax = -180. & lonmin = 180. endelse latmin = 90. & latmax = -90. endelse ; March along the 4 uv edges of the plot, keeping the range in ; geo-coordinates. for iside = 0, 3 do begin s0 = box[[sidex[iside], sidey[iside]]] ;Start of edge s1 = box[[sidex[iside+1], sidey[iside+1]]] ;End of edge LastPointGood = 0 for i = 0, del-1 do begin s2 = s0 + (s1-s0) * (float(i)/del) ;UV coordinate n2 = s2 * [!x.s[1], !y.s[1]] + [!x.s[0], !y.s[0]] ;Normalized coord p = convert_coord(n2[0], n2[1], /NORMAL, /TO_DATA) if finite(p[0]) then begin if dateline and (p[0] lt 0) then p[0] = p[0] + 360 lonmax = lonmax > p[0] lonmin = lonmin < p[0] latmax = latmax > p[1] latmin = latmin < p[1] found[iside] = 1 ;Found a legit point ; This is a special KLUDGE to detect where the dateline crosses the left ; or right edge. if (dateline + north + south eq 0) and LastPointGood and $ (iside eq 1 or iside eq 3) then begin if abs(lonlast - p[0]) gt 180 then begin ;Dateline cross? dateline = 1 ;Yes... goto, again ;And rescale... endif endif LastPointGood = 1 lonlast = p[0] endif else LastPointGood = 0 endfor ;for i endfor ;for iside if total(found) eq 4 then begin ;If we didn't values along 4 sides then punt !map.ll_box = [latmin, lonmin, latmax, lonmax] ; print, !map.ll_box, format='(4f10.2)' endif end PRO map_set_limits, limit, uvrange, CYLINDRICAL=is_cyl ; Limit = 4 or 8 point vector containing lat/lon limits. (input). ; WARNING: limit may be changed. ; uvrange = uv extent of map in form of [umin, vmin, umax, vmax] ; (output). Also, on output, limit is converted to the 8 point form. ; COMPILE_OPT hidden ; Set unit scaling for convert_coord: !x.s = [0,1] !y.s = [0,1] if n_elements(limit) eq 4 then begin ;cvt 4 pt to 8 pt limit. ; *********** limit = [latmin, lonmin, latmax, lonmax] ; Finding the UV extent of a map, given a 4 point limit is a process ; frought with perils. Previous methods didn't always work, making a ; map that didn't show the entire range of limits. This uses the ; brute force method..... ; *********** limit = [latmin, lonmin, latmax, lonmax] londel = float(limit[3] - limit[1]) ;Check grid within the lat/lon limits if londel gt 360. then message, $ 'Longitude limits exceed 360 degrees', /CONTINUE if londel lt 0 then londel = londel + 360 latdel = float(limit[2] - limit[0]) nlon = (1+ceil(londel / 10.)) > 13 ;At least 13 pnts or every 10 degs nlat = (1+ceil(latdel / 10.)) > 9 ;At least 9 pnts or every 10 degs lats = findgen(nlat) * (latdel / (nlat-1)) + limit[0] lons = findgen(nlon) * (londel / (nlon-1)) + limit[1] eps = 1.0e-4 ;Fudge factor for cylindricals ;If cylindrical and the map covers 360 degrees of longitude, don't use ;the whole extent because of the ambiguity of +/- 180. if keyword_set(is_cyl) and londel eq 360 then $ lons = findgen(nlon) * ((londel - eps) / (nlon-1)) + (limit[1]+eps/2.) uv = convert_coord(lons # replicate(1., nlat), $ ;m x n array of lons replicate(1.,nlon) # lats, $ ;m x n array of lats /TO_NORM, /DATA) u = uv[0,*]; if total(finite(u)) lt 2 then $ ;Any points good? message, 'At least two limit points must be mappable', /INFO ; Extract the range uvrange = [min(u,/NAN, MAX=umax), min(uv[1,*],/NAN, MAX=vmax), $ umax, vmax] if uvrange[0] eq uvrange[2] or uvrange[1] eq uvrange[3] then $ message, 'Invalid map limits', /INFO endif else if n_elements(limit) eq 8 then begin ;8 point case is much simpler ; 8 point limit in form of: ; [latLeft,lonLeft, latTop, lonTop, LatRt, lonRt, LatBot, LonBot] ; for i=0,3 do begin uv = convert_coord(limit[i*2+1], limit[i*2], /TO_NORM, /DATA) if finite(uv[0]) eq 0 then $ message, 'Unmappable limit point: ' + $ string(limit[i*2+1], limit[i*2]), /INFO ; Left = uvrange(0)=uv(0); Top = uvrange(3) = uv(1); Right = ; uvrange(2) = uv(0); Bot = uvrange(1) = uv(1) uvrange[([0,3,2,1])[i]] = uv[([0,1,0,1])[i]] endfor endif else message, 'Map limit must have 4 or 8 points' if uvrange[0] ge uvrange[2] or uvrange[1] ge uvrange[3] then $ message, 'Warning, MAP limits are invalid.', /CONTINUE end FUNCTION map_rotxyz, p, rx, ry, rz ;Rotate the vector p(3) counterclockwise ; about the X, Y, and Z, by the amounts rx, ry, rz, degrees, in order. ; COMPILE_OPT hidden p1 = p dtor = !dpi/ 180. if rx ne 0.0 then begin sx = sin(rx * dtor) & cx = cos(rx * dtor) t = [[1,0,0],[0,cx, sx], [0, -sx, cx]] p1 = t # p1 endif if ry ne 0.0 then begin sy = sin(ry * dtor) & cy = cos(ry * dtor) t = [[cy, 0, -sy], [0, 1, 0], [sy, 0, cy]] p1 = t # p1 endif if rz ne 0.0 then begin sz = sin(rz * dtor) & cz = cos(rz * dtor) t = [[cz, sz, 0], [ -sz, cz, 0], [0,0,1]] p1 = t # p1 endif return, p1 end FUNCTION GREAT_CIRCLE, lon0, lat0, lon1, lat1, RADIANS=radians, PRINT=print ; Return: [sin(c), cos(c), sin(az), cos(az)] from the great circle ; distance between two points. c is the great circle angular distance, ; and az is the azimuth between two points. COMPILE_OPT hidden if KEYWORD_SET(radians) eq 0 then k = !dtor else k = 1.0 coslt1 = cos(k*lat1) sinlt1 = sin(k*lat1) coslt0 = cos(k*lat0) sinlt0 = sin(k*lat0) cosl0l1 = cos(k*(lon1-lon0)) sinl0l1 = sin(k*(lon1-lon0)) cosc = sinlt0 * sinlt1 + coslt0 * coslt1 * cosl0l1 ;Cos of angle between pnts sinc = sqrt(1.0 - cosc^2) if abs(sinc) gt 1.0e-6 then begin cosaz = (coslt0 * sinlt1 - sinlt0*coslt1*cosl0l1) / sinc ;Azmuith sinaz = sinl0l1*coslt1/sinc endif else begin ;Its antipodal cosaz = 1.0 sinaz = 0.0 endelse if keyword_set(print) then begin print, 'Great circle distance: ', acos(cosc) / k print, 'Azimuth: ', atan(sinaz, cosaz)/k ; print,'sinaz, cosaz = ', sinaz, cosaz endif return, [sinc, cosc, sinaz, cosaz] end PRO map_satellite_limit, n, xr, yr ; Return n or less points, defining the limit of the limb of a satellite ; projection COMPILE_OPT hidden a = findgen(n) * (2 * !pi / (n-1)) st = !map.p[0] > 1.0 st = sqrt((st-1.)/(st+1.)) ;Map limit radius for satellite projection xr = cos(a) * st ;Taken from Snyder Pg. 175 yr = sin(a) * st a = yr * !map.cosr + xr * !map.sinr xr = xr * !map.cosr - yr * !map.sinr yr = a a = yr * !map.p[2]/(!map.p[0]-1.0) + !map.p[3] good = where(a gt 0.0, count) if count eq 0 then begin ;If nothing, return nothing trash = temporary(xr) & trash = temporary(yr) return endif if count lt n then begin ;Extract & shift so pnts on map are contiguous i0 = 0 ;Get 1st good pnt following a bad pnt for i=1,count-1 do if good[i-1]+1 ne good[i] then i0 = i good = shift(good,-i0) endif a = a[good] xr = xr[good] * !map.p[3] / a yr = yr[good]/a end PRO map_proj_info, iproj, CURRENT=current, $ UV_RANGE=r, UVRANGE=rdummy, SCALE=meters, NAME=name, CYLINDRICAL=cyl, $ PROJ_NAMES = proj_names_p, CIRCLE=is_circular, CONIC=is_conic, $ AZIMUTHAL=is_azimuthal, UV_LIMITS = uv_limits, $ LL_LIMITS = ll_limits ; Return information about maps in general, or projection iproj. COMPILE_OPT hidden ; Iproj = projection number, Input. If CURRENT is set, then the ; current map projection is returned in Iproj. ; CURRENT = if set, use current proj number, and return that index in iproj. ; UV_RANGE = UV limits of the given projection, [umin, vmin, umax, ; vmax]. Also keyword UVRANGE is accepted here for bkwds compat. ; SCALE = (output) the number of meters to one unit of UV space. The ; radius of the earth, from Clarke's 1866 ellipsoid is used (6378206). ; NAME = the name of the projection. ; CYLINDRICAL = 1 if the projection is cylindrical or pseudo-cylindrical. ; CONIC = 1 if the projection is conic. ; PROJ_NAMES = string array containing the names of all the ; available projections. ; UV_LIMITS = returned rectangle of the current map in UV coordinates, ; [Umin, Vmin, Umax, Vmax]. ; LL_LIMITS = returned rectangle of the current map in latitude/longitude ; coordinates, [Lonmin, Latmin, Lonmax, Latmax]. Note, the ; range may not always be available, and the min and max may ; both be set to 0. ; Issues for a new projection: ; UVrange, proj_scale, is_cyl, is_azimuthal, circle, name, horizon, clipping ; common map_set_common, proj_scale, proj_is_cyl, proj_is_conic, $ proj_names, circle, proj_is_azim, user_defined on_error, 2 if n_elements(proj_scale) eq 0 then begin ;Init constant arrays? p2 = 2 * !pi s8 = 2 * sqrt(8) ro = 2 * 2.628 ; St Or LC La Gn Az Sa Cy Me Mo Si Ai Ha Al Tm Mi Ro Le Go ; 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 proj_scale = $ [0, 4, 2, 4, 4, 4, p2, 1, p2, p2, s8, p2, p2, s8, 4, 0, p2,ro, 0, p2] proj_is_cyl = $ [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1] proj_is_azim = $ [0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] proj_is_conic = $ [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0] circle = $ [0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0] proj_names = $ ["", "Stereographic", "Orthographic", "LambertConic", "LambertAzimuthal", $ "Gnomic", "AzimuthalEquidistant", "Satellite", "Cylindrical", $ "Mercator", "Mollweide", "Sinusoidal", "Aitoff", "HammerAitoff", $ "AlbersEqualAreaConic", "TransverseMercator", $ "MillerCylindrical", "Robinson", "LambertConicEllipsoid", $ "GoodesHomolosine"] user_defined = bytarr(n_elements(proj_names)) endif max_proj = n_elements(proj_names)-1 if arg_present(proj_names_p) then proj_names_p = proj_names if KEYWORD_SET(current) then iproj = !map.projection ;Get current projection if ARG_PRESENT(uv_limits) then uv_limits = !map.uv_box if ARG_PRESENT(ll_limits) then ll_limits = !map.ll_box if n_elements(iproj) eq 0 then return ;Rest of params are proj specific if iproj lt 1 or iproj gt max_proj then $ message, 'Projection number must be within range of 1 to'+ strtrim(max_proj,2) ones = [-1., -1., 1., 1.] r_earth = 6378206.4d0 ;Earth equatorial radius, meters, ; Clarke 1866 ellipsoid cyl = proj_is_cyl[iproj] name = proj_names[iproj] meters = proj_scale[iproj] * r_earth ;Scale [= 0 if 1:1] is_circular = circle[iproj] is_conic = proj_is_conic[iproj] is_azimuthal = proj_is_azim[iproj] if arg_present(r) or arg_present(rdummy) or arg_present(meters) then begin case iproj of 1: r = ones ;Stereo 2: r = ones ;Ortho 3: r = 2 * ones ;Lambert Conic for sphere. 4: r = 2 * ones ;lambert 5: r = 2 * ones ;gnomic 6: r = !DPI * ones ;azimuth 7: if !map.p[1] ne 0 then begin ;Complicated satellite? map_satellite_limit, 60, xr, yr if n_elements(xr) eq 0 then r = [0,0,0,0] $ else r = [min(xr, max=maxxr), min(yr, max=maxyr), maxxr, maxyr] meters = (r[2]-r[0]) * r_earth endif else begin st = !map.p[0] > 1.0 ;For satellite r = sqrt((st-1)/(st+1)) * ones ;Simple case satellite meters = (r[2]-r[0]) * r_earth endelse 8: r = [-1., -.5, 1., .5] * !DPI ;Cylind equidistant 9: r = [-!DPI, -2.43, !DPI, 2.43] ;Mercator (+/- 80 degs of lat) ~ 98% 10: r = sqrt(2.0) * [-2.,-1.,2.,1.] ;Mollweide 11: r = [-1., -.5, 1., .5] * !DPI ;sinusoidal 12: r = [-1., -.5, 1., .5] * !DPI ;Aitoff 13: r = sqrt(2.0) * [-2., -1., 2., 1.] ; Hammer-aitoff 14: r = 2 * ones ;Albers Conic 15: r = 2*!pi * r_earth * ones ;Transverse Mercator, Units = meters. 16: r = [-!DPI, -1.832, !DPI, 1.832] ;Miller Cylind (+/- 80 degs of lat) 17: r = [-2.628, -1.349, 2.628, 1.349] ;Robinson 18: r = 2 * !pi * r_earth * ones ;Lambert Conic ellipsoid, Units = meters 19: r = [-!DPI, -1.37, !DPI, 1.37] ;Goodes Homolosine endcase rdummy = r ;For backwards compatibility if meters ne 0 then meters = meters / (r[2]-r[0]) $ ;Get scale in UV units else meters = 1.0 ;UV Units are meters. endif end PRO map_GoodeSeg, x0, y0, x1, y1, FILL=fill, _EXTRA=e ; Create the outline of a Goodes Homolosine map. If FILL is a ; non-zero scalar, store the points of the outline in FILL, otherwise ; just draw the outline. If FILL is set, at the end, the array FILL ; will contain all the points for the outline. ; COMPILE_OPT hidden nsegs = fix(18 > (abs(x1-x0) > abs(y1-y0))/5) ;Interval: 18 segs or every 5 degs x = x0 + findgen(nsegs+1) * (float(x1-x0) / nsegs) + !map.p0lon y = y0 + findgen(nsegs+1) * (float(y1-y0) / nsegs) if keyword_set(fill) eq 0 then plots, x, y, NOCLIP=0, _EXTRA=e $ else begin if n_elements(fill) eq 1 then fill = transpose([[x],[y]]) $ else fill = [[fill],[transpose([[x],[y]])]] endelse end PRO map_horizon, FILL=fill, NVERTS=n, ZVALUE=zvalue, _EXTRA=Extra ; Draw/fill the horizon (limb) of a map COMPILE_OPT hidden map_proj_info, /CURRENT, NAME=name, UV_RANGE=p, CIRCLE=is_circular, $ CONIC=is_conic r = !map.rotation ;Our rotation if is_conic then return ;Can't do conics if n_elements(n) le 0 then n = 60 ;# of vertices a = findgen(n+1) * (2 * !pi / n) ;N angles from 0 to 2 pi bFill = keyword_set(fill) if name eq "Satellite" and !map.p[1] ne 0 then begin ;Satellite proj, tilted map_satellite_limit, 60, xr, yr r = 0.0 ;No more rotation endif else if is_circular then begin ;Std circular projection xr = (p[2]-p[0])/2. * cos(a) yr = (p[3]-p[1])/2. * sin(a) endif else if name eq "Sinusoidal" then begin ;Sinusoidal ; flon = (!map.out(3)-!map.out(2))/360. flon = 1.0 xr = flon * !pi * cos(a) yr = fltarr(n+1) for i=0,n do case fix(a[i]/(!pi/2)) of 0: yr[i] = a[i] 1: yr[i] = !pi-a[i] 2: yr[i] = !pi-a[i] 3: yr[i] = a[i]-2*!pi 4: yr[i] = 0.0 ;Last pnt endcase endif else if name eq "Robinson" then begin ;Robinson lp = [2.628, 2.625, 2.616, 2.602, 2.582, 2.557, 2.532, 2.478, 2.422, $ 2.356, 2.281, 2.195, 2.099, 1.997, 1.889, 1.769, 1.633, 1.504, 1.399] lm = [0., 0.084, 0.167, 0.251, 0.334, 0.418, 0.502, 0.586, 0.669, 0.752, $ 0.833, 0.913, 0.991, 1.066, 1.138, 1.206, 1.267, 1.317, 1.349] xr = [lp, -reverse(lp), -lp, reverse(lp)] yr = [lm, reverse(lm), -lm, -reverse(lm)] endif else if name eq "GoodesHomolosine" then begin delta = 0.001 ;Fudge factor in degs for splitting N90 = 90-delta S90 = delta-90 map_GoodeSeg, -180+delta, S90, -180+delta, N90, FILL=bFill, _EXTRA=Extra map_GoodeSeg, -40-delta, N90, -40-delta, 0, FILL=bFill, _EXTRA=Extra map_GoodeSeg, -40+delta, 0, -40+delta, N90, FILL=bFill, _EXTRA=Extra map_GoodeSeg, 180-delta, N90, 180-delta, S90, FILL=bFill, _EXTRA=Extra map_GoodeSeg, 80+delta, S90, 80+delta, 0, FILL=bFill, _EXTRA=Extra map_GoodeSeg, 80-delta, 0, 80-delta, S90, FILL=bFill, _EXTRA=Extra map_GoodeSeg, -20+delta, S90, -20+delta, 0, FILL=bFill, _EXTRA=Extra map_GoodeSeg, -20-delta, 0, -20-delta, S90, FILL=bFill, _EXTRA=Extra map_GoodeSeg, -100+delta, S90, -100+delta, 0, FILL=bFill, _EXTRA=Extra map_GoodeSeg, -100-delta, 0, -100-delta, S90, FILL=bFill, _EXTRA=Extra if n_elements(bFill) gt 1 then begin ; This is a kludge because filling in map coordinates makes wrap around ; errors around the north pole. bFill = CONVERT_COORD(REVERSE(bFill,2), /DATA, /TO_NORM) polyfill, bFill, NOCLIP=0, _Extra=Extra, /NORM endif return endif else if name eq "UserDefined" then begin ;User projection endif else begin ;Rectangular xr = [p[0], p[2], p[2], p[0], p[0]] yr = [p[1], p[1], p[3], p[3], p[1]] endelse if r ne 0.0 then begin t1 = xr * !map.cosr + yr * !map.sinr ;Now rotate by - angle... yr = yr * !map.cosr - xr * !map.sinr xr = t1 endif xtsave = !x.type if n_elements(zvalue) eq 0 THEN zvalue = 0 !x.type=0 ;Plot in UV space if keyword_set(fill) then polyfill, xr, yr, _EXTRA=Extra, NOCLIP=0 $ else plots, xr, yr, zvalue, _EXTRA=Extra, NOCLIP=0 !x.type=xtsave end PRO map_set_split, ADD=add, EQUATOR = equator ; Set the splitting hemisphere. If EQUATOR is set, set the splitting ; point as if the center of projection was on the equator, i.e. use ; a latitude of 0. COMPILE_OPT hidden pole = !map.pole[4:6] ;Locn of pole if keyword_set(equator) then begin ;use lat = 0 split = [!map.p0lon, 0, -sin(!map.u0), cos(!map.u0), 0.0, 0.0] endif else begin coslat = cos(!map.v0) p0 = [ cos(!map.u0) * coslat, sin(!map.u0) * coslat, sin(!map.v0)] pln = crossp(p0, pole) split=[!map.p0lon, !map.p0lat, pln, 0.0] endelse noadd = keyword_set(add) eq 0 ;= 0 if we're adding MAP_CLIP_SET, RESET=noadd, SPLIT=split, TRANSFORM=noadd end PRO map_set_clip, PNAME=pname, IPROJ=iproj, CLIP_RADIUS = r0 COMPILE_OPT hidden ; Set up the default clipping/splitting for the given projection. ; if n_elements(pname) eq 0 then map_proj_info, iproj, NAME=pname if n_elements(iproj) eq 0 then begin map_proj_info, PROJ_NAMES=pnames pname = pnames[(where(pname eq pnames, count))[0]] if count eq 0 then message, "No such projection available: " + pname endif r = -1.0 ;Great circle clipping radius... map_proj_info, iproj, CONIC=is_conic if is_conic then begin ;Conical projection? isign = 2*(!map.p[3] ge 0.0)-1 ;Assume its in one hemisphere map_set_split ;Clip both poles ;Opposite hemisphere at 10 degs map_clip_set, CLIP_PLANE=[0,0,isign, sin(!dtor * 10)] ;pole in same hemisphere at 75 degs map_clip_set, CLIP_PLANE=[0,0,- isign, sin(!dtor * 75)] endif else if pname eq 'Mercator' or pname eq 'MillerCylindrical' then begin rm = 0.975 ;Go out 97.5% of the way to the poles map_set_split map_clip_set, CLIP_PLANE=[-!map.pole[4:6], rm] map_clip_set, CLIP_PLANE=[!map.pole[4:6], rm] endif $ ;Other cylindricals else if pname eq "Stereographic" or pname eq "Orthographic" then r = -1.0e-5 $ else if pname eq "Satellite" then r = - 1.01/!map.p[0] $ else if pname eq "Gnomic" then r = -0.5 $ else if pname eq "LambertAzimuthal" or pname eq "AzimuthalEquidistant" then $ r=0.8 $ else if pname eq "TransverseMercator" then begin ; Clip plane is perpendicular to the XY projection of the vector to the ; center of projection, and to the XY Plane. map_clip_set, CLIP_PLANE=[cos(!map.u0), sin(!map.u0), 0, 0] endif else if pname eq "GoodesHomolosine" then begin r = -1.0 ; What a screwey interrupted projection. Splits are at 180 degrees ; from the seams. splits = [-180., -40, -100, -20, 80]+180. + !map.p0lon for i=0, n_elements(splits)-1 do begin ;Add each split theta = !dtor * splits[i] MAP_CLIP_SET, SPLIT=[splits[i], 0, -sin(theta), cos(theta), 0., 0.] endfor endif else begin ;If we haven't setup a clip, add splitting map_proj_info, /CURRENT, CYLINDRICAL=is_cyl if is_cyl then map_set_split endelse ; If r is set, set clipping to points on a plane whose normal ; passes thru the center (u0,v0), and at a distance of r from the ; origin. r < 0 is the side towards (u0,v0). if r ne -1.0 then begin if n_elements(r0) ne 0 then r = r0 ;User supplied value? map_clip_set, CLIP_PLANE=[cos(!map.u0) * !map.coso, $ sin(!map.u0) * !map.coso, !map.sino, r] endif map_clip_set, /TRANSFORM end ; ***************************************************************** PRO MAP_SET, p0lat, p0lon, rot, $ ; ********** Projection keywords: PROJECTION=proj, $ ;The projection index NAME=name, $ ;The projection name as a string STEREOGRAPHIC = stereographic, $ ORTHOGRAPHIC = orthographic, $ CONIC = conic, $ LAMBERT = lambert, $ GNOMIC = gnomic, $ AZIMUTHAL = azimuth, $ SATELLITE = satellite, $ ;Also called general perspective proj CYLINDRICAL = cylindrical, $ MERCATOR = mercator, $ MILLER_CYLINDRICAL=miller, $ MOLLWEIDE = mollweide, $ SINUSOIDAL = sinusoidal, $ AITOFF = aitoff, $ ;Original Aitoff projection HAMMER = hammer, $ ;This is really the Hammer-Aitoff projection ALBERS = albers, $ ;Alber's equal-area conic TRANSVERSE_MERCATOR = utm, $ ;Transverse mercator for an ellipsoid ROBINSON = robinson, $ ;Robinson projection GOODESHOMOLOSINE = goodes, $ ;Goode's Homolosine. Origin should be (0,0). ; Also called Gauss-Krueger in Europe. ; **** Projection specific keywords: ELLIPSOID = ellips, $ ;Defines the ellipsoid for projections ;that handle the ellipsoidal case, ;currently Transverse Mercator and ;Lambert's conical. ;3 elements: a = equatorial radius (meters), e^2 ;= eccentricity squared, k0 = scale on ;central meridian. e^2 = 2*f-f^2, ;where f = 1-b/a, b = polar radius ;Default is the Clarke 1866 ellipsoid, = ; [6378206.4d0, 0.00676866d0, 0.9996d0] CENTRAL_AZIMUTH=cent_azim, $ ;Angle of central azimuth (degrees) for: ; Cylindrical, Mercator, Miller, ; Mollw, Robinson, and Sinusoidal ; projections. Default = 0. ; The pole is placed at an azimuth of ; CENTRAL_AZIMUTH degrees CCW of north. STANDARD_PARALLELS = std_p, $ ;One or two standard parallels for conics, ; One or two element array. SAT_P = Sat_p, $ ;Satellite parameters: Altitude expressed in ; units of radii, Omega, and rotation. ; Rotation may also specified by the ; rot parameter. ; ********** MAP_SET specific keywords: CLIP=clip, $ ;Default = do map specific clipping, ; CLIP=0 to disable REVERSE=reverse, $ ;Reverse X and/or Y axes. 0=none, 1 = reverse X, ; 2 = reverse Y, 3 = reverse both. SCALE=scale, $ ;Construct isotropic map with a given scale. ; Map Scale is 1:scale, otherwise fit to window ISOTROPIC = iso, $, ;If set, make X and Y scales equal LIMIT = limit, $ ;4 or 8 point lat/lon limit: ; 4 point: [latmin, lonmin, latmax, lonmax] ; 8 point: [latLeft,lonLeft, latTop, ; lonTop, LatRt, lonRt, LatBot, LonBot] ; ********** 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, $ ; E_HORIZON = structure containing ; extra keywords passed to the ; map_horizon procedure, e.g. ; E_HORIZON={FILL:1} ; ********** MAP_CONTINENTS keywords: CONTINENTS = continents, E_CONTINENTS=econt, $ ;E_CONTINENTS = structure ; containing extra keywords to be ; passed to ; map_continents, e.g. E_CONTINENTS={FILL:1} USA=usa, HIRES = hires, $ MLINESTYLE=mlinestyle, MLINETHICK=mlinethick, CON_COLOR=con_color, $ ; ********** MAP_GRID keywords: GRID=grid, E_GRID=egrid, $ ;E_GRID = extra keywords structure 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 ; ***************************************************************** ;+ ; Limit = A four or eight element vector. ; If a four element vector, [Latmin, Lonmin, Latmax, ; Lonmax] specifying the boundaries of the region ; to be mapped. (Latmin, Lonmin) and (Latmax, ; Lonmax) are the latitudes and longitudes of ; two diagonal points on the boundary with ; Latmin < Latmax and Lonmin < Lonmax. ; For maps that cross the international dateline, ; specify west longitudes as angles between ; 180 and 360. ; If an eight element vector: [ lat0, lon0, ; lat1, lon1, lat2, lon2, lat3, lon3] specify ; four points on the map which give, respectively, ; the location of a point on the left edge, ; top edge, right edge, and bottom edge of ; the map extent. ;- ; If the map range can be simply expressed as [latmin, lonmin, latmax, ; lonmax], then they are contained here. With many maps, this is ; impossible. These values are used to speed MAP_GRID and ; MAP_CONTINENTS, by not drawing elements that are clearly off the ; map. common map_set_common on_error, 2 ;Return to caller if error. del = 1.0e-6 dtor = !dpi/ 180. dpi2 = !dpi/2. !map.ll_box = 0 ;Assume no limits if n_elements(proj) then iproj = proj $ else if n_elements(name) then begin ;Name supplied, match it map_proj_info, PROJ_NAMES = names ;Available names name1 = strcompress(strupcase(name), /REMOVE_ALL) ;Ignore case & blanks iproj = (where(name1 eq strupcase(names)))[0] ;Match it if iproj lt 1 then message, 'Projection "'+name+'" not available' endif else if keyword_set(stereographic) then iproj =1 $ else if keyword_set(orthographic) then iproj =2 $ else if keyword_set(conic ) then iproj =3 $ else if keyword_set(lambert) then iproj =4 $ else if keyword_set(gnomic ) then iproj =5 $ else if keyword_set(azimuth) then iproj =6 $ else if keyword_set(satellite) then iproj =7 $ else if keyword_set(mercator) then iproj =9 $ else if keyword_set(mollweide) then iproj = 10 $ else if keyword_set(sinusoidal) then iproj = 11 $ else if keyword_set(aitoff) then iproj = 12 $ else if keyword_set(hammer) then iproj = 13 $ else if keyword_set(albers) then iproj = 14 $ else if keyword_set(utm) then iproj = 15 $ else if keyword_set(miller) then iproj = 16 $ else if keyword_set(robinson) then iproj = 17 $ else if keyword_set(goodes) then iproj = 19 $ else iproj=8 ;Assume cylindrical ;Special case for lambert's conic on an ellipsoid. if keyword_set(ellips) and iproj eq 3 then iproj = 18 !map.projection = iproj map_proj_info, /CURRENT, NAME=pname, CONIC=is_conic ;Get symbolic name ; Initial erase? if !P.multi[0] eq 0 and keyword_set(advance) then erase if (not KEYWORD_SET(noerase)) and (not KEYWORD_SET(advance)) THEN erase ; ******** Supply defaults ********** if n_params() lt 3 then rot = 0.0d0 if n_params() lt 2 then p0lon = 0d0 if n_params() lt 1 then p0lat = 0d0 if n_elements(cent_azim) le 0 then cent_azim = 0.0 if n_elements(color) eq 0 then color = !p.color ;Default color if n_elements(title) eq 0 then title = " " if n_elements(t3d) le 0 then t3d = 0 if n_elements(zvalue) eq 0 then zvalue = 0. if n_elements(charsize) eq 0 then charsize = !p.charsize if charsize le 0.0 then charsize = 1.0 if n_elements(clip) eq 0 then clip = 1 ;default for CLIP is ON if pname eq "GoodesHomolosine" and p0lat ne 0 then begin message,"Goode's Homolosine: center is always on equator.", /CONTINUE p0lat = 0.0 endif ;Hack for CONIC backwards compatibility if pname eq "LambertConic" and (rot ne 0) and n_elements(std_p) eq 0 then begin std_p = [p0lat, p0lon] ;Two standard parallels p0lon1 = rot ;and rotation is the longitude endif else p0lon1 = p0lon if abs(p0lat) gt 90.0 then $ message,'Latitude must be in range of +/- 90 degrees' if abs(p0lon) gt 360.0 then $ message,'Longitude must be in range of +/- 360 degrees' u = double(p0lon1) ;Reduce lon to +/- 180 while u lt -180.0 do u = u + 360. while u gt 180.0 do u = u - 360. !map.p0lon = u !map.p0lat = p0lat !map.u0 = u * dtor !map.v0 = p0lat * dtor !x.type = 3 ;New map code !y.type = 0 ; Simple projection? (used for (psuedo) cylindricals) !map.simple = abs(p0lat) le del and abs(cent_azim) le del !map.sino = sin(!map.v0) !map.coso = cos(!map.v0) !map.rotation = rot !map.sinr = sin(rot * dtor) !map.cosr = cos(rot * dtor) ; Compute position of Pole which is a distance of !PI/2 from ; (p0lon, p0lat), at an azimuth of rot CCW of north. ; pole = [0.0, sin(cent_azim * dtor), cos(cent_azim * dtor)] p2 = map_rotxyz(pole, 0, -p0lat, u) ;Rotate to put origin at (0,1,0) plat = asin(p2[2]) cosla = sqrt(1.0-p2[2]^2) if cosla eq 0.0 then begin ;On pole? plon = 0.0 & sinln = 0.0 & cosln = 1.0 endif else begin plon = atan(p2[1], p2[0]) sinln = p2[1] / cosla & cosln = p2[0] / cosla endelse ; lon/lat, sin(lat), cos(lat), xyz, Location of pole !map.pole = [plon, plat, p2[2], cosla, p2] if n_elements(ellips) ne 3 then $ ;Default == Clarke 1866 ellipsoid ellips1 = [6378206.4d0, 0.00676866d0, 0.9996d0] $ else ellips1 = ellips if is_conic then begin if n_elements(std_p) eq 0 then begin if p0lat ne 0 then s = [p0lat, p0lat] else s = [20., 20.] endif else if n_elements(std_p) eq 1 then s = [std_p, std_p] $ else if n_elements(std_p) eq 2 then s = std_p $ else message, 'STANDARD_PARALLELS must have 1 or 2 elements' s = s * dtor if ABS(!map.p0lat) eq 90. then $ message, 'Center of projection may not be pole for conics' if pname eq "LambertConic" then begin ; ***** if s[0] eq s[1] then n = sin(s[0]) $ else n = alog(cos(s[0])/cos(s[1])) / $ alog(tan(!dpi/4+s[1]/2)/tan(!dpi/4+s[0]/2)) if n eq 0.0 then message, $ 'Equator is illegal standard parallel for conic' F = cos(s[0]) * tan(!dpi/4 + s[0]/2)^n/n r0 = F/tan(!dpi/4+!map.v0/2.)^n !map.p = [n, F, r0, s, 0.0] endif else if pname eq "LambertConicEllipsoid" then begin ;***** if s[0] * s[1] lt 0 then $ message, 'Standard parallels must be in same hemisphere' e = sqrt(ellips1[1]) ;Compute n, F, r0 using formulae from Snyder, Page 107-108. m = cos(s) / sqrt(1.0d0- ellips1[1] * sin(s)^2) t = [!map.v0, s[0], s[1]] ;Angles for t t = tan(!dpi/4 - t/2.0) / $ ((1.0 - e * sin(t)) / (1.0 + e*sin(t)))^(e/2.) if s[0] eq s[1] then n = sin(s[0]) $ else n = (alog(m[0]) - alog(m[1])) / (alog(t[1]) - alog(t[2])) F = m[0] / (n * t[1] ^ n) r0 = ellips1[0] * F * t[0] ^ n !map.p = [n, F, r0, s, ellips1] endif else begin ;Albers conic n = (sin(s[0]) + sin(s[1]))/2. if n eq 0 then begin message,'Standard parallels should not be 0', /INFO n = 0.5 endif c = cos(s[0])^2 + 2 * n * sin(s[0]) r0 = sqrt(c - 2*n*sin(!map.v0))/n !map.p = [n, c, r0, s] endelse endif if pname eq "Satellite" then begin ;Special params for satellite projection. if n_elements(sat_p) eq 0 then !map.p[0] = 2.0 $ ;Altitude in radii else !map.p[0] = sat_p[0] if n_elements(sat_p) le 1 then omega = 0d0 else omega = dtor * sat_p[1] ;Save em. sat(1) = TRUE for simple case (Vertical perspective) !map.p[1] = omega !map.p[2] = sin(omega) ;Somega = p[1] !map.p[3] = cos(omega) ;comega if n_elements(sat_p) ge 3 then begin ; For backwards compatibility, if ; Sat_p contains three elements, interpret the 3rd element as the rotation: !map.rotation = sat_p[2] + rot !map.sinr = sin(!map.rotation * dtor) !map.cosr = cos(!map.rotation * dtor) endif ;3 elements endif map_proj_info, /CURRENT, UV_RANGE=uvrange, CYLINDRICAL = is_cyl ; Set up cylindrical projections if pname eq "TransverseMercator" then begin ;Special params for UTM ; ellips is a 3 element array containing the ellipsoid parameters: ; [a, e^2, k0] ; a = Equatorial radius, in meters. ; b = polar radius, b = a * (1-f), f = 1-b/a ; e^2 = eccentricity^2 = 2*f-f^2, where f = flattening. ; k0 = scale on central meridian, = 0.9996 for UTM. e_2 = ellips1[1]/(1.0-ellips1[1]) ; k0 e2 e_2 a b m0 m !map.p = [ellips1[2], ellips1[1], e_2, ellips1[0], ellips1[1], 0., 0.] !x.s = [0,1] & !y.s = [0,1] & !x.type = 3 ;Fake coordinate system. q = convert_coord(u, p0lat, /data, /to_norm) ;Get m0 !map.p[5] = !map.p[6] ;Set m0 from Center of projection endif else if is_cyl then begin ;other Cylindrical or pseudo-cylind proj? ; I don't know why we can't solve for this angle (the azimuth of (u0,v0) ; from (xp, yp)) using the law of sines. This is solving it the hard ; and long way..... But it works.... az = atan(!map.coso * sin(!map.u0-plon), $ cosla * !map.sino - p2[2] * !map.coso *cos(!map.u0-plon)) !map.p[0] = dpi2-az ; print, 'Pole: ', plon/dtor, plat /dtor, ', Az=',az / dtor endif uvrange_orig = uvrange ;Save original UV if is_cyl then begin ;Adjust for rotation & wrap in cylind projs ; Set the x coordinate of the edge in UV space for wrap-around detection. ; This is only revelant in (pseudo) cylindrical projections. u = uvrange[[0,0,2,2]] ;The 4 corners. v = uvrange[[1,3,1,3]] t1 = u * !map.cosr + v * !map.sinr ;Now rotate by - angle... v = v * !map.cosr - u * !map.sinr uvrange[0] = min(t1, max=t2) uvrange[2] = t2 uvrange[1] = min(v, max=t2) uvrange[3] = t2 endif if n_elements(limit) gt 0 then begin ;Get UV range. if keyword_set(scale) then $ message, 'Conflicting keywords specified: LIMIT and SCALE', /INFO lim = limit ;Save old if n_elements(lim) eq 4 then begin ; Adjust longitude range: lonmax > lonmin, and lonmin in the range of ; -180 to 180. lonmin = 1, lonmax = 3 while lim[3] le lim[1] do lim[3] = lim[3] + 360. while lim[3] gt (lim[1]+360.) do lim[3] = lim[3] - 360. t1 = [0, 360., 0, 360.] while lim[1] ge 180. do lim = lim - t1 while lim[1] lt -180 do lim = lim + t1 !map.ll_box = lim ;Save if 4 point limit endif map_set_limits, lim, uvrange, CYLINDRICAL=is_cyl endif ;***************************************************************** if keyword_set(position) then $ plot, [0,1], xsty=5, ysty=5, /NODATA, $ TITLE=title, /NOERASE, POSITION=position, CHARSIZE=charsize, COLOR=color $ else begin ;Position not set ; Use margins if n_elements(xmargin) eq 1 then xmar = [xmargin, xmargin] $ else if n_elements(xmargin) eq 2 then xmar = xmargin $ else xmar = [1,1] if n_elements(ymargin) eq 1 then ymar = [ymargin, ymargin] $ else if n_elements(ymargin) eq 2 then ymar = ymargin $ else ymar = [1,2] plot, [0,1], XSTYLE=5, YSTYLE=5, /NODATA, COLOR=color, $ XMARGIN= xmar, YMARGIN= ymar, TITLE=title, /NOERASE, $ CHARSIZE = charsize endelse !x.type = 3 ;For new maps; !map.uv_box = uvrange if keyword_set(noborder) eq 0 then begin ; fudge is used to add a bit of spacing between the border and the ; extent of the map region. fudge of 0.01 indicates that the extra ; spacing (internal map margin) should be 1% of the original map extent. ; ; this extra space is now turned off when the noborder keyword is set ; fudge = 0.01 s = (uvrange[2]-uvrange[0]) * fudge ss = (uvrange[3] - uvrange[1]) * fudge uvrange = uvrange + [-s, -ss, s, ss] endif ; Figure the size of the drawing area on the screen in normalized units x_size = !x.window[1]-!x.window[0] y_size = !y.window[1]-!y.window[0] x_cen = total(!x.window)/2. ;Midpoints in X & Y y_cen = total(!y.window)/2. ; Compute the X and Y scale factors, !x.s[1], !y.s[1]: if keyword_set(scale) then begin ;Absolute scale provided? ;Absolute size of entire map, meters / uvrange map_proj_info, /CURRENT, SCALE=meters ;Meters per UV unit s = meters / scale ;scaled ss = !d.x_size / !d.x_px_cm / 100. ;width of window in meters !x.s[1] = s / ss !y.s[1] = !x.s[1] * !d.x_size / !d.y_size ;Adjust for aspect ratio ;Adjust the UV range to include the plotting window **** t = [x_size / !x.s[1], y_size / !y.s[1]] / 2 ;half the UV range c = [uvrange[0]+uvrange[2], uvrange[1]+uvrange[3]]/2 ;UV center ; Take smallest UV box: uvrange = [c[0]-t[0] > uvrange[0], c[1]-t[1] > uvrange[1], $ c[0]+t[0] < uvrange[2], c[1]+t[1] < uvrange[3]] ;New uv range !map.uv_box = uvrange endif else if keyword_set(iso) then begin ; Scale in pixels/uvunit, correct for ISOMETRIC aspect ratio. Use smaller ; of X and Y scale factors sx = x_size * !d.x_size /(uvrange[2]-uvrange[0]) ;X scale sy = y_size * !d.y_size/(uvrange[3]-uvrange[1]) ;Y scale !x.s[1] = (sx < sy) / !d.x_size ;Set smaller !y.s[1] = !x.s[1] * !d.x_size / !d.y_size if sx gt sy then $ ;Resize window to fit map area !x.window = x_cen + sy/sx/2 * [-x_size, x_size] $ else !y.window = y_cen + sx/sy/2 * [-y_size, y_size] endif else begin ;Scale to fit WINDOW !x.s[1] = x_size/(uvrange[2]-uvrange[0]) !y.s[1] = y_size/(uvrange[3]-uvrange[1]) endelse i = n_elements(reverse) eq 0 ? 0 : reverse ;Reverse axes? if i and 1 then !x.s[1] = -!x.s[1] ;Flip left to right? if (i and 2) ne 0 then !y.s[1] = -!y.s[1] ;Top to bottom? ; Compute offsets to center UV rectangle in the window !x.s[0] = x_cen - (uvrange[0]+uvrange[2])/2. * !x.s[1] !y.s[0] = y_cen - (uvrange[1]+uvrange[3])/2. * !y.s[1] map_clip_set, /RESET ;Clear clipping pipeline. if keyword_set(clip) then begin ;Setup clipping r0 = 2 ; Special case check for polar projection where the latitude range is ; limited. if (n_elements(lim) eq 4) then begin if (abs(p0lat) eq 90) and $ ;Center at pole? (p0lat eq lim[0] or p0lat eq lim[2]) then $ ;a lat limit is at pole ;Fudge a little to so the last ;parallel isn't clipped. (the 1.0e-5) r0 = -cos((abs(lim[2]-lim[0])+ 1.0e-5) * !dtor) endif if r0 eq 2 then map_set_clip, IPROJ=iproj $ ;Set default clipping else map_set_clip, IPROJ=iproj, CLIP_RADIUS=r0 ;Set our radius if (n_elements(lim) gt 0) or keyword_set(scale) then $ ;Clip UV? MAP_CLIP_SET, CLIP_UV = uvrange ; Set up plotting clip window t1 = (uvrange[[0,2]] * !x.s[1] + !x.s[0]) * !d.x_size !p.clip[[0,2]] = [t1[0] < t1[1], t1[0] > t1[1]] t1 = (uvrange[[1,3]] * !y.s[1] + !y.s[0]) * !d.y_size !p.clip[[1,3]] = [t1[0] < t1[1], t1[0] > t1[1]] endif if clip and total(abs(!map.ll_box)) eq 0 then $ ; Try to make a lon/lat range map_set_ll_box if keyword_set(noborder) eq 0 then $ ;Draw the border. plots, !x.window[[0,1,1,0,0]], !y.window[[0,0,1,1,0]], $ COLOR=color, zvalue, /NORM, /NOCLIP, T3D=t3d ; Collect the common graphics keywords map_struct_append, egraphics, "COLOR", color map_struct_append, egraphics, "T3D", t3d map_struct_append, egraphics, "ZVALUE", zvalue ; Process MAP_HORIZON keywords: ************** if keyword_set(horizon) or keyword_set(ehorizon) then begin map_struct_merge, ehorizon, egraphics ;Add common graphics keywords MAP_HORIZON, _EXTRA=ehorizon endif ; Process MAP_CONTINENT keywords: ************** if n_elements(mlinestyle) then map_struct_append, econt, "LINESTYLE", mlinestyle if n_elements(mlinethick) then map_struct_append, econt, "THICK", mlinethick if n_elements(con_color) then map_struct_append, econt, "COLOR", con_color if n_elements(hires) then map_struct_append, econt, "HIRES", hires if n_elements(continents) then map_struct_append, econt, "CONTINENTS", continents if n_elements(usa) then map_struct_append, econt, "USA", usa if n_elements(econt) gt 0 or keyword_set(continents) then begin map_struct_merge, econt, egraphics ;Add common graphics kwrds MAP_CONTINENTS, _EXTRA=econt endif ; Process MAP_GRID keywords: ************** if n_elements(label) then map_struct_append, egrid, "LABEL", label if n_elements(latlab) then map_struct_append, egrid, "LATLAB", latlab if n_elements(lonlab) then map_struct_append, egrid, "LONLAB", lonlab if n_elements(latdel) then map_struct_append, egrid, "LATDEL", latdel if n_elements(londel) then map_struct_append, egrid, "LONDEL", londel if n_elements(latalign) then map_struct_append, egrid, "LATALIGN", latalign if n_elements(lonalign) then map_struct_append, egrid, "LONALIGN", lonalign do_grid = (keyword_set(grid) + n_elements(egrid) + n_elements(glinestyle) + $ n_elements(glinethick)) ne 0 map_struct_merge, egrid, egraphics if n_elements(glinestyle) then $ ;These keywords supersede those in egraphics map_struct_append, egrid, "LINESTYLE", glinestyle if n_elements(glinethick) then $ map_struct_append, egrid, "THICK", glinethick if do_grid then MAP_GRID, CHARSIZE=charsize, _EXTRA=egrid if KEYWORD_SET(ADVANCE) and !P.Multi[0] gt 0 THEN $ !P.Multi[0] = !P.Multi[0] - 1 $ else !p.multi[0] = !p.multi[1] * !p.multi[2] - 1 end