-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy path全球图.py
executable file
·120 lines (103 loc) · 4.33 KB
/
全球图.py
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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
# -*- coding: utf-8 -*-
"""
Created on Tue Mar 31 19:10:03 2015
@author: qqf
"""
import numpy as np
import Scientific.IO.NetCDF as S
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from scipy import stats
f=S.NetCDFFile('gpcp.precip.mon.mean.nc',mode='r')
precip=f.variables['precip'].getValue()
time=f.variables['time'].getValue()
lon=f.variables['lon'].getValue()
lat=f.variables['lat'].getValue()
f.close()
f=np.loadtxt("ersst3b.nino.mth.81-10.ascii",dtype=np.str)
Nino=f[1:,8].astype(np.float)
precip_12=precip[11:431:12,:,:]*31.0
precip_01=precip[12::12,:,:]*31.0
precip_02=precip[13::12,:,:]*28.0
tmp=(precip_12+precip_01+precip_02)/3.0
precip_djf_ave=tmp.mean(axis=0)
precip_djf_anomaly=tmp-precip_djf_ave
Nino_12=Nino[361:781:12]
Nino_01=Nino[362:781:12]
Nino_02=Nino[363:781:12]
tmp=(Nino_12+Nino_01+Nino_02)/3.0
Nino_djf_ave=tmp.mean(axis=0)
Nino_djf_anomaly=tmp-Nino_djf_ave
slope=np.empty(precip_djf_anomaly[0,:,:].shape)
intercept=np.empty(precip_djf_anomaly[0,:,:].shape)
r_value=np.empty(precip_djf_anomaly[0,:,:].shape)
p_value=np.empty(precip_djf_anomaly[0,:,:].shape)
std=np.empty(precip_djf_anomaly[0,:,:].shape)
for yy in range(len(slope[:,0])):
for xx in range(len(slope[0,:])):
slope[yy,xx],intercept[yy,xx],r_value[yy,xx],p_value[yy,xx],std[yy,xx]=stats.linregress(Nino_djf_anomaly,precip_djf_anomaly[:,yy,xx])
pre_1997=slope*Nino_djf_anomaly[1997-1979]+intercept
pre_1998=slope*Nino_djf_anomaly[1998-1979]+intercept
m=Basemap(llcrnrlon=0,llcrnrlat=-90,urcrnrlon=360,urcrnrlat=90,projection='mill')
lons,lats=np.meshgrid(lon, lat)
x,y=m(lons,lats)
fig1 = plt.figure(figsize=(16,20))
ax = fig1.add_axes([0.1,0.1,0.8,0.8])
clevs = np.linspace(-240,240,21)
CS1 = m.contour(x,y,pre_1997,clevs,linewidths=0.5,colors='k')
CS2 = m.contourf(x,y,pre_1997,clevs,cmap=plt.cm.RdBu_r)
m.drawmapboundary(fill_color='#99ffff')
m.drawcoastlines(linewidth=1.5)
m.drawparallels(np.arange(-90,90,20),labels=[1,1,0,0])
m.drawmeridians(np.arange(0.,360.,20.),labels=[0,0,0,1])
cbar = m.colorbar(CS2,location='bottom',pad="5%")
cbar.set_label('mm')
plt.title('Regression Precipitation Anomaly in 1997')
plt.savefig("./1997.png",dpi=500)
fig2 = plt.figure(figsize=(16,20))
ax = fig2.add_axes([0.1,0.1,0.8,0.8])
clevs = np.linspace(-240,240,21)
CS1 = m.contour(x,y,pre_1998,clevs,linewidths=0.5,colors='k')
CS2 = m.contourf(x,y,pre_1998,clevs,cmap=plt.cm.RdBu_r)
m.drawmapboundary(fill_color='#99ffff')
m.drawcoastlines(linewidth=1.5)
m.drawparallels(np.arange(-90,90,20),labels=[1,1,0,0])
m.drawmeridians(np.arange(0.,360.,20.),labels=[0,0,0,1])
cbar = m.colorbar(CS2,location='bottom',pad="5%")
cbar.set_label('mm')
plt.title('Regression Precipitation Anomaly in 1998')
plt.savefig("./1998.png",dpi=500)
pre_1997_test=np.ma.array(np.empty(precip_djf_anomaly[0,:,:].shape),mask=False)
pre_1998_test=np.ma.array(np.empty(precip_djf_anomaly[0,:,:].shape),mask=False)
for yy in range(len(slope[:,0])):
for xx in range(len(slope[0,:])):
pre_1997_test[yy,xx]=pre_1997[yy,xx]
pre_1998_test[yy,xx]=pre_1998[yy,xx]
if (p_value[yy,xx]>0.05): pre_1997_test.mask[yy,xx] = True
if (p_value[yy,xx]>0.05): pre_1998_test.mask[yy,xx] = True
fig3 = plt.figure(figsize=(16,20))
ax = fig3.add_axes([0.1,0.1,0.8,0.8])
clevs = np.linspace(-240,240,21)
CS1 = m.contour(x,y,pre_1997_test,clevs,linewidths=0.5,colors='k')
CS2 = m.contourf(x,y,pre_1997_test,clevs,cmap=plt.cm.RdBu_r)
m.drawmapboundary(fill_color='#99ffff')
m.drawcoastlines(linewidth=1.5)
m.drawparallels(np.arange(-90,90,20),labels=[1,1,0,0])
m.drawmeridians(np.arange(0.,360.,20.),labels=[0,0,0,1])
cbar = m.colorbar(CS2,location='bottom',pad="5%")
cbar.set_label('mm')
plt.title('Regression Precipitation Anomaly in 1997 with significance test')
plt.savefig("./mask1997.png",dpi=500)
fig4 = plt.figure(figsize=(16,20))
ax = fig4.add_axes([0.1,0.1,0.8,0.8])
clevs = np.linspace(-240,240,21)
CS1 = m.contour(x,y,pre_1998_test,clevs,linewidths=0.5,colors='k')
CS2 = m.contourf(x,y,pre_1998_test,clevs,cmap=plt.cm.RdBu_r)
m.drawmapboundary(fill_color='#99ffff')
m.drawcoastlines(linewidth=1.5)
m.drawparallels(np.arange(-90,90,20),labels=[1,1,0,0])
m.drawmeridians(np.arange(0.,360.,20.),labels=[0,0,0,1])
cbar = m.colorbar(CS2,location='bottom',pad="5%")
cbar.set_label('mm')
plt.title('Regression Precipitation Anomaly in 1998 with significance test')
plt.savefig("./mask1998.png",dpi=500)