- 积分
- 1336
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2015-8-21
- 最后登录
- 1970-1-1
|
发表于 2019-2-27 10:47:10
|
显示全部楼层
王老师,这是我的程序:
#Set data folders
st_datadir='H:/ymy/Precipitaion6h/2015/12h/06/'
st=datetime.datetime(2015,6,1,12)
et=datetime.datetime(2015,6,1,12)
fns=[]
#Read station data
while st<=et:
fn=os.path.join(st_datadir,'SURF_WEA_PRE_6HOUR-'+st.strftime('%Y%m%d%H')+'.txt')
ncol = numasciicol(fn)
nrow = numasciirow(fn)
a = asciiread(fn,shape=(nrow,ncol))
lon= a[:,6]
lat = a[:,5]
rain= a[:,8]
#griddata function - interpolate
x = arange(75, 135, 0.28125)
y = arange(10, 55, 0.28125)
prg= griddata((lon, lat), rain, xi=(x, y), method='idw', radius=0.5)
fns.append(fn)
st=st+datetime.timedelta(days=1)
#Plot
proj = projinfo(proj='lcc', lon_0=105, lat_1=25, lat_2=47)
axesm(projinfo=proj, position=[0, 0, 0.9, 1], axison=False, gridlabel=False, frameon=False)
geoshow('cn_province', edgecolor='lightgray')
bou1_layer = geoshow('cn_border', facecolor=(0,0,255))
geoshow('cn_cities', facecolor='r', size=4, labelfield='NAME', fontname=u'楷体', fontsize=16, yoffset=15)
china_layer = geoshow('china', visible=False)
levs = [0, 0.1, 10, 25, 50, 100]#降水量级
cols = [(255,255,255),(170,240,255),(120,230,240),(200,220,50),(240,220,20),(255,120,10),(255,90,10), \
(240,40,0),(180,10,0),(120,10,0)]
layer = contourfm(x, y, prg, levs, colors=cols)
masklayer(china_layer, [layer])
colorbar(layer, shrink=0.5, aspect=15)
axism([78, 130, 14, 53])#经纬度
text(95, 53, u'全国降水量实况图', fontname=u'黑体', fontsize=18)
text(95, 51, u'(2015-06-01 12:00 至 2015-06-02 12:00)', fontname=u'黑体', fontsize=16)
#Add south China Sea
sc_layer = bou1_layer.clone()
axesm(position=[0.1,0.05,0.15,0.2], axison=False, frameon=True)
geoshow(sc_layer, facecolor=(0,0,255))
xlim(106, 123)
ylim(2, 23)
#savefig('D:/Temp/test/rain_test.png', 800, 800) |
|