from pyproj import Proj
import os
from math import ceil
import json
import coorTransform

# 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 wgs84toutm_reverse(zone, utm_x, utm_y):
#     '''
#     :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")
#     lon, lat = utm_pj(utm_x, utm_y, inverse=True)  # 从经纬度转平面坐标
#     return lon, lat
#
# def tran(v):
#     if v < 0.0001:
#         return 5
#     elif v < 0.001:
#         return 4
#     elif v < 0.005:
#         return 3
#     elif v < 0.01:
#         return 2
#     elif v < 0.1:
#         return 1
#     else:
#         return 0
#
# json_path = r'D:\PycharmProjects\aermod\real_time\aermod_stastic\lankao\output\20201105\lankao_20201105\AERMOD\102_result.json'
#
# with open(json_path, 'r', encoding='utf-8') as f:
#     j_dict = json.load(f)
# print(j_dict[8][0])
# one_hour = j_dict[8]
#
# arr = []
# print(len(one_hour))
# i = 0
#
#
# for point in one_hour:
#     # p['lng'] = point['lng']
#     # p['lat'] = point['lat']
#     if point['count'] < 0.0001:
#         continue
#     zone, utm_x, utm_y = wgs84toutm(point['lng'], point['lat'])
#     utm_1, utm_2 = wgs84toutm_reverse(zone, utm_x-200, utm_y+200)
#     lng, lat = coorTransform.wgs84_to_bd09(point['lng'], point['lat'])
#     lng1, lat1 = coorTransform.wgs84_to_bd09(utm_1, utm_2)
#     p = [point['lng'], point['lat'], tran(point['count']), lng1, lat1]
#     arr.append(p)
#
# with open('aaa.json', 'w', encoding='utf-8') as f:
#     f.write(json.dumps(arr, ensure_ascii=False))

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib  # Windows平台
from pyproj import Proj
import json

matplotlib.rcParams['font.family'] = ['SimHei']  # 设置中文字体
matplotlib.rcParams['axes.unicode_minus'] = False  # 坐标轴的负号正常显示


def read_plt_grid(file_path, json_path):
    data = pd.read_csv(file_path, skiprows=8, header=None,
                       usecols=[0, 1, 2], delimiter='\s+', names=['X', 'Y', 'V'])
    print(data.head())
    data = data.values
    X = data[:, 0]
    Y = data[:, 1]
    C = data[:, 2]
    nnn = X.size
    print(nnn)
    XX = np.reshape(X, (int(nnn / 24), 24), order='F')
    YY = np.reshape(Y, (int(nnn / 24), 24), order='F')
    CC = np.reshape(C, (int(nnn / 24), 24), order='F')

    print("max:", np.max(CC))

    nnn2 = XX[:, 0].size
    ncell = int(np.sqrt(nnn2))
    print(ncell)
    X0 = np.reshape(XX[:, 10], (ncell, ncell), order='F')
    Y0 = np.reshape(YY[:, 10], (ncell, ncell), order='F')

    for i in range(1):
        i = 14
        print("hour:-------------- --------------    ", i)
        C0 = np.reshape(CC[:, i], (ncell, ncell), order='F')
        sss = 'hour ' + str(i + 1) + ' '

        f = open(json_path, "r", encoding="utf-8")
        jsonObj = json.loads(f.read())  # 打开json文件

        des_lon = 114.823391
        des_lat = 34.830386
        utm_x1, utm_y1 = relative_distance_coord(jsonObj, des_lon, des_lat)
        row_index, col_index = get_index_oushi(utm_x1, utm_y1, X0, Y0)
        print("兰考环保局站点值 ：", C0[row_index, col_index])

        des_lon = 114.835053
        des_lat = 34.8436
        utm_x2, utm_y2 = relative_distance_coord(jsonObj, des_lon, des_lat)
        row_index, col_index = get_index_oushi(utm_x2, utm_y2, X0, Y0)
        print("文化交流中心站点值 ：", C0[row_index, col_index])

        levels = np.array([0.001, 0.01, 0.05, 0.1, 1, 10., 100.])
        grid_countour_plot(X0, Y0, C0, sss, levels, [utm_x1, utm_y1], [utm_x2, utm_y2])


def grid_countour_plot(X0, Y0, C0, sss, levels, marker_point1, marker_point2):
    fig, ax1 = plt.subplots(1)
    import matplotlib.colors as colors
    import matplotlib as mpl
    cmap = mpl.colors.ListedColormap(
        ['greenyellow', 'yellow', 'orange', 'red', 'mediumvioletred', 'darkmagenta', 'maroon'])
    norm = colors.BoundaryNorm(boundaries=levels, ncolors=len(levels))
    contours1 = ax1.contourf(X0, Y0, C0, levels=levels, norm=norm, cmap=cmap)
    plt.plot(marker_point1[0], marker_point1[1], '*', color='k')
    plt.plot(marker_point2[0], marker_point2[1], 'o', color='k')
    # 去掉边框
    ax1.spines['top'].set_visible(False)
    ax1.spines['right'].set_visible(False)
    ax1.spines['bottom'].set_visible(False)
    ax1.spines['left'].set_visible(False)
    ax1.spines['top'].set_color('none')
    ax1.axis('off')

    plt.savefig('2222.png', dpi=300, bbox_inches='tight', pad_inches=0,transparent=True)
    plt.show()



def get_index_oushi(from_lon, from_lat, lon_model, lat_model):
    # 通过欧氏距离获得索引
    a = np.sqrt((from_lon - lon_model) ** 2 + (from_lat - lat_model) ** 2)
    index = np.unravel_index(a.argmin(), a.shape)
    row_index = index[0]
    col_index = index[1]
    return row_index, col_index


def wgs84toutm(lon, lat):
    '''
    :param lon:
    :param lat:
    :return:UTM投影带号和相应UTM坐标,及相对距离坐标
    '''
    zone = np.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 utm_x, utm_y


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

    return relative_x, relative_y


if __name__ == '__main__':
    # file_path = r'D:\code\aermod\real_time\aermod_stastic\chenghua\work\20201022\chenghua_20201022\AERMOD\101\all_1.plt'
    file_path = r'real_time\aermod_stastic\lankao\work\20201105\lankao_20201105\AERMOD\101\all_1.plt'

    json_path = r'real_time\aermod_stastic\lankao\work\20201022\lankao_20201022.json'
    read_plt_grid(file_path, json_path)