;******************************************************** ; Author : Ahmed lasheen * ; Email : ahmed4kernel@gmail.com * ; Date : 16 jan 2011 * ; revised: 12 jan 2012 ; 23 jan 2012 ; purpose: printing the following image at all time step* ; 1-MSL Pressure * ; 2-Temperature at surface * ; 3-wind at 1000 mb * ;NCL Version : 5.2.1,6.0.0 * ;******************************************************** 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" ;; to use time function load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" ;;to draw the ndc coordinate load "viewport_3.ncl" ;; this helps me to adjust the location of the label bar. begin ;************************************************************* package_path = getenv("HOME")+"/vis_ncl_package/" file_input=package_path+"MMOUTP_DOMAIN1.nc" MODEL="mm5" DIRECTORY="default" TITLE="yes" ;################################################################### ;##### copying the colormap needs for the drawing ################## ;################################################################### system("cp -f ahmed_red_blue.rgb $NCARG_ROOT/lib/ncarg/colormaps/.") ;################################################################### ;##### creating the folders ########################## ;################################################################### data_path = package_path+"/" ;################################################################## wrf = addfile("MMOUTP_DOMAIN1.nc","r") ;opening the grib file ;################################################################# ;################reading the data ################################ ;################################################################# u = wrf->u10(:,:,:) ;reading the U wind component at 1000 mb. v = wrf->v10(:,:,:) ;reading the V wind component at 1000 mb. tmp = wrf->t2(:,:,:) ;reading the temperature at surface . mslp = wrf->psealvlc(:,:,:) lat=wrf->latitdot lon=wrf->longidot for_time=wrf->time ini_time= "minutes since 1986-04-25 12:00:00" ;;############################################################ ;;;some reportation on the variables tmp =tmp -273.15 ;converting kelvin to celsius mslp =mslp*0.01 ;converting the pressure from Pa to hpa wind=sqrt(u(:,:,:)^2+v(:,:,:)^2) max_wind = max(wind) min_wind = min(wind) max_temp = max(tmp) min_temp = min(tmp) ;##################################################################### month_abbr = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"/) days_abbr = (/ "SUN","MON","TUE","WED","THR","FRI","SAT" /) ;;;;;;;;;;;;;;;;;Dealing with the intial time*** stc_ini_time=stringtochar(ini_time);string to character(stc) print(stc_ini_time) ;print if you want ;;for further information about coverting form string to integers see the following ;http://www.ncl.ucar.edu/Support/talk_archives/2005/0483.html iyear=stringtointeger((/stc_ini_time(14:17)/)) imonth=stringtointeger((/stc_ini_time(19:20)/)) iday=stringtointeger((/stc_ini_time(22:23)/)) ;iminute=stringtointeger((/stc_ini_time(12:13)/)) iminute=stringtointeger((/stc_ini_time(28:29)/)) ihour=stringtointeger((/stc_ini_time(25:26)/)) ;ini_time_print =sprinti(" %0.2i",iminute)+"Z "+days_abbr(day_of_week(iyear,imonth,iday))+sprinti(" %0.2i ",iday)+month_abbr(imonth)+sprinti(" %0.4i ",iyear) ini_time_print =sprinti(" %0.2i",ihour)+"Z "+days_abbr(day_of_week(iyear,imonth,iday))+sprinti(" %0.2i ",iday)+month_abbr(imonth)+sprinti(" %0.4i ",iyear) print(ini_time_print) ;print if you want ;;;;;;;;;;;;;;;;;;Dealing with the forcast time*** ;this is the calender start time to find the complete dates of the forcast time based on the ;the forcast hours. ;in other word , given the start date and the hour increaments , find the date coorsponds ;to these hours from the calender. ;for_time@units = "hours since "+iyear+"-"+imonth+"-"+iday+" "+iminute for_time@units = "hours since "+iyear+"-"+imonth+"-"+iday+" "+ihour+":"+iminute+":0.0" utc_date = ut_calendar(for_time,0) ;print(utc_date) ;print if you want fyear =floattointeger(utc_date(:,0)) fmonth =floattointeger(utc_date(:,1)) fday =floattointeger(utc_date(:,2)) fhour =floattointeger(utc_date(:,3)) fminute=floattointeger(utc_date(:,4)) fsecond=floattointeger(utc_date(:,5)) for_time_print =sprinti(" %0.2i",fhour)+"Z "+days_abbr(day_of_week(fyear,fmonth,fday))+sprinti(" %0.2i ",fday)+month_abbr(fmonth)+sprinti(" %0.4i ",fyear) print(for_time_print) ;print if you want not=dimsizes(for_time) ; number of times (not) , this count starts from 1 ;##################################################################### dir_output = "../output/"+MODEL+"/surface/" system("rm -rfv "+dir_output) system("mkdir -p " +dir_output ) ;surface folder do t=0,not-1 ;do loop over the time. print("working on the " +dir_output+sprinti("%0.3i",t)+".png map at_time " + for_time_print(t) +" Forcast Hours ") wks = gsn_open_wks("png" ,dir_output+sprinti("%0.3i",t)) ; ps,pdf,x11,ncgm,eps gsn_merge_colormaps(wks,"wind_17lev","ahmed_red_blue") ;merging the two colors ;##################################################################### ;####### map resources ############################################### ;##################################################################### resm = True ;use this resources resm@gsnMaximize = True ;maximize the picture resm@tfDoNDCOverlay = True ; do not transform data(native grid) resm@mpProjection = "CylindricalEquidistant" ;projection type resm@mpLimitMode = "Corners" ;limiting modes resm@mpLeftCornerLatF = lat(0,0) resm@mpLeftCornerLonF = lon(0,0) resm@mpRightCornerLatF = lat(dimsizes(lat)-1,0) resm@mpRightCornerLonF = lon(dimsizes(lon)-1,0) resm@pmTickMarkDisplayMode = True resm@mpOutlineOn = True ;enables the map area outline. resm@mpOutlineBoundarySets = "AllBoundaries" ;turn on country boundaries. resm@mpGeophysicalLineColor = "black" ;color of cont. outlines. resm@mpGeophysicalLineThicknessF = 1. ;double the thickness of geophysical boundaries. resm@mpNationalLineThicknessF = 1. ;double the thickness of national boundaries. resm@mpFillOn = False ;donot fill the map with a color. resm@mpDefaultFillColor = True ;##titles font## resm@gsnStringFont ="triplex_italic" ;use this font type in the subtitle ;##Left Title## resm@gsnLeftStringFontHeightF = 0.010 resm@gsnLeftStringOrthogonalPosF = 0.05 ;left string vertical pos resm@gsnLeftString = " " ;left string ;##Right Title## resm@gsnRightStringFontHeightF = 0.009 resm@gsnRightStringOrthogonalPosF = 0.01 resm@gsnRightString = "Forcast Time " +for_time_print(t)+"~C~"+ \ "Intial time "+ini_time_print +"~C~"+ \ "WRF NMM " +for_time(t)+" Hrs forcast" ;##Center Title## resm@gsnCenterStringFontHeightF = 0.013 resm@gsnCenterStringOrthogonalPosF = 0.15 ;right string vertical pos if (TITLE .eq. "yes") then resm@gsnCenterString = " Egyptian Meteorological Authority~C~" \ +"Cairo Numerical Weather prediction Center" else if(TITLE .eq. "no")then resm@gsnCenterString = " " end if end if if(MODEL .eq. "gfs") then resm@mpLimitMode = "latlon" resm@mpMinLonF =min(lon) resm@mpMaxLonF =max(lon) resm@mpMinLatF =min(lat) resm@mpMaxLatF =max(lat) end if ;##Donot draw now## resm@gsnDraw = False ; do not draw the plot resm@gsnFrame = False ; do not advance the frame ;####################################################################### ;######pressure contour resources ###################################### ;####################################################################### resc = True ;use this resources resc@cnFillOn = False ;donot shade resc@cnLinesOn = True ;draw line contour resc@cnLineThicknessF = 3 ;the thickness of the line contour resc@cnLineColor = 1 ;the color of the line ;;;;;infolabel resc@cnInfoLabelOn = False resc@cnInfoLabelFont = "triplex_italic" resc@cnInfoLabelString = "Hight contours from $CMN$ to $CMX$"; by $CIU$" resc@cnInfoLabelFontColor = "black" ;;red font color resc@cnInfoLabelPerimOn = False ;;not to draw abox arount information bar resc@cnInfoLabelParallelPosF = 1.0 ;;position in x axis resc@cnInfoLabelOrthogonalPosF = 0.1 ;;position in y axis resc@cnInfoLabelFontAspectF = 2 resc@cnInfoLabelFontHeightF = 0.01 ;;;;;LineLabel resc@cnLabelMasking = True ;; not to overlay the label and contour line ;resc@cnLineLabelAngleF = 0.0;make the labelline label horizontal resc@cnLineLabelFontColor = "black" resc@cnLineLabelPerimOn = True ;; draw box around the label resc@cnLineLabelFont = "triplex_italic" resc@cnLineLabelFontAspectF = 2 ;##titles font## resc@gsnStringFont ="triplex_italic" ;use this font type in the subtitle ;##Left Title## resc@gsnLeftStringFontHeightF = 0.010 resc@gsnLeftStringOrthogonalPosF = 0.04 ;left string vertical pos resc@gsnLeftString = "MSL Pressure Contours(mb)" ;left string ;##Right Title## resc@gsnRightStringFontHeightF = 0.010 resc@gsnRightStringOrthogonalPosF = 0.01 resc@gsnRightString =" " ;##Center Title## resc@gsnCenterStringFontHeightF = 0.015 resc@gsnCenterStringOrthogonalPosF = 0.2 ;right string vertical pos resc@gsnCenterString = " " ;;;Donot draw know resc@gsnDraw = False ; do not draw the plot resc@gsnFrame = False ; do not advance the frame ;####################################################################### ;######Temperature contour resources ###################################### ;####################################################################### rest = True ;use this resources rest@cnFillOn = True ;draw shaded rest@cnLinesOn = False ;donot draw contour line rest@gsnSpreadColors = True ;use all colors in the colormap rest@gsnSpreadColorStart = 20 rest@gsnSpreadColorEnd = 255 rest@cnLevelSelectionMode = "ManualLevels" rest@cnLevelSpacingF = 2 rest@cnMinLevelValF = min_temp rest@cnMaxLevelValF = max_temp ;##titles font## rest@gsnStringFont ="triplex_italic" ;use this font type in the subtitle ;##Left Title## rest@gsnLeftStringFontHeightF = 0.010 rest@gsnLeftStringOrthogonalPosF = 0.01 ;left string vertical pos rest@gsnLeftString = "temperature shaded(C)" ;left string ;##Right Title## rest@gsnRightStringFontHeightF = 0.010 rest@gsnRightStringOrthogonalPosF = 0.01 rest@gsnRightString =" " ;##label bar## rest@lbLabelsOn = True rest@lbLabelStride = 3 rest@pmLabelBarOrthogonalPosF = -0.2 rest@pmLabelBarParallelPosF = 0.5 rest@lbLabelFont = "triplex_italic" rest@lbLabelFontAspectF = 2 rest@lbLabelFontHeightF = 0.01 rest@lbTopMarginF = 0.7 rest@lbTitleOn = True rest@lbTitleFont = "triplex_italic" rest@lbTitleString = "Temperature in Celsius" rest@lbTitleFontAspectF = 2 rest@lbTitleFontHeightF = 0.01 ;##stop drawing now## rest@gsnDraw = False ; do not draw the plot rest@gsnFrame = False ; do not advance the frame ;####################################################################### ;####################################################################### ;######Vectors contour resources ###################################### ;####################################################################### resv = True ;use this resource ;;vector properties resv@vcGlyphStyle = "CurlyVector" resv@vcMinDistanceF = 0.022 resv@vcRefMagnitudeF = (max_wind+min_wind)/2 ;define vector ref mag resv@vcRefLengthF = 0.04 ; define length of vec ref resv@vcLevelSelectionMode ="ManualLevels" resv@vcMaxLevelValF =max_wind resv@vcMinLevelValF =min_wind ;;reference vector resv@vcRefAnnoOn = True resv@vcRefAnnoSide = "Bottom" resv@vcRefAnnoJust = "BottomLeft" resv@vcRefAnnoParallelPosF = 0.00;parallel postion resv@vcRefAnnoOrthogonalPosF = -0.34 ;vertical postion resv@vcRefAnnoString2 = "" ;the word reference vector ;resv@vcRefAnnoString1 = "" ;the mag of the wind ;##titles font## resv@gsnStringFont ="triplex_italic" ;use this font type in the subtitle ;##Left Title## resv@gsnLeftStringFontHeightF = 0.010 resv@gsnLeftStringOrthogonalPosF = 0.07 ;left string vertical pos resv@gsnLeftString = "wind at 1000mb (m/s)" ;left string ;##Right Title## resv@gsnRightStringFontHeightF = 0.012 resv@gsnRightStringOrthogonalPosF = 0.01 resv@gsnRightString =" " ;;donot draw know resv@gsnDraw = False ; do not draw the plot resv@gsnFrame = False ; do not advance the frame ;############################################################################### ;############################################################################### ;;;;; overlay all plots print("we are hereeeeeeeeeeeeeeeeeeeeeeeeeeeee") plot_m =gsn_csm_map(wks,resm) print("we are hereeeeeeeeeeeeeeeeeeeeeeeeeeeee") plot_c =gsn_csm_contour(wks,mslp(t,:,:),resc) plot_t =gsn_csm_contour(wks,tmp(t,:,:),rest) getvalues plot_t "cnLevels" : levels end getvalues rest@lbLabelStrings = sprintf("%3.0f",levels) ; Format the labels plot_t =gsn_csm_contour(wks,tmp(t,:,:),rest) plot_v =gsn_csm_vector(wks,u(t,:,:),v(t,:,:),resv) overlay(plot_m,plot_t) overlay(plot_m,plot_c) overlay(plot_m,plot_v) ;>>>>>>>>>>>>>>>>>>>>drawing the view port ,bounding box and NDC ;draw_bb_box(wks,plot_m) ;draw_vp_box(wks,plot_m) ;drawNDCGrid(wks) ; Draw helpful grid lines showing NDC square. ;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> draw(plot_m) ;*********************************************************************** ;;;;making the high and lows symbols wmsetp("ezf",1) ;Tell wmap we are using an existing map projection. ;;;High pressure wmsetp("hib",4) ;high background . wmsetp("his",1) ;high shadow Black . wmsetp("hic",5) ;high bouding circle . wmsetp("sht",0.01) ;increase the size. ;;;Low pressure wmsetp("LOB",4) ;low background color. wmsetp("LOC",1) ;low circumscribed. wmsetp("LOF",5) ;low bouding circle . wmsetp("sht",0.01) ;increase the size. ;>>>>>> NCL ways for finding the max and min indicies ;>>>>>>>this methode is computationally good xMax = max(mslp(t,:,:)) xMin = min(mslp(t,:,:)) ; the following function convert the array from two dim to one dim ;nd : means number of dims ;to : to ;one: one ; d : dim x1D = ndtooned(mslp(t,:,:)) ; the indMin has two values which are the coordinate i,j indMin = ind_resolve(ind(x1D.eq.xMin),dimsizes(mslp(t,:,:))) ; locations of min indMax = ind_resolve(ind(x1D.eq.xMax),dimsizes(mslp(t,:,:))) ; locations of max wmlabs(wks,lat(indMin(0,0)),lon(indMin(0,1)), "LO") ; Draw Low symbols. wmlabs(wks,lat(indMax(0,0)),lon(indMax(0,1)), "HI") ; Draw High symbols. ; this if condition just for if we have two high if (dimsizes(ind(x1D.eq.xMax)) .gt. 1) then wmlabs(wks,lat(indMax(1,0)),lon(indMax(1,1)), "HI") ; Draw High symbols. end if ; this if condition just for if we have two low if (dimsizes(ind(x1D.eq.xMin)) .gt. 1) then wmlabs(wks,lat(indMin(1,0)),lon(indMin(1,1)), "LO") ; Draw High symbols. end if ; we must put the two delete as sometimes we have two max or min so that ; the dim of indMin or indMax will change so that the value of the indMin ,indMax ; must be deleted each time so that it can take the new values and dim without feeling connflict ;; this is important. delete(indMin) delete(indMax) ;>>>>>>>>>> My way for finding the indicies ;>>>>>>>>>>> this way is cost computaionally. ;do i=0,dimsizes(lat)-1 ; do j=0,dimsizes(lon)-1 ;if (((.not. (ismissing(mslp(t,i,j))))) .and. ((mslp(t,i,j) .eq. max(mslp(t,:,:)) ))) then ; print("At that time there is max MSL pressure "+mslp(t,i,j)+" mb at ("+lat(i)+","+ lon(j)+")") ; wmlabs(wks,lat(i),lon(j), "HI") ; Draw High symbols. ;end if ;if (((.not. (ismissing(mslp(t,i,j))))) .and. ((mslp(t,i,j) .eq. min(mslp(t,:,:)) ))) then ; print("At that time there is min MSL pressure "+mslp(t,i,j)+" mb at ("+lat(i)+","+ lon(j)+")") ; wmlabs(wks,lat(i),lon(j), "LO") ; Draw Low symbols. ;end if ; end do ;end do ;********************************************************************************* frame(wks) system("convert -trim "+dir_output+sprinti("%0.3i",t)+".png"+" "+dir_output+sprinti("%0.3i",t)+".png") end do end