Comp_ADCP_MP.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
import warnings
from itertools import *
from shapely.geometry import mapping, LineString
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_ADCP_MP import plot_profile_3d_ADCP_MP
from scipy.interpolate import RegularGridInterpolator
import matplotlib.tri as mtri
import math

# import pykrige.kriging_tools as kt
# from pykrige.ok import OrdinaryKriging

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

x_mes = []
y_mes = []
cord_mes = open('Avl_US_DM270918sf12W1_0_010_20180927114429_0_44.303205N_4.741652E_gps_ASC.TXT').read().splitlines()
for x_l in cord_mes:
    y, x = x_l.split(',')
    if x == NODATA or y == NODATA:
        logger.info("Missing GPS values on, at least, one point from ADCP GPS file")
    else:
        x_mes.append(x)
        y_mes.append(y)
x_mes = [float(a) for a in x_mes]
y_mes = [float(a) for a in y_mes]
inProj = Proj("+init=EPSG:%i" % 4326)
outProj = Proj("+init=EPSG:%i" % 2154)
x_mes, y_mes = transform(inProj, outProj, x_mes, y_mes)
SCHEMA = {'geometry': 'LineString',
          'properties': {'nom': 'str'}}
with fiona.open('Trace_ADCP_GPS_TEST.shp', 'w', 'ESRI Shapefile', SCHEMA, crs=from_epsg(2154)) as out_shp:  # out_E.out_E
    Ltest = LineString([(x_2, y_2) for x_2, y_2 in zip(x_mes, y_mes)])
    elem = {}
    elem['geometry'] = mapping(Ltest)
    elem['properties'] = {
        'nom': 'ADCP line'}
    out_shp.write(elem)
p_raw = RawProfileObj('Avl_US_DM270918sf12W1_0_010_ASC.TXT')
processing_settings = {'proj_method': 2}
startingpoint = dict(start=Vector(0, 0))
p0 = ProcessedProfileObj(p_raw, processing_settings, startingpoint)

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

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


############
# WARNING point départ ADCP - MP
############

x_pt_depart_adcp = x_mes[0]
y_pt_depart_adcp = y_mes[0]

x_mes_physique = []
y_mes_physique = []
with open(csv_MP, '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_physique.append(X)
        y_mes_physique.append(Y)
inProj = Proj("+init=EPSG:%i" % 27563)
outProj = Proj("+init=EPSG:%i" % 2154)
x_mes_physique, y_mes_physique = transform(inProj, outProj, x_mes_physique, y_mes_physique)

distance_ADCP_MP = ((x_pt_depart_adcp-x_mes_physique[-1])**2+(y_pt_depart_adcp-y_mes_physique[-1])**2)**0.5
print(f"Attention : Le point de départ de la mesure modèle physique se situe à une distance du début de la mesure ADCP de {distance_ADCP_MP:.2f} m.")

############################

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'].iloc[-1])**2+(df_MP['Y_ADCP_RN'][i]-df_MP['Y_ADCP_RN'].iloc[-1])**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['Vu_RN'][i], df_MP['Vv_RN'][i], df_MP['Vw_RN'][i], df_MP['Vnorm_RN'][i], distance_MP[i]])

df_MP['distance_MP'] = distance_MP

triang_MP = mtri.Triangulation(df_MP['distance_MP'], df_MP['D_RN'])
interpolator_MP = mtri.LinearTriInterpolator(triang_MP, df_MP["Vnorm_RN"])

ADCP_res = []
for el in p0.ensembles:
    if el.void != True :
        for c in el.cells:
            ADCP_res.append(
            [el.position.x, el.position.y, c.z_position, c.velocity.x, c.velocity.y, c.velocity.z,
             (c.velocity.x ** 2 + c.velocity.y ** 2 + c.velocity.z ** 2)**0.5,
             (el.position.x ** 2 + el.position.y ** 2)**0.5])

grille_adcp = pd.DataFrame(ADCP_res,
                               columns=['X', 'Y', 'Z', 'U', 'V', 'W', 'Magnitude', 'Distance'])

Vitesse_MP_sur_grille_adcp = interpolator_MP(grille_adcp['Distance'],grille_adcp['Z'])

grille_adcp['Vitesse_MP_sur_grille_adcp'] = Vitesse_MP_sur_grille_adcp
grille_adcp_emprise_MP = grille_adcp.dropna() # drop les cellules ADCP qui sont hors de l'emprise du modèle physique


grille_adcp_emprise_MP["Diff_ADCP_MP"] = grille_adcp_emprise_MP["Magnitude"] - grille_adcp_emprise_MP["Vitesse_MP_sur_grille_adcp"]

plot_profile_3d_ADCP_MP(p0,cfg={'style': 'gradient'},df_MP=df_MP, df=grille_adcp_emprise_MP,shift=0)

nash = 1 - (sum((grille_adcp_emprise_MP['Magnitude'] - grille_adcp_emprise_MP['Vitesse_MP_sur_grille_adcp']) ** 2) / sum((grille_adcp_emprise_MP['Magnitude'] - grille_adcp_emprise_MP['Magnitude'].mean()) ** 2))
rmse = (((grille_adcp_emprise_MP['Magnitude'] - grille_adcp_emprise_MP['Vitesse_MP_sur_grille_adcp']) ** 2).mean()) ** 0.5
ecart_min = min(grille_adcp_emprise_MP["Diff_ADCP_MP"])
ecart_max = max(grille_adcp_emprise_MP["Diff_ADCP_MP"])
first_q = grille_adcp_emprise_MP["Diff_ADCP_MP"].quantile(0.25)
second_q = grille_adcp_emprise_MP["Diff_ADCP_MP"].quantile(0.5)
thrid_q = grille_adcp_emprise_MP["Diff_ADCP_MP"].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_ADCP_MP.csv',index=False, sep=';')