Project

General

Profile

20170922_Traitement_Viab.py

carlo carmagnola, 11/14/2018 04:12 PM

 
1
#! /usr/bin/env python
2
# -*- coding: utf-8 -*-
3
#
4
#        Viability computation - post s2m
5
#        Author: P.Spandre
6
#        Date:         2017-10-22
7

    
8
import netCDF4
9
import csv
10
import numpy as np
11

    
12
################################################################
13
##
14
##                DEFINITIONS FONCTIONS
15
##
16
################################################################
17

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

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

    
25
    code = []
26
    ind = []
27
    snow_type = []
28
    code_loc = []
29
    mp_part = []
30

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

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

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

    
49
    for loc in range(len_dim_loc):#2 4444757 667076 1200
50
      # Code de description du fichier netcdf
51
      if alti[loc] > 1000.:
52
        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))
53
      else:
54
        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))
55
      # Recherche dans le fichier de locations
56
      if code_dict.has_key(area_code):
57
        mp_part_gr[loc] = float(code_dict[area_code])
58
      # Si non défini, le mp_part reste à 0
59
    return mp_part_gr
60

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

    
70
    for loc in range(len_dim_loc):
71
      # Code de description du fichier netcdf
72
      if alti[loc] > 1000.:
73
        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))
74
      else:
75
        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))
76
      # Recherche dans le fichier de locations
77
      if code_dict.has_key(area_code):
78
        mp_part_nc[loc] = code_dict[area_code]
79
      # Si non défini, le mp_part reste à 0
80
    return mp_part_nc
81

    
82
    ##                 VIABILITE PAR JOUR PAR INDICATIF
83
def viab(SWE, list_mp, zs, list_days):
84
    viab_days = np.zeros(len(list_days))
85
    for k in range(len(list_days)):
86
      viab_loc = np.zeros(len(zs))
87
      for loc in range(len(zs)):
88
        if SWE[k,loc] >= 100:
89
          viab_loc[loc] = 1
90
      viab_days[k] = np.sum(np.multiply(list_mp,viab_loc))
91
    return viab_days
92

    
93
#####################################################
94
##        LECTURE DU FICHIER DE POIDS
95
fichier_loc = '20170616_Fichier_poids_man215.csv'
96

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

    
101
for year in [2007]:
102
  print year
103

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

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

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

    
119
    ##                VACANCES NOEL
120
    if times[i].year==year and times[i].month==12 and times[i].day==20:
121
      deb_noel = i
122
    if times[i].year==year+1 and times[i].month==1 and times[i].day==5:
123
      fin_noel = i
124

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

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

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

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

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

    
147
  list_ind = ['3802']
148

    
149
  ###################################################
150

    
151
  for indicatif in list_ind:
152
    print(indicatif)
153

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

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

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

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

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

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

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

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

    
182
    out_file.write('\n') 
183

    
184
out_file.close()