[ncl-talk] need help to understand how to change the scales in a spectrum plot

Ruksana Abedin ruksana.abedin at gmail.com
Thu Dec 8 09:43:34 MST 2016


Dear Adam and Rick,

I tried to plot the spectrums one at a time and then all three together
with the following script, and this time I could have my combined plot.
Thank you very much for your kind support and encouragement.


Best regards,
Ruksana

;*********************************
;
;**********************************
load "/usr/share/ncarg/nclscripts/csm/gsn_code.ncl"
load "/usr/share/ncarg/nclscripts/csm/gsn_csm.ncl"
load "/usr/share/ncarg/nclscripts/csm/contributed.ncl"
load "/usr/share/ncarg/nclscripts/csm/shea_util.ncl"
load "/usr/share/ncarg/nclscripts/csm/popRemap.ncl"
load "/usr/share/ncarg/nclscripts/contrib/calendar_decode2.ncl"
;************************************************
begin
;*****************************************************
;read in data
;*****************************************************
f0 = addfile
("/.../1WAH_daily_precip_b387clim_Bangladesh_100_1988_2007_land_only.nc",
"r");
f1 = addfile(".../APHRO_1988_2007_correct_daily_precip.nc", "r")
f2 = addfile("/.../GPCC_1988_2007_correct_daily_precip.nc","r")
;*****************************************************
; parameters
;*****************************************************

PPT_mod = f0->precip(:,0,:,:)
PPT_mod = PPT_mod*60*60*24
PPT_obs1 = f1->precip
PPT_obs2 = f2->precip

fldmean0 = wgt_areaave_Wrap(PPT_mod,1.0,1.0,1)
fldmean1 = wgt_areaave_Wrap(PPT_obs1,1.0,1.0,1)
fldmean2= wgt_areaave_Wrap(PPT_obs2,1.0,1.0,1)

printVarSummary(PPT_mod)
printVarSummary(PPT_obs1)
printVarSummary(PPT_obs2)

;************************************************
; set function arguments
;************************************************
  d   = 0    ; detrending opt: 0=>remove mean 1=>remove mean + detrend
  sm  = 21   ; smooth: should be at least 3 and odd
  pct = 0.10 ; percent taper: (0.0 <= pct <= 1.0) 0.10 common.
;************************************************
; calculate spectrum
;************************************************
  sdof0 = specx_anal(fldmean0,d,sm,pct)
  sdof1 = specx_anal(fldmean1,d,sm,pct)
  sdof2 = specx_anal(fldmean2,d,sm,pct)
;************************************************
; calculate confidence interval [here 5 and 95%]
; return 4 curves to be plotted
;************************************************
  splt0 = specx_ci (sdof0, 0.05, 0.95)
  splt1 = specx_ci (sdof1, 0.05, 0.95)
  splt2 = specx_ci (sdof2, 0.05, 0.95)
;************************************************
; plotting
;************************************************
   wks  = gsn_open_wks("eps","spec")              ; Opens a eps file

   r = True                                      ; plot mods desired
   r at tiMainString = "Daily Precipitation 1988-2007"
; title
   r at tiXAxisString = "Frequency (cycles/day)"  ; xaxis
   r at tiYAxisString = "Variance"                  ; yaxis
