;************************************************ ; ncl_lab_atm401.ncl ; ; Author - Philippe Papin ; Date - 13 Feb 2014 ; ; About - This program is a script that will ; plot a series of variables in order to ; make a map ; ; ; ;************************************************* ; When first starting out an ncl program, you have to read in multiple "Root" scripts that allow you to use the proper command lingo ; that makes the codes of ncl work properly. This should go at the top of your program even before you use the "begin" command to ; initiate your script. If you create ncl functions, you would also load these in here using the proper directory. ; ; Don't modify unless you know a particular function needs this script set deleted/added ;*************PROTOTYPE SCRIPTS******************* ; ; These are "prototype scripts" that are necessary to read ; in other functions for ncl. If a function is not ; running properly (you get an error) check to make sure ; its not because you didn't call the prototype script ; at the beginning of your script. ; ;************************************************* 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" load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl" ;*****************Beginning of Program*************** ; ; ; All ncl scripts start with a "begin" and finish with ; an "end". This signifies the start and stop place ; for the script to run. ; ;***************************************************** begin ;*****************Program Format********************** ; ; ; To start things off, we want to have a basic ; structure of our program. NCL traditionally is ; structured in such a way as bulleted below: ; ; 1.) Read in data ; ; 2.) Adjust data as needed so it can be use in ; the program as designed ; ; 3.) Create and modify plot options for image ; creation (if you are making an image) ; ; 4.) Match data with plot options to create your ; image as desired ; ; Not all program will obviously take this format ; (especially if you are not creating images). ; However it always a good idea to know what your ; program is expected to do so you know how your code ; should be best structured. This sample script will ; create a simple sea level pressure map of the arctic. ; It is meant to be easily configurable so you can use ; this example code and expand upon it to create many ; other images ; ; Other Important Things: ; ; Dataset: How you read in your data will often be ; determined on what dataset you are reading in. ; In some cases, reading in is relatively easy ; (e.g. most NetCDF files). With this and other ; data sources you can easily find out what data ; is contained in the file buy using this command: ; ; (note you do this in command line outside of program) ; ; ncl_filedump "name of file.nc" ; ; where data is presented in this format ; ; array type / name / (dimension name 1,dimension name 2...dimension name N) ; ; a real example looks like this: ; ; short ir ( time, lat, lon ) ; ; typically there is plenty of descriptors underneath that ; will explain that array dataset in detail ; ; This sample script will read in CFSR data which ; we have available in house in netCDF format. I'll ; provide the directory where you find this data. ; ; ;***************************************************** ;***************************************************** ;********************((STEP ONE))********************* ;***************************************************** ;*****************((READ IN THE DATA))**************** ;***************************************************** ;*****************Reading In Time********************** ; ; Time can sometimes be tricky within ncl. Depending ; on the dataset, time can be just one period in a ; dataset, or one dataset array can contain a bunch of times. ; For the CFSR dataset, the files are divided in yearly intervals ; from 1979 to 2010. We are just going to focus on a 2 day period ; in this sample script. ; ;****************************************************** ; Time units in these arrays are not done in yyyymmddhh format ; because they must go up sequentially. Thus they are put in ; values like # of days/hours/minutes since a certain date in the ; past. However, we can easily convert our yyyymmddhh time frame ; into such a time period that can be read by the array ;Time units used by CFSR timeUnits = "hours since 1800-01-01 00:00:00" ;Our start and end dates for the period we want to read in from the dataset ; ; NOTE: You can set these values as functions ahead of time and just call ; the function after ; ; ex. syyyy = 2010 (then put in syyyy) in where the year would normally go ;*******************; ; ; ; DATE SECTION ; ; ; ;*******************; syyyy = ; enter starting year smm = ; enter starting month sdd = ; enter starting day shh = ; enter starting hour eyyyy = ; enter ending year emm = ; enter ending month edd = ; enter ending day ehh = ; enter ending hour sdate = cd_inv_calendar(syyyy,smm,sdd,shh,00,00,timeUnits,0) ;puts starting time into same time units as CFSR data edate = cd_inv_calendar(eyyyy,emm,edd,ehh,00,00,timeUnits,0) ;puts ending time into same time units as CFSR data ; Next we want to know how many data periods we will have so we can properly catalogue ; each plot and assign it the proper date in the imagery. Note CFSR only has data every ; 6 hours. ntimes = toint(edate-sdate)/6+1 ;timestepfactor = 4 ;in 6hr intervals ; ; We just modify the ntimes variable here to create a new array "times" for use later on times = new( (/ntimes/), "float") times@units = timeUnits do n=0,ntimes-1 times(n) = tofloat(sdate + (6*n)) end do ; Print statements are very useful to troubleshoot your code. You can print out what ; certain variables are like the example below print("reading in files from year: "+syyyy) ; Note when you just want to print simple text not in a function use the "put text here" ; to discriminate it from traditional functions. ; OUTPUT DIRECTORY: This is the important place where you are sending your output! ; it makes sense to do this early on so that you can easily modify ; where you want your files to go ;*********************; ; ; ; OUTPUT SECTION ; ; ; ;*********************; ;imgoutdir = "/nfs/lb13/ppapin/classes/atm401/ncl_lab/output/" imgoutdir = ; FILE NAME: In addition, you also want to give your file a useful name so you can ; tell what type of output it holds (ex. "img_slp.png" for an image of sea level pressure) plotName = "500_vort_vvel" plotType = "png" ;*****************Reading In Data********************** ; ; You are finally at the stage where we are reading in ; an actual gridded dataset. Again we are using CFSR so our location ; of our input file should be noted below. ; ;****************************************************** ;*********************; ; ; ; DATA SECTION ; ; ; ;*********************; print("Reading in Variables") filename_u = "/cfsr/data/"+syyyy+"/u."+syyyy+".0p5.anl.nc" filename_v = "/cfsr/data/"+syyyy+"/v."+syyyy+".0p5.anl.nc" filename_g = "/cfsr/data/"+syyyy+"/g."+syyyy+".0p5.anl.nc" filename_w = file_u = addfile(filename_u,"r") file_v = addfile(filename_v,"r") file_g = addfile(filename_g,"r") file_w = ; Next, we need to specify what variable we are reading in. We are interested in pressure at mean sea level (pmsl) ; ; most datasets take this array structure. If you want to read in the entire array just use a colin (:) ; You must specify all arrays now. For example the pmsl variable in this CFSR data file has three dimentions (time, lat, lon) ; In this case we are going to read in our specified time (as we defined above) and all the latitudes and longitudes u_500 = file_u->u({sdate:edate},{500},:,:) ; Zonal Winds v_500 = file_v->v({sdate:edate},{500},:,:) ; Meridional Winds g_500 = file_g->g({sdate:edate},{500},:,:) ; Geopotential Height w_500 = ; Vertical Velocity ; And thats it, you just read in the file! Congrats if you have made it so far! ; This bottom piece below is a little bit of house cleaning, but ncl if not specified will display the long name and ; units of a variable on the top right and left sides of the image. In order to clear that (we will do this ourselves later) ; just clear these fields out by stating the variable (ex. u_500) and "@" with the other characteristics of that variable. ; A clear field is simply "" u_500@long_name = "" u_500@units = "" v_500@long_name = "" v_500@units = "" g_500@long_name = "" g_500@units = "" w_500@long_name = "" w_500@units = "" ;***************************************************** ;********************((STEP TWO))********************* ;***************************************************** ;*****************((ADJUST THE DATA))***************** ;***************************************************** ; This section is where you modify the data to be ; what you need it to be. For example you might have read in ; multiple fields and need to do calculations with those fields ; to get the value you want to plot or use for another program. ; In some cases you also want to multiply/divide a value to get the proper units ; In this case in order to get decameters we need to divide our g_500 values by 10 since ; They are currently in meters (m) ;************************; ; ; ; CALCULATE SECTION ; ; ; ;************************; g_500 = g_500/10 ; In some cases, there are built in ncl functions that will ; aid in adjusting the data. For example, we don't have ; vorticity files within CFSR because you can easily ; obtain the vorticity through a simple ncl function ; that is designed to calculate the value globally ; given that you have u and v winds. vor_500 = vor_500 = u_500 = v_500 = ; printMinMax(vor_500,True) ; printMinMax(g_500,True) ; printMinMax(w_500,True) ; ; for more information on the different functions ; already created within ncl visit this link: ; ; http://www.ncl.ucar.edu/Document/Functions/list_type.shtml ;***************************************************** ;********************((STEP THREE))******************* ;***************************************************** ;*******************((PLOT OPTIONS))****************** ;***************************************************** ; This is the part that might get a little overwhelming at times. ; NCL is well known for its lavish options you have when creating ; a plot. There are literally so many options its impossible to ; memorize them all. This is why its very convient that ncl has such ; a well organized guide and list of functions/plot options that you ; can look at to see how they modify the plot. In order to best learn ; how ncl can benefit you, take your time to mess around with all the ; little nuisances on how you can perfect your plot. ;*****************Resources********************** ; ; Resources are the key aspect of modifying your plot ; options. The way they work is that you are designing ; a set of options to use for each variable. In this case ; its most effective to use one set of resources with one ; variable and another set of resources with another variable. ; In addition, its better to also have a "master resource" ; block that all of the resources are then saved to for the ; plot. This master resource will do the inital work setting ; up the background of the plot and the tickmarks / latitude ; markers ect. Let the sub resources focus on the individual ; variables that you will have on your plot. ; ; Learn more about resources here: ; ; http://www.ncl.ucar.edu/Document/Graphics/Resources/ ; ;****************************************************** ;***********************; ; ; ; RESOURCE SECTION ; ; ; ;***********************; ; To start a new resource set it as true res = True ; This makes sure that the plot is not drawn or framed too early res@gsnDraw = False res@gsnFrame = False res@gsnLeftString = "" res@gsnRightString = "" ; Here is where you set the other resources equal to the "master resource" ; which in this case is "res" g_500res = res vor_500res = res w_500res = res wind_500res = res ; The master resource is not a fill map so therefore it won't need to have fill options res@mpFillOn = False ; Choose the latitude and longitude domain you want your image to be located in. Lets focus ; on the eastern United States because thats where the interesting weather is going on res@mpMinLatF = ; specify the plot domain res@mpMaxLatF = ; res@mpMinLonF = ; res@mpMaxLonF = ; ; Have fun with these options. This change the style of the tickmarks and overall ; color and style of the land borders. Again go on ncl's site to see all the ; available options for each resource res@pmTickMarkDisplayMode = "Always" res@mpOutlineOn = True ; turn the map outline on res@gsnDraw = False ; do not draw the plot res@gsnFrame = False ; do not advance the frame ;res@gsnAddCyclic = False ; longitude coordinate array is not cyclic so make false res@mpGridAndLimbOn = True ; turn on grid lines res@mpGridLatSpacingF = 10.0 ; spacing of grid lat lines res@mpGridLonSpacingF = 10.0 ; spacing of grid lon lines res@mpGridLineDashPattern = 2 ; choose dash pattern of lat/lon lines res@mpGridLineThicknessF = 2.0 ; choose thickness of grid dash lat/lon lines res@tmXBLabelFontHeightF = .01 res@mpOutlineBoundarySets = "GeophysicalAndUSStates" res@mpDataSetName = "Earth..4" res@mpGeophysicalLineThicknessF = 2.0 res@mpDataBaseVersion = 1 ;medium resolution res@mpDataResolution = 2 ;medium data resolution ;res@gsnMaximize = True ; maximize the plot res@gsnPaperOrientation = "landscape" ; when maximizing, keep the orientation as portrait res@mpGridAndLimbDrawOrder = "Postdraw" ; Set the resources for the first variable. ; Once again you need to set the resource true ; even though we set it equal to the master ; "res" resource. The resource here will only ; modify the attributes associated with the ; variable it is associated with. ;;;;; ;;;;; vor_500res@gsnPaperOrientation = "landscape" ; when maximizing, keep the orientation as portrait vor_500res@cnLevelSelectionMode = "ManualLevels" ; manually select the levels vor_500res@cnMinLevelValF = 6 ; set the minimum contour level vor_500res@cnMaxLevelValF = 36 ; set the maximum contour level vor_500res@cnLevelSpacingF = 2. ; set the contour interval vor_500res@cnFillOn = True ; turn on contour shading vor_500res@cnLinesOn = False ; don't draw the contour lines vor_500res@cnInfoLabelOn = False ; don't draw the info label vor_500res@cnLineLabelsOn = False ; don't draw the contour line labels vor_500res@tiMainFontHeightF = 0.012 ; set the main title font height vor_500res@cnFillDrawOrder = "PreDraw" vor_500res@lbAutoManage = False vor_500res@pmLabelBarHeightF = 0.05 vor_500res@pmLabelBarWidthF = 0.25 vor_500res@pmLabelBarOrthogonalPosF = 0.00 vor_500res@pmLabelBarParallelPosF = 0.30 vor_500res@lbLabelAutoStride = False vor_500res@lbLabelStride = 2 vor_500res@lbTopMarginF = 1.3 vor_500res@lbBottomMarginF = -1.3 vor_500res@lbLabelFontHeightF = 0.007 vor_500res@lbLabelOffsetF = 0.1 vor_500res@lbBoxLinesOn = False vor_500res@lbTitleString = "[10~S~-5~N~s~S~-1~N~]" vor_500res@lbTitleFontHeightF = 0.007 vor_500res@lbTitlePosition = "Right" vor_500res@lbTitleDirection = "Across" vor_500res@cnFillColors = (/0,165,175,185,195,200,205,210,215,220,225,230,235,240,245,254,255/) ;use in NHem ;;; correction for vertical labelbar vor_500res@lbOrientation = "Vertical" vor_500res@pmLabelBarHeightF = 0.25 vor_500res@pmLabelBarWidthF = 0.03 vor_500res@pmLabelBarOrthogonalPosF = 0.0 vor_500res@pmLabelBarParallelPosF = 0.75 ;;;;; ;;;;; g_500res@cnFillOn = False g_500res@cnInfoLabelOn = False g_500res@cnLineLabelsOn = True g_500res@cnLineLabelFontHeightF = .006 g_500res@cnLineLabelPlacementMode = "Constant" g_500res@cnLinesOn = True g_500res@lbLabelBarOn = True g_500res@cnLevelSelectionMode = "ManualLevels" g_500res@cnMinLevelValF = 480 g_500res@cnMaxLevelValF = 588 g_500res@cnLevelSpacingF = 6.0 g_500res@cnLineThicknessF = 3.0 g_500res@cnLineColor = "black" ;;;; ;;;; w_500res@cnFillOn = False w_500res@cnInfoLabelOn = False w_500res@cnLineLabelsOn = False w_500res@cnLinesOn = True w_500res@lbLabelBarOn = False w_500res@cnLevelSelectionMode = "ManualLevels" w_500res@cnMinLevelValF = -10 w_500res@cnMaxLevelValF = -.75 w_500res@cnLevelSpacingF = .25 w_500res@cnLineThicknessF = 2.0 w_500res@cnLineColor = ;;;; ;;;; wind_500res@vcGlyphStyle = "WindBarb" wind_500res@vcWindBarbCalmCircleSizeF = 0 wind_500res@vcWindBarbLineThicknessF = 1.8 wind_500res@vcRefLengthF = 0.018 wind_500res@vcRefMagnitudeF = 10.0 wind_500res@vcRefAnnoOn = False wind_500res@vcRefAnnoString1 = "$VMG$ m s~S~-1" wind_500res@vcRefAnnoString2On = False wind_500res@vcMinDistanceF = 0.019 wind_500res@vcRefAnnoOrthogonalPosF = -1.08 wind_500res@vcWindBarbTickLengthF = 0.37 wind_500res@vcWindBarbTickSpacingF = 0.15 wind_500res@vcWindBarbTickAngleF = 55.0 wind_500res@gsnFrame = False ;***************************************************** ;********************((STEP FOUR))******************** ;***************************************************** ;*******************((CREATE PLOTS))****************** ;***************************************************** ;**********Setting Up Image Creation****************** ; ; After setting the resources you are ready to start ; the image making process. Since we are reading in ; more than one date (48 hours worth of data at ; 6 hour intervals) we need to set up a loop so that ; each image is created for that specific time and does ; not overlap with other images ; ;****************************************************** ; That "time" function we created earlier is going to ; use here. We use a do loop to start at 0 and end at ; the number of total times you read in minus one ; (since we started at zero). do j=0,dimsizes(times)-1 ;**********Work Sheet********************************** ; ; Here we just define a special set of ; resources that will define the name of our ; output files and let us choose an appropriate size ; for our output file image. Numbers are in pixels. ; ; For more information go here: ; ; http://www.ncl.ucar.edu/Document/Graphics/Resources/wk.shtml ; ;****************************************************** wks_type = "png" wks_type@wkWidth = 1500 wks_type@wkHeight = 1500 ; Lets add this incrementing number "j" to our ; image so each image has a unique plot name wks_name = plotName+"_"+j ; The gsn_open_wks command finally lets us open up our image ; for the first time defining the image type (png) and the ; name we have given our file (wks_name) wks = gsn_open_wks(wks_type,wks_name) ;******************************************************* ; ; Define The Colormap ; ; we now we made a fill color resource option ; earlier, so we also know we want to define ; a colormap. If not selected, ncl will choose ; the default colormap, which is super ugly ; (before V6.1 beta). Thus, lets define a colormap ; that will look good for this mean sea level pressure ; data. I've picked a halfway decent one. ; ; For more color table options check out this link: ; ; http://ncl.ucar.edu/Document/Graphics/color_table_gallery.shtml ; ; NOTE: You can also CREATE your own color tables. Not ; this is actually a relatively simple process, but ; does require some patience to learn the process. ; Fortunately, the ncl website also has a nice guide ; to help you through this process which is located ; here: ; ; http://ncl.ucar.edu/Document/Graphics/create_color_table.shtml ; ; ;******************************************************** gsn_define_colormap(wks,"BlueWhiteOrangeRed") ;;*****************************************************;; ;; ;; ;; Time to finish off the plotting resources ;; ;; ;; ;;*****************************************************;; ; Remember the variable you cleared its long name ; and units before. Lets fill those spaces up with useful ; information about the image. I like to use a date and a ; title to help identify the image. ; ; There are MANY other ways to create titles and identify ; data and units, but this way is somewhat intuitive to use ; and clean, since you can easily adjust the size and values ; of these strings ; Left String is Time (the ut_string converts the time back to ; units we can understand) Read the ncl documentation for more ; information on this function res@gsnLeftString = ut_string(times(j),"%H%M UTC %d %c %Y") res@gsnLeftStringFontHeightF = .0075 ; size of text res@gsnLeftStringOrthogonalPosF = -.010 ; position of text in y direction ; The right string is just regular text to define what the variable plotted is res@gsnRightString = "CFSR 500 hPa Heights, Cyclonic Vorticity, Vertical Velocity, Winds" res@gsnRightStringFontHeightF = .0075 ; size of text res@gsnRightStringOrthogonalPosF = -.010 ; position of text in y direction ; Its always nice to have a print statement that tells you want date is being ; created currently print("Now Creating Image For: "+ut_string(times(j),"%H%M UTC %d %c %Y")) ;**********PUTTING THE PLOT TOGETHER****************** ; ; You are finally ready to create the images ; ; In this step, you put all the resource pieces together ; to create a image. There are several pieces depending ; on how many variables you are plotting. In this ; sample script, we have only the background of the plot ; (res) and the one variable slp (pmslres) ; ;****************************************************** ;*******************; ; ; ; PLOT SECTION ; ; ; ;*******************; ; Plot the initial map background plot_map = gsn_csm_map(wks,res) ; Plot the mslp variable as a contour ; ; remember its three dimensions, and the j ; do loop counter is incrementing in the proper ; time units. Remember the resource we are using for ; this variable is pmslres ; ; the functions used here are special graphics routines ; that create the images. For more information go here: ; ; http://www.ncl.ucar.edu/Document/Functions/graphics_routines.shtml ; plot_vor = gsn_csm_contour(wks,vor_500(j,:,:),vor_500res) plot_g = gsn_csm_contour(wks,g_500(j,:,:),g_500res) plot_w = plot_wind = gsn_csm_vector(wks,u_500(j,:,:),v_500(j,:,:),wind_500res) ; Now that you have plotted each element, make sure ; that they overlay on top of one another overlay(plot_map, plot_vor) overlay(plot_map, plot_g) overlay(plot_map, plot_wind) ; if you have more variables you would also ; overlay them over the initial map (plot_map) ; maximizing output makes sure image is 100% size maximize_output(wks,True) ; Now that the image is created and variable is ; overlaid, lets export it to the proper directory ; These system functions are REALLY useful. You can pretty ; much use any linux command within them. Note that the convert ; -trim and -density functions are just imagemagick commands ; that can be found here: ; ; http://www.imagemagick.org/script/index.php ; ; the other commands are regular linux commands like "mv" system("convert -trim -density 100x100 "+wks_name+".png "+wks_name+".png") system("mv "+wks_name+".png "+imgoutdir) print("Figure "+j+" is finished!") end do ; end of the do loop creating the images ; Now that all the images have been created and sent off the proper directory ; this script's job is done! All thats left to do is do an "end" command to ; end the program. print( "Figures are finished") system( "date" ) end