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

Ruksana Abedin ruksana.abedin at gmail.com
Wed Dec 7 09:07:00 MST 2016


Dear Adam,

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.

Ruksana

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/Applications/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
>>>>
>>>>
>>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20161207/97998386/attachment.html 


More information about the ncl-talk mailing list