[ncl-talk] Masking, area averaging and box whisker plots
Melissa Lazenby
M.Lazenby at sussex.ac.uk
Tue Mar 21 10:12:40 MDT 2017
Hi All
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?
My overall aim is to:
1. mask out all ocean regions
2. mask out all the wetting regions
3. area average these changes over southern africa
4. create boxplots of dP and 6 decomposed precipitation components
I feel the code does this correctly for dP, but not convinced for the other components.
Any advice would very much be appreciated regarding this issue. Below is the full code to tackle task 1 -4
Many thanks in advance!
Kindest Regards
Melissa
;*********************************************
; boxplot.ncl
;
; Concepts illustrated:
; - Drawing box plots
; - Setting the color of individual boxes in a box plot
; - Setting the width of individual boxes in a box plot
;
;*********************************************
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
;**********************************************
latS = -30.
latN = 0.
lonL = 10.
lonR = 40.
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
lsm_dPShift = landsea_mask(b->LSMASK,dPShift&lat,dPShift&lon) ; read in land sea mask
lsm_dPtwrh = landsea_mask(b->LSMASK,dPtwrh&lat,dPtwrh&lon) ; read in land sea mask
lsm_dPcross = landsea_mask(b->LSMASK,dPcross&lat,dPcross&lon) ; read in land sea mask
lsm_dPT = landsea_mask(b->LSMASK,dPT&lat,dPT&lon) ; read in land sea mask
lsm_dPweak = landsea_mask(b->LSMASK,dPweak&lat,dPweak&lon) ; read in land sea mask
lsm_dPrh = landsea_mask(b->LSMASK,dPrh&lat,dPrh&lon) ; read in land sea mask
dP = mask(dP,lsm_dP.eq.0,False) ;mask out all ocean points)
dPShift = mask(dPShift,lsm_dPShift.eq.0,False) ;mask out all ocean points)
dPtwrh = mask(dPtwrh,lsm_dPtwrh.eq.0,False) ;mask out all ocean points)
dPcross = mask(dPcross,lsm_dPcross.eq.0,False) ;mask out all ocean points)
dPT = mask(dPT,lsm_dPT.eq.0,False) ;mask out all ocean points)
dPweak = mask(dPweak,lsm_dPweak.eq.0,False) ;mask out all ocean points)
dPrh = mask(dPrh,lsm_dPrh.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
dPShift_dry = new ((/dimsizes(time),dimsizes(lat),dimsizes(lon)/), typeof(dPShift))
dPShift_dry!0="time"
dPShift_dry&time=time
dPShift_dry!1="lat"
dPShift_dry&lat=lat
dPShift_dry!2="lon"
dPShift_dry&lon=lon
dPShift_dry = mask(dPShift,dPShift.gt.0, False) ; mask out all the wetting regions i.e. any values above 0 for dP
dPtwrh_dry = new ((/dimsizes(time),dimsizes(lat),dimsizes(lon)/), typeof(dPtwrh))
dPtwrh_dry!0="time"
dPtwrh_dry&time=time
dPtwrh_dry!1="lat"
dPtwrh_dry&lat=lat
dPtwrh_dry!2="lon"
dPtwrh_dry&lon=lon
dPtwrh_dry = mask(dPtwrh,dPtwrh.gt.0, False) ; mask out all the wetting regions i.e. any values above 0 for dP
dPcross_dry = new ((/dimsizes(time),dimsizes(lat),dimsizes(lon)/), typeof(dPcross))
dPcross_dry!0="time"
dPcross_dry&time=time
dPcross_dry!1="lat"
dPcross_dry&lat=lat
dPcross_dry!2="lon"
dPcross_dry&lon=lon
dPcross_dry = mask(dPcross,dPcross.gt.0, False) ; mask out all the wetting regions i.e. any values above 0 for dP
dPT_dry = new ((/dimsizes(time),dimsizes(lat),dimsizes(lon)/), typeof(dPT))
dPT_dry!0="time"
dPT_dry&time=time
dPT_dry!1="lat"
dPT_dry&lat=lat
dPT_dry!2="lon"
dPT_dry&lon=lon
dPT_dry = mask(dPT,dPT.gt.0, False) ; mask out all the wetting regions i.e. any values above 0 for dP
dPweak_dry = new ((/dimsizes(time),dimsizes(lat),dimsizes(lon)/), typeof(dPweak))
dPweak_dry!0="time"
dPweak_dry&time=time
dPweak_dry!1="lat"
dPweak_dry&lat=lat
dPweak_dry!2="lon"
dPweak_dry&lon=lon
dPweak_dry = mask(dPweak,dPweak.gt.0, False) ; mask out all the wetting regions i.e. any values above 0 for dP
dPrh_dry = new ((/dimsizes(time),dimsizes(lat),dimsizes(lon)/), typeof(dPrh))
dPrh_dry!0="time"
dPrh_dry&time=time
dPrh_dry!1="lat"
dPrh_dry&lat=lat
dPrh_dry!2="lon"
dPrh_dry&lon=lon
dPrh_dry = mask(dPrh,dPrh.gt.0, False) ; mask out all the wetting regions i.e. any values above 0 for dP
;if ((gg.eq.0) .or. (gg.eq.9) .or. (gg.eq.10) .or. (gg.eq.11)) then
dP_fldmean = wgt_areaave(dP_dry(:,{-30:0},{10:40}),1.0, 1.0, 0)
dPShift_fldmean = wgt_areaave(dPShift_dry(:,{-30:0},{10:40}),1.0, 1.0, 0)
dPtwrh_fldmean = wgt_areaave(dPtwrh_dry(:,{-30:0},{10:40}),1.0, 1.0, 0)
dPcross_fldmean = wgt_areaave(dPcross_dry(:,{-30:0},{10:40}),1.0, 1.0, 0)
dPT_fldmean = wgt_areaave(dPT_dry(:,{-30:0},{10:40}),1.0, 1.0, 0)
dPweak_fldmean = wgt_areaave(dPweak_dry(:,{-30:0},{10:40}),1.0, 1.0, 0)
dPrh_fldmean = wgt_areaave(dPrh_dry(:,{-30:0},{10:40}),1.0, 1.0, 0)
;else
;dP_fldmean = wgt_areaave(dP_dry(:,0,{-30:0},{10:40}),1.0, 1.0, 0)
;dPShift_fldmean = wgt_areaave(dPShift_dry(:,0,{-30:0},{10:40}),1.0, 1.0, 0)
;dPtwrh_fldmean = wgt_areaave(dPtwrh_dry(:,0,{-30:0},{10:40}),1.0, 1.0, 0)
;dPcross_fldmean = wgt_areaave(dPcross_dry(:,0,{-30:0},{10:40}),1.0, 1.0, 0)
;dPT_fldmean = wgt_areaave(dPT_dry(:,{-30:0},0,{10:40}),1.0, 1.0, 0)
;dPweak_fldmean = wgt_areaave(dPweak_dry(:,0,{-30:0},{10:40}),1.0, 1.0, 0)
;dPrh_fldmean = wgt_areaave(dPrh_dry(:,0,{-30:0},{10:40}),1.0, 1.0, 0)
;end if
print(dP_fldmean)
print(dPShift_fldmean)
;Create new arrays of 20x40x240
if (gg.eq.0) then
dPfldavg = new ((/20/), typeof(dP_fldmean))
dPfldavg!0="model"
dPfldavg&model=model
end if
dPfldavg(gg) = dP_fldmean(0)
if (gg.eq.0) then
dPfldavgShift = new ((/20/), typeof(dPShift_fldmean))
dPfldavgShift!0="model"
dPfldavgShift&model=model
end if
dPfldavgShift(gg) = dPShift_fldmean(0)
if (gg.eq.0) then
dPfldavgtwrh = new ((/20/), typeof(dPtwrh_fldmean))
dPfldavgtwrh!0="model"
dPfldavgtwrh&model=model
end if
dPfldavgtwrh(gg) = dPtwrh_fldmean(0)
if (gg.eq.0) then
dPfldavgcross = new ((/20/), typeof(dPcross_fldmean))
dPfldavgcross!0="model"
dPfldavgcross&model=model
end if
dPfldavgcross(gg) = dPcross_fldmean(0)
if (gg.eq.0) then
dPfldavgT = new ((/20/), typeof(dPT_fldmean))
dPfldavgT!0="model"
dPfldavgT&model=model
end if
dPfldavgT(gg) = dPT_fldmean(0)
if (gg.eq.0) then
dPfldavgweak = new ((/20/), typeof(dPweak_fldmean))
dPfldavgweak!0="model"
dPfldavgweak&model=model
end if
dPfldavgweak(gg) = dPweak_fldmean(0)
if (gg.eq.0) then
dPfldavgrh = new ((/20/), typeof(dPrh_fldmean))
dPfldavgrh!0="model"
dPfldavgrh&model=model
end if
dPfldavgrh(gg) = dPrh_fldmean(0)
;end do
dimt = dimsizes(dPfldavg) ; should be 20
x25 = round(.25*dimt,3)-1 ; -1 to account for NCL indexing starting
x75 = round(.75*dimt,3)-1 ; at 0
qsort(dPfldavg)
print(dPfldavg) ; Print the values
dimt1 = dimsizes(dPfldavgShift) ; should be 20
x25 = round(.25*dimt1,3)-1 ; -1 to account for NCL indexing starting
x75 = round(.75*dimt1,3)-1 ; at 0
qsort(dPfldavgShift)
print(dPfldavgShift) ; Print the values
dimt2 = dimsizes(dPfldavgtwrh) ; should be 20
x25 = round(.25*dimt2,3)-1 ; -1 to account for NCL indexing starting
x75 = round(.75*dimt2,3)-1 ; at 0
qsort(dPfldavgtwrh)
print(dPfldavgtwrh) ; Print the values
dimt3 = dimsizes(dPfldavgT) ; should be 20
x25 = round(.25*dimt3,3)-1 ; -1 to account for NCL indexing starting
x75 = round(.75*dimt3,3)-1 ; at 0
qsort(dPfldavgT)
print(dPfldavgT) ; Print the values
dimt4 = dimsizes(dPfldavgcross) ; should be 20
x25 = round(.25*dimt4,3)-1 ; -1 to account for NCL indexing starting
x75 = round(.75*dimt4,3)-1 ; at 0
qsort(dPfldavgcross)
print(dPfldavgcross) ; Print the values
dimt5 = dimsizes(dPfldavgweak) ; should be 20
x25 = round(.25*dimt5,3)-1 ; -1 to account for NCL indexing starting
x75 = round(.75*dimt5,3)-1 ; at 0
qsort(dPfldavgweak)
print(dPfldavgweak) ; Print the values
dimt6 = dimsizes(dPfldavgrh) ; should be 20
x25 = round(.25*dimt6,3)-1 ; -1 to account for NCL indexing starting
x75 = round(.75*dimt6,3)-1 ; at 0
qsort(dPfldavgrh)
print(dPfldavgrh) ; Print the values
end do
iarr=new((/7,5/),"float",-999.) ; fill with minimum, 25th absile, median, 75th absile, maximum of each timeseries
iarr(0,:) = (/min(dPfldavg),dPfldavg(x25),dim_avg(dPfldavg),dPfldavg(x75),max(dPfldavg)/)
iarr(1,:) = (/min(dPfldavgShift),dPfldavgShift(x25),dim_avg(dPfldavgShift),dPfldavgShift(x75),max(dPfldavgShift)/)
iarr(2,:) = (/min(dPfldavgtwrh),dPfldavgtwrh(x25),dim_avg(dPfldavgtwrh),dPfldavgtwrh(x75),max(dPfldavgtwrh)/)
iarr(3,:) = (/min(dPfldavgcross),dPfldavgcross(x25),dim_avg(dPfldavgcross),dPfldavgcross(x75),max(dPfldavgcross)/)
iarr(4,:) = (/min(dPfldavgT),dPfldavgT(x25),dim_avg(dPfldavgT),dPfldavgT(x75),max(dPfldavgT)/)
iarr(5,:) = (/min(dPfldavgweak),dPfldavgweak(x25),dim_avg(dPfldavgweak),dPfldavgweak(x75),max(dPfldavgweak)/)
iarr(6,:) = (/min(dPfldavgrh),dPfldavgrh(x25),dim_avg(dPfldavgrh),dPfldavgrh(x75),max(dPfldavgrh)/)
print(iarr)
x = (/-9.,-6.,-3.,0.,3.,6.,9./)
;**********************************************
; create plot
;**********************************************
wks = gsn_open_wks("X11","NCL_JOC_Boxplot_Land_Drying_abs_20_OND")
;**********************************************
; resources for plot background
;**********************************************
res = True ; plot mods desired
res at tmXBLabels = (/"~F8~D~F21~P","~F8~D~F21~P~B~Shift","~F8~D~F21~P~B~twrh","~F8~D~F21~P~B~Cross","~F8~D~F21~P~B~T","~F8~D~F21~P~B~Weak","~F8~D~F21~P~B~RH"/) ; labels for each box ;
res at tiMainString = "OND Box Plot over SA Drying (Land-Only) (10-40E 0S-30S)"
res at trYMinF = -0.6 ; set minimum Y-axis value
res at trYMaxF = 0.6 ; set maximum Y-axis value
res at tiXAxisString = "Precipitation Components"
res at tiYAxisString = "Absolute Change (mm/day/K)"
res at tmXBLabelFontHeightF = 0.015
;**********************************************
; resources for polylines that draws the boxes
;**********************************************
llres = True
llres at gsLineThicknessF = 2.5 ; line thickness
;**********************************************
; resources that control color and width of boxes
;**********************************************
opti = True
opti at boxWidth = 2.0 ; Width of box (x units)
;colors = (/"red","blue","green","yellow","magenta","orange","cyan" /)
opti at boxColors = (/"red","blue","green","yellow","purple","orange","cyan" /) ; Color of box(es);
;gsn_define_colormap(wks, colors)
;***********************************************
plot = boxplot(wks,x,iarr,opti,res,llres) ; All 3 options used...
draw(wks) ; box plot does not call
frame(wks) ; these for you
;end do
end
________________________________
From: Melissa Lazenby
Sent: 21 March 2017 12:41
To: ncl-talk at ucar.edu
Subject: Masking, area averaging and box whisker plots
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/7e6b0ff9/attachment.html
More information about the ncl-talk
mailing list