;;; NCL script for creating contour plots of NARCCAP RCM data ;;; ;;; Usage: Invoke using nclwrap or via command line like so: ;;; ncl -x varname=\"var\" infile=\"filename.nc\" [etc.] plot.ncl ;;; ;;; Other options that can be specified: timestep (defaults to last ;;; step in file); title (defaults to filename, plus timestamp of ;;; timestep for time-varying data); colormap (defaults to ;;; nrl_sirkes); wkstype (defaults to ps); fillmode (defaults to ;;; raster); global (defaults to False, unless data is lat-lon) ;;; ;;; Auto-detects time-varying vs 2-D data and lat-lon vs projected ;;; data and Does The Right Thing for each case. 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" ;; open input file fin = addfile(infile, "r") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; detect input type and set defaults ;; time-varying or static? istime = any("time" .eq. getfilevarnames(fin)) if (istime) then time = fin->time ;; 'julian' is legit (no missing leap year every 100th year), ;; but NCL doesn't know it if (isatt(time, "calendar") .and. time@calendar .eq. "julian") then time@calendar = "gregorian" end if if (.not. isvar("timestep")) then timestep = dimsizes(time) - 1 end if if (.not. isvar("title")) then if (isatt(time,"calendar")) then date = ut_calendar(time(timestep),0) timestamp = ""+date(0,0)+"/"+date(0,1)+"/"+date(0,2)+" "+sprintf("%02.0f", date(0,3))+":"+sprintf("%02.0f", date(0,4)) title = systemfunc("basename "+infile)+", "+timestamp else title = systemfunc("basename "+infile)+", "+sprintf("%0.3f",time(timestep)) end if end if end if ;; overrideable plotting defautls if (.not. isvar("title")) then title = systemfunc("basename "+infile) end if if (.not. isvar("colormap")) then colormap = "nrl_sirkes" ;; amwg_blueyellowred ? end if if (.not. isvar("outfile")) then wkstype = "x11" outfile = systemfunc("basename "+infile) end if if (.not. isvar("wkstype")) then wkstype = "ps" end if if (.not. isvar("fillmode")) then fillmode = "RasterFill" end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; read in data lat = fin->lat lon = fin->lon if (istime) then data = fin->$varname$(timestep,:,:) else data = fin->$varname$ end if ;; global or projected plot? if (.not. isvar("global")) then if (isatt(data, "grid_mapping")) then global = False else global = True end if end if if (isatt(data, "grid_mapping")) then data@lat2d = lat data@lon2d = lon end if if (.not. global) nx = dimsizes(fin->xc) ny = dimsizes(fin->yc) pname = data@grid_mapping proj = fin->$pname$ end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; map projection & bounds res = True if(global) then res@mpLimitMode = "LatLon" res@mpMinLatF = min(lat) res@mpMaxLatF = max(lat) res@mpMinLonF = min(lon) res@mpMaxLonF = max(lon) res@mpCenterLonF = (max(lon) + min(lon))/2 else if (lower_case(proj@grid_mapping_name) .eq. "lambert_conformal_conic") then res@mpProjection = "LambertConformal" res@mpLambertMeridianF = proj@longitude_of_central_meridian res@mpLambertParallel1F = proj@standard_parallel(0) res@mpLambertParallel2F = proj@standard_parallel(1) end if if (lower_case(proj@grid_mapping_name) .eq. "transverse_mercator") then res@mpProjection = "Mercator" res@mpCenterLatF = proj@latitude_of_projection_origin res@mpCenterLonF = proj@longitude_of_central_meridian end if if (lower_case(pname) .eq. "polar_stereographic") then res@mpProjection= "Stereographic" res@mpCenterLonF= proj@straight_vertical_longitude_from_pole - 360 res@mpCenterLatF = proj@latitude_of_projection_origin end if if (lower_case(proj@grid_mapping_name) .eq. "rotated_latitude_longitude") then res@mpProjection = "Mercator" res@mpCenterLonF = proj@grid_north_pole_longitude - 180 res@mpCenterLatF = 90 - proj@grid_north_pole_latitude end if ;; map boundaries ;; expand slighltly to be sure we plot things on the borders res@mpLimitMode = "Corners" res@mpLeftCornerLatF = lat(0,0) ;; SW corner res@mpLeftCornerLonF = lon(0,0) ;; SW corner res@mpRightCornerLatF = lat(yc|ny-1,xc|nx-1) ;; NE corner res@mpRightCornerLonF = lon(yc|ny-1,xc|nx-1) ;; NE corner end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; general plot settings res@mpFillOn = False res@cnFillMode = fillmode res@cnFillOn = True res@cnLinesOn = False res@gsnAddCyclic = False res@gsnMaximize = True res@gsnSpreadColors = True res@lbBoxLinesOn = False res@lbLabelAutoStride = True res@tiMainString = title if (wkstype .eq. "x11") then ;; res@gsnFrame = False end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; actually create the plot wks = gsn_open_wks(wkstype, outfile) gsn_define_colormap(wks,colormap) plot = gsn_csm_contour_map(wks, data, res) exit