热带气旋【CH报文数据插值】中央气象台-台风路径数据每小时插值

对CH报文数据进行每小时插值

1.原始数据文件

2.数据

文件中的原始数据长这样：

可以看到原始数据大部分都是三小时一次的报文数据，报文内容包含时间，台风路径经纬度、气压、强度。
小部分数据是六小时一次的报文数据。

3.需求

目前咱们的需求是按小时补齐热带气旋路径信息。
效果如下：

报文中台风的起报时间和结束时间不变，中间时间是每小时一条数据。
完成该需求功能的代码如下：

插值代码

# 对ch文件插值

import pandas as pd
import datetime
import os


def interpolate_ch_one_hour (file_name):
    new_file_name=file_name.split('.')[0]+'_new.txt'

    with open(file_name,'r') as f:
        content=f.readlines()

    def write_line(line):
        if not line.endswith('\n'):
            line+='\n'
        try:
            with open(new_file_name,'r') as f:
                content=f.readlines()
                if line == content[-1]:
                    print(line+' is already in the file!')
                    return  # 重复行不写入
        except Exception as e:
            pass
        # print('开始写入数据:',line)
        with open(new_file_name,'a') as f:
            f.write(line)
    # 每个台风路径插值
    for i in range(len(content)):
        if content[i].startswith('66666'):
            write_line(content[i])  # 写入第一行
            for j in range(i+1,len(content)-1):
                if content[j+1].startswith('66666'):
                    write_line(content[j])  # 写入中间行
                    break
                else:
                    st_line=content[j].split(' ')
                    et_line=content[j+1].split(' ')
                    if st_line[2]=='':
                        st_line=st_line[:2]+st_line[3:]
                    if et_line[2]=='':
                        et_line=et_line[:2]+et_line[3:]
                    if st_line[3]=='':
                        st_line=st_line[:3]+st_line[4:]
                    if et_line[3]=='':
                        et_line=et_line[:3]+et_line[4:]
                    # 输入前后时间和经纬度，等小时间隔插值
                    st_time = st_line[0]
                    st_lat = float(st_line[2])
                    st_lon = float(st_line[3])
                    et_time = et_line[0]
                    et_lat = float(et_line[2])
                    try:

                        et_lon = float(et_line[3])
                    except Exception as e:
                        print(e)
                        print(et_line)
                        print(et_line[3])
                        print(j)

                    time_interval = (datetime.datetime.strptime(et_time,'%Y%m%d%H')-datetime.datetime.strptime(st_time,'%Y%m%d%H')).seconds/60/60
                    lat_interval = (et_lat-st_lat)/time_interval
                    lon_interval = (et_lon-st_lon)/time_interval
                    for k in range(int(time_interval)):
                        time_now = datetime.datetime.strptime(st_time,'%Y%m%d%H')+pd.Timedelta(minutes=k*60)
                        lat_now = st_lat+lat_interval*k
                        lon_now = st_lon+lon_interval*k
                        line_now = time_now.strftime('%Y%m%d%H')+' '+st_line[1]+' '+str(int(lat_now))+' '+str(int(lon_now))+' '+' '.join(st_line[4:])+''
                        write_line(line_now)

    write_line(content[-1])

if __name__ == '__main__':
    list_file=['CH2021BST.txt']  #　修改成需要插值的文件名，支持多个文件
    # list_file=['CH2023BST.txt','CH2022BST.txt','CH2021BST.txt']
    for file_name in list_file:
        interpolate_ch_one_hour(file_name)

新生成的文件名是原始文件名+‘‘_new.txt’’，例如原始文件是‘‘CH2021BST.txt’’，插值后的保存文件为“CH2021BST_new.txt”。