Comp_MP_T3D.py

#!/usr/bin/env python
"""
ADCP / TELEMAC comparison in 3D
"""
#
import matplotlib.tri as mtri
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.interpolate import griddata
import fiona
import csv
from itertools import *
from shapely.geometry import mapping, LineString, Point
import sys
from adcploader import *
from pyteltools.geom import Shapefile
from pyteltools.slf.interpolation import MeshInterpolator
from pyteltools.slf import Serafin
from pyteltools.conf import settings
from pyteltools.slf.flux import PossibleFluxComputation, FluxCalculator
from pyteltools.geom.transformation import Transformation
from pyteltools.utils.cli_base import logger, PyTelToolsArgParse
from pyproj import Proj, transform
from fiona.crs import from_epsg
from tqdm import tqdm
from quickviz_MP_T3D import plot_profile_3d_MP_T3D
from scipy.interpolate import LinearNDInterpolator
# import geopandas as gpd


NODATA = '-32768'
var_ID_x = 'U'
var_ID_y = 'V'
var_ID_z = 'W'
var_ID_2D = ["U" , "V" , "H"]


csv_file_path = 'ADCPlabo_CAL_2018.09.27_3_PT1001.csv'
output_shapefile = 'TRACE_MP_GPS.shp'
x_mes = []
y_mes = []


with open(csv_file_path, 'r') as csv_file:
    csv_reader = csv.DictReader(csv_file, delimiter=';')
    for row in csv_reader:
        X = float(row['X_ADCP_RN'])
        Y = float(row['Y_ADCP_RN'])
        x_mes.append(X)
        y_mes.append(Y)
inProj = Proj("+init=EPSG:%i" % 27563)
outProj = Proj("+init=EPSG:%i" % 2154)
x_mes, y_mes = transform(inProj, outProj, x_mes, y_mes)
with fiona.open(output_shapefile,'w',driver='ESRI Shapefile',crs='EPSG:2154', schema={'geometry': 'LineString','properties': {'name': 'str'}}) as out_shp:
    Ltest = LineString([(x_2, y_2) for x_2, y_2 in zip(x_mes, y_mes)])
    elem = {}
    elem['geometry'] = mapping(Ltest)
    elem['properties'] = {
        'name': 'ADCP line'}
    out_shp.write(elem)

lines = []
for poly in Shapefile.get_lines('Trace_MP_GPS.shp', shape_type=3):
    lines.append(poly)

with Serafin.Read('r2d_last_C42bis.slf', 'fr') as in_slf_q:
    in_slf_q.read_header()
    if not in_slf_q.header.is_2d:
        logger.error("File is not TELEMAC-2D")
    in_slf_q.get_time()
    in_slf_q.header.transform_mesh([Transformation(0, 1, 1, 838000, 6355500, 0)])
    section_names = ['Section %i' % (i + 1) for i in range(len(lines))]
    for var_ID_2d in var_ID_2D:
        if var_ID_2d not in in_slf_q.header.var_IDs:
            logger.error('Upper variable "%s" is not in Serafin file' % var_ID_2d)
    calculator_Q = FluxCalculator(FluxCalculator.DOUBLE_LINE_INTEGRAL, ['M', 'H'], in_slf_q,
                                  section_names, lines, 1)
    calculator_S = FluxCalculator(FluxCalculator.LINE_INTEGRAL, ['H'], in_slf_q,
                                  section_names, lines, 1)
    calculator_S.construct_triangles()
    calculator_S.construct_intersections()
    result_S = calculator_S.run(fmt_float='{0:.6f}')
    calculator_Q.construct_triangles(tqdm)
    intersections_Q = calculator_Q.construct_intersections()
    intersections_S = calculator_S.construct_intersections()  # DO THE SAME THING ABOUT Q
    result = []
    values_Q = []
    values_S = []
    values_V = []
    for time_index, time in enumerate(tqdm(in_slf_q.time, unit='frame')):
        i_result = [str(time)]
        valuesQ = []
        valuesS = []

        for vara_ID in calculator_Q.var_IDs:
            valuesQ.append(in_slf_q.read_var_in_frame(time_index, vara_ID))

        for varb_ID in calculator_S.var_IDs:
            valuesS.append(in_slf_q.read_var_in_frame(time_index, varb_ID))

        for j in range(len(lines)):
            intersections_QQ = calculator_Q.intersections[j]
            intersections_SS = calculator_S.intersections[j]
            flux = calculator_Q.flux_in_frame(intersections_QQ, valuesQ)
            flux_S = calculator_S.flux_in_frame(intersections_SS, valuesS)

            i_result.append(settings.FMT_FLOAT.format(flux, flux_S))

        vmoy_telemac = flux / flux_S
