Re: rainfall calculations

Hi Anil,

I bet the units error that you are getting is a result of rainc not
having any metdata. You need to use addfiles_GetVar to keep the
metadata, or you will have to assign the lat/lon/time metadata yourself.

Instead of creating a new file for use in WRF_map_c as Mary suggested, I
think you should just read in a single WRF file and use it solely in the
WRF_map_c procedure. Try running the script as it is below and see if
that works. Adam

;*************************************************
; WRF: RPRECIPITATION: Total, Cumulus and non-cumulus prc
;************************************************
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/wrf/WRF_contributed.ncl"
begin
;************************************************
; open file
; Read Cumulus (rinc) and Non-cumulus (rainnc) prc
;************************************************
diri = "/datatmp3/anilkuma/temp/"
fili = systemfunc("cd "+diri +"; ls *d02*")
fili = fili +".nc" ; explicitly add file extension

F = addfiles (diri+fili, "r") ; refers to *all* files
f = addfile (diri+fili(0), "r") ; refers to 1st file only,
                                              ; will use this later
;************************************************
; open file and read in data
; (1) Read RAINC
; (2) Read character variable Times; Convert to string for plots
;************************************************
times = chartostring(F[:]->Times) ; built-in function
ntim = dimsizes(times) ; # of time steps
print(ntim)
print(times)

print("addfiles_GetVar should be used to retain metadata")
  rainc = addfiles_GetVar(F,fili,"RAINC")
  rainnc= addfiles_GetVar(F,fili,"RAINNC")
; rainc = F[:]->RAINC ; (Time, south_north, west_east)
; rainnc = F[:]->RAINNC

;************************************************
; Use NCL operator > to make sure all values >=0.0
; Sum components and assign attributes
;************************************************
  rainc = rainc > 0.0
  rainnc = rainnc > 0.0
  rainTot = rainc ; create new array that contains metadata

  rainTot = (/ rainc + rainnc /) ; overwrite data, but keep metadata
  rainTot_at_description = "Total Precipitation"
  rainTot_at_units = rainc_at_units
;************************************************
; create plots: create colormap using named colors
; unequal contour levels
;************************************************
  wks = gsn_open_wks("ps" ,"rain_cumulus2512") ;
ps,pdf,x11,ncgm,eps
  colors = (/"white","black" \ ; {for/back}ground
            ,"white" \
            ,"white","palegreen","green", "greenyellow" \
            ,"yellow","goldenrod","orange","orangered" \
            ,"red","deeppinK", "violet","darkviolet" \
            ,"blueviolet","blue" /)
  gsn_define_colormap(wks, colors)

  res = True ; plot mods desired
;;res_at_gsnMaximize = True ; uncomment to maximize size
  res_at_gsnSpreadColors = True ; use full range of colormap
  res_at_cnFillOn = True ; color plot desired
  res_at_cnLinesOn = False ; turn off contour lines
  res_at_cnLineLabelsOn = False ; turn off contour labels

  res_at_cnLevelSelectionMode = "ExplicitLevels" ; explicit [unequal] cn
levels
  res_at_cnLevels =
(/0,10,20.5,50,70.5,100,150,250,325,437.5,550,675,700,825,950/)
  res_at_cnFillMode = "RasterFill"
  res_at_lbOrientation = "Vertical" ; default is horizontal
;************************************************
; Use WRF_contributed procedure to set map resources
;************************************************
  WRF_map_c(f, res, 0) ; reads info from file f
;************************************************
; set True for native projection (faster)
;************************************************
   res_at_tfDoNDCOverlay = True

;****************************************************************************
; create panel of different components
;****************************************************************************
  plts = new (3 , "graphic") ; 1d array to hold plots

  res_at_gsnDraw = False ; do not draw
  res_at_gsnFrame = False ; do not advance 'frame'
  res_at_lbLabelBarOn = False ; turn off individual lb's
;************************************************
; create panel: panel plots have their own set of resources
;************************************************
  resP = True ; modify the panel plot
  resP_at_gsnMaximize = True ; maximize panel area
  resP_at_gsnPanelRowSpec = True ; specify 1 top, 2 lowerlevel
  resP_at_gsnPanelLabelBar = True ; add common colorbar
  resP_at_pmLabelBarWidthF = 0.85 ; make label wider
  resP_at_lbLabelFontHeightF = 0.015 ; default 0.02

  nt = 8 ; demo only one time
do nt=0,ntim-1 ; uncomment to loop over
all times
     plts(0) = gsn_csm_contour_map(wks,rainTot(nt,:,:),res)
     plts(1) = gsn_csm_contour_map(wks,rainnc(nt,:,:),res)
     plts(2) = gsn_csm_contour_map(wks,rainc(nt,:,:),res)

     resP_at_txString = f_at_TITLE+": "+times(nt)
     gsn_panel(wks,plts,(/1,2/),resP) ; now draw as one plot
end do
end

