Hi 'NCL-talkers'
I want to plot the box and whiskers for average concentration (aerosols) vs pressure level. The average concentration is interpolated per pressure level (1000 - 100 hPa) from the WRF-Chem model output. I have done the statistical calculation for the median, max, min, IQR's 5th, etc. for each pressure level. How can I plot this data (x-y plot) for the calculated static against vertical pressure at every time interval of my model output? [plotting Pressure along the y-axis and Concentration along the x-axis]
Below is the data I have calculated:
-00001 00048 00054 00060 00066 00072
1000.0 0.07685 0.07876 0.07662 0.07625 0.07781
850.0 0.11857 0.12207 0.12101 0.11989 0.11870
700.0 0.13153 0.13172 0.13315 0.13145 0.12998
500.0 0.08571 0.08544 0.08443 0.08408 0.08450
400.0 0.08540 0.08547 0.08515 0.08509 0.08513
300.0 0.08506 0.08481 0.08445 0.08427 0.08445
250.0 0.08177 0.08155 0.08131 0.08152 0.08197
200.0 0.07668 0.07684 0.07683 0.07787 0.07859
150.0 0.07175 0.07265 0.07297 0.07368 0.07441
100.0 0.05387 0.05402 0.05366 0.05384 0.05415
The left (first) column is the pressure level, while the 2nd -6th columns are the average concentrations. The top row (except the first element '-00001') are the output time intervals. Data is output every 6 hrs (herein 48, 54, 60 & 72 hrs from the initialisation time. part of the script I am trying to use is:
----------------------------------------------------------------------------------------------------------------------
;list the species/variables to output
specieslist = (/"co"/)
nspecies = dimsizes(specieslist)
stat_list = (/"avg","stddev","min","max","median","x25","x75","x5","x95"/)
nstat = dimsizes(stat_list)
nspec = dimsizes(specieslist)
fName_list = new((/nspec, nstat/), "string")
pressure_levels = (/1000., 850., 700., 500., 400., 300., 250., 200., 150., 100./)
nlevels = dimsizes(pressure_levels)
dt = 6 ; calculate stats every 6 files
do ifil = 0, numFILES-1, dt
; do level = 0, nlevels-1
; pressure = pressure_levels(level)
do spec = 0, nspecies-1
species = specieslist(spec)
print("starting the process")
level = new(11,float)
Avg = new(11,float)
;define to filepath to where the data is
filepath_avg = datadir+"/"+species+ "/"+ "avg" +"_"+species+".txt"
filepath_stddev = datadir+"/"+species+ "/"+ "stddev" +"_"+species+".txt"
filepath_min = datadir+"/"+species+ "/"+ "min" +"_"+species+".txt"
filepath_x5 = datadir+"/"+species+ "/"+ "x5" +"_"+species+".txt"
filepath_x25 = datadir+"/"+species+ "/"+ "x25" +"_"+species+".txt"
filepath_median = datadir+"/"+species+ "/"+ "median" +"_"+species+".txt"
filepath_x75 = datadir+"/"+species+ "/"+ "x75" +"_"+species+".txt"
filepath_x95 = datadir+"/"+species+ "/"+ "x95" +"_"+species+".txt"
filepath_max = datadir+"/"+species+ "/"+ "max" +"_"+species+".txt"
;all_avg_data = asciiread("avg_co.txt", -1, "float")
all_avg_data = asciiread(filepath_avg, -1, "float")
all_stddev_data = asciiread(filepath_stddev , -1, "float")
all_min_data = asciiread(filepath_min, -1, "float")
all_x5_data = asciiread(filepath_x5, -1, "float")
all_x25_data = asciiread(filepath_x25, -1, "float")
all_median_data = asciiread(filepath_median, -1, "float")
all_x75_data = asciiread(filepath_x75, -1, "float")
all_x95_data = asciiread(filepath_x95, -1, "float")
all_max_data = asciiread(filepath_max, -1, "float")
;find out how many elements there are
data_size = dimsizes(all_avg_data)
;Assume there are two columns in each file:
columns_tot =6
rows_tot = data_size/columns_tot
all_avg_data_ =6D onedtond(all_avg_data, (/rows_tot,columns_tot/) )
all_stddev_data_6D = onedtond(all_stddev_data, (/rows_tot,columns_tot/) )
all_min_data_6D = onedtond(all_min_data, (/rows_tot,columns_tot/) )
all_x5_data_6D = onedtond(all_x5_data, (/rows_tot,columns_tot/) )
all_x25_data_6D = onedtond(all_x25_data, (/rows_tot,columns_tot/) )
all_median_data_6D = onedtond(all_median_data, (/rows_tot,columns_tot/))
all_x75_data_6D = onedtond(all_x75_data, (/rows_tot,columns_tot/) )
all_x95_data_6D = onedtond(all_x95_data, (/rows_tot,columns_tot/) )
all_max_data_6D = onedtond(all_max_data, (/rows_tot,columns_tot/) )
;level is the left column:
level = all_avg_data_6D(:, 1)
;Avg is the right:
Avg = all_avg_data_6D(:, 1)
Stddev = all_stddev_data_6D(:, 1)
Min = all_min_data_6D(:, 1)
x5 = all_x5_data_6D(:, 1)
x25 = all_x25_data_6D(:, 1)
Median = all_median_data_6D(:, 1)
x75 = all_x75_data_6D(:, 1)
x95 = all_x95_data_6D(:, 1)
Max = all_max_data_6D(:, 1)
conv = Median*1000
;------------------------------------------------------------------------------------------------------
xval = new((/11,5/),"float",-.001)
xval(:,0) = x5 ;bottom value
xval(:,1) = x25 ;bottom value of box
xval(:,2) = Median ;mid-value of box
xval(:,3) = x75 ;top-value of box
xval(:,4) = x95 ;top value
print(xval)
;xval ;denotes the X-axis values where the boxes are drawn
y = level
; end do
end do
;--------------------------------------------------------------------------------------------------
; set the plotting options for the graphs (plot background)
res = True ; plot mods desired
llres = True ; to draw boxes
opti = True ; to control color and width of boxes
tres = True ; text modes desired
mpres = True ; desired marker modes
res@tiMainString = "Time series for: " + species ; plot title
res@tiYAxisString = "level" ; y axis lable
res@tiXAxisString = "Contn (ppb)" ; time (x) axis
; res@cnLevelSpacingF = 2. ; contour interval
; res@tmYLMode = "Explicit"
res@tmYLLabels = (/"1000","850","700","500","400","300","250","200","150","100"/)
res@tmYLLabels(:) = " "
;res@tiXAxisFontHeightF = 0.02
;res@tiYAxisFontHeightF = 0.02
;res@tmXMinF = 0.0 ; min'm
;res@tmXMaxF = 2.400 ; max
;res@tmXValues = "Explicit"
res@tmXUseBottom = True
res@tmXTLabelsOn = False
res@tmXTOn = True
res@tmYROn = True
;res@tmXBTickStartF = Time(0)
;res@tmXBTickEndF = 100
;res@tmXBTickSpacingF = 6.0
llres@gsLineThicknessF = 2.5 ; box line thickness
res@tmXBLabelStride = 2
mpres@gsMarkerIndex = 3
mpres@gsMarkerSizeF = 20.
mpres@gsMarkerColor = "red"
tres@txFontHeightF = 0.024 ; change font height
opti@boxWidth = 1.8 ; Width of box (x units)
;------------------------------------------------------------------------
Firsttime = True
;make the plots
opts = res
opts@cnLineColor = "White"
opts@cnFillOn = True
opts@cnFillMode = "AreaFill"
if (species .eq. "co")
opts@cnLevelSelectionMode = "ExplicitLevels"
opti@boxColors = (/"blue"/) ; Color of box(es)
xval = (xval*1000)
xi = Avg*1000
print("ploting data for:"+ species)
delete(opts)
end if
plot = boxplot(wks,y,xval,opti,res,llres) ; All 3 options used...
; dum1 = gsn_add_polymarker(wks, y, x, xi,mpres)
draw(wks) ; box plot does not call
frame(wks) ; these for you
print("plotting the profiles")
end do
end
;----------------------------------------------------------------------------------------------------------------------
Any help would be appreciated please,
Thank you,
M. Wiston, Postgraduate
Centre for Atmospheric Science (School of Earth, Atmospheric and Environmental Science),
University of Manchester,
M13 9PL,
UK
email: modise.wiston@postgrad.manchester.ac.uk
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Jun 02 03:48:33 2014
This archive was generated by hypermail 2.1.8 : Sat Jun 07 2014 - 11:03:12 MDT