import numpy as np
import os
import pandas as pd
import time
import datetime
import math


# 读csv文件
def read_csv(csv_path):
    df = pd.read_csv(csv_path, encoding='GBK')
    return df


# 写一个OQA文件
def write_OQA_file(station_id, time_str, OQA_file_name, csv_dir, time_str_08, days, time_detal, station_height, lon, lat):
    # 组织需要参数  # TODO: 待格式确定，参数部分写入参数py中
    # BASE_PATH = r'D:\PycharmProjects\aermod'

    # station_id = 53487

    miss_value_list = ['', '-9', '99999', '99999', '999', '9999', '09999', '09999', '09999', '09999', '09999',
                       '09999', '9999', '9999', '99', '999', '99999', '999', '999', '999', '999', '999', '999']

    date_start = datetime.datetime.strptime(time_str, '%Y%m%d')
    diff_day = datetime.timedelta(days=days-1)
    date_end = date_start + diff_day
    time_end = date_end.strftime('%Y%m%d')
    # 写开头
    t_start = str(time_str[2:4]) + '/' + str(time_str[4:6]) + '/' + str(time_str[6:8])
    t_end = str(time_end[2:4]) + '/' + str(time_end[4:6]) + '/' + str(time_end[6:8])
    # 表头信息
    head_1 = '*% SURFACE'
    head_2 = '*      XDATES    ' + t_start + ' TO ' + t_end
    point_lat0 = round(lat, 3)
    point_lon0 = round(lon, 3)
    # 站点编号 经纬度
    # TODO： tadjust 参数，当地时与时间时的差别，东经为负
    # time_detal = int(np.round(-point_lon0 / 15, 0))
    head_3 = f'*@     LOCATION  {station_id}  {point_lat0}N  {point_lon0}E  {time_detal}  {station_height}'
    head_4 = '*** EOH: END OF SURFACE QA HEADERS'

    head_list = [head_1, head_2, head_3, head_4]
    str_list = [line + '\n' for line in head_list]  # 在list中加入换行符

    fo = open(OQA_file_name, "w")
    fo.writelines(str_list)
    fo.close()


    # 获取地面的数据
    csv_name = '{0}_{1}_{2}_{3}.csv'.format(time_str, 'meteo', station_id, time_str_08)
    csv_path = os.path.join(csv_dir, csv_name)
    surface_df = read_csv(csv_path)

    n = 24 * days  # 一共要写多少个小时
    # 整理开始时间和递增差
    date_str = time_str + '00'
    date_str_ = datetime.datetime.strptime(date_str, '%Y%m%d%H')
    diff_hour = datetime.timedelta(hours=1)
    for j in range(n):
        with open(OQA_file_name, 'a') as f:
            date_str = date_str_.strftime('%Y%m%d%H')

            # 降水  缺失值：-9
            PRCP = surface_df.loc[j, '降水'] if not np.isnan(surface_df.loc[j, '降水']) else -9
            # 站点气压  缺失值：99999
            PRES = surface_df.loc[j, '气压'] if not np.isnan(surface_df.loc[j, '气压']) else 99999
            # 总云量  缺失值：999
            if not np.isnan(surface_df.loc[j, '气压']):
                tem = int(surface_df.loc[j, '总云量'])
                if tem < 10:
                    if tem < 0:
                        tem = 0
                    TSKC = f'0{tem}0{tem}'
                else:
                    TSKC = f'1010'
            else:
                TSKC = str(999)
            # 干球温度（气温）  缺失值：-9990
            TMPD = surface_df.loc[j, '温度'] if not np.isnan(surface_df.loc[j, '温度']) else -9990
            # 相对湿度  缺失值：999
            RHUM = int(surface_df.loc[j, '相对湿度']) if not np.isnan(surface_df.loc[j, '相对湿度']) else 999
            # 风向  缺失值：999
            WDIR = surface_df.loc[j, '风向'] if not np.isnan(surface_df.loc[j, '风向']) else 999
            # 风速  缺失值：999
            WSPD = surface_df.loc[j, '风速'] if not np.isnan(surface_df.loc[j, '风速']) else 999

            # 控制最大最小值
            if PRCP != -9:
                PRCP = int(PRCP * 100)
                if PRCP < 0:
                    PRCP = 0
                elif PRCP > 25400:
                    PRCP = 25400
            if PRES != 99999:
                PRES = int(PRES * 10)
                if PRES < 9000:
                    PRES = 9000
                elif PRES > 10999:
                    PRES = 10999
            if TMPD != -9990:
                TMPD = int(TMPD * 10)
                if TMPD < -300:
                    TMPD = -300
                elif TMPD > 360:
                    TMPD = 360
            if RHUM != 999:
                if RHUM < 0:
                    RHUM = 0
                elif RHUM > 100:
                    RHUM = 100
            if WDIR != 999:
                WDIR = int(WDIR / 10)
                if WDIR <= 0:  # 风向不为0
                    WDIR = 1
                elif WDIR > 36:
                    WDIR = 36
            if WSPD != 999:
                WSPD = int(WSPD * 10)
                if WSPD < 0:
                    WSPD = 0
                elif WSPD > 500:
                    WSPD = 500

            f.write(' {0:>8s} {1:>5d} {2:>5s} {3:>5d} {4:>5s} {5:>5s} {6:>5s} {7:>5s} {8:>5s} {9:>5s} {10:>5s}\n'
                    '         {11:>5s} {12:>5s} {13:>5s} {14:>5s} {15:>5s} {16:>5s} {17:>5d} {18:>5s} {19:>5s}'
                    ' {20:>5d} {21:>5d} {22:>5d}  {23}\n'.
                    format(date_str[2:], PRCP, miss_value_list[2], PRES, miss_value_list[4], TSKC, miss_value_list[6],
                           miss_value_list[7], miss_value_list[8], miss_value_list[9], miss_value_list[10],
                           miss_value_list[11], miss_value_list[12], miss_value_list[13], miss_value_list[14],
                           miss_value_list[15], miss_value_list[16], TMPD, miss_value_list[18], miss_value_list[19],
                           RHUM, WDIR, WSPD, 'N'))

        # 加一个小时
        date_str_ += diff_hour


if __name__ == '__main__':
    station_id = 53487
    time_str = '20201115'
    OQA_file_name = os.path.abspath(f'real_time/aermod_stastic/lankao/weatherdata/{time_str}/{station_id}_surface.OQA')
    csv_dir = os.path.abspath(f'real_time/data/{time_str}')
    time_str_08 = '2020102108'
    days = 2
    time_detal = 0  # 参数，固定为0（若数据使用北京时间记录，时间转换因子为0；若使用世界时记录，转换因子在中国地区设置为-8)
    station_height = 1067.2
    lon = 113.25
    lat = 40.05
    # station_height = 73.3
    # lon = 114.49
    # lat = 34.51
    write_OQA_file(station_id, time_str, OQA_file_name, csv_dir, time_str_08, days, time_detal, station_height, lon, lat)