PRO polycntr1, filename, Z, X, Y, color_index = color_index, pattern = pat $ , delete_file = delfile ;+ ; NAME: POLYCNTR1 ; PURPOSE: ; Fill the contours defined by a path file created by CONTOUR. ; CATEGORY: ; Graphics. ; CALLING SEQUENCE: ; POLYCNTR1, Filename, Z [, X, Y, COLOR_INDEX = color index] ; INPUTS: ; Filename = name of file containing contour paths. This ; file must have been made using CONTOUR, PATH=Filename, ... ; ; Z = 2 dimensional data array used to generate the contours. ; ; X = vector or 2 dimensional array specifying the X coordinates ; for the contour surface. If a vector, each element specifies the ; X coordinate for a column of Z (Example: X(0) specifies the X ; coordinate for Z(0,*)). If X is a two-dimensional array, each ; element of X specifies the X coordinate of the corresponding ; point in Z ( X(i,j) specifies the X coordinate for Z(i,j)). ; ; Y = as for X above, but for the Y coordinates. ; ; ; KEYWORD PARAMETERS: ; COLOR_INDEX = color index array. Element i contains the color ; of contour level number i-1. Element 0 contains the ; background color. There must be one more color index than ; levels. ; PATTERN = optional array of patterns, dimensioned (NX, NY, NPATTERN). ; DELETE_FILE If present and non-zero, delete FILENAME once contours ; are read. ; OUTPUTS: ; The contours are filled on the display using normalized ; coordinates using the POLYFILL procedure. ; COMMON BLOCKS: ; None. ; SIDE EFFECTS: ; Contour temp file filename.tmp created in current directory. ; Display is updated. ; PROCEDURE: ; Scans file. Eliminates contours of 1 or 2 points. Open contours ; are closed clockwise and counterclockwise around the plot border ; and both contours are written to the contour temp file. Area ; of the two plots and starting bytes in the contour temp file ; are both noted. The 2 dimensional array Z is used to determine ; which of these 2 contours encloses the higher ground. Closed ; contours are written directly to the contour temp file. Areas ; are sorted and contours are reread from the contour temp file ; and drawn in order of decreasing area. ; EXAMPLE: ; Create a 10 by 10 array of random numbers, contour to get ; path file, polycontour it, then overdraw the contour lines: ; b = RANDOMU(seed, 10,10) ;Insert random numbers ; CONTOUR, b, /SPLINE, PATH = 'path.dat' ;Make path file ; POLYCNTR1, 'path.dat',b ;Fill contours ; CONTOUR, b, /SPLINE, /NOERASE ;Overplot lines & labels ; ; Note: result looks best if the plot region is filled with data. ; Setting xstyle=1, ystyle=1, on CONTOUR accomplishes this. ; Contours are only be filled for the region spanned by the data. ; MODIFICATION HISTORY: ; DMS, AH, January, 1989. ; PAS = P. A. Sterne, sterne1@llnl.gov ; PAS, LLNL, October 1991. Shades in order of contour area. ; PAS, LLNL, November 1992. Shades open contours too. ; PAS, LLNL, March 1993. Added in X,Y, and bugfixes. ; PAS, LLNL, Sept 1993. Fill background when all contours are closed. ; Fix for IDL 3.0.0 pathfile format. ; PAS, LLNL, June 1994. Bugfix to IDL-specific component. ;- header = {contour_header,$ type : 0B, $ high_low : 0B, $ level : 0, $ num : 0L, $ value : 0.0 } ; ; ******************** IMPORTANT ************************** ; ; set the variable pvoridl to 0 if you are using PVWave or early ; versions of IDL ( before 3.0.0, I think). Set it to 1 if you are ; using later versions of IDL. IDL changed the way in which contours ; are written, necessitating this change. ; pvoridl = 0 ; PVWave or early IDL versions ;pvoridl = 1 ; IDL version 3.0.0 and later ; tol = 1.e-5 asize = 100 ;# of elements in contour arrays n = 0 cval = 0 ;Contour index carea = fltarr(asize) ;Contour area cstart = lonarr(asize) ;Starting byte of record dtog0 = [!X.s(0), !Y.s(0)] dtog1 = [!X.s(1), !Y.s(1)] pos = fltarr(4) sz = size(Z) IF sz(0) NE 2 THEN BEGIN print, 'POLYCNTR1: Data array is not 2 dimensional. Returning' RETURN ENDIF regularxy = 0 ;zero for x,y vectors or not passed CASE n_params() OF 2: BEGIN ; X and Y not passed to POLYCNTR1 x1 = findgen(sz(1)) ; make X a vector y1 = findgen(sz(2)) ; make Y a vector END 3: BEGIN ; X passed as argument, but not Y szx = size(x) x1 = x IF szx(0) EQ 1 THEN y1 = findgen(sz(2)) $ ; make Y a vector ELSE IF szx(0) EQ 2 THEN BEGIN y1 = replicate(1.0, sz(1))#findgen(sz(2)) ; make Y an array regularxy = 1 ENDIF END 4: BEGIN ;X and Y both passed as arguments x1 = x y1 = y szx = size(x) szy = size(y) IF max([szx(0), szy(0)]) EQ 2 THEN BEGIN regularxy = 1 IF szx(0) EQ 1 THEN x1 = x#replicate(1.0, sz(2)) ; make X an array IF szy(0) EQ 1 THEN y1 = replicate(1.0, sz(1))#y ; make Y an array ENDIF END ENDCASE IF regularxy EQ 0 THEN BEGIN pos = [dtog0+[x1(0), y1(0)]*dtog1, dtog0+[x1(sz(1)-1), y1(sz(2)-1)]*dtog1] scalexy = [min(abs(x1-shift(x1, 1))),min(abs(y1-shift(y1,1)))]*0.01 ENDIF ELSE BEGIN isx = indgen(sz(1)) isy = sz(1)*(indgen(sz(2)-2)+1) bpts = [isx,isy+sz(1)-1,rotate(isx,2)+sz(1)*(sz(2)-1),rotate(isy,2)] nbpts = 2*(sz(1)+sz(2)-2) xy2 = [[x1(bpts)], [y1(bpts)]] z2 = z(bpts) xyn = fltarr(2, nbpts) FOR i = 0, nbpts-1 DO xyn(*, i) = xy2(i, *)*dtog1 + dtog0 scale = 0.01*min(sqrt(sum((shift(xy2, 1, 0)-xy2)^2, 1))) corners = [ 0, sz(1)-1,sz(1)+sz(2)-2,2*sz(1)+sz(2)-3] xpoly = fltarr(4) zpoly = fltarr(4) ENDELSE xyvert = pos([[0, 1], [2, 1], [2, 3], [0, 3]]) vlist = shift(pos,-1) openr,unit,filename,/get_lun,delete=KEYWORD_SET(delfile) openw,unit2,filename+'.tmp',/get_lun,/delete while (not eof(unit)) do begin ;First pass reading contours IF n EQ asize THEN BEGIN ;Expand our arrays? cstart = [cstart, cstart] carea = [carea, carea] asize = 2 * asize ENDIF a = fstat(unit2) ;File position readu, unit, header ;Read header IF (header.type EQ 0 AND header.num GT 1) THEN BEGIN ; Open contours cval = cval > fix(header.level) ; contour level xyarr = fltarr(2, header.num) ; Contour polygon hn = header.num - 1 header1 = header header2 = header header1.type = 1 ;Reset header type to closed contour. header2.type = 1 ;Reset header type to closed contour. IF (pvoridl EQ 0) THEN readu, unit, xyarr $ ELSE BEGIN xyarrp = transpose(xyarr) readu, unit, xyarrp xyarr= transpose(xyarrp) ENDELSE IF regularxy EQ 0 THEN BEGIN ; determine which side the open contour enters (side1) and exits (side2). xlst = [[1, 0, 1, 0], [2*hn+1, 2*hn, 2*hn+1, 2*hn]] sid1 = where(abs(xyarr(xlst(*, 0))-vlist) LT tol, count1) sid2 = where(abs(xyarr(xlst(*, 1))-vlist) LT tol, count2) IF count1 LE 0 AND header.num GT 2 THEN print, $ 'error open contour side1, n= ',n,' vertex=', xyarr(*, 0) IF count2 LE 0 AND header.num GT 2 THEN print, $ 'error open contour side2, n= ',n,' vertex=', xyarr(*, hn) IF (count1 LE 0) OR (count2 LE 0) THEN GOTO, RDREC ; close open contour side1 = max(sid1) ; side contour enters side2 = max(sid2) ; side contour exits vadd = shift(xyvert, 0, -side2-1) tmp = reverse(vadd, 2) ; close contour counter-clockwise (xyarr1) and clockwise (xyarr2) ivt = (-side2+side1+4) MOD 4 IF ivt EQ 0 THEN BEGIN ; same side dotdir = total((xyarr(*, hn)-xyarr(*, 0))*(vadd(*, 3)-vadd(*, 0))) IF dotdir LT 0 THEN ivt = 4 - ivt ENDIF xyarr1 = fltarr(2, hn+2+ivt) xyarr2 = fltarr(2, hn+6-ivt) header1.num = hn + 2 + ivt header2.num = hn + 6 - ivt IF ivt NE 0 THEN xyarr1(*, hn+1:hn+ivt) = vadd(*, 0:ivt-1) IF ivt NE 4 THEN xyarr2(*, hn+1:hn+4-ivt) = tmp(*, 0:3-ivt) xyarr1(*, 0:hn) = xyarr xyarr2(*, 0:hn) = xyarr xyarr1(*, header1.num-1) = xyarr(*, 0) xyarr2(*, header2.num-1) = xyarr(*, 0) ; Check which of xyarr1/2 is higher. sdxy = side2 MOD 2 ; 0 for side2=0,2; 1 for side2=1,3 sdsign = -2*(side2/2) + 1 ; +1 for side2=0,1 and -1 side2=2,3 dta = (xyarr(*, hn)-dtog0)/dtog1 IF sdxy EQ 0 THEN BEGIN vmin = min(abs(x1-dta(0)), it1) it1 = it1 > 1 < (sz(1)-2) it2 = ([-1, 0, 1]+it1) xpoly = [x1(it2), dta(0)] - dta(0) zpoly = [z(it2, (side2/2)*(sz(2)-1)), header.value] ENDIF ELSE BEGIN vmin = min(abs(y1-dta(1)), it1) it1 = it1 > 1 < (sz(2)-2) it2 = ([-1, 0, 1]+it1) xpoly = [y1(it2), dta(1)] - dta(1) zpoly = [reform(z(((sdsign+1)/2)*(sz(1)-1), it2)), header.value] ENDELSE xsplin = xpoly(sort(xpoly)) zsplin = zpoly(sort(xpoly)) - header.value dtsplin = [-1., 1.]*scalexy(sdxy) val = rotate(spline(xsplin, zsplin, dtsplin), 2-2*(side2/2)) ind = [0B, 0B] FOR i = 0, 1 DO IF val(i) GT 0. THEN ind(i) = 1b IF ind(0) EQ ind(1) THEN print, ' both sides min or max, n=', n ; Check which of xyarr1/2 is higher. Add & subtract 0.005 of the +ve side ; vector to the last xyarr point. Adding puts point in xyarr1, subtracting ; puts it in xyarr2. If xyarr1 pt > xyarr2 pt, then color = c for xyarr1 ; and c-1 for xyarr2. Reverse for xyarr1 pt < xyarr2 pt. ENDIF ELSE BEGIN ; 2 dimensional X and Y arrays ; determine the boundary points on either side of the open contour ; where it enters (ptba1) and exits (ptba2). dta1 = (xyarr(*, 0)-dtog0)/dtog1 dta2 = (xyarr(*, hn)-dtog0)/dtog1 vmin1 = min(((xy2(*,0)-dta1(0))^2+(xy2(*,1)-dta1(1))^2), it1) vmin2 = min(((xy2(*,0)-dta2(0))^2+(xy2(*,1)-dta2(1))^2), it2) itp1 = (it1 + 1) MOD nbpts itm1 = (it1 - 1 + nbpts) MOD nbpts itp2 = (it2 + 1) MOD nbpts itm2 = (it2 - 1 + nbpts) MOD nbpts grdp1 = xy2(itp1, *) - xy2(it1, *) grdm1 = xy2(itm1, *) - xy2(it1, *) grdp2 = xy2(itp2, *) - xy2(it2, *) grdm2 = xy2(itm2, *) - xy2(it2, *) diff1 = dta1 - xy2(it1, *) diff2 = dta2 - xy2(it2, *) IF total(diff1*diff1) GT tol THEN BEGIN dir1 = [total(grdp1*diff1), total(grdm1*diff1)] a1 = min(dir1/[total(grdp1*grdp1), total(grdm1*grdm1)], indir1) ptba1 = [it1-indir1, it1-indir1+1] ; bndry pt [below,above] ENDIF ELSE ptba1 = [it1, it1] indir2 = 0 IF total(diff2*diff2) GT tol THEN BEGIN dir2 = [total(grdp2*diff2), total(grdm2*diff2)] a1 = max(dir2/[total(grdp2*grdp2), total(grdm2*grdm2)], indir2) ptba2 = [it2-indir2, it2-indir2+1] ; contour boundary pt. ENDIF ELSE ptba2 = [it2, it2] ; close contour counter-clockwise (xyarr1) and clockwise (xyarr2) xy3 = shift(xyn, 0,-ptba2(1)) ptba1m = (ptba1 + nbpts - ptba2(1)) MOD nbpts xyarr1 = [[xyarr], [xy3(*, 0:ptba1m(0))], [xyarr(*, 0)]] xyarr2 = [[xyarr], [reverse(xy3(*, ptba1m(1):nbpts-1), 2)], $ [xyarr(*, 0)]] header1.num = (size(xyarr1))(2) header2.num = (size(xyarr2))(2) ; Check which of xyarr1/2 is higher. pt2 = ([it2-1, it2, it2+1]+nbpts) MOD nbpts FOR i = 0 , 3 DO IF (it2 EQ corners(i)) THEN BEGIN IF indir2 EQ 0 THEN pt2 = (pt2+1) MOD nbpts IF indir2 EQ 1 THEN pt2 = (pt2-1+nbpts) MOD nbpts ENDIF FOR i = 0, 2 DO xpoly(i) = sqrt(total((xy2(pt2(i), *)-dta2)^2)) xpoly(0) = -xpoly(0) xpoly(1) = (2*indir2-1)*xpoly(1) xpoly(3) = 0.0 zpoly(0:2) = z2(pt2) zpoly(3) = header.value xsplin = xpoly(sort(xpoly)) zsplin = zpoly(sort(xpoly)) dtsplin = [-1., 1.]*scale val = rotate(spline(xsplin, zsplin, dtsplin), 2) ind = [0B, 0B] FOR i = 0, 1 DO IF val(i) GT header.value THEN ind(i) = 1b IF ind(0) EQ ind(1) THEN print, $ ' 2-dim: both sides min or max, n=', n ENDELSE header1.high_low = ind(0) cstart(n) = a.cur_ptr ;Position carea(n) = poly_area(xyarr1(0, *), xyarr1(1, *)) writeu, unit2, header1 writeu, unit2, xyarr1 n = n + 1 a = fstat(unit2) ;File position header2.high_low = ind(1) cstart(n) = a.cur_ptr ;Position carea(n) = poly_area(xyarr2(0, *), xyarr2(1, *)) writeu, unit2, header2 writeu, unit2, xyarr2 n = n + 1 ENDIF ELSE IF (header.num LE 2) THEN BEGIN ;Skip short contours xyarr = fltarr(2, header.num) IF (pvoridl EQ 0) THEN readu, unit, xyarr $ ELSE BEGIN xyarrp = transpose(xyarr) readu, unit, xyarrp xyarr= transpose(xyarrp) ENDELSE ENDIF ELSE BEGIN ;Closed contour cval = cval > fix(header.level) ;Color to draw cstart(n) = a.cur_ptr ;Position xyarr = fltarr(2, header.num) ;Define point to skip IF (pvoridl EQ 0) THEN readu, unit, xyarr $ ELSE BEGIN xyarrp = transpose(xyarr) readu, unit, xyarrp xyarr= transpose(xyarrp) ENDELSE carea(n) = poly_area(xyarr(0, *), xyarr(1, *)) writeu, unit2, header ;Write header writeu, unit2, xyarr n = n + 1 ENDELSE RDREC: ENDWHILE ; End of contour file free_lun, unit ;print, 'cval,n=', cval,n ; max contour index ( = #contours-1) and # cstart = cstart(0:n-1) ; Truncate carea = carea(0:n-1) order = rotate(sort(carea), 2) ;Subscripts of order c_index = bindgen(cval+2)*byte(float(!D.n_colors)/(cval+1)) if keyword_set(color_index) eq 0 then color_index = c_index if max(color_index) lt 0 then color_index = c_index ;print,color_index IF regularxy EQ 0 THEN xyv1 = [[xyvert], [xyvert(*, 0)]] $ ELSE xyv1 = [[xyn], [xyn(*, 0)]] col = color_index(0) polyfill,/normal, color=col, xyv1 ;shade in background color for i=0,n-1 do begin ;Draw each contour j = order(i) ;Index of record point_lun,unit2,cstart(j) readu,unit2,header ;Reread header ;print,'n,level,high_low,value',j,header.level,header.high_low,header.value xyarr = fltarr(2, header.num) ;Define polygon array readu,unit2,xyarr ;Read polygon inds = fix(header.level) + fix(header.high_low) col = color_index(inds) ;Drawing color ;print,'color',col IF N_ELEMENTS(pat) NE 0 THEN BEGIN s = size(pat) if s(0) ne 3 then begin print,'POLYCONTR1 - pattern array not 3d' return endif polyfill,/NORMAL, pattern=pat(*,*, i mod s(3)), xyarr endif else $ polyfill, /NORMAL, color= col,xyarr ;Fill contour endfor free_lun, unit2 ;Done end