Python利用meteinfo来计算后向轨迹

Python利用meteinfo来计算后向轨迹

## 上面的基本就可以实现了,接下来来尝试算2019年所有天的轨迹
## 先做2019年的轨迹
# 先获取所有的气象数据于一个列表中
import os

train_dir=r'H:meteorologygdas_meteodata2019'
datanames = os.listdir(train_dir)

#确定要跑的那一天
#先建一个该年所有天的年月日
import pandas as pd
import datetime
days=pd.date_range('20151120','20181231',freq='1D')
days=pd.to_datetime(days)
for k in range(len(days)):
    #以最后一天为例
    test_day=days[k]
    print(test_day)
    year_flag=test_day.year
    month_flag=test_day.month
    day_flag=test_day.day
    this_year_str=str(year_flag)[-2:]

    # 需要跑的那一天
    stime = datetime.datetime(year_flag,month_flag,day_flag)
    flag_this_month=stime.strftime('%b').lower()
    
    import datetime
    from dateutil.relativedelta import relativedelta
    # 上一个月是
    last_month_day=test_day- relativedelta(months=+1)
    temp=datetime.datetime(last_month_day.year,last_month_day.month,last_month_day.day)
    flag_last_month=temp.strftime('%b').lower()
    last_year_str=str(last_month_day.year)[-2:]
    
    # 开始搜集气象数据
    aim_meteo=[]
    #本月
    for item in datanames:
        if flag_this_month+this_year_str in item:
            aim_meteo.append(item)
    #上月最后两个星期
    for item in datanames:
        if flag_last_month+last_year_str in item:
            if 'w4' in item:
                aim_meteo.append(item)
            if 'w5' in item:
                aim_meteo.append(item)
    # 整理完毕
    # 开始跑这天的轨迹：
    import calendar
    import os
    import datetime

    # Set working directory
    metDir = 'H:/meteorology/gdas_meteodata/2019'
    outDir = 'D:/CB2_traj/traj_output'
    workingDir = 'D:/meteinfo/TrajStat_Plugin_1.4.9/TrajStatworking'
    os.chdir(workingDir)
    print('Current directory: ' + os.getcwd())

    # Set parameters
    lon = '166.37377'
    lat = '-77.243138'

    hours = '-72'
    vertical = '0'
    top = '10000.0'

    # Set meteorological data files
    fns = aim_meteo

    shour_list=['00','06','12','18']
    for shour in shour_list:
        heights = ['500']
        #heights = [str(interg) for interg in np.random.randint(201,size=(3,))]
        hnum = len(heights)
        print(hnum)

        # Write ConTROL file
        ctFile = './CONTROL'
        print(stime.strftime('%Y-%m-%d ') + shour + ':00')
        ctf = open(ctFile, 'w')
        ctf.write(stime.strftime('%y %m %d ') + shour + "n")
        ctf.write(str(hnum) + 'n')
        for i in range(0,hnum):
            ctf.write(lat + ' ' + lon + ' ' + heights[i] + 'n')
        ctf.write(hours + 'n')
        ctf.write(vertical + 'n')
        ctf.write(top + 'n')
        fnnum = len(fns)
        ctf.write(str(fnnum) + 'n')
        for i in range(0,fnnum):
            ctf.write(metDir + '/' + 'n')
            ctf.write(fns[i] + 'n')
        ctf.write(outDir + '/' + 'n')
        outfn = stime.strftime('traj_%Y%m%d'+shour)
        ctf.write(outfn)
        ctf.close()
        # Calculate trajectories
        os.system('D:/meteinfo/TrajStat_Plugin_1.4.9/TrajStat/working/hyts_std.exe')

        print (shour+'Finish...')    

将轨迹进行合成，成为一个netcdf文件

################ 将所有的轨迹进行制作
import datetime
import numpy as np
import pandas as pd
import os
#************************************************************需要改
path=r'D:CB2_traj2018'
file_list=os.listdir(path)

data_list=[]
s='PRESSURE'
for file in file_list:
    print(file)
    #************************************************************需要改
    f = open(r'D:CB2_traj2018/'+file,'r').readlines()
    for i in range(len(f)):
        if s in f[i]:
            print ('line:',i+1)
            break
    #************************************************************需要改
    data = np.loadtxt(r'D:CB2_traj2018/'+file,skiprows=i+1)
    data=pd.Dataframe(data)
    data.columns=['traj_no','x1','start_year','start_mon','start_day',
                  'start_hour','x2','x3','age_hour','latitude','longitude','height','press','x4']
    for k in range(len(data)):
        data.iloc[k,1]=file+'_'+str(data.iloc[k,0])
    data_list.append(data)

ini=data_list[0]
for df_i in range(1,len(data_list)):
    print(df_i)
    ini=ini.append(data_list[df_i])
data=ini

del ini
data.shape
del data_list
data.iloc[4,:]

#新建一个ID列，并初始化为0
data['ID']=data['start_day']*0

