# -*- coding: utf-8 -*-
import datetime
import json
from pyproj import Proj
import os
from math import ceil
import shutil
import generate_receptor_grid
import time
import treat_json
import logging
import numpy as np
from osgeo import gdal
import rasterio as rio
import subprocess
import write_surface_OQA
import write_upair_OQA
import generate_json_config
import deal_result
from OSS import OSS


def wgs84toutm(lon, lat):
    '''
    :param lon:
    :param lat:
    :return:UTM投影带号和相应UTM坐标,及相对距离坐标
    '''
    zone = ceil((lon + 180) / 6.0000001)
    utm_pj = Proj("+proj=utm +zone=" + str(zone) + ", +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
    utm_x, utm_y = utm_pj(lon, lat, inverse=False)  # 从经纬度转平面坐标
    return zone, utm_x, utm_y


def relative_distance_coord(center_lon, center_lat, des_lon, des_lat):
    '''
    将全球坐标（经纬度）转为项目本地坐标（相对距离坐标）
    :param x:经度
    :param y:纬度
    :return:相对距离坐标的x,y
    '''
    center_x, center_y = wgs84toutm(center_lon, center_lat)[1:3]
    des_x, des_y = wgs84toutm(des_lon, des_lat)[1:3]
    relative_x = des_x - center_x
    relative_y = des_y - center_y
    return relative_x, relative_y


def copytree(src, dst, symlinks=False, ignore=None):
    if not os.path.exists(dst):
        os.makedirs(dst)
    for item in os.listdir(src):
        s = os.path.join(src, item)
        d = os.path.join(dst, item)
        if os.path.isdir(s):
            copytree(s, d, symlinks, ignore)
        else:
            if not os.path.exists(d) or os.stat(s).st_mtime - os.stat(d).st_mtime > 1:
                shutil.copy2(s, d)


def prepare_dirs(model_dir, work_dir, calculateID):

    work_path = work_dir + "/" + '{0}'.format(calculateID)
    work_path_aermap = work_path + "/" + "AERMAP"
    work_path_aermet = work_path + "/" + "AERMET"
    work_path_building = work_path + "/" + "building"
    work_path_aermod = work_path + "/" + "AERMOD"
    if os.path.exists(work_path):
        # 注释掉下面语句，不影响效果，方便跑第二遍调式
        # shutil.rmtree(work_path, ignore_errors=True)
        # time.sleep(1)
        # os.makedirs(work_path)
        # os.makedirs(work_path_aermap)
        # os.makedirs(work_path_aermet)
        # os.makedirs(work_path_building)
        # os.makedirs(work_path_aermod)
        pass
    else:
        os.makedirs(work_path)
        os.makedirs(work_path_aermap)
        os.makedirs(work_path_aermet)
        os.makedirs(work_path_building)
        os.makedirs(work_path_aermod)
    # TODO：要对输入文件完整性做check
    copytree(model_dir, work_path_aermap)
    copytree(model_dir, work_path_aermet)
    copytree(model_dir, work_path_building)
    # copytree(model_dir, work_path_aermod)
    shutil.copy(model_dir + "/aermod.exe", work_path_aermod)
    shutil.copy(model_dir + "/Bpipprm.exe", work_path_aermod)
    return work_path


def get_id(obj_json, id_target, key='ID'):
    """
    :param obj_json:  要遍历查询的json根节点对象
    :param id_target:  目标ID 数字
    :param key:   遍历的匹配字段，对排放源来说就为‘ID’
    :return: json根节点对象的index
    """
    nlength = len(obj_json)
    for index in range(nlength):
        if obj_json[index][key] == id_target:
            return index


def get_elevation_h(demfile, lon, lat):
    """
    从地形.tif 根据经纬度获取地形高度
    :param demfile:
    :param lon:
    :param lat:
    :return:
    """
    with rio.open(demfile) as src:
        row, col = src.index(lon, lat)
        dem = src.read(1, masked=True)  # 读取第一个band
        ele = dem[row, col]
    return ele


def get_source_area(source_details):
    if source_details['type'] in ['VOLUME', 'AREA', 'AREAPOLY', 'AREACIRC']:
        npoints = len(source_details['locationY'])
        x_list = []
        y_list = []
        for i in range(npoints):
            zone, utm_x, utm_y = wgs84toutm(source_details['locationX'][i], source_details['locationY'][i])
            x_list.append(utm_x)
            y_list.append(utm_y)
        # https://stackoverflow.com/questions/24467972/calculate-area-of-polygon-given-x-y-coordinates
        area = 0.5 * np.abs(np.dot(x_list, np.roll(y_list, 1)) - np.dot(y_list, np.roll(x_list, 1)))
    else:
        area = 0
    return area


def run_aermap(jsonObj, select_project):
    grid_RE = 'grid_location.RE'
    write_aermap_input(jsonObj, select_project, grid_RE)
    # TODO:运行aermap
    os.system('aermap.exe aermap.inp aermap.out')
    print('aermap  {0},运行完成'.format(grid_RE))


def write_aermap_input(jsonObj, select_project, grid_RE, if_discrete_point=1):
    cent_lon = jsonObj['project_location']['lng']
    cent_lat = jsonObj['project_location']['lat']
    # select_receptor_project_id = select_project['select_receptor_project_ID']
    # index = get_id(jsonObj['project_data']['AERMOD']['AERMOD_receptor'], select_receptor_project_id, key='ID')
    # discrete_point = jsonObj['project_data']['AERMOD']['AERMOD_receptor'][index]['discrete_point']
    discrete_point = jsonObj['project_data']['AERMOD']['AERMOD_receptor']['discrete_point']

    wmap = open('aermap.inp', 'w')
    # TODO：创建aermap.inp文件
    wmap.write('CO STARTING\n')
    wmap.write('   TITLEONE AERMAP-受体方案{0}\n'.format(1))
    wmap.write('   TERRHGTS EXTRACT\n')  # 受体点高度从原始DEM文件中提取还是用户指定
    wmap.write('   DATATYPE NED\n')  # 地形数据类型，NED和DEM均可
    wmap.write('   DATAFILE terrain.tif tiffdebug  METERS\n')
    zone, cent_x, cent_y = wgs84toutm(cent_lon, cent_lat)
    # TODO：评价区域西南角和东北角坐标和经度带的确定方法？可能存在跨经度带的情况,目前采用自中心点外扩30km
    wmap.write('   DOMAINXY {0:.1f} {1:.1f} {2:.1f} {3:.1f} {4:.1f} {5}\n'.
               format(cent_x - 30000, cent_y - 30000, zone, cent_x + 30000, cent_y + 30000, zone))
    wmap.write('   ANCHORXY 0.0 0.0 {0:.1f} {1:.1f} {2} 0\n'.format(cent_x, cent_y, zone))
    wmap.write('   DEBUGOPT ALL\n')  # 是否为受体点、源、山高生成额外的调试文件
    wmap.write('   RUNORNOT RUN\n')
    wmap.write('CO FINISHED\n')
    wmap.write('RE STARTING\n')
    # TODO:可直接从json文件读取受体点文件
    if if_discrete_point == 1:
        for di in discrete_point:
            rel_x, rel_y = relative_distance_coord(cent_lon, cent_lat, di[1], di[2])
            wmap.write('   DISCCART {0} {1}\n'.format(rel_x, rel_y))
    # wmap.write('   INCLUDED grid_location.RE\n')
    wmap.write('   INCLUDED ' + grid_RE + '\n')
    wmap.write('RE FINISHED\n')
    wmap.write('OU STARTING\n')
    wmap.write('   RECEPTOR AERMAP.REC\n')
    wmap.write('   DEBUGREC RECdetail.OUT RECndem.OUT RECelv.OUT\n')
    wmap.write('   DEBUGHIL hilldebug.OUT\n')
    wmap.write('OU FINISHED')
    wmap.close()


def extrat_aermap_rec(out_file):
    ## 将grid_2的AERMAP.REC 数据写入主AERMAP.REC
    f = open('AERMAP.REC', 'r')
    outfile = open(out_file, 'a')
    ii = 1
    for line in f:
        if ii > 12:  # 略过12行
            outfile.write(line)
        ii = ii + 1
    outfile.close()


def run_aermap_multi(outdir, jsonObj, select_project, n_sub=4):
    grid_REs = []
    for i in range(n_sub):
        grid_REs.append("grid_"+str(i+1)+".RE")

    for grid_RE in grid_REs:
        if grid_RE == 'grid_1.RE':
            if_discrete_point = 1
        else:
            if_discrete_point = 0
        tmp_path = grid_RE[:-3]
        if not os.path.exists(tmp_path):
            os.makedirs(tmp_path)
        os.chdir(tmp_path)
        write_aermap_input(jsonObj, select_project, grid_RE, if_discrete_point=if_discrete_point)
        os.system('copy ..\\aermap.exe .')
        os.system('copy ..\\terrain.tif .')
        os.system('copy ..\\terrain.tif .')
        shutil.copy('../'+grid_RE, grid_RE)
        subprocess.Popen(["start", "aermap.exe"], stdout=subprocess.PIPE,
                         stderr=subprocess.PIPE, shell=True, universal_newlines=True)
        time.sleep(10)
        os.chdir("..\\")

    time.sleep(60)
    #------------ check all aermap is finished ------------------
    need_wait = n_sub
    while need_wait > 0:
        need_wait = n_sub
        for grid_RE in grid_REs:
            tmp_path = grid_RE[:-3]
            if os.path.exists(tmp_path + "/" + "AERMAP.REC") and os.path.exists(tmp_path +"/" + "hilldebug.OUT"):
                if os.path.getsize(tmp_path +"/" + "AERMAP.REC") > 1 \
                        and os.path.getsize(tmp_path +"/" + "hilldebug.OUT") > 0:
                    need_wait = need_wait - 1
            if os.path.exists(tmp_path + "/" + "AERMAP.REC") and os.path.exists(tmp_path +"/" + "aermap.out"):
                if os.path.getsize(tmp_path + "/" + "aermap.out") > 0 \
                        and os.path.getsize(tmp_path +"/" + "AERMAP.REC")==0 \
                        and not os.path.exists(tmp_path +"/" + "hilldebug.OUT"):
                    msg_str = "        aermap.exe 运行失败"
                    running_monitor(outdir, msg_str)
                    touch_finished(outdir)  #结束整个程序
        if need_wait > 0:
            print("waiting for the finish of aermap.exe")
            msg_str = "        aermap.exe 正在计算"
            running_monitor(outdir, msg_str)
            time.sleep(30) #120s

    #------------  extract AERMAP.REC data ------------------------
    shutil.copy('./grid_1/AERMAP.REC', "AERMAP.REC")
    for grid_RE in grid_REs:
        tmp_path = grid_RE[:-3]
        if tmp_path != 'grid_1':
            os.chdir(tmp_path)
            extrat_aermap_rec("../AERMAP.REC")
            os.chdir("..\\")
    print("extract grid_* OK")
    for grid_RE in grid_REs:
        tmp_path = grid_RE[:-3]
        shutil.rmtree(tmp_path)


def run_aermet(jsonObj, select_project, project_path, first_day):
    # aermet_id = select_project['select_AERMET_project_ID']
    # index = get_id(jsonObj['project_data']['AERMOD']['AERMET_project'], aermet_id, key='ID')
    # select_aermet = jsonObj['project_data']['AERMOD']['AERMET_project'][index]
    select_aermet = jsonObj['project_data']['AERMOD']['AERMET_project']


    # TODO： AERMOD能采用几个气象文件？暂时默认为1个
    start_time = select_aermet['start_time']
    end_time = select_aermet['end_time']
    start_time_ = start_time[:4] + "/" + start_time[5:7] + "/" + start_time[8:10]
    end_time_ = end_time[:4] + "/" + end_time[5:7] + "/" + end_time[8:10]
    print("-------------------------------------------------")
    print(start_time_, end_time_)

    surface_file = select_aermet['select_surface_data']['surface_sites'][0]['data_file']
    upair_file = select_aermet["select_upair_data"]['upair_sites'][0]['data_file']
    # upair_file_path = project_path + "/basedata/weatherdata/weather_upair/" + upair_file
    upair_file_path = os.path.join(project_path, 'weatherdata', str(first_day), upair_file)
    # surface_file_path = project_path + "/basedata/weatherdata/weather_surface/" + surface_file
    surface_file_path = os.path.join(project_path, "weatherdata", str(first_day), surface_file)
    shutil.copy(upair_file_path, './')
    shutil.copy(surface_file_path, './')


    # TODO:生成aermet第一步合并气象数据的input文件
    wmet2 = open('met2.inp', 'w')
    wmet2.write('JOB\n')
    wmet2.write('  REPORT     MET2.RPT\n')
    wmet2.write('  MESSAGES   MET2.MSG\n')
    wmet2.write('\n')
    wmet2.write('UPPERAIR\n')
    # ss_ = '  QAOUT      ' + upair_file_path + '\n'
    ss_ = '  QAOUT      ' + upair_file + '\n'
    wmet2.write(ss_)
    wmet2.write('\n')
    wmet2.write('SURFACE\n')
    # ss_ = '  QAOUT      ' + surface_file_path + '\n'
    ss_ = '  QAOUT      ' + surface_file + '\n'
    wmet2.write(ss_)
    wmet2.write('\n')
    wmet2.write('MERGE\n')
    wmet2.write('  OUTPUT     MET2.MET\n')
    wmet2.write('  XDATES     {0} TO {1}\n'.format(start_time_, end_time_))
    wmet2.close()

    time.sleep(1)
    os.system('rename met2.inp aermet.inp')
    os.system('aermet.exe')
    time.sleep(5)
    # TODO：生成aermod所需的SFC和PFL文件
    wmet3 = open('met3.inp', 'w')
    wmet3.write('JOB\n')
    wmet3.write('  REPORT    met3.RPT\n')
    wmet3.write('  MESSAGES  met3.MSG\n')
    wmet3.write('METPREP\n')
    wmet3.write('  DATA      MET2.MET\n')
    wmet3.write('  OUTPUT    MET3.SFC\n')
    wmet3.write('  PROFILE   MET3.PFL\n')
    wmet3.write('  XDATES    {0} TO {1}\n'.format(start_time_, end_time_))
    select_site = select_aermet['select_surface_data']['surface_sites'][0]
    wmet3.write('  LOCATION  {0}  {1}  {2}  {3}\n'.format(select_site['site_id'], select_site['site_lng'],
                                                          select_site['site_lat'], -8))
    wmet3.write('  METHOD    REFLEVEL  SUBNWS\n')
    wmet3.write('  METHOD    UASELECT  SUNRISE\n')
    if select_aermet['random_windDirection']:
        wmet3.write('  METHOD    WIND_DIR  RANDOM\n')
    wmet3.write('  NWS_HGT   WIND      10.0\n')  # 界面未留接口，默认测风高度为10m
    ground_features = select_aermet['ground_features']  # 当前项目地表参数
    time_period = select_aermet['time_period']  # 当前项目生成地表参数的时间周期
    seperate_degree = select_aermet['separating_degree']  # 当前项目扇区数列
    len_sd = select_aermet['number_of_sectors']  # 当前项目扇区数
    # TODO:地面时间周期按年时的地表特征参数表
    if time_period == 'year':
        if len_sd < 2:
            wmet3.write('  FREQ_SECT ANNUAL  1\n')
            wmet3.write('  SECTOR    1    0    360\n')
            wmet3.write('  SITE_CHAR 1    1    {0}    {1}    {2}\n'.
                        format(ground_features[0][0]['Alb'], ground_features[0][0]['Bo'], ground_features[0][0]['Zo']))
            wmet3.close()
        else:
            wmet3.write('  FREQ_SECT ANNUAL  {0}\n'.format(len_sd))
            for i in range(len_sd - 1):
                wmet3.write('  SECTOR    {0}    {1}    {2}\n'.format(i + 1, seperate_degree[i], seperate_degree[i + 1]))
            wmet3.write('  SECTOR    {0}    {1}    {2}\n'.format(len_sd, seperate_degree[len_sd - 1], seperate_degree[0]))
            for i in range(len_sd):
                wmet3.write('  SITE_CHAR 1    {0}    {1}    {2}    {3}\n'.
                            format(i+1, ground_features[i][0]['Alb'], ground_features[i][0]['Bo'],
                                   ground_features[i][0]['Zo']))
            wmet3.close()
    # TODO:地面时间周期按季节时的地表特征参数表
    elif time_period == 'season':
        if len_sd < 2:
            wmet3.write('  FREQ_SECT SEASONAL  1\n')
            wmet3.write('  SECTOR    1    0    360\n')
            for j in range(4):
                wmet3.write('  SITE_CHAR {0}    1    {1}    {2}    {3}\n'.
                            format(j+1, ground_features[0][j]['Alb'], ground_features[0][j]['Bo'],
                                   ground_features[0][j]['Zo']))
            wmet3.close()
        else:
            wmet3.write('  FREQ_SECT SEASONAL  {0}\n'.format(len_sd))
            for i in range(len_sd - 1):
                wmet3.write('  SECTOR    {0}    {1}    {2}\n'.format(i + 1, seperate_degree[i], seperate_degree[i + 1]))
            wmet3.write(
                '  SECTOR    {0}    {1}    {2}\n'.format(len_sd, seperate_degree[len_sd - 1], seperate_degree[0]))
            for j in range(4):
                for i in range(len_sd):
                    wmet3.write('  SITE_CHAR {0}    {1}    {2}    {3}    {4}\n'.
                                format(j+1, i+1, ground_features[i][j]['Alb'], ground_features[i][j]['Bo'],
                                       ground_features[i][j]['Zo']))
            wmet3.close()
    # TODO:地面时间周期按月时的地表特征参数表
    else:
        if len_sd < 2:
            wmet3.write('  FREQ_SECT MONTHLY  1\n')
            wmet3.write('  SECTOR    1    0   360\n')
            for j in range(12):
                if j < 2:
                    frequency_index = 1
                elif 1 < j < 5:
                    frequency_index = 2
                elif 4 < j < 8:
                    frequency_index = 3
                elif 7 < j < 11:
                    frequency_index = 4
                else:
                    frequency_index = 1
                wmet3.write('  SITE_CHAR {0}    1    {1}    {2}    {3}\n'.
                            format(frequency_index, ground_features[0][j]['Alb'], ground_features[0][j]['Bo'],
                                   ground_features[0][j]['Zo']))
            wmet3.close()
        else:
            wmet3.write('  FREQ_SECT MONTHLY  {0}\n'.format(len_sd))
            for i in range(len_sd - 1):
                wmet3.write('  SECTOR    {0}    {1}    {2}\n'.format(i + 1, seperate_degree[i], seperate_degree[i + 1]))
            wmet3.write('  SECTOR    {0}    {1}    {2}\n'.
                        format(len_sd, seperate_degree[len_sd - 1], seperate_degree[0]))
            for j in range(12):
                if j < 2:
                    frequency_index = 1
                elif 1 < j < 5:
                    frequency_index = 2
                elif 4 < j < 8:
                    frequency_index = 3
                elif 7 < j < 11:
                    frequency_index = 4
                else:
                    frequency_index = 1
                for i in range(len_sd):
                    wmet3.write('  SITE_CHAR {0}    {1}    {2}    {3}    {4}\n'.
                                format(frequency_index, i + 1, ground_features[i][j]['Alb'],
                                       ground_features[i][j]['Bo'], ground_features[i][j]['Zo']))
            wmet3.close()
    time.sleep(1)
    os.system('rename aermet.inp met2.inp')
    time.sleep(1)
    os.system('rename met3.inp aermet.inp')
    os.system('aermet.exe')
    print('气象方案{0}, 运行完成'.format(1))


def run_bpip(jsonObj, select_project, select_sources_id, terrain_path):
    # TODO：假定每个emissionCase中的每个点源都参与建筑下洗计算
    # TODO： 在每个emissionCase运行过程中先运行Bpipprm.exe，再运行aermod.exe
    consider_downwash = select_project['consider_downwash']
    print(consider_downwash)
    if not consider_downwash:
        print("本案例不考虑建筑下洗")
        return
    else:
        print("本案例考虑建筑下洗")
        build_project_id = select_project['select_building_ID']
        index = get_id(jsonObj['project_data']['basedata']['buildings_project'], build_project_id, key='ID')
        select_bulid = jsonObj['project_data']['basedata']['buildings_project'][index]
        # TODO:创建bpipprm.inp文件
        bpi = open('buliding.inp', 'w', encoding='UTF-8')
        bpi.write("'建筑物下洗方案{0}'\n".format(build_project_id + 1))
        bpi.write("'p'\n")
        bpi.write("'METERS' 1.0\n")  # 单位和缩放比例，保持不变
        bpi.write("'UTMY' 0.00\n")
        num = len(select_bulid['bulidings'])
        bpi.write(str(num) + '\n')
        floors_num = 1
        for i, build in enumerate(select_bulid['bulidings']):
            # 取第一个点的地形高度做为建筑物整体的高度
            building_base_height = get_elevation_h(terrain_path, build['location'][0][0], build['location'][0][1])
            if jsonObj['project_data']['AERMOD']['AERMOD_project'][0]['consider_terrain'] == 1:  # 如果不考虑地形，基准高度设置为0
                building_base_height = 0

            bpi.write("'{0}' {1} {2:.1f}\n".format(build['name'], floors_num, building_base_height))
            nverts = len(build['location'])-1
            bpi.write('{0} {1:.1f}\n'.format(nverts, build['building_height']))
            for n in range(nverts):
                lonlat = build['location'][n]
                zone, z_utm_lon, z_utm_lat = wgs84toutm(lonlat[0], lonlat[1])
                bpi.write('{0:.1f} {1:.1f}\n'.format(z_utm_lon, z_utm_lat))

        # TODO: 只有点源才考虑建筑下洗; 假定每个emissionCase中的每个点源都参与建筑下洗计算
        # 确定点源个数
        point_source_num = 0
        for j, ID in enumerate(select_sources_id):
            index = get_id(jsonObj['project_data']['basedata']['sources'], ID, key='ID')
            source_detail = jsonObj['project_data']['basedata']['sources'][index]
            if source_detail['type'] == 'POINT':
                point_source_num = point_source_num + 1
        bpi.write(str(point_source_num) + '\n')

        # 写入每个点源的名称，源的地基高度，烟囱高度，相对坐标
        for j, ID in enumerate(select_sources_id):
            index = get_id(jsonObj['project_data']['basedata']['sources'], ID, key='ID')
            source_detail = jsonObj['project_data']['basedata']['sources'][index]
            if source_detail['type'] == 'POINT':
                point_source_num = point_source_num + 1
                source_x, source_y, base_height = cal_location_source(jsonObj, source_detail, 'POINT', terrain_path)
                bpi.write("'{0}' {1:.1f} {2:.1f} {3:.1f} {4:.1f}\n".
                          format(str(source_detail['ID']), base_height, source_detail['Stkhgt'], source_x, source_y))
        bpi.close()
        # TODO：运行bpipprm,提取*.out文件SO字段内容至*.so文件
        if point_source_num > 0.999:
            os.system('Bpipprm.exe buliding.inp building.out building.sum')
            bpout = open('building.out', 'r', encoding='UTF-8')
            bpso = open('build.BSO', 'w', encoding='UTF-8')
            for line in bpout.readlines():
                if 'SO ' in line:
                    bpso.write(line.lstrip())
            bpout.close()
            bpso.close()
            print('建筑下洗方案{0}, bpipprm运行完成'.format(build_project_id))
        else:
            print('建筑下洗方案{0}, 点源个数小于1'.format(build_project_id))


def aermod_deposition_gas(select_aermod_project):
    str_MODELOPT = ''
    str_GASDEPVD = ''
    str_GDSEASON = ''
    str_GDLANUSE = ''

    if select_aermod_project['Is_Direct_gasDryDepostionValue']:
        # 如果选择直接给入气体干沉降， 只设置 DDEP
        str_MODELOPT = ' ' + 'DDEP'
    else:
        if select_aermod_project['consider_TotalDepostion']:
            str_MODELOPT = str_MODELOPT + ' ' + 'DEPOS'
        if select_aermod_project['consider_DryDepostion']:
            str_MODELOPT = str_MODELOPT + ' ' + 'DDEP'
        if select_aermod_project['consider_WetDepostion']:
            str_MODELOPT = str_MODELOPT + ' ' + 'WDEP'
    if select_aermod_project['Is_Direct_gasDryDepostionValue']:
        gasDepValue = select_aermod_project['gasDepostionValue']
        str_GASDEPVD = '   GASDEPVD {0:f}'.format(gasDepValue) + '\n'
    elif select_aermod_project['consider_TotalDepostion'] or select_aermod_project['consider_DryDepostion']:
        GDSEASON = '   GDSEASON'
        for b in select_aermod_project['gdseasonMonth_value']:
            GDSEASON = GDSEASON + ' ' + str(b)
        GDSEASON += '\n'
        str_GDSEASON = GDSEASON

        GDLANUSE = '   GDLANUSE'
        for c in select_aermod_project['gdlandUseSector_type']:
            GDLANUSE = GDLANUSE + ' ' + str(c)
        GDLANUSE += '\n'
        str_GDLANUSE = GDLANUSE

    return str_MODELOPT, str_GASDEPVD, str_GDSEASON, str_GDLANUSE


def aermod_NO2(pollutant_detail, select_aermod_project):
    str_MODELOPT = ''
    str_NO2EQUIL = ''
    str_NO2STACK = ''
    str_OZONEVAL = ''
    str_OZONunit = ''
    pollutant_name = pollutant_detail['name_en']
    if pollutant_name == 'NO2' and select_aermod_project['consider_NO2']:
        str_MODELOPT = ' ' + select_aermod_project['NO2_conversion_method']
        str_POLLUTID = '   POLLUTID {0}\n'.format('NO2')
    else:
        # TODO：#暂时将不是NO2的，全部设置为OTHER，注意与AVERTIME有关系
        str_POLLUTID = '   POLLUTID {0}\n'.format('OTHER')

    if pollutant_name == 'NO2' and select_aermod_project['consider_NO2'] and \
            select_aermod_project['NO2_conversion_method'] == 'OLM':
        # TODO: NO2EQUIL 字段为可选字段，界面暂时取消，设置ambient_equil_NO2_NOx_ratio这里采用内部默认值0.1
        # str_NO2EQUIL = '   NO2EQUIL {0}\n'.format(
        #     select_aermod_project['NO2_conversion_para']['OLM_para']['ambient_equil_NO2_NOx_ratio'])
        str_NO2STACK = '   NO2STACK {0}\n'.\
            format(select_aermod_project['NO2_conversion_para']['OLM_para']['in_stack_NO2_NOx_ratio'])
        # TODO:臭氧的浓度单位写死，ug/m3
        str_OZONEVAL = '   OZONEVAL {0} UG/M3\n'.\
            format(select_aermod_project['NO2_conversion_para']['OLM_para']['background_O3'])
    elif pollutant_name == 'NO2' and select_aermod_project['consider_NO2'] and \
            select_aermod_project['NO2_conversion_method'] == 'PVMRM':
        # TODO: 设置ambient_equil_NO2_NOx_ratio
        str_NO2EQUIL = '   NO2EQUIL {0}\n'.\
            format(select_aermod_project['NO2_conversion_para']['PVMRM_para']['ambient_equil_NO2_NOx_ratio'])
        str_NO2STACK = '   NO2STACK {0}\n'.\
            format(select_aermod_project['NO2_conversion_para']['PVMRM_para']['in_stack_NO2_NOx_ratio'])
        # TODO:臭氧的浓度单位写死，ug/m3
        str_OZONEVAL = '   OZONEVAL {0} UG/M3\n'.\
            format(select_aermod_project['NO2_conversion_para']['PVMRM_para']['background_O3'])
    return str_MODELOPT, str_NO2EQUIL, str_NO2STACK, str_OZONEVAL, str_POLLUTID, str_OZONunit


def cal_location_source(jsonObj, source_detail, type, terrain_path):
    origal_lon = jsonObj['project_location']['lng']
    origal_lat = jsonObj['project_location']['lat']
    # TODO: 计算相对坐标和排放源的海拔
    if type in ['POINT', 'VOLUME', 'AREACIRC']:
        # 取这三种源的中心
        try:
            source_lon = source_detail['locationX']
            source_lat = source_detail['locationY']
        except KeyError:
            source_lon = source_detail['location'][0]  # POINT 源的json经纬度格式不同
            source_lat = source_detail['location'][1]
        source_x, source_y = relative_distance_coord(origal_lon, origal_lat, np.mean(source_lon),
                                                     np.mean(source_lat))
    elif type in ['AREA', 'AREAPOLY', 'OPENPIT']:
        # 取这三种源的第一点，作为LOCATION的输入
        source_lon = source_detail['locationX']
        source_lat = source_detail['locationY']
        source_x, source_y = relative_distance_coord(origal_lon, origal_lat, source_lon[0],
                                                     source_lat[0])
    elif type in ['LINE', 'RLINE', 'RLINEXT', 'BUOYLINE']:
        # 对线源而言，取其一个端点。界面用户提交的线源是多跟线组成的
        source_lon = source_detail['locationX']
        source_lat = source_detail['locationY']
        source_x = []
        source_y = []
        for i in range(len(source_lon)):
            source_x1, source_y1 = relative_distance_coord(origal_lon, origal_lat, source_lon[i],
                                                           source_lat[i])
            source_x.append(source_x1)
            source_y.append(source_y1)

    base_height = get_elevation_h(terrain_path, np.mean(source_lon), np.mean(source_lat))
    if jsonObj['project_data']['AERMOD']['AERMOD_project'][0]['consider_terrain'] == 1:  # 如果不考虑地形，基准高度设置为0
        base_height = 0

    return source_x, source_y, base_height


def aermod_source_input(w, jsonObj, select_aermod_project, select_sources_id, pollutant_detail, terrain_path):
    """
    :param w:
    :param jsonObj:
    :param select_sources_id:
    :param pollutant_name:
    :param str_for_deposition: 对气体而言，它为 str_GASDEPOS; 对颗粒物而言，它为一块字符
    :return:
    """
    origal_lon = jsonObj['project_location']['lng']
    origal_lat = jsonObj['project_location']['lat']
    pollutant_name = pollutant_detail['name_en']
    w.write('SO STARTING\n')
    point_source_num = 0
    for j, ID in enumerate(select_sources_id):
        index = get_id(jsonObj['project_data']['basedata']['sources'], ID, key='ID')
        source_detail = jsonObj['project_data']['basedata']['sources'][index]
        if source_detail['type'] == 'POINT':
            point_source_num = point_source_num + 1
            source_x, source_y, base_height = cal_location_source(jsonObj, source_detail, 'POINT', terrain_path)
            Vlemis = (source_detail['pollutants'][pollutant_name]) / 3.6
            w.write('   LOCATION {0} {1} {2:.1f} {3:.1f} {4:.1f}\n'.
                    format(source_detail['ID'], source_detail['type'], source_x, source_y, base_height))
            # base_height为可选项，可不填
            w.write('   SRCPARAM {0} {1:.6f} {2:.1f} {3:.1f} {4:.1f} {5:.1f}\n'.
                    format(source_detail['ID'], Vlemis, source_detail['Stkhgt'], source_detail['Stktmp'] + 273.15,
                           source_detail['Stkvel'], source_detail['Stkdia']))
            if not pollutant_detail['is_gaseous']:
                deposition_str = deposition_particle_input(select_aermod_project,pollutant_detail,str(source_detail['ID']))
            else:
                deposition_str = deposition_gas_input(select_aermod_project, pollutant_detail, str(source_detail['ID']))
            w.write(deposition_str)
        elif source_detail['type'] == 'VOLUME':
            source_x, source_y, base_height = cal_location_source(jsonObj, source_detail, 'VOLUME', terrain_path)
            Vlemis = (source_detail['pollutants'][pollutant_name])/3.6
            Relhgt = source_detail['sourceVerticalLength'] / 2.0  # 排放高度
            Szinit = source_detail['sourceVerticalLength'] / 2.15  # 垂直初始长度，暂时不考虑不在地表的elevated sources
            source_area = get_source_area(source_detail)
            Yinit = np.sqrt(source_area)  # # 侧边初始长度
            print("体源的初始侧边长度：", Yinit)
            print("----- source_type:VOLUME")
            w.write('   LOCATION {0} {1} {2:.1f} {3:.1f} {4:.1f}\n'.
                    format(source_detail['ID'], source_detail['type'], source_x, source_y, base_height))
            w.write('   SRCPARAM {0} {1:.6f} {2:.1f} {3:.1f} {4:.1f} \n'.
                    format(source_detail['ID'], Vlemis, Relhgt, Yinit, Szinit))
            if not pollutant_detail['is_gaseous']:
                deposition_str = deposition_particle_input(select_aermod_project, pollutant_detail, str(source_detail['ID']))
            else:
                deposition_str = deposition_gas_input(select_aermod_project, pollutant_detail, str(source_detail['ID']))
            w.write(deposition_str)

        elif source_detail['type'] == 'AREA':
            source_x, source_y, base_height = cal_location_source(jsonObj, source_detail, 'AREA', terrain_path)
            source_area = get_source_area(source_detail)
            Aremis = (source_detail['pollutants'][pollutant_name])/(3.6*source_area)
            w.write('   LOCATION {0} {1} {2:.1f} {3:.1f} {4:.1f}\n'.format(source_detail['ID'], source_detail['type'],
                                                                           source_x, source_y, base_height))
            w.write('   SRCPARAM {0} {1:.6f} {2:.1f} {3:.1f} {4:.1f} {5:.1f}\n'.
                    format(source_detail['ID'], Aremis, source_detail['Relhgt'], source_detail['Xinit'],
                           source_detail['Yinit'], source_detail['Angle'], source_detail['Szinit']))
            if not pollutant_detail['is_gaseous']:
                deposition_str = deposition_particle_input(select_aermod_project, pollutant_detail, str(source_detail['ID']))
            else:
                deposition_str = deposition_gas_input(select_aermod_project, pollutant_detail, str(source_detail['ID']))
            w.write(deposition_str)

        elif source_detail['type'] == 'AREAPOLY':
            source_x, source_y, base_height = cal_location_source(jsonObj, source_detail, 'AREAPOLY', terrain_path)
            source_area = get_source_area(source_detail)
            print(source_detail)
            print(pollutant_name)
            Aremis = (source_detail['pollutants'][str(pollutant_name)])/(3.6*source_area)
            nverts = len(source_detail['locationX'])-1
            w.write('   LOCATION {0} {1} {2:.1f} {3:.1f} {4:.1f}\n'.format(source_detail['ID'], source_detail['type'],
                                                                           source_x, source_y, base_height))
            if source_detail['Relhgt'] == "":
                source_detail['Relhgt'] = source_detail['Szinit']
                w.write('   SRCPARAM {0} {1:.6f} {2:.1f} {3} {4:.1f} \n'.
                        format(source_detail['ID'], Aremis, source_detail['Relhgt'], nverts, source_detail['Szinit']))
            else:
                w.write('   SRCPARAM {0} {1:.6f} {2:.1f} {3} {4:.1f} \n'.
                        format(source_detail['ID'], Aremis, source_detail['Relhgt'], nverts, source_detail['Szinit']))
            areavert = '   AREAVERT {0}'.format(source_detail['ID'])
            for point in range(nverts):
                point_x, point_y = relative_distance_coord(origal_lon, origal_lat, source_detail['locationX'][point],
                                                           source_detail['locationY'][point])
                areavert = areavert + ' ' + '{0:.1f}'.format(point_x) + ' ' + '{0:.1f}'.format(point_y)
            areavert = areavert + '\n'
            w.write(areavert)
            if not pollutant_detail['is_gaseous']:
                deposition_str = deposition_particle_input(select_aermod_project, pollutant_detail, str(source_detail['ID']))
            else:
                deposition_str = deposition_gas_input(select_aermod_project, pollutant_detail, str(source_detail['ID']))
            w.write(deposition_str)

        elif source_detail['type'] == 'AREACIRC':
            source_x, source_y, base_height = cal_location_source(jsonObj, source_detail, 'AREACIRC', terrain_path)
            source_area = get_source_area(source_detail)
            Aremis = (source_detail['pollutants'][str(pollutant_name)])/(3.6*source_area)
            w.write('   LOCATION {0} {1} {2:.1f} {3:.1f} {4:.1f}\n'.
                    format(source_detail['ID'], source_detail['type'], source_x, source_y, base_height))
            w.write('   SRCPARAM {0} {1:.6f} {2:.1f} {3:.1f} {4:.1f} \n'.
                    format(source_detail['ID'], Aremis, source_detail['Relhgt'], source_detail['Radius'],
                           source_detail['Szinit']))
            if not pollutant_detail['is_gaseous']:
                deposition_str = deposition_particle_input(select_aermod_project, pollutant_detail, str(source_detail['ID']))
            else:
                deposition_str = deposition_gas_input(select_aermod_project, pollutant_detail, str(source_detail['ID']))
            w.write(deposition_str)

        elif source_detail['type'] == 'LINE':
            # todo:多段线源以每段的一个端点x,y表示location
            # TODO: 每段线要循环写 , 注意线源排放速率单位为g/s*m2
            source_x, source_y, base_height = cal_location_source(jsonObj, source_detail, 'LINE', terrain_path)
            Aremis = (source_detail['pollutants'][str(pollutant_name)]) / 3.6
            for i in range(len(source_x) - 1):
                id_str = str(source_detail['ID']) + "_" + str(i)
                w.write('   LOCATION {0} {1} {2:.1f} {3:.1f} {4:.1f} {5:.1f} {6:.1f}\n'.
                        format(id_str, source_detail['type'], source_x[i], source_y[i], source_x[i + 1],
                               source_y[i + 1], base_height))
                w.write('   SRCPARAM {0} {1:.6f} {2:.1f} {3:.1f} {4:.1f}\n'.
                        format(id_str, Aremis, source_detail['Relhgt'], source_detail['Width'], source_detail['Szinit']))
                if not pollutant_detail['is_gaseous']:
                    deposition_str = deposition_particle_input(select_aermod_project, pollutant_detail, id_str)
                else:
                    deposition_str = deposition_gas_input(select_aermod_project, pollutant_detail, id_str)
                w.write(deposition_str)  # TODO: 每一个小线源污染物沉降的方法字段就得写一遍

        elif source_detail['type'] == 'BUOYLINE':
            # todo:多段线源以每段的一个端点x,y表示location
            # TODO: 每段线要循环写
            source_x, source_y, base_height = cal_location_source(jsonObj, source_detail, 'BUOYLINE', terrain_path)
            Blemis = (source_detail['pollutants'][str(pollutant_name)]) / 3.6
            for i in range(len(source_x) - 1):
                id_str = str(source_detail['ID']) + "_" + str(i)
                w.write('   LOCATION {0} {1} {2:.1f} {3:.1f} {4:.1f} {5:.1f} {6:.1f}\n'.
                        format(id_str, source_detail['type'], source_x[i], source_y[i], source_x[i+1], source_y[i+1],
                               base_height))
                w.write('   SRCPARAM {0} {1:.6f} {2} \n'.format(id_str, Blemis, source_detail['Relhgt']))
                if not pollutant_detail['is_gaseous']:
                    deposition_str = deposition_particle_input(select_aermod_project, pollutant_detail, id_str)
                else:
                    deposition_str = deposition_gas_input(select_aermod_project, pollutant_detail, id_str)
                w.write(deposition_str)
            w.write('   BLPINPUT {0} {1} {2} {3} {4} {5}\n'.
                    format(source_detail['Blavgblen'], source_detail['Blavgbhgt'], source_detail['Blavgbwid'],
                           source_detail['Blavglwid'], source_detail['Blavgbsep'], source_detail['Blavgfprm']))
    if select_aermod_project['consider_urban']:
        w.write('   URBANSRC ALL\n')
    if select_aermod_project['consider_downwash'] and point_source_num > 0.999:
        # TODO: 当点源个数大于等于1时，才考虑建筑
        w.write('   INCLUDED build.BSO\n')
    w.write('   SRCGROUP ALL\n')
    w.write('SO FINISHED\n')


def aermod_RE_ME_OU_input(w, select_aermod_project, aermet_project, terrain_path):
    # TODO: RE字段-------------------------------------------------------
    w.write('RE STARTING\n')
    if select_aermod_project['consider_terrain'] != 1:   # 1，表示不考虑地形
        w.write('   INCLUDED ..\\..\\AERMAP\\AERMAP.REC\n')
    else:
        w.write('   INCLUDED ..\\..\\AERMAP\\grid_location.RE\n')
    w.write('RE FINISHED\n')

    # TODO: ME字段-------------------------------------------------------
    w.write('ME STARTING\n')
    w.write('   SURFFILE MET3.SFC\n')
    w.write('   PROFFILE MET3.PFL\n')
    w.write('   SURFDATA {0} {1}\n'.format(aermet_project['select_surface_data']['surface_sites'][0]['site_id'],
                                           aermet_project['select_surface_data']['surface_sites'][0]['data_year']))
    w.write('   UAIRDATA {0} {1}\n'.format(aermet_project['select_upair_data']['upair_sites'][0]['site_id'],
                                           aermet_project['select_upair_data']['upair_sites'][0]['data_year']))
    # 目前探空OQA数据站点编号均为99999

    # TODO: 假设aermet_project下面只有一个地面站和探空站。
    site_lon = aermet_project['select_surface_data']['surface_sites'][0]['site_lng']
    site_lat = aermet_project['select_surface_data']['surface_sites'][0]['site_lat']
    site_elevaion = get_elevation_h(terrain_path, site_lon, site_lat)
    if select_aermod_project['consider_terrain'] == 1:  # 如果不考虑地形，地形高度设置为0
        site_elevaion = 0
    w.write('   PROFBASE {0} METERS\n'.format(site_elevaion))
    w.write('ME FINISHED\n')
    # TODO: OU字段-----------------------------------------------------------
    w.write('OU STARTING\n')
    RECTABLE = '   RECTABLE ALLAVE'
    # TODO: output_high_value 可能为一个整数，可能是一个字符串多个数据，如'1,8,19'
    if type(select_aermod_project['output_high_value']) == int:
        high_value = select_aermod_project['output_high_value']
        RECTABLE = RECTABLE + ' ' + str(high_value)
    else:
        high_value = select_aermod_project['output_high_value'].split(",")
        for z in high_value:
            RECTABLE = RECTABLE + ' ' + str(z)
    RECTABLE += '\n'
    # w.write(RECTABLE)
    # w.write('   DAYTABLE ALLAVE\n')
    # todo: 输出各平均时间的逐步值文件用于叠加方案计算，小时文件过大且非必须，暂不输出
    avgtime = select_aermod_project['aermod_meanTime']
    print(avgtime)
    if type(select_aermod_project['output_high_value']) == int:
        high_value = str(select_aermod_project['output_high_value']).split(",")
    else:
        high_value = select_aermod_project['output_high_value'].split(",")
    if avgtime == 'PERIOD':
        # w.write('   POSTFILE PERIOD ALL PLOT all_PERIOD.plt\n')
        # w.write('   PLOTFILE PERIOD all average_period.txt\n')
        pass
    elif type(avgtime) == int:
        if avgtime == '24':
            w.write('   POSTFILE 24 ALL PLOT all_24.plt\n')
            for h in high_value:
                w.write('   PLOTFILE 24 all {0} 24_{1}.txt\n'.format(h, h))
        else:
            for c in high_value:
                # w.write('   PLOTFILE {0} all {1} {2}_{3}.txt\n'.format(avgtime, c, avgtime, c))
                w.write('   POSTFILE 1 ALL PLOT all_1.plt\n')
    else:
        meantime = select_aermod_project['aermod_meanTime'].split(",")
        for a in meantime:
            if a == 'PERIOD':
                w.write('   POSTFILE PERIOD ALL PLOT all_PERIOD.plt\n')
                w.write('   PLOTFILE PERIOD all average_period.txt\n')
            elif a == '24':
                w.write('   POSTFILE 24 ALL PLOT all_24.plt\n')
                for h in high_value:
                    w.write('   PLOTFILE 24 all {0} 24_{1}.txt\n'.format(h, h))
            else:
                for c in high_value:
                    # w.write('   PLOTFILE {0} all {1} {2}_{3}.txt\n'.format(a, c, a, c))
                    w.write('   POSTFILE 1 ALL PLOT all_1.plt\n')
    w.write('OU FINISHED\n')


def run_aermod_gas(jsonObj, select_aermod_project, aermet_project, pollutant_detail, select_sources_id, terrain_path):
    # Todo:创建aermod.inp
    w = open('aermod.inp', 'w', encoding='UTF-8')
    # TODO: CO字段-------------------------------------------------------
    w.write('CO STARTING\n')
    w.write('   TITLEONE {0}\n'.format(select_aermod_project['ID']))
    MODELOPT = '   MODELOPT CONC'
    # TODO: 气体沉降逻辑
    str_MODELOPT, str_GASDEPVD, str_GDSEASON, str_GDLANUSE\
        = aermod_deposition_gas(select_aermod_project)
    MODELOPT = MODELOPT + str_MODELOPT
    # TODO： NO2转化逻辑
    str_MODELOPT, str_NO2EQUIL, str_NO2STACK, str_OZONEVAL, str_POLLUTID,str_OZONunit \
        = aermod_NO2(pollutant_detail, select_aermod_project)
    MODELOPT = MODELOPT + str_MODELOPT
    if select_aermod_project['consider_terrain'] == 1:  # 1,表示不考虑地形
        MODELOPT = MODELOPT + ' ' + 'FLAT'
    else:
        MODELOPT = MODELOPT + ' ' + 'ELEV'
    # TODO: ALPHA选项，小风和干湿沉降时需要;但是似乎一直保留，也没有关系
    MODELOPT = MODELOPT + ' ' + 'ALPHA'  # 启用干湿沉降时需要
    MODELOPT += '\n'
    w.write(MODELOPT)

    AVERTIME = '   AVERTIME'
    # TODO: to test  "aermod_meanTime":[1,3,8,24,"PERIOD"]
    avgtime = select_aermod_project['aermod_meanTime']
    if avgtime == 'PERIOD':
        AVERTIME = AVERTIME + ' PERIOD\n'
    elif type(avgtime) == int:
        AVERTIME = AVERTIME + ' ' + avgtime + '\n'
    else:
        avgtime = select_aermod_project['aermod_meanTime'].split(",")
        for a in avgtime:
            AVERTIME = AVERTIME + ' ' + a
        AVERTIME += '\n'
    w.write(AVERTIME)

    #TODO: 指定离地高度
    above_ground_height = select_aermod_project['above_ground_height']
    flagpole = '   FLAGPOLE' + ' ' + str(above_ground_height) + '\n'
    w.write(flagpole)


    # TODO: 默认设置城市粗糙度为0.5m，同时设定SO URBANSRC ALL 指定所有源都考虑城市效应
    # TODO： 只考虑城市区域为1个情况；后续可以完善考虑多个城市区域的情况
    # 注意：Beginning with version 19191, the urban modeling option has been incorporated as an ALPHA option for
    # the BUOYLINE, RLINE, and RLINEXT source types
    if select_aermod_project['consider_urban']:
        w.write('   URBANOPT {0}\n'.format(select_aermod_project['urban_pop'] * 10000))
        # w.write('   URBANOPT {0} {1} {2}\n'.format(select_aermod_project['urban_pop']*10000, 'urbanName', 0.5))
        # 城市名和粗糙度为可选项，version 09292开始，粗糙度默认为1m，不为1m的粗超度值均被当作非默认选项处理
    w.write(str_POLLUTID)
    if pollutant_detail['name_en'] == "SO2" or pollutant_detail['name_en']=='so2':
        half_life = 14400  #默认为4个小时
    else:
        half_life = 999999999  #默认为4个小时
    w.write('   HALFLIFE {0}\n'.format(half_life))
    w.write(str_GASDEPVD)
    w.write(str_GDSEASON)
    w.write(str_GDLANUSE)
    w.write(str_NO2EQUIL)
    w.write(str_NO2STACK)
    w.write(str_OZONEVAL)
    w.write(str_OZONunit)
    w.write('   RUNORNOT RUN\n')

    # TODO:用于事件处理的输入文件
    # w.write('   EVENTFIL event.inp\n')
    w.write('   ERRORFIL ERRORS.OUT\n')
    w.write('CO FINISHED\n')

    # TODO： SO字段------------
    aermod_source_input(w, jsonObj, select_aermod_project, select_sources_id, pollutant_detail, terrain_path)
    # TODO： RE ME OU字段------
    aermod_RE_ME_OU_input(w, select_aermod_project, aermet_project, terrain_path)

    w.close()
    subprocess.Popen(["start", "aermod.exe"], stdout=subprocess.PIPE,
                     stderr=subprocess.PIPE, shell=True, universal_newlines=True)
    # os.system('aermod.exe')


def deposition_gas_input(select_aermod_project, pollutant_detail, id_str):
    # TODO: 参考手册《Table 3-1 Summary of Deposition Options》，暂时定为只要有其中之一的沉降关键字，求给出颗粒物method字段
    # TODO：所有源涉及到气体沉降的SO 的语句也必须跟LOCATION一致
    # TODO： 干空气沉降直接设定沉降速度的逻辑在外部实现。
    str_GASDEPOS = ''
    if ((select_aermod_project['Is_Direct_gasDryDepostionValue'] == 0) and
            (select_aermod_project['consider_TotalDepostion'] or select_aermod_project['consider_DryDepostion'] or
             select_aermod_project['consider_WetDepostion'])):
        # 如果是气体，没有直接给入气体干沉降速度，同时考虑了总沉降或干沉降或湿沉降，就需要设置GASDEPOS
        a1 = pollutant_detail['AERMOD_gas_subsidence_parameters']['diffusion_in_air']
        a2 = pollutant_detail['AERMOD_gas_subsidence_parameters']['diffusion_in_water']
        a3 = pollutant_detail['AERMOD_gas_subsidence_parameters']['rebound_resistance']
        a4 = pollutant_detail['AERMOD_gas_subsidence_parameters']['Henry_constant']
        # print("diffusion_in_water:, -----------------------", a2)
        if a1==0 or a1=='':
            a1 = 0.1053
        if a2<=1 or a2=='':
            a2 = 1.285E5
        if a3==0 or a3=='':
            a3 = 745
        if a4==0 or a4=='':
            a4 = 1750
        str_GASDEPOS = str_GASDEPOS + '   GASDEPOS {0} {1} {2} {3} {4}\n'.\
            format(id_str, a1, a2, a3, a4)
    return str_GASDEPOS


def deposition_particle_input(select_aermod_project, pollutant_detail, id_str):
    # TODO: 参考手册《Table 3-1 Summary of Deposition Options》，暂时定为只要有其中之一的沉降关键字，求给出颗粒物method字段
    # TODO：所有源涉及到颗粒物沉降的SO METHOD 沉降有关的语句也必须跟LOCATION一致
    str_PM_depostions = ''
    if select_aermod_project['consider_TotalDepostion'] or select_aermod_project['consider_DryDepostion'] or \
            select_aermod_project['consider_WetDepostion']:
        if pollutant_detail['PM_grain_size'][0]['type'] == 'methodOne':
            PARTDIAM = '   PARTDIAM' + ' ' + id_str
            for a in pollutant_detail['PM_grain_size'][0]['methodOne']['D50']:
                PARTDIAM = PARTDIAM + ' ' + '{0:.2f}'.format(a)
            str_PM_depostions = str_PM_depostions + PARTDIAM + '\n'

            MASSFRAX = '   MASSFRAX' + ' ' + id_str
            for a in pollutant_detail['PM_grain_size'][0]['methodOne']['percentage']:
                MASSFRAX = MASSFRAX + ' ' + '{0:.2f}'.format(a/100)
            str_PM_depostions = str_PM_depostions + MASSFRAX + '\n'

            PARTDENS = '   PARTDENS' + ' ' + id_str
            for a in pollutant_detail['PM_grain_size'][0]['methodOne']['true_density']:
                PARTDENS = PARTDENS + ' ' + '{0:.2f}'.format(a)
            str_PM_depostions = str_PM_depostions + PARTDENS + '\n'

        else:
            str_PM_depostions = str_PM_depostions + '   METHOD_2 {0} {1} {2}\n'.\
                format(id_str, pollutant_detail['PM_grain_size'][0]['methodTwo']['mass_percent'],
                       pollutant_detail['PM_grain_size'][0]['methodTwo']['average_size'])
    return str_PM_depostions


def run_aermod_particle(jsonObj,select_aermod_project,aermet_project,pollutant_detail,select_sources_id,terrain_path):
    # Todo:创建aermod.inp
    w = open('aermod.inp', 'w', encoding='UTF-8')
    # TODO: CO字段-------------------------------------------------------
    w.write('CO STARTING\n')
    w.write('   TITLEONE {0}\n'.format(select_aermod_project['ID']))
    MODELOPT = '   MODELOPT CONC'
    # TODO: 颗粒物沉降逻辑
    str_MODELOPT = ''
    if select_aermod_project['consider_TotalDepostion']:
        str_MODELOPT = str_MODELOPT + ' ' + 'DEPOS'
    if select_aermod_project['consider_DryDepostion']:
        str_MODELOPT = str_MODELOPT + ' ' + 'DDEP'
    if select_aermod_project['consider_WetDepostion']:
        str_MODELOPT = str_MODELOPT + ' ' + 'WDEP'
    MODELOPT = MODELOPT + str_MODELOPT
    if select_aermod_project['consider_terrain'] == 1:
        MODELOPT = MODELOPT + ' ' + 'FLAT'
    else:
        MODELOPT = MODELOPT + ' ' + 'ELEV'
    # TODO: ALPHA选项，小风和干湿沉降时需要;但是似乎一直保留，也没有关系
    MODELOPT = MODELOPT + ' ' + 'ALPHA'  # 启用干湿沉降时需要
    MODELOPT += '\n'
    w.write(MODELOPT)

    AVERTIME = '   AVERTIME'
    # TODO: to test  "aermod_meanTime":[1,3,8,24,"PERIOD"]
    avgtime = select_aermod_project['aermod_meanTime']
    if avgtime == 'PERIOD':
        AVERTIME = AVERTIME + ' PERIOD\n'
    elif type(avgtime) == int:
        AVERTIME = AVERTIME + ' ' + avgtime + '\n'
    else:
        for a in avgtime:
            AVERTIME = AVERTIME + ' ' + a
        AVERTIME += '\n'
    w.write(AVERTIME)

    # TODO: 指定离地高度
    above_ground_height = select_aermod_project['above_ground_height']
    flagpole = '   FLAGPOLE' + ' ' + str(above_ground_height) + '\n'
    w.write(flagpole)

    # TODO: 默认设置城市粗糙度为0.5m，同时设定SO URBANSRC ALL 指定所有源都考虑城市效应
    # TODO： 只考虑城市区域为1个情况；后续可以完善考虑多个城市区域的情况
    # 注意：Beginning with version 19191, the urban modeling option has been incorporated as an ALPHA option for
    # the BUOYLINE, RLINE, and RLINEXT source types
    if select_aermod_project['consider_urban']:
        w.write('   URBANOPT {0}\n'.format(select_aermod_project['urban_pop'] * 10000))
        # w.write('   URBANOPT {0} {1} {2}\n'.format(select_aermod_project['urban_pop']*10000,
        #                                            'urbanName', 0.5))
        # 城市名和粗糙度为可选项，粗糙度默认为1m，version 09292开始，不为1m的粗超度值均被当作非默认选项处理

    w.write('   POLLUTID {0}\n'.format('OTHER'))
    # w.write('   HALFLIFE {0}\n'.format(pollutant_detail['half_life']))
    w.write('   RUNORNOT RUN\n')

    # TODO:用于事件处理的输入文件
    # w.write('   EVENTFIL event.inp\n')
    w.write('   ERRORFIL ERRORS.OUT\n')
    w.write('CO FINISHED\n')

    # TODO： SO字段------------
    aermod_source_input(w, jsonObj, select_aermod_project, select_sources_id, pollutant_detail, terrain_path)
    # TODO： RE ME OU字段------
    aermod_RE_ME_OU_input(w, select_aermod_project, aermet_project, terrain_path)

    w.close()
    # os.system('aermod.exe')
    subprocess.Popen(["start", "aermod.exe"], stdout=subprocess.PIPE,
                     stderr=subprocess.PIPE, shell=True, universal_newlines=True)


def run_aermod(jsonObj, select_aermod_project, aermet_project, pollutant_detail, select_sources_id, terrain_path):
    """
    :param jsonObj:
    :param select_aermod_project: 多个排放情景共享一套aermod参数设置
    :param aermet_project:
    :param pollutant_detail: 污染物对象
    :param select_sources_id:  该排放情景涉及到的多个污染源的id
    :return:
    """
    if pollutant_detail['is_gaseous']:
        print("is gaseous, ---------------------------------------",pollutant_detail['is_gaseous'])
        run_aermod_gas(jsonObj, select_aermod_project, aermet_project, pollutant_detail, select_sources_id, terrain_path)
    else:
        print("is particle, ---------------------------------------",pollutant_detail['is_gaseous'])
        run_aermod_particle(jsonObj, select_aermod_project, aermet_project, pollutant_detail, select_sources_id,
                            terrain_path)


def running_monitor(outdir, msg_str):
    filename = outdir + "/" + "running_monitor.txt"
    rnm = open(filename, "a", encoding='utf-8')
    rnm.write(msg_str)
    rnm.write("\n")
    rnm.close()


def touch_finished(outdir):
    filename = outdir + "/" + "_finished_"
    rnm = open(filename, "w")
    rnm.close()
    quit()


def run_aermod_batch(jsonObj, select_project, terrain_path, outdir):
    emissionCases = select_project['select_emissionCase']
    # select_aermet_project_id = select_project['select_AERMET_project_ID']
    # print(select_aermet_project_id)
    # print(type(select_aermet_project_id))
    # print(jsonObj['project_data']['AERMOD']['AERMET_project']['ID'])
    # print(type(jsonObj['project_data']['AERMOD']['AERMET_project']['ID']))
    # index = get_id(jsonObj['project_data']['AERMOD']['AERMET_project'], select_aermet_project_id, key='ID')
    # select_aermet_project = jsonObj['project_data']['AERMOD']['AERMET_project'][index]
    select_aermet_project = jsonObj['project_data']['AERMOD']['AERMET_project']

    for i in range(len(emissionCases)):
        emission_case_id = emissionCases[i]['emissionCase_ID']
        index = get_id(jsonObj['project_data']['basedata']['emissionCase'], emission_case_id, key='ID')
        emission_case_detail = jsonObj['project_data']['basedata']['emissionCase'][index]
        select_pollutant_id = emission_case_detail['pollutant_ID']
        print("select_pollutant_id", select_pollutant_id)
        index = get_id(jsonObj['project_data']['basedata']['pollutants'], select_pollutant_id, key='pollutants_ID')
        pollutant_detail = jsonObj['project_data']['basedata']['pollutants'][index]
        select_sources_id = emission_case_detail['select_Source_ID']

        aermod_path = '{0}'.format(emission_case_id)
        if os.path.exists(aermod_path):
            # shutil.rmtree(aermod_path, ignore_errors=True)
            # time.sleep(0.5)
            # os.makedirs(aermod_path)
            os.chdir(aermod_path)
            os.system('copy ..\\aermod.exe .')
            os.system('copy ..\\Bpipprm.exe .')
            # os.system('copy ..\\..\\AERMAP\\AERMAP.REC .')

            # os.system('copy ..\\..\\AERMET\\MET3.PFL .')
            # os.system('copy ..\\..\\AERMET\\MET3.SFC .')
        else:
            os.makedirs(aermod_path)
            os.chdir(aermod_path)
            os.system('copy ..\\aermod.exe .')
            os.system('copy ..\\Bpipprm.exe .')
            # os.system('copy ..\\..\\AERMAP\\AERMAP.REC .')
            os.system('copy ..\\..\\AERMET\\MET3.PFL .')
            os.system('copy ..\\..\\AERMET\\MET3.SFC .')

        try:
            if not os.path.exists("_aermod_ok_"):
                # TODO：运行建筑下洗前处理模块，产生build.BSO
                run_bpip(jsonObj, select_project, select_sources_id, terrain_path)

                # TODO: 运行一次aermod.exe,
                # TODO： 注意传递进去的是 select_sources_index
                run_aermod(jsonObj, select_project, select_aermet_project, pollutant_detail, select_sources_id, terrain_path)
                print('///////////////////////////////////////////////////////////////////////////////////////')
                print('方案{0}, 排放情景{1},运行完成'.format(select_project['ID'], emission_case_id))
                print('///////////////////////////////////////////////////////////////////////////////////////')
                msg_str = '    排放情景{0} aermod正在运行'.format(emission_case_id)
                running_monitor(outdir, msg_str)
            else:
                print("本次aermod已经运算成功")
                msg_str = '排放情景{0} aermod结果已经存在'.format(emission_case_id)
                running_monitor(outdir, msg_str)
            pass
        except Exception as ex:
            msg_str = '    排放情景{0} aermod运行失败'.format(emission_case_id)
            running_monitor(outdir, msg_str)
            template = "        An error of type {0} occurred. Arguments:\n{1!r}"
            msg_str = template.format(type(ex).__name__, ex.args)
            running_monitor(outdir, msg_str)
            logging.exception(ex)
            rnm = open("_aermod_fail_", "w")
            rnm.close()
        else:
            pass
            # rnm = open("_aermod_ok_", "w")
            # rnm.close()
        time.sleep(1)
        # os.remove('aermod.out')
        os.chdir('../')


def run_control(new_project_json_path, json_path, model_dir, work_dir, terrain_path):
    # TODO: 实际运行时根据项目的 MD5建立工作文件名，不会首先获取json信息
    f = open(new_project_json_path, "r", encoding="utf-8")
    jsonObj = json.loads(f.read())  # 打开json文件

    distance = jsonObj['project_data']['AERMOD']['AERMOD_receptor']['max_distance']
    os.chdir(work_dir)

    outdir, paraName = os.path.split(new_project_json_path)
    project_path, projectName = os.path.split(new_project_json_path)
    project_path = json_path
    # TODO: 默认一个json文件只含有一个项目
    # select_project = jsonObj2
    select_project = jsonObj['project_data']['AERMOD']['AERMOD_project'][0]
    first_day = select_project['first_day']
    calculateID = first_day + "/" + select_project['ID']
    cal_path = work_dir + "/" + calculateID

    # TODO: 生成运行文件，copy *exe和输入文件
    prepare_dirs(model_dir, work_dir, calculateID)
    os.chdir(cal_path)  # cal_path 为计算ID目录
    print("current directory is : " + os.getcwd())

    # TODO： 生成 grid_recetor.RE for aermap
    os.chdir('AERMAP')
    n_sub = 2
    # generate_receptor_grid.gen_grid_xy_auto(distance, n_sub=n_sub)
    generate_receptor_grid.gen_grid_xy_simple(distance, n_sub=n_sub)
    os.remove('grids.RE')
    os.chdir(cal_path)
    #
    # # TODO: 运行aermap
    # select_project['consider_terrain'] = 1
    print("select_project['consider_terrain']: ",select_project['consider_terrain'])
    if select_project['consider_terrain'] != 1:  # 1, 表示不考虑地形
        t1 = time.perf_counter()
        os.chdir('AERMAP')
        try:
            if not os.path.exists("_aermap_ok_"):
                # run_aermap(jsonObj, select_project)
                run_aermap_multi(outdir, jsonObj, select_project, n_sub=n_sub)
            else:
                print("aermap file already exits!")
                msg_str = "aermap files already exits!"
                running_monitor(outdir, msg_str)
            pass
        except Exception as ex:
            msg_str = "        aermap 运行失败"
            running_monitor(outdir, msg_str)
            template = "        An error of type {0} occurred. Arguments:\n{1!r}"
            msg_str = template.format(type(ex).__name__, ex.args)
            running_monitor(outdir, msg_str)
            logging.exception(ex)
        else:
            if os.path.exists("AERMAP.REC"):
                if os.path.getsize("AERMAP.REC") > 1:
                    rnm = open("_aermap_ok_", "w")
                    rnm.close()
                    msg_str = "        aermap 运行成功"
                    running_monitor(outdir, msg_str)

        os.chdir(cal_path)
        t2 = time.perf_counter()
        msg_str = "        aermap运行时间（s）：{0:.2f}".format(t2 - t1)
        print(msg_str)
        running_monitor(outdir, msg_str)
    else:
        print('本次计算不考虑地形影响')
        msg_str = "        本次计算不考虑地形影响"
        running_monitor(outdir, msg_str)

    # TODO：运行aermet
    t1 = time.perf_counter()
    os.chdir('AERMET')
    try:
        if not os.path.exists("_aermet_ok_"):
            run_aermet(jsonObj, select_project, project_path, first_day)
        else:
            print("aermet files already exits!")
            msg_str = "aermet files already exits!"
            running_monitor(outdir, msg_str)
        pass
    except Exception as ex:
        msg_str = "        aermet 运行失败"
        running_monitor(outdir, msg_str)
        template = "        An error of type {0} occurred. Arguments:\n{1!r}"
        msg_str = template.format(type(ex).__name__, ex.args)
        running_monitor(outdir, msg_str)
        logging.exception(ex)
    else:
        if os.path.exists("MET3.SFC"):
            if os.path.getsize("MET3.SFC") > 1:
                msg_str = "        aermet 运行成功"
                running_monitor(outdir, msg_str)
                rnm = open("_aermet_ok_", "w")
                rnm.close()

    os.chdir(cal_path)
    t2 = time.perf_counter()
    msg_str = "        aermet运行时间（s）：{0:.2f}".format(t2 - t1)
    print(msg_str)
    running_monitor(outdir, msg_str)

    # TODO: 对select_emissionCase循环，依次运行aermod
    # t1 = time.perf_counter()
    os.chdir('AERMOD')
    run_aermod_batch(jsonObj, select_project, terrain_path, outdir)
    os.chdir(cal_path)
    # t2 = time.perf_counter()
    # msg_str = "        aermod运行时间（s）：{0:.2f}".format(t2 - t1)
    # print(msg_str)
    # running_monitor(outdir, msg_str)


def main0():

    time_str = datetime.datetime.today().strftime('%Y%m%d')
    # time_str = '20201115'
    levels = [1000, 950, 925, 900, 850, 700, 600]
    download_levels = ['meteo', '1000', '950', '925', '900', '850', '700', '600', '500']  # 为了OSS下载遍历文件

    # 兰考
    station_id = 57093
    region = 'lankao'
    station_height = 73.3  # lankao
    lon = 114.835208  # lankao
    lat = 34.798364  # lankao

    # # 成华
    # station_id = 56294
    # region = 'chenghua'
    # station_height = 100  # chenghua
    # lon = 104.176  # chenghua
    # lat = 30.6872  # chenghua

    # 大同
    # station_id = 53487
    # region = 'datong'
    # station_height = 1067.2  # datong
    # lon = 113.25  # datong
    # lat = 40.05  # datong

    if station_id == 53487:
        # 如果是大同的这个检测站点
        levels = [850, 700, 600]

    surface_OQA_file_name = os.path.abspath(f'real_time/aermod_stastic/{region}/weatherdata/{time_str}/{station_id}_surface.OQA')
    upair_OQA_file_name = f'real_time/aermod_stastic/{region}/weatherdata/{time_str}/{station_id}_upair.OQA'
    csv_dir = os.path.abspath(f'real_time/data/{time_str}')

    days = 4  # 运行几天，与写OQA文件天数一致
    time_detal = 0  # 参数，固定为0（若数据使用北京时间记录，时间转换因子为0；若使用世界时记录，转换因子在中国地区设置为-8)
    base_path = os.path.dirname(os.path.abspath('real_time'))
    json_path = os.path.abspath(f"real_time/aermod_stastic/{region}")
    model_dir = os.path.abspath(f'real_time/aermod_stastic/{region}/model')
    work_dir = os.path.abspath(f'real_time/aermod_stastic/{region}/work')
    output_dir = os.path.abspath(f'real_time/aermod_stastic/{region}/output')
    data_dir = os.path.abspath(f'real_time/data')
    print(work_dir)

    terrain_path = model_dir + "/" + "terrain.tif"
    today = datetime.datetime.strptime(time_str, '%Y%m%d')

    # 从OSS上获取需要日期的数据  # TODO: 正式使用时使用，避免每次调试都下载一遍文件
    bucket_name = 'zhonglan-meteorology'  # 远程数据的容器的名字，涉及具体项目时需要联系管理员获取
    bucket = OSS(bucket_name)
    # 如果下载失败返回False
    # time_str_end = '2020121820'
    time_str_end = bucket.download_csv(time_str, data_dir, station_id, download_levels)
    if time_str_end:
        print('文件全部下载成功，正常运行')
    else:
        print('文件下载失败，请查看错误原因')
        return
    # 写探空和地面OQA
    OQA_file_dir = os.path.dirname(upair_OQA_file_name)
    if not os.path.exists(OQA_file_dir):
        os.makedirs(OQA_file_dir)
    write_upair_OQA.write_OQA_file(levels, station_id, time_str, csv_dir, upair_OQA_file_name, time_str_end, days, lon, lat)
    write_surface_OQA.write_OQA_file(station_id, time_str, surface_OQA_file_name, csv_dir, time_str_end, days, time_detal, station_height, lon, lat)
    for day_i in range(days):
        try:  # TODO: main process
            # TODO：aermod运算
            # 这里搬运并修改json
            offset = datetime.timedelta(days=day_i)
            this_day = today + offset
            this_time_str = this_day.strftime('%Y%m%d')
            base_json_path = os.path.join(json_path, f'{region}.json')
            # 运行过程中，有目录切换操作，每次运行，切换回初始路径
            os.chdir(base_path)
            new_project_json_dir = os.path.abspath(f'real_time/aermod_stastic/{region}/work/{time_str}')
            # 创建每日的项目目录
            if not os.path.exists(new_project_json_dir):
                os.makedirs(new_project_json_dir)

            # 执行一个模型运行之前，把上次项目删除
            project_path = os.path.join(new_project_json_dir, f'{region}_{this_time_str}')
            if os.path.exists(project_path):
                print(project_path)
                shutil.rmtree(project_path)
            # 配置json的路径, 后续运算的运行都依赖这个
            new_project_json_path = os.path.join(new_project_json_dir, f'{region}_{this_time_str}.json')
            # 生成配置json
            generate_json_config.generate_json_config(base_json_path, new_project_json_path, time_str, this_time_str, region)

            try:
                # 运行模型
                run_control(new_project_json_path, json_path, model_dir, work_dir, terrain_path)
                # 完成后，做后处理， 将plt文件转为png存到固定位置
                is_success = deal_result.deal_result(new_project_json_path, output_dir, lon, lat)
                if is_success:
                    pass
                else:
                    raise Exception('aermod运算失败')
            except Exception as ex:
                msg_str = "aermod运算失败"
                running_monitor(json_path, msg_str)
                template = "        An error of type {0} occurred. Arguments:\n{1!r}"
                msg_str = template.format(type(ex).__name__, ex.args)
                logging.exception(ex)
                running_monitor(json_path, msg_str)
                touch_finished(json_path)
            else:
                msg_str = "aermod运算成功"
                running_monitor(json_path, msg_str)
    #
            time.sleep(1)
        except Exception as ex:
            print(ex)
            touch_finished(json_path)


        # # TODO: extra任务，删除几天前的中间文件夹，防止占用过多存储
        # ndays_ago = 2
        # del_working_oldFiles.del_old_dir(work_dir, ndays_ago)

        # python -m nuitka --module run_aermod.py


if __name__ == '__main__':
    main0()