;***********************************************
; Generate log plot showing "red noise" confidence bounds
; (a) log scaling and (b) the Band Width
;***********************************************
   r at trYLog              = True                 ; log scaling
   r at trYMinF             = 0.0010                 ; manually set lower limit
   r at trYMaxF             = 360.0                 ;   "          upper
   r at gsnFrame            = False                ; do not advance frame
   r at gsnDraw = False
   r at gsnFrame = False

   r at xyLineColor = "black"
   plot0  = gsn_csm_xy(wks,sdof0 at frq, splt0,r)
   r at xyLineColor = "blue"
   plot1  = gsn_csm_xy(wks,sdof1 at frq, splt1,r)
   r at xyLineColor = "red"
   plot2  = gsn_csm_xy(wks,sdof2 at frq, splt2,r)

   xf0   = (/0.50, 0.50+sdof0 at bw/)                ; Create band width line
   ys0   = 0.75*max(sdof0 at spcx)                   ; 75% up Y axis
   yv0   = (/ys0,ys0/)
   rpl0  = True                                  ; resources for polyline
   rpl0 at gsLineThicknessF  = 2                    ; Define line thickness
   gsn_polyline(wks,plot0,xf0,yv0,rpl0)             ; Draw BandWidth

   xf1   = (/0.50, 0.50+sdof1 at bw/)                ; Create band width line
   ys1   = 0.75*max(sdof1 at spcx)                   ; 75% up Y axis
   yv1   = (/ys1,ys1/)
   rpl1  = True                                  ; resources for polyline
   rpl1 at gsLineThicknessF  = 2                    ; Define line thickness
   gsn_polyline(wks,plot1,xf1,yv1,rpl1)             ; Draw BandWidth

   xf2   = (/0.50, 0.50+sdof2 at bw/)                ; Create band width line
   ys2   = 0.75*max(sdof2 at spcx)                   ; 75% up Y axis
   yv2   = (/ys2,ys2/)
   rpl2  = True                                  ; resources for polyline
   rpl2 at gsLineThicknessF  = 2                    ; Define line thickness
   gsn_polyline(wks,plot2,xf2,yv2,rpl2)             ; Draw BandWidth

   overlay(plot0,plot1)
   overlay(plot0,plot2)

   draw(plot0)

   frame (wks)
end




On Thu, Dec 8, 2016 at 12:44 PM, Ruksana Abedin <ruksana.abedin at gmail.com>
wrote:

