[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 05:44:39 MST 2016


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_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
>>>>>>>>> ;************************************************
>>>>>>>>>   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/e8ff25f4/attachment-0001.html 
-------------- next part --------------
;*********************************
;
;**********************************
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
   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)
end

fatal:syntax error: line 121 in file spectrum.ncl before or near 
plot0 = gsn_csm_xy(wks,spec0 at frq,spec0 at spcx0,res)
-------------- next part --------------
A non-text attachment was scrubbed...
Name: spectrum figures.docx
Type: application/vnd.openxmlformats-officedocument.wordprocessingml.document
Size: 739234 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20161208/e8ff25f4/attachment-0001.bin 


More information about the ncl-talk mailing list