func cgd /* DOCUMENT cgd these routines are some that I think are of general use to the climate and global dynamics division. Current contents: find_nice_contours: returns a nice set of contour levels, given an array and optionally a desired number of contour levels conint: returns a set of contour levels spaced a prescribed interval apart; find_ctlabs: return a nice set of labels for color tables, given an array and number of labels. Differs from find_nice_contours because the number of labels is not adjusted. This is necessary because the number of labels should have a specific and precise relationship to the number of colors in the colortable. fillc: a routine to determine filled areas and the corresponding contours plfc: uses results from fillc to plot the filled contours (including above and below regions) plfr: uses results from fillc to plot the filled contours (excluding above and below regions) plc_labels: uses results from fillc to plot the labels for contours and optionally the contours themselves plsys_set, plsys_get set and retrieve the plot system coordinate. get_hfile: open a netcdf version of a ccm history file */ { } extern plsys_coord; plsys_coord = 1; func span_round (start, stop, n) /* DOCUMENT span_round same as span, but rounds anything very small compared to the endpoints back to zero */ { x = span(start, stop, n) list = where(abs(x/max(abs(x))) < 1.e-8); if (numberof(list)> 0) x(list) = 0; return x; } func conset(inar, conint) { /*DOCUMENT conset(inar, conint); returns contour levels that are even multiples of conint that span the bounds of array inar if conint is nil or zero, calls find_nice_contours */ if (is_void(conint) || conint == 0. ) { levs = find_nice_contours(inar); } else { mini = min(inar); maxi = max(inar); mini2 = floor(mini/conint)*conint; maxi2 = ceil(maxi/conint)*conint; print, mini, maxi, mini2, maxi2; nint = int((maxi2-mini2)/conint)+1; levs = span(mini2,maxi2,nint); } return levs; } func plsys_set (n) /* DOCUMENT plfsys_set(n) store the coordinate system away for later access with plsys_get then call plsys. This allows one to retrieve the current coordinate system so you can reset it after working in a particular system. SEE ALSO: plsys, plsys_get */ { extern plsys_coord; plsys_coord = n; plsys,n; } func plsys_get(n) /* DOCUMENT plsys_get return the coordinate system last set with plsys_set This allows one to retrieve the current coordinate system so you can reset it after working in a particular coordinate system. SEE ALSO: plsys, plsys_set */ { extern plsys_coord; return plsys_coord; } struct fillc_struc { pointer zp, yp, xp, np, zc, yc, xc, nc; } func fillc(z, levs, x, y, contours=, miss_value=, miss_rect=, xperiodic=, nested=, maxseg=, maxpoly=, maxnode= ) /* DOCUMENT fillc( z, levs, x, y, contours=, miss_value=, miss_rect=, xperiodic=, nested=, maxseg=, maxpoly=, maxnode= ) FILLC locates contours of a function z(x,y), and returns a list of coordinates of polygons to be filled to accomplish contour-filled shading. (function plfp may be used to fill the polygons; e.g., f= fillc(z, levs); plfp, *f.zp, *f.yp, *f.xp, *f.np) Input: z -- 2-D array of values to be contoured levs -- contour levels array, with one or more values x, y -- 1-D or 2-D mesh array of X, Y coordinate values [If the X/Y coordinate values are omitted, a "logical" coordinate system (0,1,2,...) is used, consistent with the two dimensions of Z] (KEYWORD PARAMETERS: 0=> false (default), 1=> true) contours -- if true, return list contour nodes as well miss_value -- missing data value (0 => no missing values, the default) miss_rect -- if true, mark missing data values using rectangles, rather than using triangles xperiodic -- if true, data is assumed to be periodic in X direction nested -- if true, shade contours in a nested fashion (a bit slower, but reduces the number of shaded polygons) maxseg -- maximum number of contour segments (default 16384) maxpoly -- maximum number of contour polygons (default 16384) maxnode -- maximum number of contour nodes (default 262144) Output: A structure of the following form is returned by FILLC: struct fillc_struc { pointer zp, yp, xp, np, zc, yc, xc, nc; } ZP: pointer a 1-D integer array containing the "level" values associated with each filled polygon. The level values range from 0 to NUMBEROF(LEVS)+1: level value 0 denotes "below range" values, level values between 1 and NUMBEROF(LEVS)-1 denote within range values, level value NUMBEROF(LEVS) denotes "above range" values, and level value NUMBEROF(LEVS)+1 denotes missing value regions. YP,XP: pointer to 1-D arrays of the form ARRAY(FLOAT,NCORNER) containing the Y/X coordinates of NCORNER corners of polygons to be filled NP: number of corners for each successive polygon in the list [i.e., NCORNER == SUM(NP)] The following are values are returned only if parameters CONTOURS is true; otherwise, they are set nil values. ZC: pointer a 1-D integer array containing the "level" values associated with each contour. The level values range from' 0 to NUMBEROF(LEVS): level values between 0 and NUMBEROF(LEVS)-1 denote contour levels, and level value NUMBEROF(LEVS) denotes boundary of missing value regions. YC,XC: pointer to 1-D arrays of the form ARRAY(FLOAT,NNODE) containing the Y/X coordinates of NNODE contour nodes to be joined. NC: number of noded for each successive contour in the list; closed contours will have the first node repeated as the last node. [i.e., NNODE == SUM(NC)] NOTE: 1. If NESTED is true, polygons should be shaded in the order provided, with missing value regions being shaded last. 2. Shading should be completed prior to drawing the contour segments. * Last modified: 7 March 1997, R. Saravanan * * SEE ALSO: trishade, flexplot */ { // Set default values for input parameters if (is_void(contours)) contours= 0; if (is_void(miss_value)) miss_value= 0; if (is_void(maxseg)) maxseg= 16384; if (is_void(maxpoly)) maxpoly= 16384; if (is_void(maxnode)) maxnode= 262144; if (is_void(xperiodic)) xperiodic= 0; if (is_void(miss_rect)) miss_rect= 0; if (is_void(nested)) nested= 0; if (is_void(levs)) error, "No contour levels specified"; // Determine dimensions of data array zdims= dimsof(z); if (zdims(1) != 2) error, "Field not a 2-dimensional array" nx= zdims(2); ny= zdims(3); nlev= numberof(levs); // Generate "logical" coordinate for missing X/Y coordinates if (is_void(x)) x= span(0., zdims(2)-1., zdims(2)); if (is_void(y)) y= span(0., zdims(3)-1., zdims(3)); // Broadcast 1-D coordinate arrays to 2-D to generate a mesh if (dimsof(x)(1) == 1) x= x(,-:1:ny); if (dimsof(y)(1) == 1) y= y(-:1:nx,); // Check data consistency xdims= dimsof(x); ydims= dimsof(y); if ((xdims(1) != 2) || (xdims(2) != nx) || (xdims(3) != ny)) error, "X coordinate dimensions are not conformable with data dimensions"; if ((ydims(1) != 2) || (ydims(2) != nx) || (ydims(3) != ny)) error, "Y coordinate dimensions are not conformable with data dimensions"; if (xperiodic) if (max(abs(z(1,) - z(nx,))) != 0.) error, "Data not periodic in X direction" // Call RECCONTOUR to contour data polar= 0; reccontour, z, x, y, nx, nx, ny, xperiodic, polar, nested, miss_rect, levs, nlev, miss_value, maxseg, maxpoly, maxnode, nseg, iseg, jseg, lseg, npoly, ipoly, jpoly, lpoly, nnode, vnode; // Truncate polygon arrays xp= vnode(1,1:nnode); yp= vnode(2,1:nnode); zp= lpoly(1:npoly); np= jpoly(1:npoly); if (!contours || (nseg == 0)) { // Return output structure for polygon fills only return fillc_struc( zp=&zp, yp=&yp, xp=&xp, np=&np ); } else { // Contour details requested; allocate coordinate array for contour nodes vcont= array(float, 2, sum(jseg)); zc= array(long, nseg); nc= array(long, nseg); ncont= 0; ncnode= 0; ncnode_prev= 0; // Loop over the contour segments for (k=1; k<=nseg; k++) { // Check if current segment is linked to the next one if ( (k < nseg) && ((lseg(k) == lseg(k+1)) && ( vnode(1,iseg(k)+jseg(k)) == vnode(1,iseg(k+1)+1) ) && ( vnode(2,iseg(k)+jseg(k)) == vnode(2,iseg(k+1)+1) ) ) ) { // It is linked; copy current contour segment, omitting last node if (jseg(k) > 1) { vcont(,ncnode+1:ncnode+jseg(k)-1)= vnode(,iseg(k)+1:iseg(k)+jseg(k)-1); ncnode += jseg(k)-1; } // Go on to next contour segment, skipping rest of this iteration continue; } // Copy current segment vcont(,ncnode+1:ncnode+jseg(k))= vnode(,iseg(k)+1:iseg(k)+jseg(k)); ncnode += jseg(k); // Copy contour level value, node count ncont++; zc(ncont)= lseg(k); nc(ncont)= ncnode - ncnode_prev; ncnode_prev= ncnode; } // Truncate contour arrays xc= vcont(1,1:ncnode); yc= vcont(2,1:ncnode); zc= zc(1:ncont); nc= nc(1:ncont); // print, ncont, ncnode, nseg, sum(jseg(1:nseg)) // Return output structure for polygon fills and contours return fillc_struc( zp=&zp, yp=&yp, xp=&xp, np=&np, zc=&zc, yc=&yc, xc=&xc, nc=&nc ); } } /* *-------------------------------------------------------------------- * For GNU Emacs: * $Msg-digest-checksum: e2e2dc91670ad493e17dad6fbf246464 $ *-------------------------------------------------------------------- * Local Variables: * mode: C * msg-digest-active: t * End: *-------------------------------------------------------------------- */ func plfr (z, y, x, ireg, levs=, colors=, region=, abscolors=) /* DOCUMENT plfr, z,y,x,ireg, levs=levels, colors=colors A simplified interface to routines to determing color filled regions and plot them EXCLUDES the regions above and below the extremes of the contours: plot filled contours of the point centered values Z on the 2D mesh Y,X,IREG. The arguments and levs= keyword are as for the plc function, but plfc fills the regions between contours instead of drawing contour lines. The colors= keyword is a list of color indices; it should be one element SHORTER than levs. the first color is above levs(1) the last color is below levs(0) If colors= is not specified, you get numberof(levs)-1 colors in the current palette. If abscolors is specified, then we assume you want the regions to map directly to color indices. Otherwise you will get an even sampling of colors from the current palette KEYWORDS: region, SEE ALSO: fillc, plc, plfc */ { dz = dimsof(z); nx = dz(2); ny = dz(3); if (!is_void(colors) && (numberof(levs)-1 != numberof(colors)) ) { emsg = " numberof levs and colors mismatch " + pr1(numberof(levs)) + " " + pr1(numberof(colors)); error, emsg; } //maxseg= 1000; //maxpoly= 1000; //maxnode= 16317; if (!is_void(ireg) && !is_void(region)) { // print, "using the region version"; miss_value = max(abs(z))*1.1; spval = miss_value; // assign a missing value so no contours out of the region of interest; z2 = z; list = where (ireg != region); if (numberof(list) > 0) { z2(list) = spval; } pz = &z2; } else { pz = &z; } // info, pz; // read, prompt="debug",dummy; fillout = fillc(*pz, levs, x, y, maxseg=maxseg, maxpoly=maxpoly, maxnode=maxnode, contours=1,miss_value=miss_value); // exclude indices for special vals and outside of contour levels; // index 0 for values < levs(1) // index 1 for values between levs(1) and levs(2) // ... // index numberof(levs)-1 for values between levs(-1) and levs(0) // index numberof(levs) for values greater than levs(0) list = where(*fillout.zp < numberof(levs) & *fillout.zp > 0 ); np = []; zp = []; yp = []; xp = []; if (numberof(list) > 0) { culist = (*fillout.np)(cum)+1; for (j=1;j<=numberof(list);j++) { i = list(j); ib = culist(i); ie = ib + (*fillout.np)(i) - 1; grow, np, (*fillout.np)(i); grow, zp, (*fillout.zp)(i); grow, yp, (*fillout.yp)(ib:ie); grow, xp, (*fillout.xp)(ib:ie); } // print, np; // print, zp; // read,prompt="debug2",dummy; } if (!is_void(abscolors)) { zp = char(zp); } if (is_void(colors)) { plfp, zp, yp, xp, np; } else { plfp, char(colors(zp)), yp, xp, np; } // read,prompt="debugit",dummy; } func plfcs (z, y, x, ireg, levs=, colors=, region=, spv=) /* DOCUMENT plfcs, z,y,x,ireg, levs=levels, colors=colors A simplified interface to routines to determing color filled regions and plot them plot filled contours of the point centered values Z on the 2D mesh Y,X,IREG. The arguments and levs= keyword are as for the plc function, but plfc fills the regions between contours instead of drawing contour lines. The colors= keyword is a list of color indices; it should be one element longer than levs. the first color is for below levs(1) the last color is for above levs(0) If colors= is not specified, you get numberof(levs)+1 colors in the current palette. KEYWORDS: region, SEE ALSO: fillc, plc, plfc */ { dz = dimsof(z); nx = dz(2); ny = dz(3); //maxseg= 1000; //maxpoly= 1000; //maxnode= 16317; if (!is_void(ireg) && !is_void(region)) { // print, "using the region version"; z2 = z; // make a copy we can mess with miss_value = 2*max(z2); // assign a missing value so no contours out of the region of interest list = where (ireg != region); if (numberof(list) > 0) { z2(list) = miss_value; } fillout = fillc(z2, levs, x, y, maxseg=maxseg, maxpoly=maxpoly, maxnode=maxnode, contours=1,miss_value=miss_value); /* np = *fillout.np; zp = *fillout.zp; yp = *fillout.yp; xp = *fillout.xp; */ list = where(*fillout.zp != numberof(levs)+1); // indices where there are no special values // list = where(*fillout.zp != 1.e7); // indices where there are no special values ib = 1; np = []; zp = []; yp = []; xp = []; for (i=1;i<=numberof(list);i++) { ie = ib + (*fillout.np)(i) - 1; grow, np, (*fillout.np)(i); grow, zp, (*fillout.zp)(i); grow, yp, (*fillout.yp)(ib:ie); grow, xp, (*fillout.xp)(ib:ie); ib = ie+1; } // read,prompt="debug",dummy; } else { // print, "normal version"; fillout = fillc(z, levs, x, y, maxseg=maxseg, maxpoly=maxpoly, maxnode=maxnode, contours=1, miss_value=spv); np = *fillout.np; zp = *fillout.zp; yp = *fillout.yp; xp = *fillout.xp; } if (is_void(colors)) { plfp, zp, yp, xp, np; } else { plfp, char(colors(zp+1)), yp, xp, np; } // read,prompt="debugit",dummy; } func plc_labels (z, y, x, ireg, levs=, region=, condraw=, type=, ngradmax=, maxlabs=, labcol=, opaque=, width=, spv=) /* DOCUMENT plc_labels, z, y, x, ireg, levs=levels, type=, ngradmax=, maxlabs=, labcol=, opaque=, width=; draw labels (and optionally draw the contours themselves) of the point centered values Z on the 2D mesh Y,X,IREG. by using the contour vertices. Labels will not be drawn if gradients are too strong. contour vertices are currently returned by svn's fillc routine zone indices used by the routine are currently returned by another call to fillc this functionality might later be done more efficiently by direct calls to yorick routines (ie mesh_loc) named arguments levs= an array of contour levels to label (not optional) region= the region number to find labels in (need ireg array) condraw= if defined, then draw the contours themselves type= type of line if condraw defined ngradmax= a measure of the non-dimensional max gradient below which labels can be drawn (larger implies more labels) maxlabs= approximate max number labels to put on contour line segment (default 1) labcol= color to draw the label with (default black) opaque= if defined, put white down behind the label KEYWORDS: region, SEE ALSO: fillc, plc */ { local dz, nx, ny, xmax, xmin, ymax, ymin, xc, yc, zc; local xcl, ycl, zcl; local xlog, ylog; dz = dimsof(z); nx = dz(2); ny = dz(3); if (is_void(maxlabs)) maxlabs = 1; if (is_void(type)) type = "solid"; if (is_void(width)) width = 1; if (is_void(ngradmax)) ngradmax = 20; if (is_void(labcol)) labcol="black"; if (is_void(levs)) { error,"the levs keyword has not been set"; } if (sum(dimsof(z) == dimsof(y) == dimsof(x))) { error,"dimsof x, y, and z arrays are inconsistent"; } //maxseg= 1000; //maxpoly= 1000; //maxnode= 16317; xmax = max(x); xmin = min(x); ymax = max(y); ymin = min(y); small = (xmax-xmin)*1.e-7; if (!is_void(ireg) && !is_void(region)) { z2 = z // make a copy we can mess with; miss_value = 2*max(z2); // assign a missing value so no contours out of the region of interest; list = where (ireg != region); if (numberof(list) > 0) { z2(list) = miss_value; } // call fillc again using a logical coordinate system to determine the zone indices for // each contour node. These are used to get an estimate of the gradient of z at the contour // vertices xlog = span(1.,nx,nx); ylog = span(1.,ny,ny); fillout = fillc(z2, levs, x , y , maxseg=maxseg, maxpoly=maxpoly, maxnode=maxnode, contours=1, miss_value=miss_value, nested=1); filloutl = fillc(z2, levs, xlog, ylog, maxseg=maxseg, maxpoly=maxpoly, maxnode=maxnode, contours=1, miss_value=miss_value, nested=1); list = where(*fillout.zc != numberof(levs)) // indices where there are no special values; // move the returned values into arrays for the region of interest. ib = 1; nc = []; zc = []; yc = []; xc = []; zcl = []; ycl = []; xcl = []; for (i=1;i<=numberof(list);i++) { // info, xcl; // info, ycl; // info, zcl; ie = ib + (*fillout.nc)(i) - 1; grow, nc, (*fillout.nc)(i); grow, zc, (*fillout.zc)(i); grow, zcl, (*filloutl.zc)(i); grow, yc, (*fillout.yc)(ib:ie); grow, ycl, (*filloutl.yc)(ib:ie); grow, xc, (*fillout.xc)(ib:ie); grow, xcl,(*filloutl.xc)(ib:ie); ib = ie+1; } numnodes_contour_log = nc; numnodes_contour = nc; // read,prompt="debug",dummy; } else { fillout = fillc(z, levs, x, y, maxseg=maxseg, maxpoly=maxpoly, maxnode=maxnode, contours=1,nested=1,miss_value=spv); zc = *fillout.zc; yc = *fillout.yc; xc = *fillout.xc; numnodes_contour = *fillout.nc; // call fillc again using a logical coordinate system to determine the zone indices for ; // each contour node. These are used to get an estimate of the gradient of z at the contour; // vertices xlog = span(1.,nx,nx); ylog = span(1.,ny,ny); filloutl = fillc(z, levs, xlog, ylog, maxseg=maxseg, maxpoly=maxpoly, maxnode=maxnode, contours=1, nested=1); zcl = *filloutl.zc; ycl = *filloutl.yc; xcl = *filloutl.xc; numnodes_contour_log = *filloutl.nc; } n = numberof(numnodes_contour_log); // a measure of the gradient, normalized by the dimensions of the plot; grada = (xmax-xmin)*(abs(z(dif,)/(x(dif,)+small)))(pcen,) + (ymax-ymin)*(abs(z(,dif)/(y(,dif)+small)))(,pcen) ; // print, "mingrada, max ", min(grada), max(grada); // read,prompt="debug",dummy; mgrad = max(grada); grada(1:3,) = grada(-2:0,) = 1000.*mgrad; grada(,1:3) = grada(,-2:0) = 1000.*mgrad; /* // a measure of the curvature (2nd derivative) normalized by the dimensions of the plot; curva = array(0.,dimsof(z)); xx = x(zcen,) yy = y(,zcen) zx = ((z(dif,)/x(dif,))(dif,))/xx(dif,); zy = ((z(,dif)/y(,dif))(,dif))/yy(,dif); curva(2:-1,) = abs(zx) + curva(2:-1,); curva(,2:-1) = abs(zy) + curva(,2:-1); curva(1:3,) = curva(-2:0,) = 1000.*max(curva) curva(,1:3) = curva(,-2:0) = 1000.*max(curva) // read,prompt="debug curva",dummy; */ labx = []; laby = []; lablev = []; nlabs = 0; kb = 1; for (k=1; k<=n; k++) { ke = kb + numnodes_contour(k) - 1; izone = int(xcl(kb:ke)+small); jzone = int(ycl(kb:ke)+small); ndm = izone + nx*(jzone-1); xa = xc(kb:ke); ya = yc(kb:ke); cz = zc(k)+1; // print, " plotting contour, level number, val, nnodes ", k, cz, levs(cz), numnodes_contour(k); if (!is_void(condraw)) plg, yc(kb:ke), xc(kb:ke), marks=0, type=type,width=width; // a measure of the wiggliness of the contour; wig = array(0.,dimsof(xa)); if (numnodes_contour(k) > 2) { dr = sqrt(xa(dif)^2 + ya(dif)^2)(cum); drd = dr(dif)+small; drz = dr(zcen); drzd = drz(dif)+small; // xx = xa(zcen); // yy = ya(zcen); y2 = ((ya(dif)/drd)(dif))/drzd; x2 = ((xa(dif)/drd)(dif))/drzd; wig(2:-1) = sqrt(y2*y2 + x2*x2) + wig(2:-1,); // print, "y2", y2; // print, "x2", x2; // print, "wig bef smooth", wig; wig = wig(zcen)(pcen); // print, "wig aft smooth", wig; wig = wig*10/(avg(wig)+1.e-20); wig(1:3,) = wig(-2:0,) = 1000.*max(wig); // print, "renormalized wig", wig; // read, prompt="wig debug",dummy; } // print, " at wig1, number nodes ", numberof(ndm); // a measure of the number of contours in the vicinity of this point; czp1 = min(numberof(levs),cz+1); czm1 = max(1, cz-1); // print, " czp1, czm1, levs, normal ", czp1, czm1, levs(czp1), levs(czm1), // (czp1-czm1)/(levs(czp1)-levs(czm1)); // print, "min, max of grada ", min(grada(ndm)), max(grada(ndm)); if (czp1 != czm1) { // ngrad will be a measure of the number of contours which would ; // span the x/y domain, if the gradient were constant at ; // a value of grada(ndm1) everywhere; // so a number greater than about 10 implies there are too many contours around to be good; // ie the gradient is too steep ngrad = grada(ndm)*(czp1-czm1)/(levs(czp1)-levs(czm1)); list = where( ngrad < ngradmax); } else { list = array(1,numberof(ndm)); } if (numberof(list) == 0) { // print, "ndm", ndm; // print, " the gradient is too steep for this contour, skipping ", ngrad; // read, prompt="steep debug", dummy; kb = ke + 1; continue; } ndm = ndm(list); xa = xa(list); ya = ya(list); wig = wig(list); // print, " at wig, number nodes ", numberof(ndm); // dont even consider any contours which are very small; extent = (max(xa)-min(xa))/(xmax-xmin) + (max(ya)-min(ya))/(ymax-ymin) ; // print, " extent is ", extent; numplt = 0; if (extent > 0.1) { while (numberof(ndm) > 0) { // now check the potential label sites against previous labels; // and eliminate any that are too close; if (nlabs > 0) { dist = abs(xa-labx(-,))/(xmax-xmin) + abs(ya-laby(-,))/(ymax-ymin); mdist = dist(,min); cmax = array(.10,numberof(xa)); // read, prompt="debug", dummy; // list = where((dist - cmax)(,min) > 0.); list = where((mdist > cmax)); // print, " reducing numberof nodes closeness from ", numberof(xa), "to ", numberof(list); ndm = ndm(list); xa = xa(list); ya = ya(list); wig = wig(list); mdist = mdist(list) } if (numberof(ndm) > 0) { // now pick out the minimum gradient for the node to actually plot; // print, " reduced to ", numberof(ndm), " points"; if (czp1 != czm1) { ngrad = grada(ndm)*(czp1-czm1)/(levs(czp1)-levs(czm1)); } else { ngrad = array(1.,numberof(ndm)); } // print, " ngrad ", ngrad; // print, " wig ", wig; // a normalized version of the curvature with a mean value of order 1; // ncurva = curva(ndm)*10/avg(curva(ndm)); // print, " ncurva ", ncurva; if (nlabs == 0) { //works penalty = ngrad+ncurva; penalty = ngrad+wig; ndist = array(0.,dimsof(wig)); } else { ndist = 1./mdist; // print, " penalty for being near another label ", ndist; //works penalty = ngrad+ncurva+ndist; penalty = ngrad+wig+ndist; } // print, " total penalty ", penalty; list =penalty(mnx); // print, " min at ", list, ngrad(list), wig(list), ndist(list), penalty(list); xl = xa(list); yl = ya(list); // wig = wig(list); // print, "plotting a label", levs(cz); plt, pr1(levs(cz)), xl, yl, tosys=1, justify="CH",color=labcol, opaque=opaque; // read,prompt=" just plotted", dummy grow, labx, xl; grow, laby, yl; nlabs = nlabs + 1; numplt = numplt + 1; if (numplt >= maxlabs) { // print, " numplt >= maxlabs "; break; } } else { // print, "ndm is zero"; } } // read, prompt="debug labels...", dummy; kb = ke + 1; } else { // print, " contour too short, skipping"; kb = ke + 1; } } } func find_nice_contours(inar, nlevs=) /* DOCUMENT find_nice_contours (inar, nlevs) return a nice set of approximately nlevs contours given data array inar if you want precisely nlevs contours use find_ctlabs */ { if (is_void(nlevs)) nlevs= 6; zmax = max(inar); zmin = min(inar); if (zmax == zmin) { if (zmax != 0) { return [0.9*zmax,zmax,1.1*zmax]; } else { return [-0.1,0.0, 0.1]; } } zrange = zmax/(zmin+1.e-20); if (zrange > 100 && zrange == -1) { print, " find_nice_contours log section not written "; return []; } else { okints = [0.1, 0.2, 0.5,1,2,5,10,20]; // acceptable values for contour intervals zinc = (zmax-zmin)/nlevs; pow = int(log10(zinc)); cints = okints*(10.^pow); nlevsout = (zmax-zmin)/cints; nlevbest = abs(nlevsout-nlevs)(mnx); zminout = short(zmin/cints(nlevbest))*cints(nlevbest); zmaxout = short(zmax/cints(nlevbest))*cints(nlevbest); ninc = long((zmaxout-zminout)/cints(nlevbest) + 1.00000001); xlabs = span(zminout, zmaxout, ninc); list = where(abs(xlabs/max(abs(xlabs))) < 1.e-5); if (numberof(list)> 0) xlabs(list) = 0; return xlabs; } } func find_ctlabs (z=, nlabs=, zmin=, zmax=) /* DOCUMENT find_ctlabs(z=, nlabs=, zmin=, zmax=) find nice colortable labels input: z input array nlabs number of labels to produce zmin zmax output: returns nlabs numbers which span zmin and zmax and are nice for contours or colortable labels */ { okints = [0.1,0.2, 0.25, 0.5, 1,2,2.5, 5,10,20] // good intervals ; if (is_void(zmin)) zmin = min(z); if (is_void(zmax)) zmax = max(z); if (is_void(nlabs)) nlabs = 11; zinc = (zmax-zmin)/nlabs; pow = int(log10(zinc)); cints = okints*(10.^pow); nlabsout = (zmax-zmin)/cints; nlevbest = where(nlabsout-nlabs < 0)(1); zminout = (short(zmin/cints(nlevbest))-1)*cints(nlevbest); zmaxout = zminout +(nlabs-1)*cints(nlevbest); ninc = long((zmaxout-zminout)/cints(nlevbest) + 1.00000001); xlabs = span(zminout, zmaxout, ninc); list = where(abs(xlabs/max(abs(xlabs))) < 1.e-5); if (numberof(list)> 0) xlabs(list) = 0; return xlabs; } func get_hfile ( fname, &fhandle, &nchandle,force=,vars=,rmafiles=, rmhista= ) /* DOCUMENT get_hfile (fname, fhandle, nchandle, force=, vars=, rmafiles=, rmhista= ) Open a file and return the filehandle, and extra netcdf info The file name (fname) is expected to be a ccm history tape. The subroutine will actually open a netcdf file with a similar name (ncfilename), if it exists. If ncfilename does not exist, get_hfile will create it by a) acquiring the file with ccm2nc b) converting it from fname to ncfilename c) if vars is non-void, only vars will be translated from fname to ncfilename If ncfilename exists, and rmafiles is not void, both fname and ncfilename will first be removed, re-acquired, and reconverted. If force is not void, ncfilename will be removed, and ccm2nc will be called again to reconvert the netcdf file. (this allows different variables to be put onto ncfilename) If rmhista is not void, the original history tape after the conversion will be removed to save room. SEE ALSO: ccm2nc (not a yorick routine) NOTE:!!!!! ccm2nc must be on your path for this routine to work */ { locname = strpart(fname,2:) + ".nc"; flocname = strpart(fname,2:) // name used by ccm2nc for history tape; print, "trying to open local file " + locname; if (!is_void(vars)) { vars = "-v " + vars + " "; } else { vars = " "; } if (!is_void(rmafiles)) { print, " rmafiles set, removing " + flocname + " and " + locname; remove, flocname; remove, locname; force = 1; } if (open(locname,"",1) && is_void(force)) { print, locname + " exists"; } else { // print, "converting ccm format " + fname + // " to netcdf file " + locname; xxx = "ccm2nc " + vars + fname + " " + locname; print, xxx; system, xxx; if (!is_void(rmhista)) { print, "conversion complete, removing history file ", flocname; remove, flocname; } } fhandle = openb(locname); nchandle = nc_file; // print, " opened ", locname; // info, fhandle; // info, nchandle; } func set_red_blue_ct { red = [255,0,0]; white = [255,255,255]; blue = [0,0,255]; blue = [100,100,255]; ct3 = char([red,white,blue]); ct200 = interp(ct3, span(0.,1.,3), span(0.,1.,200),2); palette, ct200(1,), ct200(2,), ct200(3,); } func plot_hbar (levs, colors, offset=) /* DOCUMENT plot_hbar, levs, colors, offset= plot a horizontal color bar below the viewport levs defines the levels seperating the region colors define the index into the color bar to be used between each level offset= the fraction of the frame height to offset the color bar by (0.3) if (is_void(levs) || is_void(colors)) error," both levs and colors must be defined"; */ { // coords of viewport port= viewport(); xd = port(2)-port(1); yd = port(4)-port(3); if (is_void(offset)) offset = 0.3; //print, "offset = ", offset; xl = port(1); xr = xl + xd; yb = port(3)-yd*offset; yt = yb + yd*offset/3.; ylabs = yb - 0.02; nlevs = numberof(levs); nc = nlevs-1; zcol = array(short,1,nc); zcol(1,) = indgen(1:nc:1); //plot colors lastsys = plsys_get(); plsys, 0 if (!is_void(colors)) { pli, char(transpose(colors(zcol))), xl, yb, xr, yt; } else { pli, transpose(zcol), xl, yb, xr, yt; } //plot labels; if (nc > 10) { // we want no more than 10 labels num_colors_per_label = nc/10. nlabs = int(nc/(int(num_colors_per_label)+1)); nskip = nlevs/nlabs; // number of levels to skip; //read, prompt="debug",dummy; } else { nlabs = nc; nskip = 1; } //print, "nlevs, nlabs, nskip", nlevs, nlabs, nskip; xndc = span(xl,xr,nlabs); yndc = ylabs; for (i=1;i<=nlevs;i=i+nskip) { // xx = xl + (levs(i)-levs(1))/(levs(0)-levs(1))*xd; xx = xl + (i-1)*xd/(nlevs-1) plt,pr1(levs(i)),xx,ylabs,height=12,justify="CH"; // print, "xl, xx, xr", xl, xx, xr // print, "i, nlevs ", i, nlevs; } pldj, xl, yb, xr, yb; pldj, xl, yt, xr, yt; pldj, xl, yb, xl, yt; pldj, xr, yb, xr, yt; plsys_set, lastsys; } func tick_labs_off /* DOCUMENT turn off ticks and labels on the current viewport */ { require, "style.i"; cpl = plsys_get(0); get_style, landscape, systems, legends, clegends; systems(cpl).ticks.horiz.flags=0x000 // turn off ticks and labels systems(cpl).ticks.vert.flags=0x000 // turn off ticks and labels set_style, landscape, systems, legends, clegends; plsys_set(cpl) }