## 上面的基本就可以实现了,接下来来尝试算2019年所有天的轨迹
## 先做2019年的轨迹
# 先获取所有的气象数据于一个列表中
import os
train_dir=r'H:\meteorology\gdas_meteodata\2019'
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/TrajStat\working'
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_traj\2018'
file_list=os.listdir(path)
data_list=[]
s='PRESSURE'
for file in file_list:
print(file)
#************************************************************需要改
f = open(r'D:\CB2_traj\2018/'+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_traj\2018/'+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()