import netCDF4 as nc4
import matplotlib.pyplot as plt
import numpy as np
import datetime
import os
import cartopy.crs as ccrs

path='E://' #文件路徑
files= os.listdir(path) #得到文件夾下的所有文件名稱

fig=plt.figure(figsize=(20,12)) #設(shè)置畫布大小
parallels = np.arange(0.,90.,3.) 
meridians = np.arange(0.0,360.,3.)   
ax = plt.axes(projection=ccrs.Robinson())   #設(shè)置投影方式
# Set figure extent & ticks
ax.set_extent([100, 150, 5, 50])  #設(shè)置緯度范圍
# plt.grid(linestyle=':',color='y')#
for file in files:  #按照順序在 files 里面進(jìn)行每一個(gè)文件的 數(shù)據(jù)名稱 循環(huán)讀取
    f = open(path+file,'r')   # 打開第一個(gè) dat 文件
    records = f.readlines()   # 讀取這個(gè)文件里面的所有數(shù)據(jù)
    f.close()                 # 關(guān)閉這個(gè)dat文件  
   
    date_t = ''               # 設(shè)置一個(gè)用來表示空格
    btk_lat = []              # 設(shè)置一個(gè)空的 待傳入數(shù)據(jù)的緯度
    btk_lon = []              # 設(shè)置一個(gè)空的 待傳入數(shù)據(jù)的經(jīng)度
    btk_vmax = []             # 風(fēng)速最大值 Maximum sustained wind speed in knots: 0 - 300 kts.
    btk_time = []             # 時(shí)間
    btk_rmw = []              # 最大風(fēng)速半徑 radius of max winds, 0 - 999 n mi.
    btk_name = []             # 臺(tái)風(fēng)名稱

    for rcd in records:       # 對(duì)這個(gè)dat文件里面，已經(jīng)讀取的每一行數(shù)據(jù)進(jìn)行循環(huán)處理
        
        strs = rcd.split(',') #將每一個(gè)','分隔開
        if(len(strs)<21):     #判斷語句，如果這個(gè)被分割開的字符 長(zhǎng)度<21 ,繼續(xù)進(jìn)行處理
            continue
        date_str = strs[2].strip(' ') #將strs這個(gè)list的索引為2的值賦給data_str,既年月日時(shí)
        if date_str == date_t:#判讀如果是一個(gè)空格值，賦給data——str
            continue
        dt = datetime.datetime(int(date_str[0:4]),int(date_str[4:6]),int(date_str[6:8]),\
        int(date_str[8:]),0,0,0)
        btk_time.append(nc4.date2num(dt,units='second since 1970-1-1 00:00:00'))#計(jì)算距離給的時(shí)間有多少秒，并從后往前排列
        #處理緯度
        lat_str = strs[6].strip()
        #判斷南北緯
        if lat_str[-1] == 'N':
            lat_t = float(lat_str[0:-1])*0.1
        else:
            lat_t = float(lat_str[0:-1])*-0.1
        btk_lat.append(lat_t)
        #處理經(jīng)度
        lon_str = strs[7].strip()
        #判斷 東西經(jīng)
        if lon_str[-1] == 'E':
            lon_t = float(lon_str[0:-1])*0.1
        else:
            lon_t = float(lon_str[0:-1])*-0.1
        btk_lon.append(lon_t)        
        #處理最大風(fēng)速
        vmax = strs[8].strip()
        btk_vmax.append(float(vmax))#轉(zhuǎn)換為單浮點(diǎn)型，（帶小數(shù)點(diǎn)）
        #時(shí)間
        date_t = date_str
        #最大風(fēng)速半徑
        rmw = strs[19].strip()
        btk_rmw.append(float(rmw))
        #處理臺(tái)風(fēng)名稱
        if(len(strs) < 27):
            btk_name.append('noname')
        else:
            name = strs[27].strip()
            btk_name.append(name)
#==============================================================================
    btk_lat = np.array(btk_lat) #將得到的list 值轉(zhuǎn)換為數(shù)組型的值，為了便于繪圖。因?yàn)槔L圖的橫縱坐標(biāo)都是數(shù)組排列
    btk_lon = np.array(btk_lon)%360 #因?yàn)樵冀?jīng)度為-180 - 0 -180 ，出現(xiàn)斷隔，為解決問題，化為 0-360
    btk_time = np.array(btk_time)  #時(shí)間轉(zhuǎn)換
    btk_vmax = np.array(btk_vmax)*0.5144 #風(fēng)速換算公式
    btk_rmw = np.array(btk_rmw)*1.852 #
    #判斷，如果數(shù)組緯度的值是0，則為nan值，既無法計(jì)算的值（無窮大，，），否則即為臺(tái)風(fēng)的名稱
    if(len(btk_lat) == 0): 
        tc_name = 'noname'
    else:
        index = btk_vmax.argmax()
        tc_name = btk_name[index]
    
    #進(jìn)行繪圖，經(jīng)度、緯度曲線
    ax.plot(btk_lon,btk_lat,color='k',linewidth=0.5,transform=ccrs.PlateCarree())
    #散點(diǎn)圖繪制，經(jīng)度、緯度、最大風(fēng)速，
    cb = ax.scatter(btk_lon,btk_lat,c=btk_vmax,s=10.0,transform=ccrs.PlateCarree()
                    ,vmin=10,vmax=60)
ax.coastlines()
ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)

plt.colorbar(cb,label='Vmax (m/s)',pad=0.07,orientation='vertical',shrink=1)
plt.title('  path')
# 保存繪制圖片 ，注意保存路徑不能放在dat文件夾中     
#fig.savefig(path2+'tester.tiff',format='tiff',dpi=100)