#根据研究区域画分网格，这里用0.5的分辨率
#纬度：16.85-67.25，可以画分为101个左右
#经度：32.73-126.98，可以划分为189个左右
dpi=1.0
LON1 = int(data['longitude'].min())-1
LON2 = int(data['longitude'].max())+1
LAT1 = int(data['latitude'].min())-1
LAT2 = int(data['latitude'].max())+1
N_lon=int((LON2-LON1)/dpi)
N_lat=int((LAT2-LAT1)/dpi)

def generalID(lon,lat):
    # 若在范围外的点，返回-1
    if lon <= LON1 or lon >= LON2 or lat <= LAT1 or lat >= LAT2:
        return -1
    # 把经度范围根据列数等分切割
    column = (LON2 - LON1)/N_lon
    # 把纬度范围根据行数数等分切割
    row = (LAT2 - LAT1)/N_lat
    # 二维矩阵坐标索引转换为一维ID，即： （列坐标区域（向下取整）+ 1） + （行坐标区域 * 列数）
    return int((lon-LON1)/column)+ 1 + int((lat-LAT1)/row) * N_lon

data.shape[0]
#对每个点计算对应网格ID
for i in range(data.shape[0]):
    print(i)
    data.iloc[i,-1]=generalID(data.iloc[i,-5],data.iloc[i,-6])
data=data.set_index(data['ID'])
# 计算每个网格中所有的轨迹数目
# 计算每个网格中所有的点数，代表总的停留时间
# 平均每条轨迹的停留时间=总的停留时间/网格中所有的轨迹数目
#找出所有不同索引值
ID=data['ID'].value_counts()
ID_arr=np.array(ID.index)
# 新建一个存放每个ID数据的数据表
grid_para=pd.Dataframe(np.arange(len(ID))+1)
grid_para['ID_pass']=ID_arr
grid_para=grid_para.set_index(grid_para['ID_pass'])
grid_para['traj_number']=grid_para['ID_pass']*0
grid_para['traj_point']=grid_para['ID_pass']*0
grid_para['residual time']=grid_para['ID_pass']*0

for i in ID_arr:
    print(i)
    temp=data[(data.index==i)]
    # 点的数目
    point_num=len(temp)
    # 轨迹的数目
    traj_num=len(temp['x1'].value_counts())
    # 平均的点数or小时数
    mean_hours=point_num/traj_num
    
    grid_para.loc[i,'traj_number']=point_num
    grid_para.loc[i,'traj_point']=traj_num
    grid_para.loc[i,'residual time']=mean_hours*3600
    
    
#下面进行格点和经纬度对应
#人为设定经纬度网格的中心点
#总共的ID/网格数目是N_lon*N_lat
LON=np.linspace(LON1+dpi/2,LON2-dpi/2,N_lon)
LAT=np.linspace(LAT1+dpi/2,LAT2-dpi/2,N_lat)
ID=np.arange(N_lon*N_lat)+1

Total_range_data=pd.Dataframe(ID,columns=(['NUM']))
Total_range_data['ID']=Total_range_data['NUM']*0
Total_range_data['lon']=Total_range_data['NUM']*0
Total_range_data['lat']=Total_range_data['NUM']*0

for i in range(N_lon):
    for j in range(N_lat):
        Total_range_data.iloc[i+j*N_lon,1]=i+j*N_lon+1
        Total_range_data.iloc[i+j*N_lon,2]=i*dpi+LON1+dpi/2
        Total_range_data.iloc[i+j*N_lon,3]=j*dpi+LAT1+dpi/2
del Total_range_data['NUM']

#将ID列设为索引值
Total_range_data=Total_range_data.set_index('ID')
Total_data=Total_range_data.join(grid_para)

import netCDF4 as nc
## create NetCDF file
newfile = nc.Dataset('newfile2018.nc', 'w', format='NETCDF4')

## define dimesions
long = newfile.createDimension('longitude', size=N_lon)
lati = newfile.createDimension('latitude', size=N_lat)
## define variables for storing data
lon = newfile.createVariable('lon', 'f4', dimensions='longitude')
lat = newfile.createVariable('lat', 'f4', dimensions='latitude')
traj_number = newfile.createVariable('traj_number', 'f4', dimensions=('longitude', 'latitude'))
traj_point = newfile.createVariable('traj_point', 'f4', dimensions=('longitude', 'latitude'))
residual_time = newfile.createVariable('residual time', 'f4', dimensions=('longitude', 'latitude'))

## add data to variables
lon[:] = LON
lat[:] = LAT
# lon_trans=np.array(Total_data['lon'])
# lon_trans=lon_trans.reshape((N_lat,N_lon)),还需要进一步转置一下
h1=np.array(Total_data['traj_number'])
h1=h1.reshape((N_lat,N_lon)).T
traj_number[:,:]=h1

h2=np.array(Total_data['traj_point'])
h2=h2.reshape((N_lat,N_lon)).T
traj_point[:,:]=h2

h3=np.array(Total_data['residual time'])
h3=h3.reshape((N_lat,N_lon)).T
residual_time[:,:]=h3

newfile.close()