关于台风移动路径频数的脚本，请王老师指点。

qxsw2016

根据王老师给的参考脚本，利用CMA热带气旋最佳路径数据集，计算1981-2017年7月和8月生成的台风经过北纬0-60度，东经100-180度，网格距2.5X2.5，的移动路径频数。想仅绘制频数等值线（填色图），直接的想法是，用topology模块的 intersects 函数判断两个几何图形是否相交，来得到区域网格的频数，存入文件，然后读取数据绘图。下面的脚本，如何去除台风路径，请老师指教。
脚本如下：
# Read typhoon data file
fn = r'F:\DataSet\BestTrack\CH%dBST.txt'
lons = []
lats = []
prss = []
for yy in range(1980,2018):
    datafn = fn%yy
    #print datafn
    tf = open(datafn)
    for line in tf:
        data = line.split()
        pn = int(data[2])

        #读取每个台风数据第一行
        line = tf.readline()
        data = line.split()         

        t = data[0]
        mm = t[4:6]

        if(mm == '07'or mm == '08'):
            #print data
            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)        
            for i in range(pn-1): #读取台风余下数据

                line = tf.readline()
                data = line.split()
                #print data        
                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)
        else:
            for j in range(pn-1):
                line = tf.readline()   

        lons.append(nan)
        lats.append(nan)
        prss.append(nan)

# Plot
axesm()
lworld = shaperead('F:/MeteoInfo/MeteoInfo/Map/country1.shp')
geoshow(lworld, facecolor=[200,200,200])
plotm(lons, lats, linewidth=1)
layer = plotm(lons, lats, linewidth=1)
xlim(100, 210)
ylim(0, 60)
#title('Typhoon pathway')

xstart = 98.75
ystart = -1.25
deltad = 2.5

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)

numsarray = zeros((yn,xn))
lat = zeros(5)
lon = zeros(5)
for j in arange(yn):
    for i in arange(xn):
        #Add a rectangle polygon
        x0 = xgrid
        x1 = x0+delta
        y0 = ygrid[j]
        y1 = y0+delta
        #lat = array([20, 20, 10, 10, 20])
        lat[:] =y1
        lat[2] = y0
        lat[3] = y0
        #lon = array([180, 190, 190, 180, 180])
        lon[:] = x0
        lon[1] = x1
        lon[2] = x1

        g1 = geoshow(lat, lon, displaytype='polygon', color=[0,0,0,0], edgecolor='w', size=0)
        #Count the number of typhoons cross the rectangle
        n = 0
        for s in layer.shapes():
            if tp.intersects(s, g1):
                n += 1
            numsarray[j,i] = n
#print yn,xn
for j in arange(yn):
    for i in arange(xn):
        if(numsarray[j,i]<10):
            numsarray[j,i] = nan
#print numsarray
# Plot
xgrid[:] = xgrid[:] +deltad/2
ygrid[:] = ygrid[:] +deltad/2
#layer1 = contour(xgrid,ygrid,numsarray)
#colorlabel(layer1)
layer2 = contourf(xgrid,ygrid,numsarray)
colorbar(layer2)

fey

你好，请问你是如何计算台风经过某个区域，网格距5X5，的移动路径频数呢？

dieam

请问楼主解决了么

polointer

感谢分享{:eb502:}{:eb502:}