;=========================================================== ; IDL Program to plot stuff ;=========================================================== PRO Animation ;----------------------------------------------------------- ; Set-up 16 colours + 0=black, 1=white ;----------------------------------------------------------- r = [0,255,0,0,0,0,66,184,245,255,255,255,230,191,148,105,161,184] g = [0,255,84,199,255,255,255,255,255,209,135,28,0,0,0,0,0,0] b = [0,255,255,255,255,178,0,0,0,0,0,0,0,0,0,20,135,186] TVLCT, r,g,b ;----------------------------------------------------------- ; Loop over each input data file ;----------------------------------------------------------- data3d = fltarr(80,40,241) dataname = 'tracer' zeroline = 0 plusminus = 0 times = 176 corefile = 'ForeStep' FOR time = 1, times DO BEGIN IF (time LT 10) THEN BEGIN file = corefile+'00'+STRING(time,FORMAT='(I1)')+'.nc' ENDIF ELSE BEGIN IF (time LT 100) THEN BEGIN file = corefile+'0'+STRING(time,FORMAT='(I2)')+'.nc' ENDIF ELSE BEGIN file = corefile+STRING(time,FORMAT='(I3)')+'.nc' ENDELSE ENDELSE PRINT, 'Reading in ',file long = 'longitude' lat = 'latitude' ;----------------------------------------------------------- ; Read data from netCDF file ;----------------------------------------------------------- fileid = NCDF_OPEN(file) varidx = NCDF_VARID(fileid, long) varidy = NCDF_VARID(fileid, lat) varid = NCDF_VARID(fileid, dataname) NCDF_VARGET, fileid, varidx, x_degrees NCDF_VARGET, fileid, varidy, y_degrees NCDF_VARGET, fileid, varid, data NCDF_CLOSE, fileid data3d(*,*,time-1) = data(*,*) ENDFOR ;----------------------------------------------------------- ; Set-up the map ; /MERCATOR, /CYLINDRICAL, /SATELLITE, etc. ; /NOBORDER ;----------------------------------------------------------- Px1 = 0.03 Px2 = Px1 + 0.7 Py1 = 0.1 Py2 = 0.9 Lx = Px2-Px1 Ly = Py2-Py1 xinc = x_degrees[1] - x_degrees[0] yinc = y_degrees[1] - y_degrees[0] minx = x_degrees[0] maxx = x_degrees[79] midx = (maxx+minx)/2.0 miny = y_degrees[0] maxy = y_degrees[39] midy = (maxy+miny)/2.0 ;----------------------------------------------------------- ; Set-up postscript file ;----------------------------------------------------------- SET_PLOT, 'ps' ;----------------------------------------------------------- ; Set-up colour grades array ;----------------------------------------------------------- minval = MIN(data3d(*,*,*)) maxval = MAX(data3d(*,*,*)) IF (plusminus EQ 1) THEN BEGIN ;We expect the data to be roughly balanced in +/- extreme values IF (ABS(minval) GT ABS(maxval)) THEN BEGIN maxval = ABS(minval) ENDIF ELSE BEGIN minval = -ABS(maxval) ENDELSE ENDIF range = maxval-minval grades = FINDGEN(16)*range/16 + minval ;----------------------------------------------------------- ; Do the plotting ;----------------------------------------------------------- FOR time = 1, times DO BEGIN ; Set-up the postscript filename IF (time LT 10) THEN BEGIN psfile = 'Testps/State000'+STRING(time,FORMAT='(I1)')+'.ps' ENDIF ELSE BEGIN IF (time LT 100) THEN BEGIN psfile = 'Testps/State00'+STRING(time,FORMAT='(I2)')+'.ps' ENDIF ELSE BEGIN psfile = 'Testps/State0'+STRING(time,FORMAT='(I3)')+'.ps' ENDELSE ENDELSE PRINT, 'Plotting ', psfile DEVICE, FILE=psfile, /PORTRAIT, /COLOR, XSIZE=15.0, YSIZE=10.0, $ XOFFSET=1.0, YOFFSET=1.0, FONT_SIZE=8, /ENCAPSULATED !P.POSITION = [Px1, Py1, Px2, Py2] ; Develop a new title that includes other information (max, min) TitleA = dataname + "t = " + STRING(time) + " " + $ '!C Min = ' + STRING(minval) + $ '!C Max = ' + STRING(maxval) MAP_SET, 0.0, 180.0, LIMIT=[miny,minx,maxy,maxx], TITLE=TitleA, $ /CYLINDRICAL, /NOERASE ;----------------------------------------------------------- ;Fill cells with colour ;----------------------------------------------------------- FOR ix = 0, 79 DO BEGIN FOR iy = 0, 39 DO BEGIN ;Calculate the colour of the box boxcolour = MAX(WHERE(data3d(ix,iy,time-1) GE grades)) ;calculate the limits of the element to plot xmin = x_degrees[ix] - xinc/2.0 xmax = x_degrees[ix] + xinc/2.0 ymin = y_degrees[iy] - yinc/2.0 ymax = y_degrees[iy] + yinc/2.0 myxcoords = [xmin,xmax,xmax,xmin,xmin] myycoords = [ymin,ymin,ymax,ymax,ymin] ;correct limits if they go over the bounds pts = WHERE(myycoords GT maxy,count) IF (count GT 0) THEN myycoords(pts) = maxy pts = WHERE(myycoords LT miny,count) IF (count GT 0) THEN myycoords(pts) = miny pts = WHERE(myxcoords GT maxx,count) IF (count GT 0) THEN myxcoords(pts) = maxx pts = WHERE(myxcoords LT minx,count) IF (count GT 0) THEN myxcoords(pts) = minx POLYFILL, myxcoords, myycoords, COLOR=boxcolour+2 ENDFOR ENDFOR ;----------------------------------------------------------- ; Add zero contour ;----------------------------------------------------------- IF (zeroline EQ 1) THEN BEGIN zerolevel = FINDGEN(1) & zerolevel(0) = 0.0 linestyle = INDGEN(1) & linestyle(0) = 1 ;linecolours = INDGEN(16) & linecolours = 0 ;All contours black ;linecolours = INDGEN(16)+2 ;Contours many colours CONTOUR, data3d(*,*,time-1), x_degrees, y_degrees, LEVELS=zerolevel, $ xstyle=1, ystyle=1, nlevels=1, $ /OVERPLOT, C_COLORS=linecolours, C_LINESTYLE=linestyle ENDIF ;----------------------------------------------------------- ; Add continents ;----------------------------------------------------------- MAP_CONTINENTS, /NOERASE ;----------------------------------------------------------- ; Add axes ;----------------------------------------------------------- blank = STRARR(10) blank(*) = ' ' AXIS, COLOR=0, XAXIS=0, TICKLEN=-0.015, XRANGE=[minx,maxx], XTITLE=' ' AXIS, COLOR=0, XAXIS=1, TICKLEN=-0.015, XRANGE=[minx,maxx], XTICKNAME=blank AXIS, COLOR=0, YAXIS=0, TICKLEN=-0.015, YRANGE=[miny,maxy], YTITLE=' ' AXIS, COLOR=0, YAXIS=1, TICKLEN=-0.015, YRANGE=[miny,maxy], YTICKNAME=blank ;----------------------------------------------------------- ; Add colourbar ;----------------------------------------------------------- COLBAR, [Px2+0.1, Py1, Px2+0.15, Py2], grades, ' ' DEVICE, /CLOSE ENDFOR END