[ncl-talk] need help to understand how to change the scales in a spectrum plot
Rick Brownrigg
brownrig at ucar.edu
Wed Dec 7 08:49:25 MST 2016
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/f49856f0/attachment.html
More information about the ncl-talk
mailing list