<div dir="ltr"><div><div><div><div><div><div><div><div><div><div>A bit involved for a simpleton such as myself.<br><br></div>[1] Code clarity:<br> <br>You repeatedly create an array and add assorted meta data. This just adds 'code clutter.' A practical fact is that the longer a code is ... the less likely people are going to invest the time to address the desired issues. <br><br></div>[2] Use a procedure or function to isolate repeated tasks. You do the following task multiple times. <br><br> dPShift_dry = new ((/dimsizes(time),dimsizes(<wbr>lat),dimsizes(lon)/), typeof(dPShift))<br>
dPShift_dry!0="time"<br>
dPShift_dry&time=time<br>
dPShift_dry!1="lat"<br>
dPShift_dry&lat=lat<br>
dPShift_dry!2="lon"<br>
dPShift_dry&lon=lon<br><div><br></div><br><br>[3] I suggest a function that isolates the task.<br><br>undef("createVarMeta")<br>function createVarMeta(time[*]:numeric,lat[*]:numeric, lon[*]:numeric \<br> ,varType[1]:string,long_name[1]:string, units[1]:string)<br><div>local x<br></div>begin<br> x = new ((/dimsizes(time),dimsizes(<wbr>lat),dimsizes(lon)/), varType)<br>
x!0="time"<br>
x&time=time<br>
x!1="lat"<br>
x&lat=lat<br>
x!2="lon"<br>
x&lon=lon<br><br> if (<a href="http://long_name.ne">long_name.ne</a>."") the<br> var@long_name = long_name<br> end if<br> if (<a href="http://units.ne">units.ne</a>."") the<br> var@units =units<br> end if<br><br> return(x)<br>end<br><div><br><div>[4] Use it like:<br><br> dPShift_dry = createVarMeta(time, lat, lon, typeof(...), "dummy" , "whatever")<br></div><div><br><div><div><div> <br></div><div>DO this for all variable you are creating.<br></div><div>The function should be placed at the top of the script <br></div><div><div><div><div><div><br><br></div></div></div></div></div></div></div></div></div><br></div></div></div></div></div></div></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Mar 21, 2017 at 10:12 AM, Melissa Lazenby <span dir="ltr"><<a href="mailto:M.Lazenby@sussex.ac.uk" target="_blank">M.Lazenby@sussex.ac.uk</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div>
<div style="direction:ltr;font-family:Tahoma;color:#000000;font-size:10pt">Hi All
<br>
<br>
Please ignore previous email below. I have managed to figure the area-averaging out but have some concerns regarding the overall code. I am happy that the code does what I would like for dP but not for the other components?<br>
<br>
My overall aim is to: <br>
<br>
1. mask out all ocean regions<br>
2. mask out all the wetting regions<br>
3. area average these changes over southern africa<br>
4. create boxplots of dP and 6 decomposed precipitation components<br>
<br>
I feel the code does this correctly for dP, but not convinced for the other components.<br>
<br>
Any advice would very much be appreciated regarding this issue. Below is the full code to tackle task 1 -4<br>
<br>
Many thanks in advance!<br>
<br>
Kindest Regards<br>
Melissa<br>
<br>
;*****************************<wbr>****************<br>
; boxplot.ncl<br>
;<br>
; Concepts illustrated:<br>
; - Drawing box plots<br>
; - Setting the color of individual boxes in a box plot<br>
; - Setting the width of individual boxes in a box plot<br>
;<span class=""><br>
;*****************************<wbr>****************<br>
load "$NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/gsn_code.ncl"<br>
load "$NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/gsn_csm.ncl"<br>
load "$NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/contributed.<wbr>ncl"<br>
load "$NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/shea_util.ncl"<br>
;*****************************<wbr>****************<br>
begin<br>
;*****************************<wbr>*****************<br>
; Add data in text format and sort in ascending order<br>
;*****************************<wbr>*****************<br>
<br></span>
latS = -30. <br>
latN = 0. <br>
lonL = 10. <br>
lonR = 40.<div><div class="h5"><br>
<br>
model = (/"bcc-csm1-1-m","BNU-ESM","<wbr>CanESM2","CCSM4","CESM1-BGC","<wbr>CESM1-CAM5","CSIRO-Mk3-6-0","<wbr>FIO-ESM","GFDL-CM3","GFDL-<wbr>ESM2G","GFDL-ESM2M","GISS-E2-<wbr>H","HadGEM2-CC","HadGEM2-ES","<wbr>IPSL-CM5A-LR","IPSL-CM5A-MR","<wbr>MIROC5","MRI-CGCM3","NorESM1-<wbr>M","NorESM1-ME"/)<br>
<br>
<br>
do gg = 0,dimsizes(model)-1<br>
<br>
;dP<br>
in1=addfile("/research/geog/<wbr>data2/DATA/mlazenby/Chadwick_<wbr>Analysis/ChadprfutureOND/<wbr>Delta_pr_"+model(gg)+".nc","r"<wbr>)
<br>
<br>
;dPShift<br>
in2=addfile("/research/geog/<wbr>data2/DATA/mlazenby/Chadwick_<wbr>Analysis/<wbr>ChadDeltaPComponentsOND/<wbr>DeltaPSpat_"+model(gg)+"_<wbr>safrica_PSpat.nc","r")<br>
;dPT<br>
in3=addfile("/research/geog/<wbr>data2/DATA/mlazenby/Chadwick_<wbr>Analysis/<wbr>ChadDeltaPComponentsOND/<wbr>DeltaPT_"+model(gg)+"_safrica_<wbr>PT.nc","r")<br>
;dPWeak<br>
in4=addfile("/research/geog/<wbr>data2/DATA/mlazenby/Chadwick_<wbr>Analysis/<wbr>ChadDeltaPComponentsOND/<wbr>DeltaPDiv_"+model(gg)+"_<wbr>safrica_PDiv.nc","r")<br>
;dPTWRH<br>
in5=addfile("/research/geog/<wbr>data2/DATA/mlazenby/Chadwick_<wbr>Analysis/<wbr>ChadDeltaPComponentsOND/<wbr>DeltaPDiv+PT+PRH_"+model(gg)+"<wbr>.nc","r")<br>
;dPCross<br>
in6=addfile("/research/geog/<wbr>data2/DATA/mlazenby/Chadwick_<wbr>Analysis/<wbr>ChadDeltaPComponentsOND/<wbr>DeltaPNL_"+model(gg)+"_<wbr>safrica_PNL.nc","r")<br>
;dPRH<br>
in7=addfile("/research/geog/<wbr>data2/DATA/mlazenby/Chadwick_<wbr>Analysis/<wbr>ChadDeltaPComponentsOND/<wbr>DeltaPRH_"+model(gg)+"_<wbr>safrica_PRH.nc","r")<br>
<br>
<br>
lat = in1->lat <br>
lon = in1->lon <br>
time = in1->time<br>
<br>
if ((gg.eq.0) .or. (gg.eq.9) .or. (gg.eq.10) .or. (gg.eq.11)) then<br>
dP = dble2flt(in1->pr(:,:,:))<br>
dPShift = in2->huss(:,:,:)<br>
dPtwrh = in5->huss(:,:,:)<br>
dPcross = in6->huss(:,:,:)<br>
dPT = in3->es2(:,:,:)<br>
dPweak = in4->huss(:,:,:)<br>
dPrh = in7->huss(:,:,:)<br>
<br>
else<br>
<br>
dP = in1->pr(:,0,:,:)<br>
;dPShift = in2->huss(:,0,:,:)<br>
;dPtwrh = in5->huss(:,0,:,:)<br>
;dPcross = in6->huss(:,0,:,:)<br>
;dPT = in3->es2(:,0,:,:)<br>
;dPweak = in4->huss(:,0,:,:)<br>
;dPrh = in7->huss(:,0,:,:)<br>
<br>
end if<br>
<br>
;Land mask<br>
<br>
b=addfile("$NCARG_ROOT/lib/<wbr>ncarg/data/cdf/<a href="http://landsea.nc" target="_blank">landsea.nc</a>","r"<wbr>)<br>
lsm_dP = landsea_mask(b->LSMASK,dP&lat,<wbr>dP&lon) ; read in land sea mask<br></div></div>
lsm_dPShift = landsea_mask(b->LSMASK,<wbr>dPShift&lat,dPShift&lon) ; read in land sea mask<br>
lsm_dPtwrh = landsea_mask(b->LSMASK,dPtwrh&<wbr>lat,dPtwrh&lon) ; read in land sea mask<br>
lsm_dPcross = landsea_mask(b->LSMASK,<wbr>dPcross&lat,dPcross&lon) ; read in land sea mask<br>
lsm_dPT = landsea_mask(b->LSMASK,dPT&<wbr>lat,dPT&lon) ; read in land sea mask<br>
lsm_dPweak = landsea_mask(b->LSMASK,dPweak&<wbr>lat,dPweak&lon) ; read in land sea mask<br>
lsm_dPrh = landsea_mask(b->LSMASK,dPrh&<wbr>lat,dPrh&lon) ; read in land sea mask<span class=""><br>
<br>
dP = mask(dP,lsm_dP.eq.0,False) ;mask out all ocean points)<br></span>
dPShift = mask(dPShift,lsm_dPShift.eq.0,<wbr>False) ;mask out all ocean points)<br>
dPtwrh = mask(dPtwrh,lsm_dPtwrh.eq.0,<wbr>False) ;mask out all ocean points)<br>
dPcross = mask(dPcross,lsm_dPcross.eq.0,<wbr>False) ;mask out all ocean points)<br>
dPT = mask(dPT,lsm_dPT.eq.0,False) <wbr> ;mask out all ocean points)<br>
dPweak = mask(dPweak,lsm_dPweak.eq.0,<wbr>False) ;mask out all ocean points)<br>
dPrh = mask(dPrh,lsm_dPrh.eq.0,False)<wbr> ;mask out all ocean points)<span class=""><br>
<br>
<br>
;Drymask<br>
dP_dry = new ((/dimsizes(time),dimsizes(<wbr>lat),dimsizes(lon)/), typeof(dP))<br>
dP_dry!0="time"<br>
dP_dry&time=time<br>
dP_dry!1="lat"<br>
dP_dry&lat=lat<br>
dP_dry!2="lon"<br>
dP_dry&lon=lon<br>
<br>
dP_dry = mask(dP,dP.gt.0, False) ; mask out all the wetting regions i.e. any values above 0 for dP<br>
<br></span>
dPShift_dry = new ((/dimsizes(time),dimsizes(<wbr>lat),dimsizes(lon)/), typeof(dPShift))<br>
dPShift_dry!0="time"<br>
dPShift_dry&time=time<br>
dPShift_dry!1="lat"<br>
dPShift_dry&lat=lat<br>
dPShift_dry!2="lon"<br>
dPShift_dry&lon=lon<br>
<br>
dPShift_dry = mask(dPShift,dPShift.gt.0, False) ; mask out all the wetting regions i.e. any values above 0 for dP<br>
<br>
dPtwrh_dry = new ((/dimsizes(time),dimsizes(<wbr>lat),dimsizes(lon)/), typeof(dPtwrh))<br>
dPtwrh_dry!0="time"<br>
dPtwrh_dry&time=time<br>
dPtwrh_dry!1="lat"<br>
dPtwrh_dry&lat=lat<br>
dPtwrh_dry!2="lon"<br>
dPtwrh_dry&lon=lon<br>
<br>
dPtwrh_dry = mask(dPtwrh,dPtwrh.gt.0, False) ; mask out all the wetting regions i.e. any values above 0 for dP<br>
<br>
dPcross_dry = new ((/dimsizes(time),dimsizes(<wbr>lat),dimsizes(lon)/), typeof(dPcross))<br>
dPcross_dry!0="time"<br>
dPcross_dry&time=time<br>
dPcross_dry!1="lat"<br>
dPcross_dry&lat=lat<br>
dPcross_dry!2="lon"<br>
dPcross_dry&lon=lon<br>
<br>
dPcross_dry = mask(dPcross,dPcross.gt.0, False) ; mask out all the wetting regions i.e. any values above 0 for dP<br>
<br>
dPT_dry = new ((/dimsizes(time),dimsizes(<wbr>lat),dimsizes(lon)/), typeof(dPT))<br>
dPT_dry!0="time"<br>
dPT_dry&time=time<br>
dPT_dry!1="lat"<br>
dPT_dry&lat=lat<br>
dPT_dry!2="lon"<br>
dPT_dry&lon=lon<br>
<br>
dPT_dry = mask(dPT,dPT.gt.0, False) ; mask out all the wetting regions i.e. any values above 0 for dP<br>
<br>
dPweak_dry = new ((/dimsizes(time),dimsizes(<wbr>lat),dimsizes(lon)/), typeof(dPweak))<br>
dPweak_dry!0="time"<br>
dPweak_dry&time=time<br>
dPweak_dry!1="lat"<br>
dPweak_dry&lat=lat<br>
dPweak_dry!2="lon"<br>
dPweak_dry&lon=lon<br>
<br>
dPweak_dry = mask(dPweak,dPweak.gt.0, False) ; mask out all the wetting regions i.e. any values above 0 for dP<br>
<br>
dPrh_dry = new ((/dimsizes(time),dimsizes(<wbr>lat),dimsizes(lon)/), typeof(dPrh))<br>
dPrh_dry!0="time"<br>
dPrh_dry&time=time<br>
dPrh_dry!1="lat"<br>
dPrh_dry&lat=lat<br>
dPrh_dry!2="lon"<br>
dPrh_dry&lon=lon<br>
<br>
dPrh_dry = mask(dPrh,dPrh.gt.0, False) ; mask out all the wetting regions i.e. any values above 0 for dP<br>
<br>
<br>
;if ((gg.eq.0) .or. (gg.eq.9) .or. (gg.eq.10) .or. (gg.eq.11)) then<br>
<br>
dP_fldmean = wgt_areaave(dP_dry(:,{-30:0},{<wbr>10:40}),1.0, 1.0, 0)<br>
dPShift_fldmean = wgt_areaave(dPShift_dry(:,{-<wbr>30:0},{10:40}),1.0, 1.0, 0)<br>
dPtwrh_fldmean = wgt_areaave(dPtwrh_dry(:,{-30:<wbr>0},{10:40}),1.0, 1.0, 0)<br>
dPcross_fldmean = wgt_areaave(dPcross_dry(:,{-<wbr>30:0},{10:40}),1.0, 1.0, 0)<br>
dPT_fldmean = wgt_areaave(dPT_dry(:,{-30:0},<wbr>{10:40}),1.0, 1.0, 0)<br>
dPweak_fldmean = wgt_areaave(dPweak_dry(:,{-30:<wbr>0},{10:40}),1.0, 1.0, 0)<br>
dPrh_fldmean = wgt_areaave(dPrh_dry(:,{-30:0}<wbr>,{10:40}),1.0, 1.0, 0)<br>
<br>
;else<br>
<br>
;dP_fldmean = wgt_areaave(dP_dry(:,0,{-30:0}<wbr>,{10:40}),1.0, 1.0, 0)<br>
;dPShift_fldmean = wgt_areaave(dPShift_dry(:,0,{-<wbr>30:0},{10:40}),1.0, 1.0, 0)<br>
;dPtwrh_fldmean = wgt_areaave(dPtwrh_dry(:,0,{-<wbr>30:0},{10:40}),1.0, 1.0, 0)<br>
;dPcross_fldmean = wgt_areaave(dPcross_dry(:,0,{-<wbr>30:0},{10:40}),1.0, 1.0, 0)<br>
;dPT_fldmean = wgt_areaave(dPT_dry(:,{-30:0},<wbr>0,{10:40}),1.0, 1.0, 0)<br>
;dPweak_fldmean = wgt_areaave(dPweak_dry(:,0,{-<wbr>30:0},{10:40}),1.0, 1.0, 0)<br>
;dPrh_fldmean = wgt_areaave(dPrh_dry(:,0,{-30:<wbr>0},{10:40}),1.0, 1.0, 0)<br>
<br>
<br>
;end if <br>
<br>
print(dP_fldmean)<br>
print(dPShift_fldmean)<br>
<br>
<br>
;Create new arrays of 20x40x240<br>
<br>
if (gg.eq.0) then<br>
dPfldavg = new ((/20/), typeof(dP_fldmean))<br>
<br>
dPfldavg!0="model"<br>
dPfldavg&model=model<br>
<br>
end if<br>
<br>
dPfldavg(gg) = dP_fldmean(0)<br>
<br>
if (gg.eq.0) then<br>
dPfldavgShift = new ((/20/), typeof(dPShift_fldmean))<br>
<br>
dPfldavgShift!0="model"<br>
dPfldavgShift&model=model<br>
<br>
end if<br>
<br>
dPfldavgShift(gg) = dPShift_fldmean(0)<br>
<br>
if (gg.eq.0) then<br>
dPfldavgtwrh = new ((/20/), typeof(dPtwrh_fldmean))<br>
<br>
dPfldavgtwrh!0="model"<br>
dPfldavgtwrh&model=model<br>
<br>
end if<br>
<br>
dPfldavgtwrh(gg) = dPtwrh_fldmean(0)<br>
<br>
if (gg.eq.0) then<br>
dPfldavgcross = new ((/20/), typeof(dPcross_fldmean))<br>
<br>
dPfldavgcross!0="model"<br>
dPfldavgcross&model=model<br>
<br>
end if<br>
<br>
dPfldavgcross(gg) = dPcross_fldmean(0)<br>
<br>
if (gg.eq.0) then<br>
dPfldavgT = new ((/20/), typeof(dPT_fldmean))<br>
<br>
dPfldavgT!0="model"<br>
dPfldavgT&model=model<br>
<br>
end if<br>
<br>
dPfldavgT(gg) = dPT_fldmean(0)<br>
<br>
if (gg.eq.0) then<br>
dPfldavgweak = new ((/20/), typeof(dPweak_fldmean))<br>
<br>
dPfldavgweak!0="model"<br>
dPfldavgweak&model=model<br>
<br>
end if<br>
<br>
dPfldavgweak(gg) = dPweak_fldmean(0)<br>
<br>
if (gg.eq.0) then<br>
dPfldavgrh = new ((/20/), typeof(dPrh_fldmean))<br>
<br>
dPfldavgrh!0="model"<br>
dPfldavgrh&model=model<br>
<br>
end if<br>
<br>
dPfldavgrh(gg) = dPrh_fldmean(0)<br>
<br>
;end do<br>
<br>
<br>
dimt = dimsizes(dPfldavg) ; should be 20<br>
x25 = round(.25*dimt,3)-1 ; -1 to account for NCL indexing starting<br>
x75 = round(.75*dimt,3)-1 ; at 0 <br>
qsort(dPfldavg)<br>
print(dPfldavg) ; Print the values<br>
<br>
dimt1 = dimsizes(dPfldavgShift) ; should be 20<br>
x25 = round(.25*dimt1,3)-1 ; -1 to account for NCL indexing starting<br>
x75 = round(.75*dimt1,3)-1 ; at 0 <br>
qsort(dPfldavgShift)<br>
print(dPfldavgShift) ; Print the values<br>
<br>
dimt2 = dimsizes(dPfldavgtwrh) ; should be 20<br>
x25 = round(.25*dimt2,3)-1 ; -1 to account for NCL indexing starting<br>
x75 = round(.75*dimt2,3)-1 ; at 0 <br>
qsort(dPfldavgtwrh)<br>
print(dPfldavgtwrh) ; Print the values<br>
<br>
dimt3 = dimsizes(dPfldavgT) ; should be 20<br>
x25 = round(.25*dimt3,3)-1 ; -1 to account for NCL indexing starting<br>
x75 = round(.75*dimt3,3)-1 ; at 0 <br>
qsort(dPfldavgT)<br>
print(dPfldavgT) ; Print the values<br>
<br>
dimt4 = dimsizes(dPfldavgcross) ; should be 20<br>
x25 = round(.25*dimt4,3)-1 ; -1 to account for NCL indexing starting<br>
x75 = round(.75*dimt4,3)-1 ; at 0 <br>
qsort(dPfldavgcross)<br>
print(dPfldavgcross) ; Print the values<br>
<br>
dimt5 = dimsizes(dPfldavgweak) ; should be 20<br>
x25 = round(.25*dimt5,3)-1 ; -1 to account for NCL indexing starting<br>
x75 = round(.75*dimt5,3)-1 ; at 0 <br>
qsort(dPfldavgweak)<br>
print(dPfldavgweak) ; Print the values<br>
<br>
dimt6 = dimsizes(dPfldavgrh) ; should be 20<br>
x25 = round(.25*dimt6,3)-1 ; -1 to account for NCL indexing starting<br>
x75 = round(.75*dimt6,3)-1 ; at 0 <br>
qsort(dPfldavgrh)<br>
print(dPfldavgrh) ; Print the values<br>
<br>
end do<br>
<br>
iarr=new((/7,5/),"float",-<wbr>999.) ; fill with minimum, 25th absile, median, 75th absile, maximum of each timeseries<br>
iarr(0,:) = (/min(dPfldavg),dPfldavg(x25),<wbr>dim_avg(dPfldavg),dPfldavg(<wbr>x75),max(dPfldavg)/)
<br>
iarr(1,:) = (/min(dPfldavgShift),<wbr>dPfldavgShift(x25),dim_avg(<wbr>dPfldavgShift),dPfldavgShift(<wbr>x75),max(dPfldavgShift)/)<br>
iarr(2,:) = (/min(dPfldavgtwrh),<wbr>dPfldavgtwrh(x25),dim_avg(<wbr>dPfldavgtwrh),dPfldavgtwrh(<wbr>x75),max(dPfldavgtwrh)/)
<br>
iarr(3,:) = (/min(dPfldavgcross),<wbr>dPfldavgcross(x25),dim_avg(<wbr>dPfldavgcross),dPfldavgcross(<wbr>x75),max(dPfldavgcross)/)<br>
iarr(4,:) = (/min(dPfldavgT),dPfldavgT(<wbr>x25),dim_avg(dPfldavgT),<wbr>dPfldavgT(x75),max(dPfldavgT)/<wbr>)<br>
iarr(5,:) = (/min(dPfldavgweak),<wbr>dPfldavgweak(x25),dim_avg(<wbr>dPfldavgweak),dPfldavgweak(<wbr>x75),max(dPfldavgweak)/)
<br>
iarr(6,:) = (/min(dPfldavgrh),dPfldavgrh(<wbr>x25),dim_avg(dPfldavgrh),<wbr>dPfldavgrh(x75),max(<wbr>dPfldavgrh)/)<br>
<br>
print(iarr)<br>
<br>
x = (/-9.,-6.,-3.,0.,3.,6.,9./)<br>
<br>
<br>
;*****************************<wbr>*****************<br>
; create plot<br>
;*****************************<wbr>*****************<br>
wks = gsn_open_wks("X11","NCL_JOC_<wbr>Boxplot_Land_Drying_abs_20_<wbr>OND")<br>
;*****************************<wbr>*****************<br>
; resources for plot background<br>
;*****************************<wbr>*****************<br>
res = True ; plot mods desired<br>
res@tmXBLabels = (/"~F8~D~F21~P","~F8~D~F21~P~<wbr>B~Shift","~F8~D~F21~P~B~twrh",<wbr>"~F8~D~F21~P~B~Cross","~F8~D~<wbr>F21~P~B~T","~F8~D~F21~P~B~<wbr>Weak","~F8~D~F21~P~B~RH"/) ; labels for each box ;<br>
res@tiMainString = "OND Box Plot over SA Drying (Land-Only) (10-40E 0S-30S)"<br>
res@trYMinF = -0.6 ; set minimum Y-axis value<br>
res@trYMaxF = 0.6 ; set maximum Y-axis value<br>
res@tiXAxisString = "Precipitation Components" <br>
res@tiYAxisString = "Absolute Change (mm/day/K)"<br>
res@tmXBLabelFontHeightF = 0.015 <br>
<br>
;*****************************<wbr>*****************<br>
; resources for polylines that draws the boxes<br>
;*****************************<wbr>***************** <br>
llres = True <br>
llres@gsLineThicknessF = 2.5 ; line thickness <br>
;*****************************<wbr>*****************<br>
; resources that control color and width of boxes<br>
;*****************************<wbr>***************** <br>
opti = True <br>
opti@boxWidth = 2.0 ; Width of box (x units)<br>
;colors = (/"red","blue","green","<wbr>yellow","magenta","orange","<wbr>cyan" /)<br>
opti@boxColors = (/"red","blue","green","<wbr>yellow","purple","orange","<wbr>cyan" /) ; Color of box(es);<br>
;gsn_define_colormap(wks, colors) <br>
;*****************************<wbr>******************<br>
plot = boxplot(wks,x,iarr,opti,res,<wbr>llres) ; All 3 options used...<br>
draw(wks) <wbr> ; box plot does not call<br>
frame(wks) <wbr> ; these for you<br>
<br>
;end do<br>
<br>
<br>
end<br>
<br>
<div style="font-family:Times New Roman;color:#000000;font-size:16px">
<hr>
<div style="direction:ltr" id="m_4712411266843353398divRpF244188"><font color="#000000" face="Tahoma" size="2"><b>From:</b> Melissa Lazenby<br>
<b>Sent:</b> 21 March 2017 12:41<br>
<b>To:</b> <a href="mailto:ncl-talk@ucar.edu" target="_blank">ncl-talk@ucar.edu</a><br>
<b>Subject:</b> Masking, area averaging and box whisker plots<br>
</font><br>
</div><div><div class="h5">
<div></div>
<div>
<div style="direction:ltr;font-family:Tahoma;color:#000000;font-size:10pt">
<div style="direction:ltr;font-family:Tahoma;color:#000000;font-size:10pt">Dear NCL Users<span class="m_4712411266843353398Apple-tab-span" style="white-space:pre-wrap">
</span>
<div><br>
</div>
<div>I am currently trying to create box whisker plots of drying regions over land regions.</div>
<div><br>
</div>
<div>I have tried creating the masks and field averages in CDO however there seems to be a bug in the CDO code I am running so I am trying to do all the masking and area averaging in an NCL script.</div>
<div><br>
</div>
<div>Essentially what I would like to do is:</div>
<div><br>
</div>
<div>1. Mask out the land regions (DONE)<br>
</div>
<div>2. Mask out the wetting regions (i.e. only have the drying regions) (DONE)<br>
</div>
<div>3. Area average the change in precipitation over a particular domain over the southern African continent (NOT SURE?? area average or dim_avg?)<br>
</div>
<div>4. Create a box whisker plot of the change in precipitation over that region for 20 different CMIP5 models to determine the robustness and spread.</div>
<div><br>
</div>
<div>Below is my code in NCL attempting to do 1-4<br>
<br>
I have managed to do steps 1 and 2 without problems, however I am struggling to complete step 3 i.e. create an field mean using the area average function.
<br>
Essentially I would like 20 output values from my 20 models of the area averages and currently I am only obtaining 1 value not 20.<br>
<br>
;*****************************<wbr>****************<br>
; boxplot.ncl<br>
;<br>
; ;*****************************<wbr>****************<br>
load "$NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/gsn_code.ncl"<br>
load "$NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/gsn_csm.ncl"<br>
load "$NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/contributed.<wbr>ncl"<br>
load "$NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/shea_util.ncl"<br>
;*****************************<wbr>****************<br>
begin<br>
;*****************************<wbr>*****************<br>
; Add data in text format and sort in ascending order<br>
;*****************************<wbr>*****************<br>
<br>
<br>
model = (/"bcc-csm1-1-m","BNU-ESM","<wbr>CanESM2","CCSM4","CESM1-BGC","<wbr>CESM1-CAM5","CSIRO-Mk3-6-0","<wbr>FIO-ESM","GFDL-CM3","GFDL-<wbr>ESM2G","GFDL-ESM2M","GISS-E2-<wbr>H","HadGEM2-CC","HadGEM2-ES","<wbr>IPSL-CM5A-LR","IPSL-CM5A-MR","<wbr>MIROC5","MRI-CGCM3","NorESM1-<wbr>M","NorESM1-ME"/)<br>
<br>
<br>
do gg = 0,dimsizes(model)-1<br>
;dP<br>
in1=addfile("/research/geog/<wbr>data2/DATA/mlazenby/Chadwick_<wbr>Analysis/ChadprfutureOND/<wbr>Delta_pr_"+model(gg)+".nc","r"<wbr>)
<br>
<br>
;dPShift<br>
in2=addfile("/research/geog/<wbr>data2/DATA/mlazenby/Chadwick_<wbr>Analysis/<wbr>ChadDeltaPComponentsOND/<wbr>DeltaPSpat_"+model(gg)+"_<wbr>safrica_PSpat.nc","r")<br>
;dPT<br>
in3=addfile("/research/geog/<wbr>data2/DATA/mlazenby/Chadwick_<wbr>Analysis/<wbr>ChadDeltaPComponentsOND/<wbr>DeltaPT_"+model(gg)+"_safrica_<wbr>PT.nc","r")<br>
;dPWeak<br>
in4=addfile("/research/geog/<wbr>data2/DATA/mlazenby/Chadwick_<wbr>Analysis/<wbr>ChadDeltaPComponentsOND/<wbr>DeltaPDiv_"+model(gg)+"_<wbr>safrica_PDiv.nc","r")<br>
;dPTWRH<br>
in5=addfile("/research/geog/<wbr>data2/DATA/mlazenby/Chadwick_<wbr>Analysis/<wbr>ChadDeltaPComponentsOND/<wbr>DeltaPDiv+PT+PRH_"+model(gg)+"<wbr>.nc","r")<br>
;dPCross<br>
in6=addfile("/research/geog/<wbr>data2/DATA/mlazenby/Chadwick_<wbr>Analysis/<wbr>ChadDeltaPComponentsOND/<wbr>DeltaPNL_"+model(gg)+"_<wbr>safrica_PNL.nc","r")<br>
;dPRH<br>
in7=addfile("/research/geog/<wbr>data2/DATA/mlazenby/Chadwick_<wbr>Analysis/<wbr>ChadDeltaPComponentsOND/<wbr>DeltaPRH_"+model(gg)+"_<wbr>safrica_PRH.nc","r")<br>
<br>
<br>
lat = in1->lat <br>
lon = in1->lon <br>
time = in1->time<br>
<br>
if ((gg.eq.0) .or. (gg.eq.9) .or. (gg.eq.10) .or. (gg.eq.11)) then<br>
dP = dble2flt(in1->pr(:,:,:))<br>
dPShift = in2->huss(:,:,:)<br>
dPtwrh = in5->huss(:,:,:)<br>
dPcross = in6->huss(:,:,:)<br>
dPT = in3->es2(:,:,:)<br>
dPweak = in4->huss(:,:,:)<br>
dPrh = in7->huss(:,:,:)<br>
<br>
else<br>
<br>
dP = in1->pr(:,0,:,:)<br>
dPShift = in2->huss(:,0,:,:)<br>
dPtwrh = in5->huss(:,0,:,:)<br>
dPcross = in6->huss(:,0,:,:)<br>
dPT = in3->es2(:,0,:,:)<br>
dPweak = in4->huss(:,0,:,:)<br>
dPrh = in7->huss(:,0,:,:)<br>
<br>
end if<br>
<br>
;Land mask<br>
<br>
b=addfile("$NCARG_ROOT/lib/<wbr>ncarg/data/cdf/<a href="http://landsea.nc" target="_blank">landsea.nc</a>","r"<wbr>)<br>
lsm_dP = landsea_mask(b->LSMASK,dP&lat,<wbr>dP&lon) ; read in land sea mask<br>
dP = mask(dP,lsm_dP.eq.0,False) ;mask out all ocean points)<br>
<br>
;Drymask<br>
<br>
dP_dry = new ((/dimsizes(time),dimsizes(<wbr>lat),dimsizes(lon)/), typeof(dP))<br>
dP_dry!0="time"<br>
dP_dry&time=time<br>
dP_dry!1="lat"<br>
dP_dry&lat=lat<br>
dP_dry!2="lon"<br>
dP_dry&lon=lon<br>
<br>
dP_dry = mask(dP,dP.gt.0, False) ; mask out all the wetting regions i.e. any values above 0 for dP<br>
<br>
<br>
printVarSummary(dP)<br>
printVarSummary(dP_dry)<br>
<br>
dP_fldmean = wgt_areaave(dP_dry((gg),{-30:<wbr>0},{10:40}),1.0, 1.0, 0)<br>
;dP_fldmean = dim_avg(dP_dry(0,{-30:0},{10:<wbr>40}))<br>
<br>
print(dP_dry)<br>
printVarSummary(dP_fldmean)<br>
print(dP_fldmean)<br>
<br>
end do <br>
<br>
end <br>
</div>
<div><br>
</div>
<div>Any advice regarding this code would very much be appreciated. </div>
<div><br>
Kindest Regards<br>
Melissa<br>
</div>
<div><br>
</div>
</div>
</div>
</div>
</div></div></div>
</div>
</div>
<br>______________________________<wbr>_________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@ucar.edu">ncl-talk@ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">http://mailman.ucar.edu/<wbr>mailman/listinfo/ncl-talk</a><br>
<br></blockquote></div><br></div>