print('débit telemac =', flux)


print('vitessee moyenne t3d-PNU:',vmoy_telemac)

with Serafin.Read('r3d_last_C42bis.slf', 'fr') as in_slf:
    in_slf.read_header()
    if in_slf.header.is_2d:
        logger.error("File is not TELEMAC-3D")

    in_slf.get_time()
    index_time = len(in_slf.time) - 1
    planes = in_slf.header.nb_planes
    nplan = planes
    nb_2Dnodes = in_slf.header.nb_nodes_2d
    output_header = in_slf.header.copy()
    output_header.transform_mesh([Transformation(0, 1, 1, 838000, 6355500, 0)])
    mesh = MeshInterpolator(output_header, True)
    if "Z" not in in_slf.header.var_IDs:
        logger.error('Upper variable Z is not in Serafin file')
    z = in_slf.read_var_in_frame(index_time, 'Z').reshape((nplan, in_slf.header.nb_nodes_2d)).T
    z_surface_libre = in_slf.read_var_in_frame(index_time, 'Z')
    print('z surface libre =', z_surface_libre.mean())
    values_x = in_slf.read_var_in_frame(index_time, var_ID_x).reshape((nplan, in_slf.header.nb_nodes_2d)).T
    values_y = in_slf.read_var_in_frame(index_time, var_ID_y).reshape((nplan, in_slf.header.nb_nodes_2d)).T
    values_z = in_slf.read_var_in_frame(index_time, var_ID_z).reshape((nplan, in_slf.header.nb_nodes_2d)).T
    values = (values_x**2 + values_y**2 + values_z**2)**0.5
    nb_nonempty, indices_nonempty, lines_interpolators, line_interpolators_internal = mesh.get_line_interpolators(
        lines)
    res = mesh.interpolate_along_lines(in_slf, 'Z', list(range(len(in_slf.time))), indices_nonempty,
                                       lines_interpolators, '{:.6e}')

    Z = in_slf.read_var_in_frame_as_3d(index_time, 'Z')
    h=0
    for i in range(0,nplan-1):
        for j in range(0,nb_2Dnodes):
            if abs(Z[i,j]-Z[i+1,j])>h:
                h = abs(Z[i,j]-Z[i+1,j])

line_interpolator, distances = line_interpolators_internal[0]
npt = len(distances)
point_y = np.empty((npt, nplan))
point_values = np.empty((npt, nplan))
point_values_x = np.empty((npt,nplan))
point_values_y = np.empty((npt,nplan))
point_values_z = np.empty((npt,nplan))
for i_pt, ((x, y, (i, j, k), interpolator), distance) in enumerate(zip(line_interpolator, distances)):
    point_y[i_pt] = interpolator.dot(z[[i, j, k]])
    point_values_x[i_pt] = interpolator.dot(values_x[[i, j, k]])
    point_values_y[i_pt] = interpolator.dot(values_y[[i, j, k]])
    point_values_z[i_pt] = interpolator.dot(values_z[[i, j, k]])
point_values = (point_values_x**2+point_values_y**2+point_values_z**2)**0.5
point_x = np.array([[distances[i]] * nplan for i in range(npt)])
triangles = [((ipt - 1) * nplan + iplan - 1, (ipt - 1) * nplan + iplan, ipt * nplan + iplan - 1)
             for ipt in range(1, npt) for iplan in range(1, nplan)] + \
            [((ipt - 1) * nplan + iplan, ipt * nplan + iplan - 1, ipt * nplan + iplan)
             for ipt in range(1, npt) for iplan in range(1, nplan)]
triang = mtri.Triangulation(point_x.flatten(), point_y.flatten(), triangles)

z_t3d = max(point_y.flatten())

csv_MP = 'ADCPlabo_CAL_2018.09.27_3_PT1001.csv'
df_MP = pd.read_csv(csv_MP, sep=';')

