登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 静言_GRMC 于 2017-5-16 09:21 编辑
参考了版主发的脚本示例(http://bbs.06climate.com/forum.p ... 6422&extra=page%3D1),版主的是用micaps资料,我也依葫芦画瓢做了一个,数据存在txt(经、纬度、PM2.5浓度),分享给大家,有代码多余或者不规范的请多多指教~
代码入下:
- #Get file names
- fn = 'f:/meteoinfoJS/AQI.txt'
- ncol = numasciicol(fn)
- nrow = numasciirow(fn)
- a = asciiread(fn,shape=(nrow,ncol))
- lon = a[:,0]
- lat = a[:,1]
- pm= a[:,2]
- #To grid data
- x = arange(112.9, 114.3, 0.2)
- y = arange(22.5, 24.2, 0.2)
- gtemp,gx,gy = griddata((lon, lat), pm, xi=(x, y), method='idw', radius=0.8)
- #Plot
- axesm()
- bou1_layer = shaperead('f:/MeteoInfo/map/guangzhou_county_new.shp')
- mlayer = shaperead('f:/MeteoInfo/map/guangzhou_county_new.shp')
- geoshow(bou1_layer, edgecolor='lightgray')
- geoshow(mlayer, visible=False)
- levs = [0, 25, 35, 50, 75, 95, 115, 130, 150,200]
- cols = [(255,255,255),(0,255,0),(127,255,0),(255,255,0),(255,215,0),(255,128,0),(255,97,0), \
- (255,0,0),(176,23,31),(135,38,87),(255,0,255)]
- #layer = contourfm(x, y, gtemp,20)
-
- layer = contourfm(x, y, gtemp,levs,colors=cols)
- #slayer = scatterm(lon, lat,pm,colors=['k'], size=10)
- slayer = scatterm(lon, lat,pm,levs,colors=cols, size=8)
- masklayer(mlayer, [layer])
- xlim(112.9, 114.1)
- ylim(22.5, 24)
- title(u'广州市空气质量站点PM2.5浓度',fontname=u'黑体',fontsize=20,bold=False,color='blue')
- text(113.3, 23.9, u'2017年5月6日13时', fontname=u'黑体', fontsize=16)
- colorbar(layer)
复制代码
|