[ncl-talk] Masking, area averaging and box whisker plots

[Melissa Lazenby]

Dear NCL Users

I am currently trying to create box whisker plots of drying regions over land regions.

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.

Essentially what I would like to do is:

1. Mask out the land regions  (DONE)
2. Mask out the wetting regions (i.e. only have the drying regions) (DONE)
3. Area average the change in precipitation over a particular domain over the southern African continent (NOT SURE?? area average or dim_avg?)
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.

Below is my code in NCL attempting to do 1-4

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.
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.

;*********************************************
; boxplot.ncl
;
; ;*********************************************
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/csm/shea_util.ncl"
;*********************************************
begin
;**********************************************
; Add data in text format and sort in ascending order
;**********************************************


 model = (/"bcc-csm1-1-m","BNU-ESM","CanESM2","CCSM4","CESM1-BGC","CESM1-CAM5","CSIRO-Mk3-6-0","FIO-ESM","GFDL-CM3","GFDL-ESM2G","GFDL-ESM2M","GISS-E2-H","HadGEM2-CC","HadGEM2-ES","IPSL-CM5A-LR","IPSL-CM5A-MR","MIROC5","MRI-CGCM3","NorESM1-M","NorESM1-ME"/)


do gg = 0,dimsizes(model)-1
;dP
  in1=addfile("/research/geog/data2/DATA/mlazenby/Chadwick_Analysis/ChadprfutureOND/Delta_pr_"+model(gg)+".nc","r")

;dPShift
in2=addfile("/research/geog/data2/DATA/mlazenby/Chadwick_Analysis/ChadDeltaPComponentsOND/DeltaPSpat_"+model(gg)+"_safrica_PSpat.nc","r")
;dPT
in3=addfile("/research/geog/data2/DATA/mlazenby/Chadwick_Analysis/ChadDeltaPComponentsOND/DeltaPT_"+model(gg)+"_safrica_PT.nc","r")
;dPWeak
in4=addfile("/research/geog/data2/DATA/mlazenby/Chadwick_Analysis/ChadDeltaPComponentsOND/DeltaPDiv_"+model(gg)+"_safrica_PDiv.nc","r")
;dPTWRH
in5=addfile("/research/geog/data2/DATA/mlazenby/Chadwick_Analysis/ChadDeltaPComponentsOND/DeltaPDiv+PT+PRH_"+model(gg)+".nc","r")
;dPCross
in6=addfile("/research/geog/data2/DATA/mlazenby/Chadwick_Analysis/ChadDeltaPComponentsOND/DeltaPNL_"+model(gg)+"_safrica_PNL.nc","r")
;dPRH
in7=addfile("/research/geog/data2/DATA/mlazenby/Chadwick_Analysis/ChadDeltaPComponentsOND/DeltaPRH_"+model(gg)+"_safrica_PRH.nc","r")


  lat    = in1->lat
  lon    = in1->lon
  time   = in1->time

if ((gg.eq.0) .or. (gg.eq.9) .or. (gg.eq.10) .or. (gg.eq.11)) then
dP = dble2flt(in1->pr(:,:,:))
dPShift = in2->huss(:,:,:)
dPtwrh = in5->huss(:,:,:)
dPcross = in6->huss(:,:,:)
dPT = in3->es2(:,:,:)
dPweak = in4->huss(:,:,:)
dPrh = in7->huss(:,:,:)

else

dP = in1->pr(:,0,:,:)
dPShift = in2->huss(:,0,:,:)
dPtwrh = in5->huss(:,0,:,:)
dPcross = in6->huss(:,0,:,:)
dPT = in3->es2(:,0,:,:)
dPweak = in4->huss(:,0,:,:)
dPrh = in7->huss(:,0,:,:)

end if

    ;Land mask

    b=addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r")
    lsm_dP = landsea_mask(b->LSMASK,dP&lat,dP&lon)    ; read in land sea mask
    dP = mask(dP,lsm_dP.eq.0,False)    ;mask out all ocean points)

    ;Drymask

        dP_dry = new ((/dimsizes(time),dimsizes(lat),dimsizes(lon)/), typeof(dP))
       dP_dry!0="time"
        dP_dry&time=time
       dP_dry!1="lat"
        dP_dry&lat=lat
        dP_dry!2="lon"
        dP_dry&lon=lon

    dP_dry = mask(dP,dP.gt.0, False)  ; mask out all the wetting regions i.e. any values above 0 for dP


printVarSummary(dP)
printVarSummary(dP_dry)

dP_fldmean = wgt_areaave(dP_dry((gg),{-30:0},{10:40}),1.0, 1.0, 0)
 ;dP_fldmean = dim_avg(dP_dry(0,{-30:0},{10:40}))

print(dP_dry)
printVarSummary(dP_fldmean)
print(dP_fldmean)

end do

end

Any advice regarding this code would very much be appreciated.

Kindest Regards
Melissa

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170321/45555359/attachment.html