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

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:}

qxsw2016

热带气旋最佳路径数据集，文件头记录含热带气旋的英文名称, 名称后加 “(-1)n” 表示副中心及其序号，把这个部分记录不列入统计，整年路径频数统计，在Lab新版本代码如下，但有两点有疑惑。一是plot函数有的版本是忽略nan值的，二是出图为什么有陆地遮盖效果。脚本如下:

import mipylib.geolib.topology as tp
import re

# Read typhoon data file
fn = 'D:/Temp/test/ty/CH%dBST.txt'
lons = [
lats = [
prss = [
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:
            #读取每个台风数据第一行
            line = tf.readline()
            data = line.split()  
        
            #t = data[0]
            #mm = t[4:6]        
            #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(rown-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)
           
        lons.append(nan)
        lats.append(nan)
        prss.append(nan)

# Plot
axesm()
geoshow('country', facecolor=[200,200,200)
layer = plot(lons, lats, linewidth=1)
xlim(100, 210)
ylim(0, 60)

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[i
        x1 = x0+deltad
        y0 = ygrid[j
        y1 = y0+deltad
        #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='c', 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
        clf()
#print yn,xn
#numsarray[numsarray<10] = nan
#print numsarray
# Plot
xgrid[: = xgrid[: +deltad/2
ygrid[: = ygrid[: +deltad/2

axesm()
geoshow('country', facecolor=[200,200,200)
layer2 = contourfm(xgrid,ygrid,numsarray)
colorbar(layer2)