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