| 1 | import numpy as N |
|---|
| 2 | from functions import solarang,time_zone,residuals |
|---|
| 3 | from constantes import * |
|---|
| 4 | from genutil import statistics |
|---|
| 5 | from scipy.optimize import leastsq |
|---|
| 6 | #------------------------------------------------------------------------------------------- |
|---|
| 7 | # function gap_fill_func gapfills the meteorlogical data |
|---|
| 8 | # argument1 : the weather dataset (weather) |
|---|
| 9 | # argument2 : the climatology dataset (clim) |
|---|
| 10 | # argument3 : the ratio between clim and weather time steps (diff_clim_weather) |
|---|
| 11 | |
|---|
| 12 | # returns a vector that contains the weather dataset gapfilled |
|---|
| 13 | def gap_fill_func(weather,clim,diff_clim_weather,weather_period,climato_period,julian,year_length,lon,lat,gapmax,avg,climatoshift,timeshift): |
|---|
| 14 | weather_clim_period=[] |
|---|
| 15 | weather_clim_period_test=[] |
|---|
| 16 | |
|---|
| 17 | weather_clim_period_nogap=[] |
|---|
| 18 | clim_nogap=[] |
|---|
| 19 | weather_clim_period_gapfill=[] |
|---|
| 20 | weather_gapfill=[] |
|---|
| 21 | |
|---|
| 22 | if(timeshift==-9999): |
|---|
| 23 | timezone=time_zone(0,lon) |
|---|
| 24 | # east of Greenwich => timeshift>0 |
|---|
| 25 | # west of Greenwich => timeshift<0 |
|---|
| 26 | if(timezone[1]<13): |
|---|
| 27 | timeshift=timezone[1] |
|---|
| 28 | else: |
|---|
| 29 | timeshift=timezone[1]-24 |
|---|
| 30 | |
|---|
| 31 | stat_vec=[] |
|---|
| 32 | |
|---|
| 33 | print 'shift to UTC time =',timeshift,' hours' |
|---|
| 34 | |
|---|
| 35 | |
|---|
| 36 | for k in range(len(weather)): |
|---|
| 37 | |
|---|
| 38 | |
|---|
| 39 | weather_clim_period.append(N.zeros(len(weather[k])/diff_clim_weather, N.float32, 0)) |
|---|
| 40 | if(k==id_precip): |
|---|
| 41 | freq_precip=N.zeros(len(weather[k])/diff_clim_weather, N.float32, 0) |
|---|
| 42 | weather_clim_period_test.append(N.zeros(len(weather[k])/diff_clim_weather, N.float32, 0)) |
|---|
| 43 | |
|---|
| 44 | |
|---|
| 45 | # climatoshift indicates to which time step or time period |
|---|
| 46 | # unit is fraction of a time period (between 2 consecutive time steps) |
|---|
| 47 | # a climatic field value corresponds to |
|---|
| 48 | # for field that is a mean value (avg=1) |
|---|
| 49 | # climatoshift=0 when the mean value is calculated from one time step to the next one |
|---|
| 50 | # climatoshift=-0.5 when the mean value is centered on one time step |
|---|
| 51 | # for field that is an instantaneous value (avg=0) |
|---|
| 52 | # climatoshift=0 when the value correspond to the current time step |
|---|
| 53 | # climatoshift=1 when the value corresponds to the next time step |
|---|
| 54 | totalshift=timeshift+climatoshift[k]*climato_period |
|---|
| 55 | print '\tVar=',label_fig[k] |
|---|
| 56 | for t in range(len(weather[k])): |
|---|
| 57 | tshift=t-(float)(totalshift)/weather_period |
|---|
| 58 | if((tshift>=0)and(tshift<len(weather[k]))): |
|---|
| 59 | # in case of a mean value calculation, we sum all weather element within each element of weather_clim_period |
|---|
| 60 | # if one weather element equals -9999, the related weather_clim_period is equal to -9999 |
|---|
| 61 | if (avg[k]==1): |
|---|
| 62 | if((weather[k][t] != -9999) and (weather_clim_period[k][tshift/diff_clim_weather] != -9999)): |
|---|
| 63 | weather_clim_period[k][tshift/diff_clim_weather]+=weather[k][t]/diff_clim_weather |
|---|
| 64 | else: |
|---|
| 65 | weather_clim_period[k][tshift/diff_clim_weather]=-9999 |
|---|
| 66 | weather_clim_period_test[k][tshift/diff_clim_weather]=1 |
|---|
| 67 | |
|---|
| 68 | if(k==id_precip): |
|---|
| 69 | if((weather[k][t] != -9999) and (freq_precip[tshift/diff_clim_weather] != -9999)): |
|---|
| 70 | if(weather[k][t]>0): |
|---|
| 71 | freq_precip[tshift/diff_clim_weather]+=1./diff_clim_weather |
|---|
| 72 | else: |
|---|
| 73 | freq_precip[tshift/diff_clim_weather]=-9999 |
|---|
| 74 | |
|---|
| 75 | # in case of a instantaneaous calculation, each weather_clim_period element corresponds |
|---|
| 76 | # to the first weather element associated to this weather_clim_period element |
|---|
| 77 | else: |
|---|
| 78 | if(weather_clim_period_test[k][tshift/diff_clim_weather] == 0): |
|---|
| 79 | if((weather[k][t] != -9999)): |
|---|
| 80 | weather_clim_period[k][tshift/diff_clim_weather]=weather[k][t] |
|---|
| 81 | else: |
|---|
| 82 | weather_clim_period[k][tshift/diff_clim_weather]=-9999 |
|---|
| 83 | weather_clim_period_test[k][tshift/diff_clim_weather]=1 |
|---|
| 84 | # in case of a mean value calculcation (avg==1) |
|---|
| 85 | # elements of weather_clim_period that have been partly filled (at the beginning or at the end) |
|---|
| 86 | # are set to -9999 |
|---|
| 87 | if (avg[k]==1): |
|---|
| 88 | if(totalshift%climato_period !=0): |
|---|
| 89 | if(totalshift>0): |
|---|
| 90 | weather_clim_period[k][(len(weather[k])-1-(float)(totalshift)/weather_period)/diff_clim_weather]=-9999 |
|---|
| 91 | else: |
|---|
| 92 | if(totalshift<0): |
|---|
| 93 | weather_clim_period[k][(-(float)(totalshift)/weather_period)/diff_clim_weather]=-9999 |
|---|
| 94 | # elements of weather_clim_period that have NOT been filled (at the beginning or at the end) |
|---|
| 95 | # are set to -9999 |
|---|
| 96 | weather_clim_period[k]=N.where(weather_clim_period_test[k]==1,weather_clim_period[k],-9999) |
|---|
| 97 | |
|---|
| 98 | |
|---|
| 99 | weather_clim_period_nogap.append([]) |
|---|
| 100 | clim_nogap.append([]) |
|---|
| 101 | for t in range(len(weather_clim_period[k])): |
|---|
| 102 | if(weather_clim_period[k][t] != -9999): |
|---|
| 103 | weather_clim_period_nogap[k].append(weather_clim_period[k][t]) |
|---|
| 104 | clim_nogap[k].append(clim[k][t]) |
|---|
| 105 | |
|---|
| 106 | |
|---|
| 107 | weather_clim_period_nogap[k]=N.array(weather_clim_period_nogap[k],N.float) |
|---|
| 108 | clim_nogap[k]=N.array(clim_nogap[k],N.float) |
|---|
| 109 | |
|---|
| 110 | |
|---|
| 111 | |
|---|
| 112 | # Evaluate the correlation |
|---|
| 113 | # between clim_nogap and weather_clim_period_nogap |
|---|
| 114 | # intercept and slope of the relation |
|---|
| 115 | # will be used for correcting clim data when filling gaps |
|---|
| 116 | a=clim_nogap[k] |
|---|
| 117 | b=weather_clim_period_nogap[k] |
|---|
| 118 | |
|---|
| 119 | if(k==id_swdown): |
|---|
| 120 | stat=leastsq(residuals,1.,args=(b,a)) |
|---|
| 121 | stat=N.array(stat) |
|---|
| 122 | stat[1]=0 |
|---|
| 123 | elif(k==id_precip): |
|---|
| 124 | stat=N.zeros(2) |
|---|
| 125 | stat[0]=sum(b)/sum(a) |
|---|
| 126 | stat[1]=0 |
|---|
| 127 | else: |
|---|
| 128 | stat=statistics.linearregression(b,x=a) |
|---|
| 129 | |
|---|
| 130 | if(str(stat[0])!='nan'): |
|---|
| 131 | slope=stat[0] |
|---|
| 132 | intercept=stat[1] |
|---|
| 133 | else: |
|---|
| 134 | slope=1 |
|---|
| 135 | intercept=0 |
|---|
| 136 | |
|---|
| 137 | stat_vec.append([slope,intercept]) |
|---|
| 138 | print '\t\tslope of the linear relation in-situ VS reanalysis=',slope |
|---|
| 139 | print '\t\tintercept of the linear relation in-situ VS reanalysis=',intercept |
|---|
| 140 | if(str(stat[0])!='nan'): |
|---|
| 141 | print '\t\tRMSE without correction=',statistics.rms(a,b) |
|---|
| 142 | print '\t\tRMSE with correction=',statistics.rms(a*slope+intercept,b) |
|---|
| 143 | |
|---|
| 144 | |
|---|
| 145 | weather_clim_period_gapfill.append([]) |
|---|
| 146 | weather_clim_period_gapfill[k]=N.where(weather_clim_period[k]==-9999, slope*clim[k]+intercept, weather_clim_period[k]) |
|---|
| 147 | |
|---|
| 148 | weather_gapfill.append([]) |
|---|
| 149 | mean_csang=0. |
|---|
| 150 | if(k==id_precip): |
|---|
| 151 | freq_precip_nogap=N.ma.masked_values(freq_precip,-9999) |
|---|
| 152 | freq_precip_nogap_nonull=N.ma.masked_values(freq_precip_nogap,0.) |
|---|
| 153 | freq_precip_scalar=freq_precip_nogap_nonull.mean() |
|---|
| 154 | print '\t\tfreq_precip_scalar=',freq_precip_scalar |
|---|
| 155 | number_precip_per_diff_clim_weather=round(1./freq_precip_scalar) |
|---|
| 156 | |
|---|
| 157 | for t in range(len(weather[k])): |
|---|
| 158 | tshift=t-(float)(totalshift)/weather_period |
|---|
| 159 | |
|---|
| 160 | if(k==id_swdown): |
|---|
| 161 | if((t-(float)(timeshift)/weather_period)<(len(weather[k])-(diff_clim_weather-1))): |
|---|
| 162 | if(tshift%diff_clim_weather==0): |
|---|
| 163 | mean_csang=0. |
|---|
| 164 | for l in range(diff_clim_weather): |
|---|
| 165 | mean_csang+=solarang(julian[(t-(float)(timeshift)/weather_period)+l],0,lon,lat,year_length[tshift])/diff_clim_weather |
|---|
| 166 | if(mean_csang == 0.): |
|---|
| 167 | ratio=0. |
|---|
| 168 | else: |
|---|
| 169 | ratio=solarang(julian[(t-(float)(timeshift)/weather_period)],0,lon,lat,year_length[(t-(float)(timeshift)/weather_period)])/mean_csang |
|---|
| 170 | else: |
|---|
| 171 | ratio=1 |
|---|
| 172 | if(k==id_lwdown): |
|---|
| 173 | ratio=1 |
|---|
| 174 | if(k==id_precip): |
|---|
| 175 | if(t%number_precip_per_diff_clim_weather==0): |
|---|
| 176 | ratio=number_precip_per_diff_clim_weather |
|---|
| 177 | else: |
|---|
| 178 | ratio=0. |
|---|
| 179 | |
|---|
| 180 | # if the current value is -9999, we have to fill the gap |
|---|
| 181 | if(weather[k][t] == -9999): |
|---|
| 182 | # if gapmax is > 0 |
|---|
| 183 | # for short gap, we will try to interpolate |
|---|
| 184 | # between the last and the next |
|---|
| 185 | # defined values in the weather dataset |
|---|
| 186 | active_linear=1 |
|---|
| 187 | if(gapmax>0): |
|---|
| 188 | tlow=t-1 |
|---|
| 189 | while((tlow>=0)and(weather[k][tlow]==-9999)and((t-tlow)<=gapmax)): |
|---|
| 190 | tlow=tlow-1 |
|---|
| 191 | tup=t+1 |
|---|
| 192 | while((tup<len(weather[k]))and(weather[k][tup]==-9999)and((tup-t)<=gapmax)): |
|---|
| 193 | tup=tup+1 |
|---|
| 194 | # if the last and the next defined values are not too far |
|---|
| 195 | # we linearly interpolate with the weather dataset |
|---|
| 196 | if(tup-tlow<=gapmax+1): |
|---|
| 197 | if(tlow<0): |
|---|
| 198 | weather_gapfill[k].append(weather[k][tup]) |
|---|
| 199 | elif(tup>=len(weather[k])): |
|---|
| 200 | weather_gapfill[k].append(weather[k][tlow]) |
|---|
| 201 | else: |
|---|
| 202 | step=(float)(t-tlow)/(tup-tlow) |
|---|
| 203 | weather_gapfill[k].append(weather[k][tlow]*(1-step)+weather[k][tup]*step) |
|---|
| 204 | else: |
|---|
| 205 | active_linear=0 |
|---|
| 206 | # if gapmax == 0 or |
|---|
| 207 | # if the last and the next defined values are too far each other |
|---|
| 208 | # we use the clim dataset for filling the gap |
|---|
| 209 | if((gapmax==0)or(not active_linear)): |
|---|
| 210 | if (avg[k]==0): |
|---|
| 211 | step=(float)(tshift%diff_clim_weather)/diff_clim_weather |
|---|
| 212 | if(tshift<0.): |
|---|
| 213 | weather_gapfill[k].append(weather_clim_period_gapfill[k][0]) |
|---|
| 214 | elif((tshift/diff_clim_weather+1)>=(len(weather_clim_period_gapfill[k]))): |
|---|
| 215 | weather_gapfill[k].append(weather_clim_period_gapfill[k][len(weather_clim_period_gapfill[k])-1]) |
|---|
| 216 | else: |
|---|
| 217 | weather_gapfill[k].append(weather_clim_period_gapfill[k][tshift/diff_clim_weather]*(1-step)+weather_clim_period_gapfill[k][tshift/diff_clim_weather+1]*step) |
|---|
| 218 | else: |
|---|
| 219 | if(tshift<0.): |
|---|
| 220 | weather_gapfill[k].append(weather_clim_period_gapfill[k][0]*ratio) |
|---|
| 221 | elif((tshift/diff_clim_weather)>=(len(weather_clim_period_gapfill[k]))): |
|---|
| 222 | weather_gapfill[k].append(weather_clim_period_gapfill[k][len(weather_clim_period_gapfill[k])-1]*ratio) |
|---|
| 223 | else: |
|---|
| 224 | weather_gapfill[k].append(weather_clim_period_gapfill[k][tshift/diff_clim_weather]*ratio) |
|---|
| 225 | else: |
|---|
| 226 | weather_gapfill[k].append(weather[k][t]) |
|---|
| 227 | |
|---|
| 228 | return weather_gapfill,weather_clim_period_gapfill,weather_clim_period_nogap,weather_clim_period,clim_nogap,stat_vec |
|---|
| 229 | # end function |
|---|
| 230 | #------------------------------------------------------------------------------------------- |
|---|
