;----------------------------------------------------------------------------- ; ; png-demo.ncl -- Demo direct PNG output, raster file format, from NCL. ; ; 2009-jan-19 Adapted from plot-smc.0114.ncl. ; By Dave Allured, NOAA/ESRL/PSD/CAB, CU/CIRES/CDC. ; Tested with pre-release ncl.ppc.5.0.1, generated on ; or before 2008-may-07. Platform = Mac OS 10.6 PPC. ; ; Command summary: ; ; ncl png-demo.ncl [switches] ; ; Output switches: pdf=1, png=1, ps=1, etc. Default is X11. See table below. ; ; lambert=1 for Lambert Conformal Conic projection. Default is rectangular. ; ; contour=1 to plot regular shaded contours. Default is grid cells (boxes). ; ; ti=nnn: Time index, zero based, to select the desired time step. ; Use nctime utility to look up the time index for the desired date. ; If missing, defaults to time step 160 (1998-jun-10), for demo purposes. ; ;----------------------------------------------------------------------------- ; ; Post processing notes: ; ; 1. Make GIF from PNG, if desired, with ImageMagick "convert". ; Our current version appears to retain full quality with lossless ; conversion from PNG to GIF. ImageMagick version 6.4.1 01/13/09 ; Q16, http://www.imagemagick.org. ; ; convert infile.png outfile.gif ; ; 2. For web applications etc., you can strip off excess white ; borders with a simple trim operator in convert: ; ; convert -trim infile.png outfile.png ; ; (1) and (2) can be used in the same command line, with no loss ; in quality. ; ;----------------------------------------------------------------------------- load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" begin ; Get command line parameters. ; Caveat: Some switches are not found here, only used lower down ; by direct reference: e.g. -- if (isdefined ("lambert")) then ... devices = (/ "ps", "eps", "pdf", "png", "x11" /) output = "x11" ; default = X11 window do i = 0, dimsizes (devices) - 1 ; find output selector in list; if (isdefined (devices(i))) then ; specify only one output = devices(i) break end if end do if (isdefined ("contour")) then fill_mode = "AreaFill" ; plot contours, not grid cells else fill_mode = "CellFill" ; plot grid cells (boxes) end if ; Run parameters. ; Data file may be downloaded by anonymous FTP here: ; ftp://ftp.cdc.noaa.gov/Datasets/cpc_us_precip/precip.1998.nc ; Base URL is subject to change some time after 2009-jan-20. infile = "/Datasets/cpc_us_precip/precip.1998.nc" varname = "precip" outfile = "demo" ; omit file extension, NCL will add it monthnames = (/ "Jan", "Feb", "Mar", "Apr", "May", "Jun", \ "Jul", "Aug", "Sep", "Oct", "Nov", "Dec" /) cmap = "WhViBlGrYeOrRe" ; specify color map ; Increase PNG resolution. Ref. ncl-talk message, 2008-aug-17, by David Brown. if (output .eq. "png") then ; must be square; deviations output@wkWidth = 3000 ; either way are wasted output@wkHeight = 3000 end if ;----------------------------------------------------------------------------- ; Open gridded input file. print ("") print ("Open input file.") print (infile+"") in = addfile (infile, "r") dims = getfilevardimsizes (in, varname) ; get data array dimensions ntimes = dims(0) ny = dims(1) nx = dims(2) print ("nx = " + nx) print ("ny = " + ny) print ("ntimes = " + ntimes) ; Check specified time index. print ("") ti_last = ntimes - 1 ; last time step in file ;; ti_default = ti_last ; set default time step ti_default = 160 ; specific time step for demo if (.not. isdefined ("ti")) then ; ti command arg is missing: print ("*** ti is not specified; use default time step.") ti = ti_default end if if (ti .lt. 0 .or. ti .gt. ti_last) then ; ti is out of range: print ("*** ti is out of range; use default time step.") ti = ti_default end if ; Get the date for this time index. time_coord = in->time(ti) ; read time coordinate value opt_ymd = 0 ymd = ut_calendar (time_coord, opt_ymd) y = floattoint (ymd(0,0)) ; convert to integer date/time m = floattoint (ymd(0,1)) d = floattoint (ymd(0,2)) date_str = y + "-" + monthnames(m-1) + sprinti ("-%02i", d) ; formatted nicely for titles date_str2 = y + sprinti ("-%02i", m) + sprinti ("-%02i", d) ; fixed format YYYY-MM-DD for file and menu item names print (sprinti ("Specified time index = %i", ti) + ", date = " + date_str) ; Read data grid at the specified time step. grid = short2flt (in->$varname$(ti,:,:)) ; read and unpack one grid print ("Min, max data for this grid (unpacked) = " \ + min (grid) + ", " + max (grid)) ;----------------------------------------------------------------------------- ; Create plot. print ("") print ("Create plot.") wks = gsn_open_wks (output, outfile) ; open output file or window gsn_define_colormap (wks, cmap) ; assign color map bkgnd = 0 ; index for background color (usually white) res = True res@mpDataBaseVersion = "Ncarg4_1" ; higher res data base res@mpOutlineOn = True ; turn on map outlines res@gsnPaperOrientation = "portrait" ; only works for ps, pdf res@gsnMaximize = True ; expand to fit view frame ; Specify color fill and contour intervals. res@cnFillOn = True ; enable shaded color res@cnLinesOn = False ; turn off contour lines res@cnLineLabelsOn = False ; turn off contour line labels res@cnFillMode = fill_mode ; contours or grid cells, above res@cnLevelSelectionMode ="ManualLevels" ; manual contour levels res@cnMinLevelValF = 0.10 ; min contour level res@cnMaxLevelValF = 1.40 ; max contour level res@cnLevelSpacingF = 0.10 ; contour spacing res@gsnSpreadColors = True ; use full range of colors ; res@gsnSpreadColorStart = 10 ; set lowest color ; res@gsnSpreadColorEnd = -1 ; set highest color, -1 = last ; Titles. res@tiMainString = "CPC Daily Precipitation for " + date_str res@gsnLeftString = "0.25 degree data" res@tiMainFontHeightF = 0.016 ; font height of title res@gsnStringFontHeightF = 0.010 ; font height of subtitles ; Control rectangular plot border. res@pmTickMarkDisplayMode = "Always" ; this gob is needed because res@tmXTOn = False ; it seems to be the only way res@tmXBOn = False ; to make the plot outline res@tmYLOn = False ; come out right with Lambert. res@tmYROn = False ;; res@tmBorderThicknesssF = 5. ; doesn't work ; Specify map projection. res@mpLimitMode = "LatLon" ; common resources res@gsnAddCyclic = False ; mandatory for subset plots if (isdefined ("lambert")) then res@mpProjection = "LambertConformal" res@mpLambertParallel1F = 29.5 res@mpLambertParallel2F = 45.5 res@mpLambertMeridianF = -96 res@mpMinLatF = 25 ; for lambert conic projection res@mpMaxLatF = 50 res@mpMinLonF = -120 res@mpMaxLonF = -73 res@tiDeltaF = 2.0 ; must adjust title position else res@mpProjection = "CylindricalEquidistant" res@mpMinLatF = 24 ; for rectangular projection res@mpMaxLatF = 50 res@mpMinLonF = -126 res@mpMaxLonF = -66 res@tiDeltaF = 0 ; must adjust title position end if print ("Map projection = " + res@mpProjection) print ("Contour mode = " + fill_mode) print ("Color map = " + cmap) ; Apply mask for US only. ; This is actually a negative mask. First, the data grid is plotted ; everywhere. Then the whole map is filled with the background color ; white, except within the mask specifier boundaries. ; ; Note that drawing country and state outlines is controlled ; explicitly, rather than by drawing and filling. print ("Apply US area mask.") res@mpMaskAreaSpecifiers = "united states" ; enable masking; specify area ; Ncarg4_1 database required res@mpFillAreaSpecifiers = (/ "water","land" /) ; fill all except above res@mpSpecifiedFillColors = (/ bkgnd, bkgnd /) ; fills=background color res@mpOutlineBoundarySets = "USStates" ; draw US and state outlines res@cnFillDrawOrder = "Predraw" ; contours first; then fills ; blank out unwanted areas ; Make plot. plot = gsn_csm_contour_map (wks, grid, res) print ("Done.") end