| 1 | #!/usr/bin/env cdat |
|---|
| 2 | # -*- coding: iso-8859-15 -*- |
|---|
| 3 | |
|---|
| 4 | import cdms, cdutil, genutil, os, sys, getopt |
|---|
| 5 | from cdms import MV |
|---|
| 6 | from dataproc.utils import cdat |
|---|
| 7 | from regrid import Regridder |
|---|
| 8 | |
|---|
| 9 | from Numeric import ones |
|---|
| 10 | |
|---|
| 11 | def usage(): |
|---|
| 12 | """Display usage information.""" |
|---|
| 13 | |
|---|
| 14 | man = """ |
|---|
| 15 | This script will get 2 SECHIBA or STOMATE restarts for ORCHIDEE code and will |
|---|
| 16 | produce a regridded restart with the grid of the first restart (and other non gridded variables) |
|---|
| 17 | and the values of the second one. |
|---|
| 18 | |
|---|
| 19 | Usage : |
|---|
| 20 | prompt> restart_interpol.py restart_outgrid.nc restart_invalues.nc' |
|---|
| 21 | This script regrids values in file 2 (namely restart_invlaues.nc) with |
|---|
| 22 | grid in restart_outgrid.nc. |
|---|
| 23 | |
|---|
| 24 | Examples : |
|---|
| 25 | > restart_interpol.py CF_CM1414A_20001231_sechiba_rest.nc CF_CM5PIRC11_29901231_sechiba_rest.nc |
|---|
| 26 | > restart_interpol.py CF_CM1414A_20001231_stomate_rest.nc CF_CM5PIRC11_29901231_stomate_rest.nc |
|---|
| 27 | > restart_interpol.py CF_CM5PIRC8_26601231_sechiba_rest.nc CF_OOL-CTRLPI-1_9_5_SPIN_v2ORC_43_01001231_sechiba_rest.nc |
|---|
| 28 | """ |
|---|
| 29 | print man |
|---|
| 30 | |
|---|
| 31 | options = [ |
|---|
| 32 | "help", |
|---|
| 33 | "version" ] |
|---|
| 34 | |
|---|
| 35 | # Extract arguments from the command line |
|---|
| 36 | try: |
|---|
| 37 | myopts, myargs = getopt.getopt(sys.argv[1:], 'Vh', longopts = options) |
|---|
| 38 | except getopt.GetoptError: |
|---|
| 39 | print 'getopt error' |
|---|
| 40 | usage() |
|---|
| 41 | sys.exit(1) |
|---|
| 42 | |
|---|
| 43 | variable="" |
|---|
| 44 | if (len(myargs) < 2): |
|---|
| 45 | print 'I need at least 2 files.' |
|---|
| 46 | usage() |
|---|
| 47 | sys.exit(1) |
|---|
| 48 | else: |
|---|
| 49 | file = myargs[:] |
|---|
| 50 | lf = len(file) |
|---|
| 51 | for i in range(2): |
|---|
| 52 | if (os.path.exists(file[i])==False): |
|---|
| 53 | print 'file number ',i+1,' : ',file[i],\ |
|---|
| 54 | ' doesn''t exist !' |
|---|
| 55 | sys.exit(1) |
|---|
| 56 | |
|---|
| 57 | f1=cdms.open( file[0],'r') |
|---|
| 58 | f2=cdms.open( file[1],'r') |
|---|
| 59 | listv=f1.listvariables() |
|---|
| 60 | listv2=f2.listvariables() |
|---|
| 61 | #listv.sort() |
|---|
| 62 | #listv2.sort() |
|---|
| 63 | print listv |
|---|
| 64 | print listv2 |
|---|
| 65 | |
|---|
| 66 | |
|---|
| 67 | f3=cdms.open( "regrid_"+file[1],'w') |
|---|
| 68 | |
|---|
| 69 | onet=None |
|---|
| 70 | |
|---|
| 71 | missing_value = 1e+20 |
|---|
| 72 | |
|---|
| 73 | for var in listv: |
|---|
| 74 | print var |
|---|
| 75 | if (var in ["nav_lon","nav_lat","nav_lev","time_steps","time", |
|---|
| 76 | "lat","bounds_lat","lon","bounds_lon","lev","day_counter", |
|---|
| 77 | "dt_days", "routingcounter", "veget_year", "date" ]): |
|---|
| 78 | var1=f1(var) |
|---|
| 79 | f3.write(var1) |
|---|
| 80 | |
|---|
| 81 | print "remplacé." |
|---|
| 82 | continue |
|---|
| 83 | |
|---|
| 84 | var1=f1(var) |
|---|
| 85 | try: |
|---|
| 86 | var2=f2(var) |
|---|
| 87 | except: |
|---|
| 88 | f3.write(var1) |
|---|
| 89 | print "remplacé." |
|---|
| 90 | continue |
|---|
| 91 | |
|---|
| 92 | # var2.info() |
|---|
| 93 | if (len(var1.shape) > 3): |
|---|
| 94 | var_1=f1(var, squeeze=1) |
|---|
| 95 | tmp=cdms.createVariable(var1) |
|---|
| 96 | var_2=f2(var, squeeze=1) |
|---|
| 97 | # print var1.id,var |
|---|
| 98 | |
|---|
| 99 | # if (var1.id == 'day_counter'): |
|---|
| 100 | # var1[0,0,0,0]=var_2 |
|---|
| 101 | # continue |
|---|
| 102 | |
|---|
| 103 | # print tmp.shape,var_2,type(var_2) |
|---|
| 104 | # axlist2=var_2.getAxisList() |
|---|
| 105 | # print tmp.shape,var_2,len(axlist2) |
|---|
| 106 | # if (len(axlist2) > 0 ): |
|---|
| 107 | |
|---|
| 108 | if (len(var_2)!=1) : |
|---|
| 109 | if (len(var_2.shape) > 2): |
|---|
| 110 | # print tmp.shape, var_2.shape,var2.shape |
|---|
| 111 | if ( len(tmp.getAxis(0)) == 1 ): |
|---|
| 112 | # print tmp.getAxis(0),var2.getAxis(0) |
|---|
| 113 | |
|---|
| 114 | if ( len(var_2.getAxis(0)) > 0 ): |
|---|
| 115 | # print var_2.getAxis(0)[:] |
|---|
| 116 | shift=( var_2.getAxis(0)[0] > 0 ) |
|---|
| 117 | outgrid=var_1[0,:,:].getGrid() |
|---|
| 118 | ingrid=var_2[0,:,:].getGrid() |
|---|
| 119 | regridfunc = Regridder(ingrid, outgrid) |
|---|
| 120 | sumpft=var_1[0,:,:]*0. |
|---|
| 121 | for fpft in var_2.getAxis(0): |
|---|
| 122 | if ( shift ): |
|---|
| 123 | pft=int(fpft) - int(var_2.getAxis(0)[0]) |
|---|
| 124 | else: |
|---|
| 125 | pft=int(fpft) |
|---|
| 126 | |
|---|
| 127 | # print pft |
|---|
| 128 | tmp_ = tmp[0,pft,:,:] |
|---|
| 129 | v2 = var_2[pft,:,:] |
|---|
| 130 | v1 = var_1[pft,:,:] |
|---|
| 131 | # mask_undef = MV.where(v2 > missing_value, 1, 0) |
|---|
| 132 | # v2_mask=MV.masked_array(v2, mask_undef) |
|---|
| 133 | # varmax2=max(max(abs(v2_mask))) |
|---|
| 134 | # v2=v2_mask/varmax2 |
|---|
| 135 | tmp__ = regridfunc(v2) |
|---|
| 136 | # tmp[0,pft,:,:]=tmp__*varmax2 |
|---|
| 137 | tmp[0,pft,:,:]=tmp__ |
|---|
| 138 | if (var in ["veget_max","veget"]): |
|---|
| 139 | sumpft[:,:]=sumpft[:,:]+tmp__ |
|---|
| 140 | |
|---|
| 141 | if (var in ["veget_max"]): |
|---|
| 142 | # On veut : |
|---|
| 143 | #somm(tmp[0,pft,:,:])+frac_nobio=1 |
|---|
| 144 | vegnorm=1-frac_nobio |
|---|
| 145 | veget_max=tmp[0,:,:,:] |
|---|
| 146 | for fpft in var_2.getAxis(0): |
|---|
| 147 | if ( shift ): |
|---|
| 148 | pft=int(fpft) - int(var_2.getAxis(0)[0]) |
|---|
| 149 | else: |
|---|
| 150 | pft=int(fpft) |
|---|
| 151 | veget_max[pft,:,:]=MV.where(sumpft[:,:] == 0., 0., veget_max[pft,:,:] * vegnorm / sumpft[:,:]) |
|---|
| 152 | tmp[0,pft,:,:]= veget_max[pft,:,:] |
|---|
| 153 | |
|---|
| 154 | if (var in ["veget"]): |
|---|
| 155 | for fpft in var_2.getAxis(0): |
|---|
| 156 | if ( shift ): |
|---|
| 157 | pft=int(fpft) - int(var_2.getAxis(0)[0]) |
|---|
| 158 | else: |
|---|
| 159 | pft=int(fpft) |
|---|
| 160 | |
|---|
| 161 | tmp[0,pft,:,:]= MV.where(tmp[0,pft,:,:] > veget_max[pft,:,:], veget_max[pft,:,:], tmp[0,pft,:,:]) |
|---|
| 162 | |
|---|
| 163 | else: |
|---|
| 164 | tmp[0,0,:,:]=var_2.regrid( var_1.getGrid() ) |
|---|
| 165 | else: |
|---|
| 166 | print "Erreur ! " |
|---|
| 167 | print tmp.getAxis(0),var2.getAxis(0) |
|---|
| 168 | # print var1.getAxisList() |
|---|
| 169 | |
|---|
| 170 | if ( len(var_2.getAxis(0)) > 0 ): |
|---|
| 171 | # print var_2.getAxis(0)[:] |
|---|
| 172 | shift=( var_2.getAxis(0)[0] > 0 ) |
|---|
| 173 | for fpft in var_2.getAxis(0): |
|---|
| 174 | if ( shift ): |
|---|
| 175 | pft=int(fpft) - int(var_2.getAxis(0)[0]) |
|---|
| 176 | else: |
|---|
| 177 | pft=int(fpft) |
|---|
| 178 | |
|---|
| 179 | # print pft |
|---|
| 180 | tmp_ = tmp[0,pft,:,:] |
|---|
| 181 | v2 = var_2[pft,:,:] |
|---|
| 182 | v1 = var_1[0,pft,:,:] |
|---|
| 183 | # mask_undef = MV.where(v2 > missing_value, 1, 0) |
|---|
| 184 | # v2_mask=MV.masked_array(v2, mask_undef) |
|---|
| 185 | # varmax2=max(abs(v2_mask)) |
|---|
| 186 | # v2=v2_mask/varmax2 |
|---|
| 187 | tmp__ = v2.regrid( v1.getGrid() ) |
|---|
| 188 | # tmp[0,pft,:,:]=tmp__*varmax2 |
|---|
| 189 | tmp[0,pft,:,:]=tmp__ |
|---|
| 190 | |
|---|
| 191 | else: |
|---|
| 192 | tmp[0,0,:,:]=var_2.regrid( var_1.getGrid() ) |
|---|
| 193 | if (var in ["frac_nobio"]): |
|---|
| 194 | frac_nobio=tmp[0,0,:,:] |
|---|
| 195 | else: |
|---|
| 196 | tmp[0,0,0,0]=var_2 |
|---|
| 197 | |
|---|
| 198 | var1=tmp |
|---|
| 199 | |
|---|
| 200 | # var1.info() |
|---|
| 201 | t=var1.getTime() |
|---|
| 202 | print 'len of t:',t |
|---|
| 203 | if t is not None: |
|---|
| 204 | if len(t)!=1: |
|---|
| 205 | raise 'Error too much time' |
|---|
| 206 | if onet is None: |
|---|
| 207 | onet=t.clone() |
|---|
| 208 | else: |
|---|
| 209 | if t[0]!=onet[0]: |
|---|
| 210 | raise 'different time' |
|---|
| 211 | |
|---|
| 212 | f3.write(var1) |
|---|
| 213 | else: |
|---|
| 214 | print "impossible ?" |
|---|
| 215 | # var2=f2(var) |
|---|
| 216 | # f3.write(var2) |
|---|
