|
#! /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()
|
