forked from gavin971/NCL_meteorology_libs
-
Notifications
You must be signed in to change notification settings - Fork 0
/
plot-scatter-AO.ncl
83 lines (57 loc) · 1.93 KB
/
plot-scatter-AO.ncl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
begin
itimes = 197901
itimee = 201001
siglvl = 0.01
nyear = (itimee - itimes)/100+1
xaxis = fspan(-2,2,nyear) ; x轴 与y轴的范围
yaxis = xaxis
;;;;read air data
f = addfile("./data/h300-197901-201412.nc", "r")
var := short2flt(f->hgt(:,0,{65:90},:))
time := f->time
YYYYMM := cd_calendar(time,-1)
ist = ind(itimes.eq.YYYYMM)
ied = ind(itimee.eq.YYYYMM)
hgt = var(ist:ied:12,:,:)
rad = 3.1415926/180 ;4.0*atan(1.0)/180.0
lat = hgt&lat
clat = cos(lat*rad)
vor = hgt(:,0,0)
vor = wgt_areaave(hgt(:,{65:90},:),clat,1.0,0)
vor = dim_standardize(vor,0)
;;;;read AOI data
ncfile = addfile("SLP-PC1-jan-34yr.nc","r")
eof_ts = ncfile->AOI
AOI = eof_ts({1979:2010}) ; 截取部分时间段,以保证与air数据一致
;;直线拟合
; y = mx+b
; m is the slope: returned from regline
; b is the y intercept: rc@yave attribute of rc returned from regline
;; 检验
;;**plot**************************************
wks = gsn_open_wks("eps","plot-scatter-AO") ; specifies a ps plot
res = True ; plot mods desired
res@gsnMaximize = True ; maximize plot in frame
res@gsnDraw = False
res@gsnFrame = False
if(prob.le.siglvl )
else
res@gsnCenterString = "r is " + sprintf("%5.2f",m) + "(insig. @"+siglvl+")"
end if
res@trYMinF = min(yaxis)
res@trYMaxF = max(yaxis)
res@trXMinF = min(xaxis)
res@trXMaxF = max(xaxis)
res@tiXAxisString = "AO"
res@tiYAxisString = "Vor."
plot = gsn_csm_xy(wks,xaxis,data,res) ; create plot
pmres = True
pmres@gsMarkerIndex = 12
pmres@gsMarkerSizeF = 0.02
pmres@gsMarkerColor = "black"
xl = eof_ts
yl = vor
dum = gsn_add_polymarker(wks,plot,xl,yl,pmres)
draw(plot)
frame(wks)
end