#! /usr/bin/env python
# -*- coding: utf-8 -*-
#
#	Viability computation - post s2m
#	Author: P.Spandre
#	Date: 	2017-10-22

import netCDF4
import csv
import numpy as np

################################################################
##
##		DEFINITIONS FONCTIONS
##
################################################################

## 		FICHIER DE POIDS
# Lecture du fichier de discrétisation des locations avec poids

def dict_loc_poids(fichier_poids, indicatif, type_man):
    code_dict =  dict()
    disc_file = list(csv.reader(open(fichier_poids,'r'),delimiter=';'))

    code = []
    ind = []
    snow_type = []
    code_loc = []
    mp_part = []

    for k in range(1,len(disc_file)):
      ind = disc_file[k][0][0:4]
      snow_type = disc_file[k][0][4:6]
      mp_part = float(disc_file[k][1])

      code_loc = disc_file[k][0][6:]
      if (ind == indicatif) and (snow_type == type_man):
	code_dict[code_loc] = mp_part
    return code_dict

## 		PART MP en DAMAGE
def mp_part_gr(code_dict, lat, lon, zs, aspect, slope):
    # Dictionnaire de correspondance des expositions (degrés / code)
    dict_aspect = {-1:0,0:1,45:2,90:3,135:4,180:5,225:6,270:7,315:8}
    
    len_dim_loc=len(zs)
    mp_part_gr = np.zeros(len_dim_loc)

    for loc in range(len_dim_loc):#2 4444757 667076 1200
      # Code de description du fichier netcdf
      if alti[loc] > 1000.:
	area_code = str(2) + str(int(lat[loc]*100000)) + str(int(lon[loc]*100000)) + str(int(alti[loc]/100)) + str(int(dict_aspect[aspect[loc]])) + str(int(slope[loc]/10))
      else:
	area_code = str(2) + str(int(lat[loc]*100000)) + str(int(lon[loc]*100000)) + str(9) + str(int(alti[loc]/100)) + str(int(dict_aspect[aspect[loc]])) + str(int(slope[loc]/10))
      # Recherche dans le fichier de locations
      if code_dict.has_key(area_code):
	mp_part_gr[loc] = float(code_dict[area_code])
      # Si non défini, le mp_part reste à 0
    return mp_part_gr

    ## 		PART MP en DAMAGE + PRODUCTION
def mp_part_nc(code_dict, lat, lon, zs, aspect, slope):
    
    # Dictionnaire de correspondance des expositions (degrés / code)
    dict_aspect = {-1:0,0:1,45:2,90:3,135:4,180:5,225:6,270:7,315:8}
    
    len_dim_loc=len(zs)
    mp_part_nc = np.zeros(len_dim_loc)

    for loc in range(len_dim_loc):
      # Code de description du fichier netcdf
      if alti[loc] > 1000.:
	area_code = str(1) + str(int(lat[loc]*100000)) + str(int(lon[loc]*100000)) + str(int(alti[loc]/100)) + str(int(dict_aspect[aspect[loc]])) + str(int(slope[loc]/10))
      else:
	area_code = str(1) + str(int(lat[loc]*100000)) + str(int(lon[loc]*100000)) + str(9) + str(int(alti[loc]/100)) + str(int(dict_aspect[aspect[loc]])) + str(int(slope[loc]/10))
      # Recherche dans le fichier de locations
      if code_dict.has_key(area_code):
	mp_part_nc[loc] = code_dict[area_code]
      # Si non défini, le mp_part reste à 0
    return mp_part_nc

    ## 		VIABILITE PAR JOUR PAR INDICATIF
def viab(SWE, list_mp, zs, list_days):
    viab_days = np.zeros(len(list_days))
    for k in range(len(list_days)):
      viab_loc = np.zeros(len(zs))
      for loc in range(len(zs)):
	if SWE[k,loc] >= 100:
	  viab_loc[loc] = 1
      viab_days[k] = np.sum(np.multiply(list_mp,viab_loc))
    return viab_days

