以下全文代码均已发布至和鲸社区，点击下面链接前往，可一键fork跑通：https://www.heywhale.com/mw/project/6195cf2bb2739500178eef9f我国东南沿海地区常年受到台风的影响，为了评估台风对当地的影响，基础的工作之一就是提取出台风影响到该地区时的资料。为了提取影响台风的数据，常用的方法是以站点为中心，给定一个影响半径，确定一个圆。根据台风7级风影响半径的统计，一般影响半径取300 km。
台风历史资料来源于中国气象局上海台风所（CMA-STI）的最佳路径数据集，数据的获取地址为：http://tcdata.typhoon.org.cn。台风历史数据格式为文本格式，时间分辨率通常为6小时，一些台风生命史过程中部分时段的时间分辨率可加密至3小时。现行最新版本的台风最佳路径数据集提供了1949以来年西北太平洋海域（含南海，赤道以北，东经180°以西）的热带气旋每6小时（部分3小时）的位置和强度。数据存储形式为每年一个txt文件，每个txt文件包含当年的所有热带气旋。CMA-STI台风数据集资料，2018年2号台风“SANBA”的部分路径数据以珠海为例，提取1949-2018年影响珠海的台风数据，最终提取到219个影响珠海的台风，平均每年3.13个台风。import numpy as npimport csvfrom math import sin,radians,cos,asin,sqrtimport datetimebegYear = 1949endYear = 2018totalYear = endYear-begYear+1radiusInflu = 300 # influence radius ,unit:KMsiteName = 'Zhuhai' #珠海latSite = 22.28lonSite = 113.58
输入两个点经纬度(degree)，输出两点距离(km)
def SphereDistance(lon1, lat1, lon2, lat2): radius = 6371.0 lon1, lat1,lon2, lat2 = map(radians,[lon1, lat1,lon2, lat2]) dlon = lon2 -lon1 dlat = lat2 -lat1 arg = sin(dlat*0.5)**2 + \ cos(lat1)*cos(lat2)*sin(dlon*0.5)**2 dist = 2.0 * radius * asin(sqrt(arg)) return dist
对每个台风的经纬度信息，计算与珠海站的距离，判断是否在影响珠海的范围内，如果是则保存相应数据。fileName = []for i in range(begYear,endYear+1): iFile = 'CH'+str(i)+'BST.txt' fileName.append(iFile) allInfluTC = [] allInfluLine = [] for i in range(0,totalYear): filePath = 'CMABSTdata1949_2018/' + fileName[i] print("Processing File : %s" %filePath ) fileRead = open(filePath) while True: line = fileRead.readline() if not line: break if line[0:5] == "66666": numberTy = line[20:25] continue newLine = ['numberCN', 'yyyymmddhh', 'latRec', 'lonRec'
, 'presRec','gradeRec'] if line[8:10] in ["00","06","12","18"]: numberTy = numberTy yyymmddhh = line[0:10] latRec = float(line[13:16]) * 0.1 lonRec = float(line[17:21]) * 0.1 presRec = line[22:26] gradeRec = line[11:12] newLine[0] = numberTy newLine[1] = yyymmddhh newLine[2] = str(latRec) newLine[3] = str(lonRec) newLine[4] = presRec newLine[5] = gradeRec if int(numberTy) == 0: continue distTy2Site = SphereDistance(float(lonRec),float(latRec),lonSite,latSite) if distTy2Site > radiusInflu: continue else: allInfluLine.append(newLine) if numberTy not in allInfluTC: allInfluTC.append(numberTy) fileRead.closeallNewLineArr = np.array(allInfluLine)outFileName = siteName+"_"+str(begYear)+"-"+str(endYear)+".csv"print("output data : ",outFileName)np.savetxt(outFileName, allNewLineArr, delimiter = ',', fmt='%s')print("end program :",datetime.datetime.now())print(allInfluTC)print(len(allInfluTC))
以上测试代码和数据，可在气海同途公众号回复关键字珠海台风获取。
