;=========================================================== ; IDL Program to plot stuff ;=========================================================== PRO BasicPlot ;----------------------------------------------------------- ; 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 ;----------------------------------------------------------- data = fltarr(80,40) name = 'tracer' zeroline = 1 plusminus = 0 calcbounds = 0 file = 'Background.nc' 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, name) NCDF_VARGET, fileid, varidx, x_degrees NCDF_VARGET, fileid, varidy, y_degrees NCDF_VARGET, fileid, varid, data NCDF_CLOSE, fileid ;----------------------------------------------------------- ; 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 ;----------------------------------------------------------- minvalactual = MIN(data(*,*)) maxvalactual = MAX(data(*,*)) IF (calcbounds EQ 1) THEN BEGIN minval = minvalactual maxval = maxvalactual ENDIF ELSE BEGIN ; minval = -7.3e-5 ; maxval = 0.00024 minval = 935.0 maxval = 1095.0 ENDELSE 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 ;----------------------------------------------------------- ; Set-up the postscript filename psfile = 'Plot.ps' 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 = name + $ '!C Min = ' + STRING(minvalactual) + $ '!C Max = ' + STRING(maxvalactual) 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(data(ix,iy) 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, data(*,*), 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 END