尝试绘制台风生成源地密度分布

qxsw2016

有文献中绘制了台风生成源地密度分布图（单位是个/(pi.250km.250km)）,尝试利用台风最佳路径数据集，台风生成位置取第一个风速大于等于17.2米/秒的位置。脚本如下：
import mipylib.geolib.topology as tp
import re
import math

# Read typhoon data file
fn = 'D:/Temp/test/ty/CH%dBST.txt'
fn = r'E:\MeteoInfo\Dataset\typhoon\CH%dBST.txt'
orglons = [
orglats =
pattern =re.compile(r'\(-\)')
for yy in range(1949,2024):
    datafn = fn%yy
    #print datafn
    tf = open(datafn)
    for line in tf:
        data = line.split()
        rown = int(data[2)
        match =  re.search(pattern,line)
        
        if(match):            
            for j in range(rown):
                line = tf.readline()
        else:
            #读取每个台风数据
            lons = [
            lats = [
            prss = [
            wss = [
            
            for i in range(rown):
                line = tf.readline()
                data = line.split()
                lat = float(data[2)        
                lats.append(lat * 0.1)
                lon = float(data[3)
                lons.append(lon * 0.1)
                prs = float(data[4)
                prss.append(prs)
                ws = int(data[5)
                wss.append(ws)

            idx = -1
            for i in range(rown):
                idx=idx+1
                if(wss[i>=17.2):
                    break
           
            if(idx>=0):
                orglons.append(lons[idx)
                orglats.append(lats[idx)
    tf.close()


num = len(orglats)
xstart = 98.75
ystart = -1.25
deltad = 2.5
R=6371.0

xgrid = arange(xstart,182.5,deltad)#98.75E-181.25E
ygrid = arange(ystart,62.5,deltad) #1.25S - 61.25N

xn = len(xgrid)
yn = len(ygrid)

print num,xn,yn
numsarray = zeros((yn,xn))
for j in range(yn):
    for i in range(xn):
        n = 0
        for k in range(num):
            x1 = xgrid[i
            y1 = ygrid[j
            x2= orglons[k
            y2 = orglats[k

            if(abs(x2-x1)<3 and abs(y2-y1)<3):
                x1=x1*3.1415926/180
                y1=y1*3.1415926/180
                x2=x2*3.1415926/180
                y2=y2*3.1415926/180
                d=R*math.acos(cos(y1)*cos(y2)*cos(x1-x2)+sin(y1)*sin(y2))
                if(d<=250.0):
                    n = n+1
        numsarray[j,i = n
print numsarray

axesm()
lworld = shaperead('E:/MeteoInfo/MeteoInfo/Map/country.shp')
geoshow(lworld, facecolor=[200,200,200)
# Plot
#layer1 = contour(xgrid,ygrid,numsarray)
#colorlabel(layer1)
layer2 = contourf(xgrid,ygrid,numsarray)
colorbar(layer2)