> Dear Adam and Rick,
>
>
> Sorry for being late to respond.
>
>
> Rick, in the addfiles commands, I have only removed the actual path
> details and put ... the 3 "."s.
>
>
> Here I am giving the variable summary:
>
>
> *Variable: PPT_mod*
> Type: float
> Total Size: 25056000 bytes
>             6264000 values
> Number of Dimensions: 3
> Dimensions and sizes:   [time | 7200] x [lat | 30] x [lon | 29]
> Coordinates:
>             time: [52950.5..60149.5]
>             lat: [15.25..29.75]
>             lon: [  83..  97]
> Number Of Attributes: 10
>   z1 :   0
>   standard_name :       precipitation_flux
>   long_name :   TOTAL PRECIPITATION RATE KG/M2/S
>   units :       kg m-2 s-1
>   _FillValue :  1e+20
>   missing_value :       1e+20
>   stash_item :  216
>   stash_section :       5
>   field_code :  90
>   cell_method : time: mean
>
> *Variable: PPT_obs1*
> Type: float
> Total Size: 25056000 bytes
>             6264000 values
> Number of Dimensions: 3
> Dimensions and sizes:   [time | 7200] x [lat | 30] x [lon | 29]
> Coordinates:
>             time: [19880101.5..20071230.5]
>             lat: [15.25..29.75]
>             lon: [  83..  97]
> Number Of Attributes: 5
>   long_name :   daily precipitation analysis interpolated onto 0.5deg
> grids [mm/day]
>   units :       mm/day
>   _FillValue :  -99.9
>   missing_value :       -99.9
>   level_description :   Earth surface
>
>
> *Variable: PPT_obs2*Type: float
> Total Size: 25056000 bytes
>             6264000 values
> Number of Dimensions: 3
> Dimensions and sizes:   [time | 7200] x [lat | 30] x [lon | 29]
> Coordinates:
>             time: [19880101.5..20071230.5]
>             lat: [15.25..29.75]
>             lon: [  83..  97]
> Number Of Attributes: 5
>   long_name :   full data daily product version 1 precipitation per grid
>   units :       mm/day
>   code :        20
>   _FillValue :  -99999.99
>   missing_value :       -99999.99
>
>
>
> The separate spectrums that I could get earlier and the script are
> attached.
>
> Thank you very much for your kind efforts in helping me.
>
> Best regards,
>
> Ruksana
>
>
>
>
> On Thu, Dec 8, 2016 at 3:58 AM, Adam Phillips <asphilli at ucar.edu> wrote:
>
>> Hi Ruksana,
>> I strongly recommend that you use print and printVarSummary calls to help
>> you debug what is going on with your script. (You should always attempt to
>> debug your own script before emailing/replying to ncl-talk.) Fatal errors
>> can sometimes result from a reference to a variable that is not present.
>>
>> Please add the following line to your script right before the gsn_csm_xy
>> call:
>> printVarSummary(spec0)
>> plot0 = gsn_csm_xy(wks,spec0 at frq,spec0 at spcx,res)
>>
>> and send ncl-talk the output of the printVarSummary statement along with
>> any other error messages you are receiving. Also, please attach your
>> current script to your reply. As always, please respond to ncl-talk and not
>> to me personally.
>> Adam
>>
>>
>>
>> On Wed, Dec 7, 2016 at 9:08 AM, Ruksana Abedin <ruksana.abedin at gmail.com>
>> wrote:
>>
>>> Dear Rick,
>>>
>>> Added the 'end' at the end of the script (sorry for this silly mistake)
>>> but nothing changes. Still getting just the same error message. No other
>>> error message. Feeling confused now.
>>>
>>> On Wed, Dec 7, 2016 at 3:49 PM, Rick Brownrigg <brownrig at ucar.edu>
>>> wrote:
>>>
>>>> Hmmm....line 121 is not the line cited...in fact its the last line of
>>>> the script.  There's an opening "begin" at the top of the script -- perhaps
>>>> match it with a closing "end" at the end of the script?
>>>>
>>>>
>>>>
>>>> On Wed, Dec 7, 2016 at 8:29 AM, Ruksana Abedin <
>>>> ruksana.abedin at gmail.com> wrote:
>>>>
>>>>> Dear Rick,
>>>>>
>>>>> Thank you for your response. I have changed accordingly but getting
>>>>> the same error message:
>>>>>
>>>>> *f*
>>>>> *atal:syntax error: line 121 in file spectrum.ncl before or near plot0
>>>>> = gsn_csm_xy(wks,spec0 at frq,spec0 at spcx,res)*
>>>>>
>>>>> Any other ideas?
>>>>>
>>>>> Best regards,
>>>>> Ruksana
>>>>>
>>>>> On Wed, Dec 7, 2016 at 2:44 PM, Rick Brownrigg <brownrig at ucar.edu>
>>>>> wrote:
>>>>>
>>>>>> Hi Ruksana,
>>>>>>
>>>>>> I'm not entirely sure what the issue is, but by my reading of the
>>>>>> docs on the function specx_anal(), it returns a variable with certain
>>>>>> attributes attached, one of which is an array named "spcx". Thus, I think
>>>>>> all your references such as:
>>>>>>
>>>>>> spec0 at spcx0
>>>>>> spec1 at spcx1
>>>>>> spec2 at spcx2
>>>>>>
>>>>>> should be:
>>>>>>
>>>>>> spec0 at spcx
>>>>>> spec1 at spcx
>>>>>> spec2 at spcx
>>>>>>
>>>>>> There are a number of such references -- I would think there would
>>>>>> have been other warnings/errors elsewhere in your script (?)
>>>>>>
>>>>>> Hope that helps...
>>>>>> Rick
>>>>>>
>>>>>>
>>>>>> On Wed, Dec 7, 2016 at 5:36 AM, Ruksana Abedin <
>>>>>> ruksana.abedin at gmail.com> wrote:
>>>>>>
>>>>>>> Dear Adam,
>>>>>>>
>>>>>>> Thank you very much for suggesting about the way to have three
>>>>>>> spectrums in one plot. I have tried the following script but getting error.
>>>>>>>
>>>>>>> ;*********************************
>>>>>>> ;
>>>>>>> ;**********************************
>>>>>>> load "/usr/share/ncarg/nclscripts/csm/gsn_code.ncl"
>>>>>>> load "/usr/share/ncarg/nclscripts/csm/gsn_csm.ncl"
>>>>>>> load "/usr/share/ncarg/nclscripts/csm/contributed.ncl"
>>>>>>> load "/usr/share/ncarg/nclscripts/csm/shea_util.ncl"
>>>>>>> load "/usr/share/ncarg/nclscripts/csm/popRemap.ncl"
>>>>>>> load "/usr/share/ncarg/nclscripts/contrib/calendar_decode2.ncl"
>>>>>>> ;************************************************
>>>>>>> begin
>>>>>>> ;*****************************************************
>>>>>>> ;read in data
>>>>>>> ;*****************************************************
>>>>>>> f0 = addfile ("/../1WAH_daily_precip_b387clim_Bangladesh_100_1988_
>>>>>>> 2007_land_only.nc", "r");
>>>>>>> f1 = addfile("/.../APHRO_1988_2007_correct_daily_precip.nc", "r")
>>>>>>> f2 = addfile("/.../GPCC_1988_2007_correct_daily_precip.nc","r")
>>>>>>> ;*****************************************************
>>>>>>> ; parameters
>>>>>>> ;*****************************************************
>>>>>>> PPT_mod = f0->precip(:,0,:,:)
>>>>>>> PPT_mod = PPT_mod*60*60*24
>>>>>>> PPT_obs1 = f1->precip
>>>>>>> PPT_obs2 = f2->precip
>>>>>>> fldmean0 = wgt_areaave_Wrap(PPT_mod,1.0,1.0,1)
>>>>>>> fldmean1 = wgt_areaave_Wrap(PPT_obs1,1.0,1.0,1)
>>>>>>> fldmean2= wgt_areaave_Wrap(PPT_obs2,1.0,1.0,1)
>>>>>>> ;printVarSummary(time)
>>>>>>> printVarSummary(PPT_mod)
>>>>>>> printVarSummary(PPT_obs1)
>>>>>>> printVarSummary(PPT_obs2)
>>>>>>> ;************************************************
>>>>>>> ; set function arguments
>>>>>>> ;************************************************
>>>>>>> ; detrending opt: 0=>remove mean 1=>remove mean and detrend
>>>>>>>   d = 0
>>>>>>> ; smoothing periodogram: (0 <= sm <= ??.) should be at least 3 and
>>>>>>> odd
>>>>>>>   sm = 7
>>>>>>> ; percent tapered: (0.0 <= pct <= 1.0) 0.10 common.
>>>>>>>   pct = 0.10
>>>>>>> ;************************************************
>>>>>>> ; calculate spectrum
>>>>>>> ;************************************************
>>>>>>>   spec0 = specx_anal(fldmean0,d,sm,pct)
>>>>>>>   spec1 = specx_anal(fldmean1,d,sm,pct)
>>>>>>>   spec2 = specx_anal(fldmean2,d,sm,pct)
>>>>>>> ;************************************************
>>>>>>> ; calculate confidence interval [here 5 and 95%]
>>>>>>> ; return 4 curves to be plotted
>>>>>>> ;************************************************
>>>>>>>   splt0 = specx_ci (spec0, 0.05, 0.95)
>>>>>>>   splt1 = specx_ci (spec1, 0.05, 0.95)
>>>>>>>   splt2 = specx_ci (spec2, 0.05, 0.95)
>>>>>>> ;************************************************
>>>>>>> ;************************************************
>>>>>>> ; plotting
>>>>>>> ;************************************************
>>>>>>>    wks  = gsn_open_wks("eps","spec_combined")                ;
>>>>>>> Opens a eps file
>>>>>>>    res = True        ; no plot mods desired
>>>>>>>    res at gsnDraw = False
>>>>>>>    res at gsnFrame = False
>>>>>>>    res at tiMainString = "Precipitation"             ; title
>>>>>>>    res at tiXAxisString = "Frequency (cycles/day)"  ; xaxis
>>>>>>>    res at tiYAxisString = "Variance"                  ; yaxis
>>>>>>>
>>>>>>> ;***********************************************
>>>>>>> ; Generate log plot showing "red noise" confidence bounds
>>>>>>> ; (a) log scaling and (b) the Band Width
>>>>>>> ;***********************************************
>>>>>>>    res at trYLog              = True                 ; log scaling
>>>>>>>    res at trYMinF             = 0.00010                 ; manually set
>>>>>>> lower limit
>>>>>>>    res at trYMaxF             = 350.0                ;   "    upper
>>>>>>>    res at gsnFrame            = False                ; do not advance
>>>>>>> frame
>>>>>>>
>>>>>>>    xf0   = (/0.75, 0.75+spec0 at bw/)                ; Create band
>>>>>>> width line
>>>>>>>    ys0   = 0.75*max(spec0 at spcx0)                   ; 75% up Y axis
>>>>>>>    yv0   = (/ys0,ys0/)
>>>>>>>    respl  = True                                  ; resources for
>>>>>>> polyline
>>>>>>>    respl at gsLineThicknessF  = 2                    ; Define line
>>>>>>> thickness
>>>>>>>
>>>>>>>    txres= True                                  ; label BW line
>>>>>>>    txres at txFontHeightF = 0.015                  ; font height
>>>>>>>    txres at txJust        = "CenterLeft"           ; Set lable location
>>>>>>>    gsn_text(wks,plot,"BW",0.41+spec0 at bw,ys0,txres); Label
>>>>>>>
>>>>>>>    xf1   = (/0.75, 0.75+spec1 at bw/)                ; Create band
>>>>>>> width line
>>>>>>>    ys1   = 0.75*max(spec1 at spcx1)                   ; 75% up Y axis
>>>>>>>    yv1   = (/ys1,ys1/)
>>>>>>>    respl  = True                                  ; resources for
>>>>>>> polyline
>>>>>>>    respl at gsLineThicknessF  = 2                    ; Define line
>>>>>>> thickness
>>>>>>>
>>>>>>>    txres= True                                  ; label BW line
>>>>>>>    txres at txFontHeightF = 0.015                  ; font height
>>>>>>>    txres at txJust        = "CenterLeft"           ; Set lable location
>>>>>>>    gsn_text(wks,plot,"BW",0.41+spec1 at bw,ys1,txres); Label
>>>>>>>
>>>>>>>    xf2   = (/0.75, 0.75+spec2 at bw/)                ; Create band
>>>>>>> width line
>>>>>>>    ys2   = 0.75*max(spec2 at spcx2)                   ; 75% up Y axis
>>>>>>>    yv2   = (/ys2,ys2/)
>>>>>>>    respl  = True                                  ; resources for
>>>>>>> polyline
>>>>>>>    respl at gsLineThicknessF  = 2                    ; Define line
>>>>>>> thickness
>>>>>>>
>>>>>>>    txres= True                                  ; label BW line
>>>>>>>    txres at txFontHeightF = 0.015                  ; font height
>>>>>>>    txres at txJust        = "CenterLeft"           ; Set lable location
>>>>>>>    gsn_text(wks,plot,"BW",0.41+spec2 at bw,ys2,txres); Label
>>>>>>>
>>>>>>> res at gsnDraw = False
>>>>>>> res at gsnFrame = False
>>>>>>> res at xyLineColor = "black"
>>>>>>> plot0 = gsn_csm_xy(wks,spec0 at frq,spec0 at spcx0,res)     ; create
>>>>>>> plot
>>>>>>> gsn_polyline(wks,plot,xf0,yv0,respl)             ; Draw BandWidth
>>>>>>> res at xyLineColor = "blue"
>>>>>>> plot1 = gsn_csm_xy(wks,spec1 at frq,spec1 at spcx1,res)     ; create plot
>>>>>>>
>>>>>>> gsn_polyline(wks,plot,xf1,yv1,respl)             ; Draw BandWidth
>>>>>>> res at xyLineColor = "red"
>>>>>>> plot2 = gsn_csm_xy(wks,spec2 at frq,spec2 at spcx2,res)     ; create plot
>>>>>>> gsn_polyline(wks,plot,xf2,yv2,respl)             ; Draw BandWidth
>>>>>>> overlay(plot0,plot1)
>>>>>>> overlay(plot0,plot2)
>>>>>>> draw(plot)
>>>>>>> frame(wks)
>>>>>>>
>>>>>>>
>>>>>>> *fatal:syntax error: line 121 in file spectrum.ncl before or near
>>>>>>> plot0 = gsn_csm_xy(wks,spec0 at frq,spec0 at spcx0,res)*
>>>>>>>
>>>>>>> *Can you please help me to identify where my mistake is?*
>>>>>>>
>>>>>>> *Best regards,*
>>>>>>> *Ruksana*
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> On Mon, Dec 5, 2016 at 7:16 PM, Adam Phillips <asphilli at ucar.edu>
>>>>>>> wrote:
>>>>>>>
>>>>>>>> Hi Ruksana,
>>>>>>>> If all you want is 3 different spectra overlaid on the same plot
>>>>>>>> you could use the overlay procedure. You will need to tell NCL to not draw
>>>>>>>> the plot nor advance the frame before creating the 3 spectra. You may also
>>>>>>>> need to set trYMaxF if portions of one spectra go off the plot.
>>>>>>>>
>>>>>>>> Here's an example:
>>>>>>>> res at gsnDraw = False
>>>>>>>> res at gsnFrame = False
>>>>>>>> res at xyLineColor = "green"
>>>>>>>> plot=gsn_csm_xy(wks,spec at frq,spec at spcx,res)     ; create plot
>>>>>>>> res at xyLineColor = "black"
>>>>>>>> plot2 = gsn_csm_xy(wks,spec2 at frq,spec2 at spcx,res)     ; create plot
>>>>>>>>
>>>>>>>> res at xyLineColor = "blue"
>>>>>>>> plot3 = gsn_csm_xy(wks,spec3 at frq,spec3 at spcx,res)     ; create plot
>>>>>>>> overlay(plot,plot2)
>>>>>>>> overlay(plot,plot3)
>>>>>>>> draw(plot)
>>>>>>>> frame(wks)
>>>>>>>>
>>>>>>>> Hope that helps. If not, please respond to the ncl-talk email list.
>>>>>>>> Adam
>>>>>>>>
>>>>>>>> On Mon, Dec 5, 2016 at 12:09 PM, Ruksana Abedin <
>>>>>>>> ruksana.abedin at gmail.com> wrote:
>>>>>>>>
>>>>>>>>> Hi,
>>>>>>>>>
>>>>>>>>> I could figure the required process out. This time followed the
>>>>>>>>> example 4 at http://www.ncl.ucar.edu/Applic
>>>>>>>>> ations/Scripts/spec_4.ncl. However, I would still like to know
>>>>>>>>> how can I plot three spectrums of three colours in one single plot. Your
>>>>>>>>> help is always useful.
>>>>>>>>>
>>>>>>>>> Thank you for putting up wonderful resources.
>>>>>>>>>
>>>>>>>>> Best regards,
>>>>>>>>> Ruksana
>>>>>>>>>
>>>>>>>>> On Mon, Dec 5, 2016 at 5:30 PM, Ruksana Abedin <
>>>>>>>>> ruksana.abedin at gmail.com> wrote:
>>>>>>>>>
>>>>>>>>>> Hi,
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> I have followed the http://www.ncl.ucar.edu/Applic
>>>>>>>>>> ations/Scripts/spec_1.ncl script to plot a spectrum for my data
>>>>>>>>>> but the specs are very difficult to visualize clearly in the plot. Data is
>>>>>>>>>> daily precipitation from model and two observations for 1988-2007.
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> 1. How to change the x and y axis scales to make it more visible.
>>>>>>>>>>
>>>>>>>>>> 2. Can I have the cycles per day instead of months? How to change
>>>>>>>>>> this aspect?
>>>>>>>>>>
>>>>>>>>>> 3. Can I plot three spectrums for three different files in a
>>>>>>>>>> panel together?
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> The variable summary, script and the resultant plot are attached.
>>>>>>>>>> Your time and efforts to help me out is highly appreciated.
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> Variable: PPT_mod
>>>>>>>>>> Type: float
>>>>>>>>>> Total Size: 25056000 bytes
>>>>>>>>>>             6264000 values
>>>>>>>>>> Number of Dimensions: 3
>>>>>>>>>> Dimensions and sizes:   [time | 7200] x [lat | 30] x [lon | 29]
>>>>>>>>>> Coordinates:
>>>>>>>>>>             time: [52950.5..60149.5]
>>>>>>>>>>             lat: [15.25..29.75]
>>>>>>>>>>             lon: [  83..  97]
>>>>>>>>>> Number Of Attributes: 10
>>>>>>>>>>   z1 :   0
>>>>>>>>>>   standard_name :       precipitation_flux
>>>>>>>>>>   long_name :   TOTAL PRECIPITATION RATE KG/M2/S
>>>>>>>>>>   units :       kg m-2 s-1
>>>>>>>>>>   _FillValue :  1e+20
>>>>>>>>>>   missing_value :       1e+20
>>>>>>>>>>   stash_item :  216
>>>>>>>>>>   stash_section :       5
>>>>>>>>>>   field_code :  90
>>>>>>>>>>   cell_method : time: mean
>>>>>>>>>>
>>>>>>>>>> Variable: PPT_obs1
>>>>>>>>>> Type: float
>>>>>>>>>> Total Size: 25056000 bytes
>>>>>>>>>>             6264000 values
>>>>>>>>>> Number of Dimensions: 3
>>>>>>>>>> Dimensions and sizes:   [time | 7200] x [lat | 30] x [lon | 29]
>>>>>>>>>> Coordinates:
>>>>>>>>>>             time: [19880101.5..20071230.5]
>>>>>>>>>>             lat: [15.25..29.75]
>>>>>>>>>>             lon: [  83..  97]
>>>>>>>>>> Number Of Attributes: 5
>>>>>>>>>>   long_name :   daily precipitation analysis interpolated onto
>>>>>>>>>> 0.5deg grids [mm/day]
>>>>>>>>>>   units :       mm/day
>>>>>>>>>>   _FillValue :  -99.9
>>>>>>>>>>   missing_value :       -99.9
>>>>>>>>>>   level_description :   Earth surface
>>>>>>>>>>
>>>>>>>>>> Variable: PPT_obs2
>>>>>>>>>> Type: float
>>>>>>>>>> Total Size: 25056000 bytes
>>>>>>>>>>             6264000 values
>>>>>>>>>> Number of Dimensions: 3
>>>>>>>>>> Dimensions and sizes:   [time | 7200] x [lat | 30] x [lon | 29]
>>>>>>>>>> Coordinates:
>>>>>>>>>>             time: [19880101.5..20071230.5]
>>>>>>>>>>             lat: [15.25..29.75]
>>>>>>>>>>             lon: [  83..  97]
>>>>>>>>>> Number Of Attributes: 5
>>>>>>>>>>   long_name :   full data daily product version 1 precipitation
>>>>>>>>>> per grid
>>>>>>>>>>   units :       mm/day
>>>>>>>>>>   code :        20
>>>>>>>>>>   _FillValue :  -99999.99
>>>>>>>>>>   missing_value :       -99999.99
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> ;************************************************
>>>>>>>>>> begin
>>>>>>>>>> ;*****************************************************
>>>>>>>>>> ;read in data
>>>>>>>>>> ;*****************************************************
>>>>>>>>>> f0 = addfile ("/../1WAH_daily_precip_b387cl
>>>>>>>>>> im_Bangladesh_100_1988_2007_land_only.nc", "r");
>>>>>>>>>> f1 = addfile("/../APHRO_1988_2007_correct_daily_precip.nc", "r")
>>>>>>>>>> f2 = addfile("/../GPCC_1988_2007_correct_daily_precip.nc","r")
>>>>>>>>>>
>>>>>>>>>> ;*****************************************************
>>>>>>>>>> ; parameters
>>>>>>>>>> ;*****************************************************
>>>>>>>>>>
>>>>>>>>>> PPT_mod = f0->precip(:,0,:,:)
>>>>>>>>>> PPT_mod = PPT_mod*60*60*24
>>>>>>>>>> PPT_obs1 = f1->precip
>>>>>>>>>> PPT_obs2 = f2->precip
>>>>>>>>>>
>>>>>>>>>> fldmean0 = wgt_areaave_Wrap(PPT_mod,1.0,1.0,1)
>>>>>>>>>> fldmean1 = wgt_areaave_Wrap(PPT_obs1,1.0,1.0,1)
>>>>>>>>>> fldmean2= wgt_areaave_Wrap(PPT_obs2,1.0,1.0,1)
>>>>>>>>>>
>>>>>>>>>> ;printVarSummary(time)
>>>>>>>>>>
>>>>>>>>>> printVarSummary(PPT_mod)
>>>>>>>>>> printVarSummary(PPT_obs1)
>>>>>>>>>> printVarSummary(PPT_obs2)
>>>>>>>>>>
>>>>>>>>>> ;************************************************
>>>>>>>>>> ; set function arguments
>>>>>>>>>> ;************************************************
>>>>>>>>>> ; detrending opt: 0=>remove mean 1=>remove mean and detrend
>>>>>>>>>>   d = 0
>>>>>>>>>> ; smoothing periodogram: (0 <= sm <= ??.) should be at least 3
>>>>>>>>>> and odd
>>>>>>>>>>   sm = 7
>>>>>>>>>> ; percent tapered: (0.0 <= pct <= 1.0) 0.10 common.
>>>>>>>>>>   pct = 0.10
>>>>>>>>>> ;************************************************
>>>>>>>>>> ; calculate spectrum
>>>>>>>>>> ;************************************************
>>>>>>>>>>   spec = specx_anal(fldmean2,d,sm,pct)
>>>>>>>>>> ;************************************************
>>>>>>>>>> ; plotting
>>>>>>>>>> ;************************************************
>>>>>>>>>>    wks  = gsn_open_wks("eps","spec")                ; Opens a
>>>>>>>>>> eps file
>>>>>>>>>>
>>>>>>>>>>    res = True                       ; no plot mods desired
>>>>>>>>>>    res at tiMainString = "Precipitation"                   ; title
>>>>>>>>>>    res at tiXAxisString = "Frequency (cycles/month)"  ; xaxis
>>>>>>>>>>    res at tiYAxisString = "Variance"                  ; yaxis
>>>>>>>>>>
>>>>>>>>>>    plot=gsn_csm_xy(wks,spec at frq,spec at spcx,res)     ; create
>>>>>>>>>> plot
>>>>>>>>>> ;***********************************************
>>>>>>>>>> end
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>> _______________________________________________
>>>>>>>>> ncl-talk mailing list
>>>>>>>>> ncl-talk at ucar.edu
>>>>>>>>> List instructions, subscriber options, unsubscribe:
>>>>>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>>>>>>
>>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> --
>>>>>>>> Adam Phillips
>>>>>>>> Associate Scientist,  Climate and Global Dynamics Laboratory, NCAR
>>>>>>>> www.cgd.ucar.edu/staff/asphilli/   303-497-1726
>>>>>>>>
>>>>>>>> <http://www.cgd.ucar.edu/staff/asphilli>
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> ncl-talk mailing list
>>>>>>> ncl-talk at ucar.edu
>>>>>>> List instructions, subscriber options, unsubscribe:
>>>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>>>>
>>>>>>>
>>>>>>
>>>>>
>>>>
>>>
>>
>>
>> --
>> Adam Phillips
>> Associate Scientist,  Climate and Global Dynamics Laboratory, NCAR
>> www.cgd.ucar.edu/staff/asphilli/   303-497-1726
>>
>> <http://www.cgd.ucar.edu/staff/asphilli>
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20161208/400674c8/attachment-0001.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: spec.jpg
Type: image/jpeg
Size: 422959 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20161208/400674c8/attachment-0001.jpg 


More information about the ncl-talk mailing list