#####################################################
##	LECTURE DU FICHIER DE POIDS
fichier_loc = '20170616_Fichier_poids_man215.csv'

##	FICHIER DE SORTIE
out_file = open('Exemple.txt','w')
out_file.write('year;ind;viab_gr;viab_100nc;viab_45nc;viab_30nc;viab_15nc' + '\n')

for year in [2007]:
  print year

  pro_file='PRO_'+str(year)+'080106_'+str(year+1)+'080106'
  pro = netCDF4.Dataset('20160629_GRO/' + pro_file + '.nc','a')
  pro_nc = netCDF4.Dataset('20160629_SM2/' + pro_file + '.nc','a')
  lat = pro.variables["LAT"]
  lon = pro.variables["LON"]
  alti = pro.variables["ZS"]
  aspect = pro.variables["aspect"]
  slope = pro.variables["slope"]

  time = pro.variables['time']
  times = netCDF4.num2date(time,time.units)

  ## 	DETERMINATION INDEX NOEL, FEVRIER
  for i in range(len(time)):

    ##		VACANCES NOEL
    if times[i].year==year and times[i].month==12 and times[i].day==20:
      deb_noel = i
    if times[i].year==year+1 and times[i].month==1 and times[i].day==5:
      fin_noel = i

    ##		VACANCES FEVRIER
    if times[i].year==year+1 and times[i].month==2 and times[i].day==5:
      deb = i
    if times[i].year==year+1 and times[i].month==3 and times[i].day==5:
      fin = i
  days_christmas = np.arange(deb_noel,fin_noel+1)
  days_february = np.arange(deb,fin+1)
  #nb: +1 pour inclure le dernier jour

  #	Ensemble de jrs pr calcul viabilite quotidienne
  list_days = np.concatenate((days_christmas,days_february),axis=1)

  ##	SIMULATIONS CONDITIONS DAMAGE
  SD_GRO = pro.variables["DSN_T_ISBA"][list_days,0,:]
  SWE_GRO = pro.variables["WSN_T_ISBA"][list_days,0,:]

  ##	SIMULATIONS CONDITIONS DAMAGE + PRODUCTION
  SD_NC = pro_nc.variables["DSN_T_ISBA"][list_days,0,:]
  SWE_NC = pro_nc.variables["WSN_T_ISBA"][list_days,0,:]

  man = ['gr','n4','n3','n1']

  list_ind = ['3802']

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

  for indicatif in list_ind:
    print(indicatif)

    out_file.write(str(year) + ';' + str(indicatif))

    for type_man in man:
      code_dict = dict_loc_poids(fichier_loc, indicatif, type_man)

      #Part MP par location
      list_mp_gr = mp_part_gr(code_dict, lat, lon, alti, aspect, slope)
      list_mp_nc = mp_part_nc(code_dict, lat, lon, alti, aspect, slope)

      # Viabilite quotidienne
      viab_days = viab(SWE_GRO, list_mp_gr, alti, list_days) + viab(SWE_NC, list_mp_nc, alti, list_days)

      # Viabilite moyenne sur periodes noel/fevrier
      viab_chr = np.mean(viab_days[0:len(days_christmas)])
      viab_feb = np.mean(viab_days[len(days_christmas):len(viab_days)])
      coeff = 0.05
      viab_indicatif = viab_chr*coeff + viab_feb*(1-coeff)
      out_file.write(';' + str(round(viab_indicatif,2)))

      if type_man=='gr':	# Simulation du cas 100%
	viab_days = viab(SWE_NC, list_mp_gr, alti, list_days)

	viab_chr = np.mean(viab_days[0:len(days_christmas)])
	viab_feb = np.mean(viab_days[len(days_christmas):len(viab_days)])
	viab_indicatif = viab_chr*coeff + viab_feb*(1-coeff)

	out_file.write(';' + str(round(viab_indicatif,2)))

    out_file.write('\n') 

out_file.close()