MP_res = []
distance_MP = [0]*len(df_MP)
z_mp = max(df_MP['Z_RN'])
vmoy_mp = sum(df_MP['Vnorm_RN'])/len(df_MP['Vnorm_RN'])
for i in range(len(df_MP)):
    if i == 0:
        distance_MP[i] = 0
    else:
        distance_MP[i] = ((df_MP['X_ADCP_RN'][i]-df_MP['X_ADCP_RN'][0])**2+(df_MP['Y_ADCP_RN'][i]-df_MP['Y_ADCP_RN'][0])**2)**0.5
    MP_res.append([df_MP['X_ADCP_RN'][i], df_MP['Y_ADCP_RN'][i], df_MP['Z_RN'][i], df_MP['Vnorm_RN'][i], distance_MP[i]])

TELEMAC_res = []
for dist, z, values, values_x, values_y, values_z in zip(list(chain(*point_x.tolist())), list(chain(*point_y.tolist())),
                           list(chain(*point_values.tolist())), list(chain(*point_values_x.tolist())), list(chain(*point_values_y.tolist())), list(chain(*point_values_z.tolist()))):
    TELEMAC_res.append([dist, z, values])
L = len(TELEMAC_res)
taille_max = max(len(MP_res), len(TELEMAC_res))
while len(MP_res) < taille_max:
    MP_res.append([0] * len(MP_res[0]))
while len(TELEMAC_res) < taille_max:
    TELEMAC_res.append([0] * len(TELEMAC_res[0]))
grille_commune = []
for i in range(taille_max):
    grille_commune.append(MP_res[i] + TELEMAC_res[i])



grille_mp_t3d = pd.DataFrame(grille_commune,
                               columns=['X_RN', 'Y_RN', 'Z_RN', 'Vnorm_RN', 'distance_RN', 'Dt3d',
                                        'Zt3d', 'Mt3d'])


####################
# Methode 2: TIN avec nos triangles --> Méthode choisie pour le moment
####################

triang = mtri.Triangulation(point_x.flatten(), point_y.flatten(), triangles)
interpolator_T3D = mtri.LinearTriInterpolator(triang, point_values.flatten())

Vitesse_T3D_sur_grille_mp = interpolator_T3D(grille_mp_t3d['distance_RN'],grille_mp_t3d['Z_RN'])

grille_mp_t3d['Vitesse_T3D_sur_grille_mp'] = Vitesse_T3D_sur_grille_mp
grille_mp_emprise_T3D = grille_mp_t3d.dropna()

grille_mp_emprise_T3D["Diff_MP_T3D"] = grille_mp_emprise_T3D['Vnorm_RN'] - grille_mp_emprise_T3D["Vitesse_T3D_sur_grille_mp"]

plot_profile_3d_MP_T3D(r3d=(point_x, point_y, point_values, triang), df=grille_mp_emprise_T3D, shift=0, vmoy_telemac=vmoy_telemac)

nash = 1 - (sum((grille_mp_emprise_T3D['Vnorm_RN'] - grille_mp_emprise_T3D['Vitesse_T3D_sur_grille_mp']) ** 2) / sum((grille_mp_emprise_T3D['Vnorm_RN'] - grille_mp_emprise_T3D['Vnorm_RN'].mean()) ** 2))
rmse = (((grille_mp_emprise_T3D['Vnorm_RN'] - grille_mp_emprise_T3D['Vitesse_T3D_sur_grille_mp']) ** 2).mean()) ** 0.5
ecart_min = min(grille_mp_emprise_T3D["Diff_MP_T3D"])
ecart_max = max(grille_mp_emprise_T3D["Diff_MP_T3D"])
first_q = grille_mp_emprise_T3D["Diff_MP_T3D"].quantile(0.25)
second_q = grille_mp_emprise_T3D["Diff_MP_T3D"].quantile(0.5)
thrid_q = grille_mp_emprise_T3D["Diff_MP_T3D"].quantile(0.75)
stat_indic = []
stat_indic.append([nash, rmse, ecart_min, ecart_max, first_q, second_q, thrid_q])
df_stat_indic = pd.DataFrame(stat_indic,
                             columns=["Nash", "RMSE [m/s]", "Ecart_min [m/s]", "Ecart_max [m/s]", "1_quantile ecart [m/s]", "Medianne ecart[m/s]",
                                      "3_quantile ecart [m/s]"])
df_stat_indic.to_csv('Indicateurs_statistiques_MP_T3D.csv',index=False, sep=';')