func plc_svn (z, y, x, ireg, region=, spv=, levs=, type=, width=, concol=, lab=, ngradmax=, maxlabs=, labcol=, opaque=, height=) /* DOCUMENT plc_svn, z, y, x, ireg, region=, spv=, levs=, type=, width=, concol=, ngradmax=, maxlabs=, labcol=, opaque=, height= Draw contours of the point centered values Z on the 2D mesh Y,X,IREG. by using the contour vertices returned by svn's fillc routine. Y,X,IREG are optional. named arguments: region= the region number to contour within spv= special value for missing data levs= an array of contour levels type= type of line to draw (default solid) width= width of line to draw (default 1) concol= color of line to draw (default black) lab= label the contours (no labels if zero) (default=1) ngradmax= a measure of the non-dimensional max gradient below which labels can be drawn (larger implies more labels), default=20 maxlabs= approximate max number labels for contour line segment (default 1) labcol= color of label (default black) opaque= if defined, put white down behind the label height= size of characters KEYWORDS: region, SEE ALSO: fillc, plc */ { local dz, nx, ny, xc, yc, zc; local xmin, xmax, ymin, ymax, small, id, jd, grada; local labx, laby, lablev, nlabs; // set void input arguments to default if (is_void(levs)) levs = set_contours(z,nlevs=10,spvspv); if (is_void(type)) type = "solid"; if (is_void(width)) width = 1; if (is_void(concol)) concol="black"; if (is_void(lab)) lab = 1; if (is_void(ngradmax)) ngradmax = 20; if (is_void(maxlabs)) maxlabs = 1; if (is_void(labcol)) labcol="black"; // check dimensions and set X,Y arrays if necessary dz = dimsof(z); nx = dz(2); ny = dz(3); if (is_void(x)) { x = indgen(dz(2))(,-:1:dz(3)); print, " plc_labels: x array not set, using indices ", dz(2); } if (is_void(y)) { y = indgen(dz(3))(-:1:dz(2),); print, " plc_labels: y array not set, using indices ", dz(2); } if (sum(dimsof(z) == dimsof(y) == dimsof(x))) { error,"plc_labels: dimsof x, y, and z arrays are inconsistent"; } // deal with regions as missing values but do not outline them if (!is_void(ireg) && !is_void(region)) { // use missing values to mimic a masked region z2 = z; // make a copy we can mess with // set z values outside specified region to missing value if (is_void(spv)) spv = 2*max(z2); list = where (ireg != region); if (numberof(list) > 0) { z2(list) = spv; } fillout = fillc(z2, levs, x, y, contours=1, nested=1, miss_value=spv, maxseg=maxseg, maxpoly=maxpoly, maxnode=maxnode); // list of contours, not region (or special value) boundaries list = where(*fillout.zc != numberof(levs)); // store the nodes for contours in the list ib = 1; nc = zc = yc = xc = []; for (i=1;i<=numberof(list);i++) { ie = ib + (*fillout.nc)(i) - 1; grow, nc, (*fillout.nc)(i); grow, zc, (*fillout.zc)(i); grow, yc, (*fillout.yc)(ib:ie); grow, xc, (*fillout.xc)(ib:ie); ib = ie+1; } numnodes_contour = nc; // read,prompt="debug",dummy; } else { fillout = fillc(z, levs, x, y, contours=1, nested=1, miss_value=spv, maxseg=maxseg, maxpoly=maxpoly, maxnode=maxnode); zc = *fillout.zc; yc = *fillout.yc; xc = *fillout.xc; numnodes_contour = *fillout.nc; } // check for nonzero number of nodes to draw if (numberof(numnodes_contour) == 0) { //print,"early return"; return; } // initialize contour labelling if (lab != 0) plc_svn_li,xmin, xmax, ymin, ymax, small, id, jd, grada, labx, laby, lablev, nlabs; // draw the contours kb = 1; for (k=1; k<=numberof(numnodes_contour); k++) { ke = kb + numnodes_contour(k) - 1; cz = zc(k)+1; // print, " plotting contour, level number, val, nnodes ", // k, cz, levs(cz), numnodes_contour(k); plg, yc(kb:ke), xc(kb:ke), marks=0, type=type,width=width,color=concol; // label the contour if (lab != 0) plc_svn_lc; kb = ke + 1; } } func plc_svn_li (&xmin, &xmax, &ymin, &ymax, &small, &id, &jd, &grada, &labx, &laby, &lablev, &nlabs) { require, "digit2.i"; xmax = max(x); xmin = min(x); ymax = max(y); ymin = min(y); small = (xmax-xmin)*1.e-7; // find the zones for each node in the contours nd = digit2(yc, xc, y, x); list = where(nd <= 0); if (numberof(list) > 0) nd(list) = 1; // anything outside the zones is just point in first zone id = nd%(nx-1); jd = (nd - id)/(nx-1) + 1; // 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; } func plc_svn_lc { izone = id(kb:ke); jzone = jd(kb:ke); ndm = izone + nx*(jzone-1); xa = xc(kb:ke); ya = yc(kb:ke); cz = zc(k)+1; // 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, "grada", grada(ndm); // print, "czp1, czm1, levs ", czp1, czm1, levs(czp1), levs(czm1); // print, " the gradient is too steep for this contour, skipping ", ngrad; // read, prompt="steep debug", dummy; return; } 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) { if (numplt >= maxlabs) { // print, " numplt >= maxlabs "; break; } // 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, height=height; // read,prompt=" just plotted", dummy grow, labx, xl; grow, laby, yl; nlabs = nlabs + 1; numplt = numplt + 1; } else { // print, "ndm is zero"; } } } }