爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 13332|回复: 7

[源代码] 用python制作维基路径图(附源代码文件及底图下载)

[复制链接]

新浪微博达人勋

发表于 2015-3-29 11:17:09 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 yangzhe1997 于 2015-3-29 11:19 编辑

转自台风论坛
python初学者,前几日了解到了强大的matplotlib模块以及其旗下的Basemap模块,花三个小时写了制作维基路径图的程序,通过导入B-Deck文件(如JTWC的BST文件)自动导出

                               
登录/注册后可看大图
欢迎交流!效果如下:
1.png

代码
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from os import startfile
         
#由风速返回SSHS级别
def category(wind):
    if wind > 137:
        cat = 'CAT5'
    elif wind > 114:
        cat = 'CAT4'
    elif wind > 96:
        cat = 'CAT3'
    elif wind > 83:
        cat = 'CAT2'
    elif wind > 64:
        cat = 'CAT1'
    elif wind > 30:
        cat = 'TS'
    else:
        cat = 'TD'
    return cat
         
#根据路径经纬度范围计算图片经纬度范围,宽长比固定为0.618
def bestscale(latmax,latmin,lonmax,lonmin):
    latmid = latmax/2 + latmin/2
    lonmid = lonmax/2 + lonmin/2
    deltalat = latmax - latmin
    deltalon = lonmax - lonmin
    if deltalat / deltalon > 0.618:
        deltalon = deltalat / 0.618
        lonmax = lonmid + deltalon / 2
        lonmin = lonmid - deltalon / 2
    elif deltalat / deltalon < 0.618:
        deltalat = deltalon * 0.618
        latmax = latmid + deltalat / 2
        latmin = latmid - deltalat / 2
    return latmax, latmin, lonmax, lonmin
         
#由级别返回颜色值,均从维基路径图上取色以保持一致
def fitcolor(cat):
    if cat == 'CAT5':
        return '#ff5e5e'
    elif cat == 'CAT4':
        return '#ff8f1d'
    elif cat == 'CAT3':
        return '#ffbe49'
    elif cat == 'CAT2':
        return '#ffe473'
    elif cat == 'CAT1':
        return '#ffffc7'
    elif cat == 'TS':
        return '#03f8f9'
    else:
        return '#59bbff'
         
#由级别返回符号标记,'^'代表三角形,'o'代表圆形
def fitmarker(cat):
    if cat == 'EX' or cat == 'DB':
        return '^'
    else:
        return 'o'
         
#主程序
def plot():
    stormnum = input('输入台风六位编号,如201201:')
    filepath = 'C:\\Users\\crazyapril\\Desktop\\BST\\%d\\WP%02d.txt' % (stormnum/100,stormnum%100)
    #这个路径是楼主电脑里BST存放路径,可按情况更改
    filehandle = open(filepath,'r')
    bstline = filehandle.readline() #打开BST文件并逐行读取,现读取第一行
    latmax, latmin, lonmax, lonmin = 0, 90, 0, 180
    latlist, lonlist, catlist = [], [], [] #纬度列表,经度列表和强度级别列表
    while bstline: #当此行不为空,执行循环
        lat, lon = int(bstline[35:38])/10.0, int(bstline[41:45])/10.0
        wind = int(bstline[48:51])
        level = bstline[59:61] #按老J BST文件格式读取经纬度、风力和级别
        latlist.append(lat)
        lonlist.append(lon) #将读取到的经纬度添加到列表中
        if wind < 64 and stormnum > 200300:
            catlist.append(level) #2003年后的老J BST文件有level信息,包含扰动和温气阶段,将其存入级别列表
        else:
            catlist.append(category(wind)) #其他均按category函数转换得到级别,存入列表
        if lat > latmax:
            latmax = lat
        if lat < latmin:
            latmin = lat
        if lon > lonmax:
            lonmax = lon
        if lon < lonmin:
            lonmin = lon #获取路径范围(最大经度、最小经度、最大纬度、最小纬度)
        bstline = filehandle.readline() #读取下一行
    filehandle.close() #文件读取完毕,关闭文件
    latmax, latmin, lonmax, lonmin = bestscale(latmax+5,latmin-5,lonmax+5,lonmin-5)
    #获取图片经纬度范围,+5、-5意为在图片中预留5纬度空间
    m = Basemap(projection='merc',resolution=None,
                llcrnrlat=latmin,urcrnrlat=latmax,
                llcrnrlon=lonmin,urcrnrlon=lonmax)
    #建立地图,采用墨卡托投影,地图四角经纬度使用上面获得的经纬度
    m.warpimage(image='E:\\nasa.jpg',scale=0.6) #调用背景图片,压缩率为0.6
    xlist, ylist = m(lonlist,latlist) #将经纬度转换为图片所用坐标
    m.plot(xlist,ylist,color='w') #首先画白线
    for i in range(len(lonlist)): #再逐个描点
        x, y = m(lonlist<i>, latlist<i>)
        m.plot(x,y,color=fitcolor(catlist<i>),marker=fitmarker(catlist<i>),markersize=6,markeredgecolor='none')
    figpath = 'C:\\Users\\crazyapril\\Desktop\\BST\\%d\\WP%02d.png' % (stormnum/100,stormnum%100) #图片保存地址
    plt.savefig(figpath,dpi=200,bbox_inches='tight',edgecolor='none',pad_inches=0) #将图片保存,dpi为200
    print '保存成功。保存地址:%s' % (figpath)
    startfile(figpath) #调用系统默认程序打开路径图片
         
plot()</i></i></i></i>


评分

参与人数 1金钱 +20 收起 理由
po_po1 + 20

查看全部评分

密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2015-3-29 11:21:12 | 显示全部楼层
另外这个安装所需的库matplotlib和basemap ,这两个库还需要另外一些库的支持,在上面这个链接里继续找到numpy、dateutil、six、pyparsing、Pillow这五个库安装即可 。
底图(代码中的nasa.jpg)在这里

                               
登录/注册后可看大图
点Original File右键保存。
JTWC Best Track:http://www.usno.navy.mil/NOOC/nmfc-ph/RSS/jtwc/best_tracks/
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-5-1 21:45:26 | 显示全部楼层
谢谢分享  
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-6-3 19:06:00 | 显示全部楼层
谢谢分享  
青春版
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-8-21 15:55:04 | 显示全部楼层
感谢分享,就是有点不太明白从经纬度到图片位置转化的那段代码。。。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-9-2 13:42:22 | 显示全部楼层
学写了,谢谢
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-10-25 21:28:57 | 显示全部楼层
好东西,学习了。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-5-3 23:16:57 | 显示全部楼层
好东西,学习了
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表