anil rohilla wrote:
> Dear lists i am now able to read the times in my files but still i am
> geeting this error
>
> ----------------------------oupte of ncl----------------------
> Variable: F
> Type: list
> Total Size: 4 bytes
> 1 values
> Number of Dimensions: 1
> Dimensions and sizes: [1]
> Coordinates:
>
>
> Variable: ntim
> Type: integer
> Total Size: 4 bytes
> 1 values
> Number of Dimensions: 1
> Dimensions and sizes: [1]
> Coordinates:
> (0) 8
>
>
> Variable: times
> Type: string
> Total Size: 32 bytes
> 8 values
> Number of Dimensions: 1
> Dimensions and sizes: [8]
> Coordinates:
> (0) 2006-07-01_00:00:00
> (1) 2006-07-01_03:00:00
> (2) 2006-07-01_06:00:00
> (3) 2006-07-01_09:00:00
> (4) 2006-07-01_12:00:00
> (5) 2006-07-01_15:00:00
> (6) 2006-07-01_18:00:00
> (7) 2006-07-01_21:00:00
> warning:Attempt to reference attribute (units) which is undefined
> fatal:Argument type mismatch on argument (0) of (WRF_map_c) can not coerce
> fatal:Execute: Error occurred at or near line 77 in file myplot.ncl
> -----------------------------------------------------------------------------
>
>
> my script is
> ------------------------------
>
> ;*************************************************
> ; WRF: RPRECIPITATION: Total, Cumulus and non-cumulus prc
> ;************************************************
> 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/wrf/WRF_contributed.ncl"
> begin
> ;************************************************
> ; open file
> ; Read Cumulus (rinc) and Non-cumulus (rainnc) prc
> ;************************************************
> diri = "/datatmp3/anilkuma/temp/"
> fili = systemfunc("cd "+diri +"; ls *d02*")
> fili = fili +".nc" ; explicitly add file
> extension
>
> F = addfiles (diri+fili, "r") ; refers to *all* files
> ; f = addfile (diri+fili(0), "r") ; refers to 1st file only
> print(F)
> ;************************************************
> ; open file and read in data
> ; (1) Read RAINC
> ; (2) Read character variable Times; Convert to string for plots
> ;************************************************
> times = chartostring(F[:]->Times) ; built-in function
> ntim = dimsizes(times) ; # of time steps
> print(ntim)
> print(times)
>
>
>
>
>
>
> ; f =
> addfile("/datatmp3/anilkuma/temp/wrfout_d02_2006-07-01_12:00:00.nc",
> "r")
> rainc = F[:]->RAINC ; (Time, south_north, west_east)
> rainnc = F[:]->RAINNC
>
> times = chartostring(F[:]->Times) ; convert to type string [plot]
> ntim = dimsizes(times) ; # time steps
> ;************************************************
> ; Use NCL operator > to make sure all values >=0.0
> ; Sum components and assign attributes
> ;************************************************
> rainc = rainc > 0.0
> rainnc = rainnc > 0.0
> rainTot = rainc + rainnc
> rainTot_at_description = "Total Precipitation"
> rainTot_at_units = rainc_at_units
> ;************************************************
> ; create plots: create colormap using named colors
> ; unequal contour levels
> ;************************************************
> wks = gsn_open_wks("ps" ,"rain_cumulus2512") ;
> ps,pdf,x11,ncgm,eps
> colors = (/"white","black" \ ; {for/back}ground
> ,"white" \
> ,"white","palegreen","green", "greenyellow" \
> ,"yellow","goldenrod","orange","orangered" \
> ,"red","deeppinK", "violet","darkviolet" \
> ,"blueviolet","blue" /)
> gsn_define_colormap(wks, colors)
>
> res = True ; plot mods desired
> ;;res_at_gsnMaximize = True ; uncomment to maximize size
> res_at_gsnSpreadColors = True ; use full range of colormap
> res_at_cnFillOn = True ; color plot desired
> res_at_cnLinesOn = False ; turn off contour lines
> res_at_cnLineLabelsOn = False ; turn off contour labels
>
> res_at_cnLevelSelectionMode = "ExplicitLevels" ; explicit [unequal] cn
> levels
> res_at_cnLevels =
> (/0,10,20.5,50,70.5,100,150,250,325,437.5,550,675,700,825,950/)
> res_at_cnFillMode = "RasterFill"
> res_at_lbOrientation = "Vertical" ; default is horizontal
> ;************************************************
> ; Use WRF_contributed procedure to set map resources
> ;************************************************
> WRF_map_c(F, res, 0) ; reads info from file
> ;************************************************
> ; set True for native projection (faster)
> ;************************************************
> res_at_tfDoNDCOverlay = True
>
> ;****************************************************************************
>
> ; create panel of different components
> ;****************************************************************************
>
> plts = new (3 , "graphic") ; 1d array to hold plots
>
> res_at_gsnDraw = False ; (a) do not draw
> res_at_gsnFrame = False ; (b) do not advance 'frame'
> res_at_lbLabelBarOn = False ; (c) turn off individual
> lb's
> ;************************************************
> ; create panel: panel plots have their own set of resources
> ;************************************************
> resP = True ; modify the panel plot
> resP_at_gsnMaximize = True ; maximize panel area
> resP_at_gsnPanelRowSpec = True ; specify 1 top, 2 lower
> level
> resP_at_gsnPanelLabelBar = True ; add common colorbar
> resP_at_pmLabelBarWidthF = 0.85 ; make label wider
> resP_at_lbLabelFontHeightF = 0.015 ; default 0.02 [demo
> make smaller]
>
> nt = 8 ; demo only one time
> do nt=0,ntim-1 ; uncomment to loop over
> all times
> plts(0) = gsn_csm_contour_map(wks,rainTot(nt,:,:),res)
> plts(1) = gsn_csm_contour_map(wks,rainnc(nt,:,:),res)
> plts(2) = gsn_csm_contour_map(wks,rainc(nt,:,:),res)
>
> resP_at_txString = f_at_TITLE+": "+times(nt)
> gsn_panel(wks,plts,(/1,2/),resP) ; now draw as one plot
> end do
> end
> --------------------------------------------------------
>
> Any help in this regards
>
>
>
>
>

-- 
--------------------------------------------------------------
Adam Phillips			             asphilli_at_ucar.edu
National Center for Atmospheric Research   tel: (303) 497-1726
ESSL/CGD/CAS                               fax: (303) 497-1333
P.O. Box 3000				
Boulder, CO 80307-3000	  http://www.cgd.ucar.edu/cas/asphilli
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk