source: tags/ORCHIDEE_1_9_5/ORCHIDEE_OL/Utilitaire/Regrid_restart/fonctions_orchidee.py @ 8

Last change on this file since 8 was 8, checked in by orchidee, 14 years ago

import first tag equivalent to CVS orchidee_1_9_5 + OOL_1_9_5
File size: 11.6 KB
Line 
1# -*- coding: iso-8859-15 -*-
2import os
3import Numeric
4
5import cdms
6import MV
7import vcs
8
9# filtre un flichier créé par cdms avec un ncdump, puis un ncgen.
10def filtre_dump (fichier):
11    os.spawnlp(os.P_WAIT,"filtre_dump","filtre_dump",fichier,fichier)
12
13
14# calcul le tableau d'index d'un axe de points de terre
15# (ie kindex2D[i,j] = land[kij]),
16# ainsi que les indices inversés
17# (ie ( iland[ kindex2D[i,j]],jland[kindex2D[i,j]] ) = ( i,j )).
18# l'algorithme est celui utilisé dans ORCHIDEE.
19def calcul_kindex(kindex2D,iland,jland,axisLand,iim,jjm):
20    kij=0
21    for ki in axisLand:
22        j = int((ki-1)/iim)
23        i = int(((ki-1) - j*iim))
24        kindex2D[i,j] = int(kij)
25        iland[kij] = i
26        jland[kij] = j
27        #kindex2Dt[jjm-1-j,i] = kij
28        kij = kij+1
29    return
30
31# calcul le vecteur des points de terre à partir du tableau des index
32# (ie land[kij] = kindex2D[i,j])
33# l'algorithme est celui utilisé dans ORCHIDEE.
34def calcul_axisLand(kindex2D,iim,jjm,missing_value):
35    axisLand = []
36    kij=0
37    for j in range(jjm):
38        for i in range(iim):
39            if (kindex2D[i,j] < missing_value):
40                ki = int(j*iim + i + 1)
41                axisLand.append(ki)
42                # vérification :
43                if (j != (ki-1)/iim or i != ((ki-1) - j*iim)):
44                    print "Erreur calcul_axisLand au ",kij,"-ième point de terre pour (i,j) = ",i,j,"  et land_i = ",ki
45                kij = kij+1
46    return(MV.array(axisLand))
47
48
49# Converti une variable d'axes [t,z,land] (avec land les points de terre)
50# en une variable 2D (lat,lon)
51def conversion2D(VarIn, kindex, nbindex, iim, jjm, missing_value):
52    VarOut = missing_value*MV.ones((jjm,iim))
53    for j in range(jjm):
54        for i in range(iim):
55            if (kindex[i,j] < missing_value):
56                #print i,j,kindex[j,i],(jjm-1)-j
57                VarOut[j,i] = VarIn[kindex[i,j]]
58    return(VarOut)
59
60# Converti une variable 2D (lat,lon) en
61# une variable d'axe land (les points de terre).
62def conversion2Dinv(VarIn, kindex, nbindex, iim, jjm, missing_value):
63    VarOut = missing_value*MV.ones(nbindex)
64    ii=0
65    for j in range(jjm):
66        for i in range(iim):
67                if (kindex[i,j] < missing_value):
68                    VarOut[ii] = VarIn[j,i]
69                    ii=ii+1
70    return(VarOut)
71
72# Pour les graphiques : même fonction, avec la latitude inversée.
73def conversion2DGr(VarIn, kindex, nbindex, iim, jjm, missing_value):
74    VarOutGr = missing_value*Numeric.ones((jjm,iim))
75    for j in range(jjm):
76        for i in range(iim):
77                if (kindex[i,j] < missing_value):
78                        #print i,j,kindex[j,i],(jjm-1)-j
79                        VarOutGr[(jjm-1)-j,i] = VarIn[kindex[i,j]]
80    return(VarOutGr)
81
82# Pour les graphiques : inverse la latitude d'une variable 2D
83def inverse2DGr(VarIn, kindex, nbindex, iim, jjm, missing_value):
84    VarOutGr = missing_value*Numeric.ones((jjm,iim))
85    for j in range(jjm):
86        for i in range(iim):
87                if (kindex[i,j] < missing_value):
88                    VarOutGr[(jjm-1)-j,i] = VarIn[j,i]
89    return(VarOutGr)
90
91# Calcul l'échelle iso d'une variable Var, dans un graphique grd, avec une valeur missing_value.
92# On décale les grandeurs min/max de decal %
93def errlevGr(Var, grd, decal, missing_value):   
94    minv = cdms.MV.minimum(Var)
95    maxv = cdms.MV.maximum(Var)
96    if (minv > 0):
97        minv = minv*(1-decal/100.)
98    else:
99        minv = minv*(1+decal/100.)
100    if (maxv > 0):
101        maxv = maxv*(1+decal/100.)
102    else:
103        maxv = maxv*(1-decal/100.)
104    diffv = (maxv - minv)/10. 
105    errlevs = vcs.mkevenlevels(minv,maxv,12)
106    errlevs.append(missing_value)
107    errlevs.insert(0,missing_value)
108    cols = vcs.getcolors(errlevs)
109
110    return(errlevs, cols)
111
112
113
114# The input VarIn is a land variable.
115# The output variable is 2D ( [jjm,iim] or [jmin:(jmax+1),imin:(imax+1)] )
116#            and if calcland==1, it returns too :
117#                zoomkindex, zoomiim, zoomjjm, nlandz,
118#                invindex, ifenetre =(imin,imax,jmin,jmax)
119def zoom_land(VarIn, kindex, nbindex, lat, lon, iim, jjm, limits, missing_value, calcland=1):
120# where limits is in the order : [limit_W, limit_E, limit_N, limit_S]
121    VarOut = missing_value*Numeric.ones((jjm,iim))
122    # inverse of kindex
123    invindex = Numeric.zeros((nbindex,2))
124    imin = iim+1
125    imax = -1
126    jmin = jjm+1
127    jmax = -1
128    for j in range(jjm):
129        for i in range(iim):
130                if (kindex[i,j] < missing_value):
131                        if ( limits[0] <= lon[j,i] and lon[j,i] <= limits[1]
132                             and limits[3] <= lat[j,i] and lat[j,i] <= limits[2] ):
133                                imin = min(imin, i)
134                                imax = max(imax, i)
135                                jmin = min(jmin, j)
136                                jmax = max(jmax, j)
137                                invindex[kindex[i,j],0]=j
138                                invindex[kindex[i,j],1]=i
139                                VarOut[j,i] = VarIn[kindex[i,j]]
140    zoomiim = imax - imin +1
141    zoomjjm = jmax - jmin +1
142    #print imin, imax, jmin, jmax
143    #print zoomiim, zoomjjm
144    # pourquoi le +1 ???? (python ??) => nécessaire en tout cas.
145    # à cause des lat/lon qui doivent anglober la zone zoomée
146    VarOutZoom = VarOut[jmin:(jmax+1),imin:(imax+1)]
147    if (calcland==1):
148        zoomkindex = Numeric.ones((zoomiim,zoomjjm))*int(missing_value)
149        nlandz=0
150        for j in range(jjm):
151            for i in range(iim):
152                if (kindex[i,j] < missing_value):
153                    if ( limits[0] <= lon[j,i] and lon[j,i] <= limits[1]
154                         and limits[3] <= lat[j,i] and lat[j,i] <= limits[2] ):
155                        zoomkindex[i-imin, j-jmin] = kindex[i,j]
156                        nlandz=nlandz+1
157        print "zoom : dimension", zoomiim, zoomjjm, nlandz
158        # print "zoom : nouvel index ",zoomkindex
159        return(VarOutZoom, zoomkindex, zoomiim, zoomjjm, nlandz,
160               invindex, (imin, imax, jmin, jmax))
161    else:
162        return(VarOutZoom)
163
164
165# The input VarIn is a land variable.
166# The output variable is 2D ( [jjm,iim] or [jmin:(jmax+1),imin:(imax+1)] )
167# Simplified version of precedent one.
168def zoom_land_simpl(VarIn, kindex, nbindex, zoomiim, zoomjjm, ifenetre, missing_value):
169# where limits is in the order : [limit_W, limit_E, limit_N, limit_S]
170    (imin, imax, jmin, jmax) = ifenetre
171    VarOutZoom = missing_value*Numeric.ones((zoomjjm,zoomiim))
172    for i in range(zoomiim):
173        for j in range(zoomjjm):
174            ki=kindex[imin+i,jmin+j]
175            if (ki < missing_value):
176                VarOutZoom[j,i] = VarIn[ki]
177    return(VarOutZoom)
178
179
180# The input VarIn is a 2D variable.
181# the output variable is 2D ( [jjm,iim] or [jmin:(jmax+1),imin:(imax+1)] )
182#            and if calcland==1, it returns too :
183#                zoomkindex, zoomiim, zoomjjm, nlandz,
184#                invindex, ifenetre =(imin,imax,jmin,jmax)
185def zoom(VarIn, kindex, nbindex, lat, lon, iim, jjm, limits, missing_value, calcland=1):
186# where limits is in the order : [limit_W, limit_E, limit_N, limit_S]
187    imin = iim+1
188    imax = -1
189    jmin = jjm+1
190    jmax = -1
191    for j in range(jjm):
192        for i in range(iim):
193                if (kindex[i,j] < missing_value):
194                        if ( limits[0] <= lon[j,i] and lon[j,i] <= limits[1]
195                             and limits[3] <= lat[j,i] and lat[j,i] <= limits[2] ):
196                                imin = min(imin, i)
197                                imax = max(imax, i)
198                                jmin = min(jmin, j)
199                                jmax = max(jmax, j)
200                                invindex[kindex[i,j],0]=j
201                                invindex[kindex[i,j],1]=i
202    zoomiim = imax - imin +1
203    zoomjjm = jmax - jmin +1
204    # pourquoi le +1 ???? (python ??) => nécessaire en tout cas.
205    # à cause des lat/lon qui doivent anglober la zone zoomée
206    VarOutZoom = VarIn[jmin:(jmax+1),imin:(imax+1)]
207    if (calcland == 1):
208        zoomkindex = Numeric.ones((zoomiim,zoomjjm))*int(missing_value)
209        nlandz=0
210        for j in range(jjm):
211            for i in range(iim):
212                if (kindex[i,j] < missing_value):
213                    if ( limits[0] <= lon[j,i] and lon[j,i] <= limits[1]
214                         and limits[3] <= lat[j,i] and lat[j,i] <= limits[2] ):
215                        zoomkindex[i-imin, j-jmin] = kindex[i,j]
216                        nlandz=nlandz+1
217        print "zoom : dimension",zoomiim, zoomjjm, nlandz
218        # print "zoom : nouvel index ",zoomkindex
219        return(VarOutZoom, zoomkindex, zoomiim, zoomjjm, nlandz,
220               invindex, (imin, imax, jmin, jmax))
221    else:
222        return(VarOutZoom)
223
224# The input VarIn is a 2D variable.
225# the output variable is 2D [jmin:(jmax+1),imin:(imax+1)]
226def zoom2D(VarIn, lat, lon, iim, jjm, limits, missing_value):
227# where limits is in the order : [limit_W, limit_E, limit_N, limit_S]
228    imin = iim+1
229    imax = -1
230    jmin = jjm+1
231    jmax = -1
232    for j in range(jjm):
233        for i in range(iim):
234            if ( limits[0] <= lon[j,i] and lon[j,i] <= limits[1]
235                 and limits[3] <= lat[j,i] and lat[j,i] <= limits[2] ):
236                imin = min(imin, i)
237                imax = max(imax, i)
238                jmin = min(jmin, j)
239                jmax = max(jmax, j)
240    zoomiim = imax - imin +1
241    zoomjjm = jmax - jmin +1
242    VarOutZoom = VarIn[jmin:(jmax+1),imin:(imax+1)]
243    return(VarOutZoom)
244
245# The input VarIn is a 2D variable.
246# the output variable is 2D ( [jjm,iim] or [jmin:(jmax+1),imin:(imax+1)] )
247# Simplified version of precedent one.
248def zoom_simpl(VarIn, zoomiim, zoomjjm, ifenetre, missing_value):
249# where limits is in the order : [limit_W, limit_E, limit_N, limit_S]
250    (imin, imax, jmin, jmax) = ifenetre
251    VarOutZoom = missing_value*Numeric.ones((zoomjjm,zoomiim))
252    VarOutZoom = VarIn[jmin:(jmax+1),imin:(imax+1)]
253    return(VarOutZoom)
254
255
256tab_mois=[31,28,31,30,31,30,31,31,30,31,30,31]
257
258# Conversion Month/Day => Num of Day in year
259def monthday_to_numday(month,day):
260#Tableau definissant la duree des mois dans un calendrier de 365 jours
261    numday=sum(tab_mois[0:(month-1)])+day
262    return(numday)
263   
264# Rajoute dans le fichier "fichier" les tableaux bounds_lat et bounds_lon
265# pour le rendre compatible avec la convention CF.
266
267# def rajoutBounds(fichier,lat,lon,nland,land,kindex,nbindex,missing_value)
268#     bounds_lon = Numeric.zeros((nland,,4))
269#     bounds_lat = Numeric.zeros((nland,,4))
270#     for j in range(jjm):
271#       for i in range(iim):
272#               if (kindex[i,j] < missing_vlaue):
273#                         #print i,j,kindex[i,j],(jjm-1)-j
274#                       VarOut[i,j] = VarIn[kindex[i,j]]
275#     return(VarOut)
276
277
278# Construction des tableaux de coordonnées 2D nav_lat et nav_lon
279def Construit_nav_l(ax_lat, ax_lon, missing_v):
280    Mones=Numeric.ones((len(ax_lat),1))
281    Slon=Numeric.array([ax_lon])
282    nvlon=MV.dot(Mones,Slon)
283    nav_lon=cdms.createVariable(nvlon, 
284                                typecode = cdms.MV.Float16,
285                                fill_value = missing_v,
286                                attributes={'name': 'nav_lon',
287                                            'valid_min': -180.,
288                                            'valid_max': 180.,
289                                            'long_name': 'Longitude',
290                                            'units': 'degrees_east' },
291                                axes = [ ax_lat, ax_lon ],
292                                id = 'nav_lon')
293
294    Mones=Numeric.swapaxes(Numeric.ones((len(ax_lon),1)),0,1)
295    tmplat=Numeric.swapaxes(Numeric.array([ax_lat]),0,1)
296    nvlat=MV.dot(tmplat, Mones)
297    nav_lat=cdms.createVariable(nvlat, 
298                                typecode = cdms.MV.Float16,
299                                fill_value = missing_v,
300                                attributes={'name': 'nav_lat',
301                                            'valid_min': -90.,
302                                            'valid_max': 90.,
303                                            'long_name': 'Latitude',
304                                            'units': 'degrees_north' },
305                                axes = [ ax_lat, ax_lon ],
306                                id = 'nav_lat')
307    return (nav_lat, nav_lon)
Note: See TracBrowser for help on using the repository browser.