[ncl-talk] Masking, area averaging and box whisker plots
Dennis Shea
shea at ucar.edu
Wed Mar 22 12:03:29 MDT 2017
A bit involved for a simpleton such as myself.
[1] Code clarity:
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.
[2] Use a procedure or function to isolate repeated tasks. You do the
following task multiple times.
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
[3] I suggest a function that isolates the task.
undef("createVarMeta")
function createVarMeta(time[*]:numeric,lat[*]:numeric, lon[*]:numeric \
,varType[1]:string,long_name[1]:string, units[1]:string)
local x
begin
x = new ((/dimsizes(time),dimsizes(lat),dimsizes(lon)/), varType)
x!0="time"
x&time=time
x!1="lat"
x&lat=lat
x!2="lon"
x&lon=lon
if (long_name.ne."") the
var at long_name = long_name
end if
if (units.ne."") the
var at units =units
end if
return(x)
end
[4] Use it like:
dPShift_dry = createVarMeta(time, lat, lon, typeof(...), "dummy" ,
"whatever")
DO this for all variable you are creating.
The function should be placed at the top of the script
On Tue, Mar 21, 2017 at 10:12 AM, Melissa Lazenby <M.Lazenby at sussex.ac.uk>
wrote:
> 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
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170322/c1d90972/attachment.html
More information about the ncl-talk
mailing list