[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 05:36:56 MST 2016
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>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20161207/e57cf14f/attachment.html
More information about the ncl-talk
mailing list