[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 08:29:54 MST 2016
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/9f84ebef/attachment-0001.html
More information about the ncl-talk
mailing list