; read GEFS data from multiple directories ; requires (or cancel) a file that has a list of times ; contact barry.h.lynn@gmail.com if you have comments, improvements to suggest ; thank you to ncl-talk participants and the NCL on-line help ; thank you to a couple of folk who posted shell scripts for reading NCEP files ; for example: "www.meteo-blog.net" 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/contrib/cd_string.ncl" external GET_IJ "./get_ij.so" ;*************************************** begin plot = new(2,graphic) ; create a plot array ; define number of time variables n_t_vars=5 n_l_times = 0 ; ascii file a_asc = (/"gefs_time.file"/) ;---Read the data in as a small1D array values = asciiread(a_asc,-1,"float") ; set how many time steps n_l_times = 39 nvars = 5 nxy = 1 nblk = nxy * nvars ; temporary ny=1 nx=1 ; create variables to hold the ascii file variables. ; n_t_vars=5 time_array = new((/n_t_vars/),integer) year = new((/n_l_times,ny,nx/),integer) month = new((/n_l_times,ny,nx/),integer) day = new((/n_l_times,ny,nx/),integer) hour = new((/n_l_times,ny,nx/),integer) minute = new((/n_l_times,ny,nx/),integer) do it=0,n_l_times-1 nstart = it*(nblk+1)+1 nend = nstart+nblk-1 year(it,:,:) = toint(reshape(values(nstart:nend:nvars),(/ny,nx/))) month(it,:,:) = toint(reshape(values(nstart+1:nend:nvars),(/ny,nx/))) day(it,:,:) = toint(reshape(values(nstart+2:nend:nvars),(/ny,nx/))) hour(it,:,:) = toint(reshape(values(nstart+3:nend:nvars),(/ny,nx/))) minute(it,:,:) = toint(reshape(values(nstart+4:nend:nvars),(/ny,nx/))) end do ; print("year = " + year(:,0,0)) ; print("day = " + day(:,0,0)) DATADir = "/home/cust100021_vol1/barry/cust100021_vol2/GEFS" ; find directories with this date, etc cmd = "find " + DATADir + "/GEFS* -type d -print" dirWRF = systemfunc (cmd) print(dirWRF) ; this is the total number of directories within the listed directory header ; it can be just the number of, for example, "GEFS_XX" directories or also all subdirectories within n_dirWRF = dimsizes(dirWRF) print("n_dirWRF" + n_dirWRF) ; set name of output file outfile="gefs_700mb_maps" wks = gsn_open_wks("pdf",outfile) ; open wk station ; http://www.ncl.ucar.edu/Document/Functions/Built-in/systemfunc.shtml ; to retrieve a list of files from another directory diri = dirWRF(0) + "/" all_files = systemfunc("ls " + diri + "ge*pgrb2b*grb") n_files = dimsizes(all_files) print("n_files = " + n_files) ;;; Now that we have identified the number of files and number of directories ; loop over number of files within each directory do n=0,n_files -1,1 ;do n=0,1,1 print("n times steps = " + n) ; loop over directories do i_dir = 0,n_dirWRF-1 ; do i_dir = 0,1,1 print("i_dir = " + i_dir) diri = dirWRF(i_dir) + "/" ; define individual file read all_files = systemfunc("ls " + diri + "ge*pgrb2b*grb") print("all_files(n) = " + all_files(n)) filename = all_files(n) print ("filename = " + filename) a = addfile(filename,"r") fin = a ; print(fin) ;************************************************* ; template panel_13.ncl ; http://www.ncl.ucar.edu/Applications/Scripts/panel_13.ncl ; ; Concepts illustrated: ; - Using "overlay" to overlay contours and vectors on separate maps ; - Paneling two plots vertically on a page ; ;************************************************ lv0 = fin->lv_ISBL0 print("lv0 = " + lv0) lv8 = fin->lv_ISBL8 print("lv8 = " + lv8) lv10 = fin->lv_ISBL10 print("lv10 = " + lv10) t = fin->TMP_P1_L100_GLL0(17,:,:) rh = fin->RH_P1_L100_GLL0(17,:,:) dims2d = dimsizes(t) print("dims2d(0) = " + dims2d(0)) print("dims2d(1) = " + dims2d(1)) z_up = fin->HGT_P1_L100_GLL0(16,:,:) z_dn = fin->HGT_P1_L100_GLL0(17,:,:) z =0.5*(z_up+z_dn) print("z ave = "+ z(20,20)) copy_VarCoords(z_up,z) ; copy coord vars to speed u_up = fin->UGRD_P1_L100_GLL0(13,:,:) u_dn = fin->UGRD_P1_L100_GLL0(14,:,:) v_up = fin->VGRD_P1_L100_GLL0(13,:,:) v_dn = fin->VGRD_P1_L100_GLL0(14,:,:) u =0.5*(u_up+u_dn) v =0.5*(v_up+v_dn) copy_VarCoords(u_up,u) ; copy coord vars to speed copy_VarCoords(v_up,v) ; copy coord vars to speed lon = fin->lon_0(:) lat = fin->lat_0(:) dims2d = dimsizes(lat) print("dim2d lat = " + dims2d(0)) print("dim2d lat = " + dims2d(1)) dims2d = dimsizes(lon) print("dim2d lon = " + dims2d(0)) print("dim2d lon = " + dims2d(1)) if (i_dir.eq.0)then t_ave = t/n_dirWRF rh_ave = rh/n_dirWRF z_ave = z/n_dirWRF u_ave = u/n_dirWRF v_ave = v/n_dirWRF end if if (i_dir.ne.0)then t_ave = t/(n_dirWRF) + t_ave rh_ave = rh/n_dirWRF + rh_ave z_ave = z/n_dirWRF + z_ave u_ave = u/n_dirWRF + u_ave v_ave = v/n_dirWRF + v_ave end if ; make plot if (i_dir.eq.n_dirWRF-1)then t_ave=t_ave-273.16 ; convert to C z_ave=z_ave/10. ; convert to dam printMinMax(z_ave,False) u = u * 3.6 v = v * 3.6 ; speed = sqrt(u_ave^2+v_ave^2) ; km/h from m/s ; speed = speed * 3.6 copy_VarCoords(u,u_ave) ; copy coord vars to speed copy_VarCoords(v,v_ave) ; copy coord vars to speed copy_VarCoords(z,z_ave) ; copy coord vars to speed copy_VarCoords(t,t_ave) ; copy coord vars to speed copy_VarCoords(rh,rh_ave) ; copy coord vars to speed ;************************************************ ; create plots ;************************************************ ;************************************************ ; wks = gsn_open_wks("png","panel") ; send graphics to PNG file ;************************************************ ; create plots gsn_define_colormap(wks,"BkBlAqGrYeOrReViWh200") ; map perameters res1 = True res1@mpDataBaseVersion = "Ncarg4_1" ; choose more recent database res1@mpDataSetName = "Earth..4" ; high resolution res1@mpOutlineBoundarySets = "National" ; National borders res1@mpNationalLineColor = "Black" ; National borders color res1@mpNationalLineThicknessF = 3.0 res1@mpGeophysicalLineThicknessF = 2.0 ; time parameters ; first calculate time new_time_units = "hours since 1800-01-01 00:00" new_time_units = new_time_units new_time_array = cd_inv_calendar( year(n,0,0), month(n,0,0), day(n,0,0), hour(n,0,0), minute(n,0,0), 0, new_time_units,0) print("new_time_array = " + new_time_array) format = "" ;; defaults to "%H%M IST %d %c %Y" format = "%H%M IST %d %c %Y" new_time_array = new_time_array + 2.0 time_var=cd_string(new_time_array,format) res1@tiMainString = time_var ;************************************************ ; temperature res1@gsnDraw = False ; don't draw res1@gsnFrame = False ; don't advance frame res1@cnFillOn = True ; turn on color ; res1@cnFillPalette = "gui_default" ; set color map res1@cnLinesOn = False ; no contour lines ; res1@gsnLeftString = "700 mb T" ; change left string res1@gsnLeftString = "700 mb Wind Vectors" ; change left string ; res1@gsnRightString = u@units ; assign right string res1@gsnRightString = "km/h" ; assign right string res1@mpFillOn = False ; no map fill res1@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels ; res1@cnLevels = (/ -26.,-24.,-22.,-20.,-18.,-16.,-14.,\ ; -12.,-10.,-8.,-6.,-4.,-2.,0.,\ ; 2.,4.,6.,8.,10.,12.,14.,16.,18.,\ ; 20.,22.,24.,26.,28.,30.,32.,34.,\ ; 36.,38.,40./) ; set levels res1@cnLevels = (/ -34.,-32.,-30.,-28.,-26.,\ -24.,-22.,-20.,-18.,-16.,-14.,\ -12.,-10.,-8.,-6.,-4.,-2.,0.,\ 2.,4.,6.,8.,10.,12.,14.,16.,18.,\ 20.,22.,24.,26.,28.,30.,32./) res1@cnFillColors = (/ 10,12,14,15,17,19,22,25,31,35,39,43,48,53,\ 61,64,67,71,89,95,99,104,108,111,115,123,\ 130,135,143,148,150,155,159,163,165/) ; set the colors to be used ; res1 label bar res1@cnInfoLabelOn = False ; turn off contour info label res1@lbAutoManage = False ; control label bar res1@pmLabelBarDisplayMode = "Always" ; turns on label bar res1@lbOrientation = "Vertical" res1@pmLabelBarSide = "Right" res1@lbLabelAutoStride = True res1@pmLabelBarWidthF = 0.20 res1@pmLabelBarHeightF = 0.7 res1@lbLabelFontHeightF = .025 res1@lbPerimOn = False res1@lbLabelFont = "Helvetica" ; label font res1@lbTitleOn = True ; turn on title res1@lbTitleString = "Temp [~S~o~N~C]" ; title string res1@lbTitlePosition = "Top" ; title position res1@lbTitleFontHeightF= .020 ; make title smaller res1@lbTitleDirection = "Across" ; title direction ; doesn't seem to help res1@lbTitleOffsetF = 0.002 ; label bar title position ; res1@lbTopMarginF = 0.001 ; rh res2 = True ; map parameters res2 = True res2@mpDataBaseVersion = "Ncarg4_1" ; choose more recent database res2@mpDataSetName = "Earth..4" ; high resolution res2@mpOutlineBoundarySets = "National" ; National borders res2@mpNationalLineColor = "Black" ; National borders color res2@mpNationalLineThicknessF = 3.0 res2@mpGeophysicalLineThicknessF = 2.0 ; time string res2@tiMainString = time_var ; Shading res2@gsnDraw = False ; don't draw res2@gsnFrame = False ; don't advance frame res2@cnFillOn = True ; turn on color res2@cnFillPalette = "gui_default" ; set color map res2@cnLinesOn = False ; no contour lines ; res2@gsnLeftString = "700 mb RH" ; change left string res2@gsnLeftString = "" ; change left string ; res2@gsnRightString = rh@units ; assign right string res2@gsnRightString = "" ; assign right string res2@mpFillOn = False ; no map fill res2@lbTitleString = "RH [%]" ; title string ; res2 label bar res2@cnInfoLabelOn = False ; turn off contour info label res2@lbAutoManage = False ; control label bar res2@pmLabelBarDisplayMode = "Always" ; turns on label bar res2@lbOrientation = "Vertical" res2@pmLabelBarSide = "Right" res2@lbLabelAutoStride = True res2@pmLabelBarWidthF = 0.20 res2@pmLabelBarHeightF = 0.7 res2@lbLabelFontHeightF = .025 res2@lbPerimOn = False res2@lbLabelFont = "Helvetica" ; label font res2@lbTitleOn = True ; turn on title res2@lbTitleString = "RH [%]" ; title string res2@lbTitlePosition = "Top" ; title position res2@lbTitleFontHeightF= .020 ; make title smaller res2@lbTitleDirection = "Across" ; title direction ; doesn't seem to help res2@lbTitleOffsetF = 0.002 ; label bar title position ; res2@lbTopMarginF = 0.001 ; z res3 = True res3@gsnDraw = False ; do not draw the plot res3@gsnFrame = False ; do not advance the frame res3@gsnLeftString = "700 mb Geopotential Height" ; change left string ; res3@gsnRightString = z@units ; assign right string res3@gsnRightString = "dam" ; assign right string res3@cnConstFLabelPerimOn = False res3@cnLevelSelectionMode = "ManualLevels" res3@cnMinLevelValF = 250 res3@cnMaxLevelValF = 350 res3@cnLevelSpacingF = 5 res3@cnLineLabelAngleF = 0.0 res3@cnLineLabelFontHeightF = 0.01 res3@cnLineLabelBackgroundColor = -1 res3@gsnContourLineThicknessesScale = 1.5 res3@tfDoNDCOverlay = True vecres = True ; vector only resources vecres@gsnDraw = False ; don't draw vecres@gsnFrame = False ; don't advance frame vecres@vcGlyphStyle = "CurlyVector" ; curly vectors vecres@vcRefMagnitudeF = 20 ; define vector ref mag vecres@vcRefLengthF = 0.045 ; define length of vec ref vecres@gsnRightString = " " ; turn off right string vecres@gsnLeftString = " " ; turn off left string vecres@tiXAxisString = " " ; turn off axis label vecres@vcRefAnnoOrthogonalPosF = 0.01 ; move ref vector into plot ;; Map Projection (min and max Lon and Lat) lat_beg = 20 lat_end = 50 lon_beg = 20 lon_end = 50 res1@gsnAddCyclic = True ; data already has cyclic point ; this must also be set for any zoom ; res1@mpProjection = "LambertConformal" ; The GEFS is not! res1@mpMinLatF = lat_beg res1@mpMaxLatF = lat_end res1@mpMinLonF = lon_beg res1@mpMaxLonF = lon_end res1@mpLimitMode="LatLon" res2@mpMinLatF = lat_beg res2@mpMaxLatF = lat_end res2@mpMinLonF = lon_beg res2@mpMaxLonF = lon_end res2@mpLimitMode="LatLon" plotA = gsn_csm_contour_map(wks,t_ave,res1) ; must manuall excerpt variable dims2d = dimsizes(z) nj = dims2d(0) ; must manuall excerpt variable dims2d = dimsizes(z) nj = dims2d(0) ni = dims2d(1) lat2d = z lon2d = z ; print("lat = " + lat) ; print("lon = " + lon) dims2d = dimsizes(lat2d) print("diml2d = " + dims2d(0)) print("diml2d = " + dims2d(1)) dims2d = dimsizes(lon2d) print("diml2d = " + dims2d(0)) print("diml2d = " + dims2d(1)) copy_VarCoords(z,lat2d) copy_VarCoords(z,lon2d) do j = 0,nj-1 do i = 0,ni-1 lat2d(j,i)=lat(j) lon2d(j,i)=lon(i) end do end do print("dim2d = " + dims2d(0)) print("dim2d = " + dims2d(1)) dims2d = dimsizes(lat) print("dim2d = " + dims2d(0)) print("dim2d = " + dims2d(1)) i_loc = 1 j_loc = 1 GET_IJ::get_ij(lat2d,lon2d,lat_beg,lon_beg,i_loc,j_loc,nj,ni) x_start = i_loc y_start = j_loc print("x_start = " + x_start) print("y_start = " + y_start) GET_IJ::get_ij(lat2d,lon2d,lat_end,lon_end,i_loc,j_loc,nj,ni) x_end = i_loc y_end = j_loc print("x_end = " + x_end) print("y_end = " + y_end) u_ave_new = u_ave(y_start:y_end,x_start:x_end) v_ave_new = v_ave(y_start:y_end,x_start:x_end) lat2d_new = lat2d(y_start:y_end,x_start:x_end) lon2d_new = lon2d(y_start:y_end,x_start:x_end) ; plotB = gsn_csm_vector(wks,u_ave_new,v_ave_new,vecres) ;;;opy_VarCoords(dumb_ave,u_ave_new) ; copy_VarCoords(dumb_ave,v_ave_new) u_ave_new!0 = "lat" u_ave_new!1 = "lon" u_ave_new&lat= lat(y_start:y_end) u_ave_new&lon= lon(x_start:x_end) v_ave_new!0 = "lat" v_ave_new!1 = "lon" v_ave_new&lat= lat(y_start:y_end) v_ave_new&lon= lon(x_start:x_end) vecres@gsnAddCyclic = False ; we will manually zoom in. plotB = gsn_csm_vector(wks,u_ave_new(::4,::4),v_ave_new(::4,::4),vecres) overlay(plotA,plotB) ; result will be plotA plot(0) = plotA ; now assign plotA to array plotC = gsn_csm_contour_map(wks,rh_ave,res2) res3@gsnAddCyclic = False ; we will manually zoom in. ; z_ave_new = z z_ave_new = z_ave(y_start:y_end,x_start:x_end) printMinMax(z_ave_new,False) z_ave_new!0 = "lat" z_ave_new!1 = "lon" z_ave_new&lat= lat(y_start:y_end) z_ave_new&lon= lon(x_start:x_end) plotD = gsn_csm_contour(wks,z_ave_new,res3) overlay(plotC,plotD) ; result is plotC plot(1) = plotC ; now assign plotC to array resP = True ; panel only resources resP@gsnMaximize = True ; maximize plots gsn_panel(wks,plot,(/2,1/),resP) ; now draw as one plot end if end do end do end