| 1 | # Reconstructing ORCHIDEE's forcing files as matrices for a gien area |
|---|
| 2 | # L. Fita, CIMA. December 2017 |
|---|
| 3 | # |
|---|
| 4 | import numpy as np |
|---|
| 5 | from netCDF4 import Dataset as NetCDFFile |
|---|
| 6 | import os |
|---|
| 7 | import re |
|---|
| 8 | import numpy.ma as ma |
|---|
| 9 | # Importing generic tools file 'generic_tools.py' |
|---|
| 10 | import generic_tools as gen |
|---|
| 11 | import nc_var_tools as ncvar |
|---|
| 12 | import subprocess as sub |
|---|
| 13 | import module_ForSci as Sci |
|---|
| 14 | from optparse import OptionParser |
|---|
| 15 | |
|---|
| 16 | parser = OptionParser() |
|---|
| 17 | parser.add_option("-y", "--year", dest="year", help="year to process", \ |
|---|
| 18 | metavar="VALUE") |
|---|
| 19 | |
|---|
| 20 | (opts, args) = parser.parse_args() |
|---|
| 21 | |
|---|
| 22 | ####### ###### ##### #### ### ## # |
|---|
| 23 | |
|---|
| 24 | filen = 'cruncep_halfdeg_' + opts.year + '.nc' |
|---|
| 25 | |
|---|
| 26 | # Variable whcih provides the indices of a 1D vector from the dimy, dimx space |
|---|
| 27 | indvar = 'land' |
|---|
| 28 | |
|---|
| 29 | # 2D longitude, latitude matrices |
|---|
| 30 | lonvar = 'nav_lon' |
|---|
| 31 | latvar = 'nav_lat' |
|---|
| 32 | |
|---|
| 33 | # Range to retrieve |
|---|
| 34 | Xmin=-90.25 |
|---|
| 35 | #Xmin='all' |
|---|
| 36 | Xmax=-33.25 |
|---|
| 37 | Ymin=-67.25 |
|---|
| 38 | Ymax=15.25 |
|---|
| 39 | |
|---|
| 40 | # Variables to reconstruct |
|---|
| 41 | #variable = 'all' |
|---|
| 42 | variable = 'time,LWdown,Qair,PSurf,SWdown,Tair,Wind_E,Wind_N,Snowf,Rainf' |
|---|
| 43 | |
|---|
| 44 | # Resolution |
|---|
| 45 | resX = 0.5 |
|---|
| 46 | resY = 0.5 |
|---|
| 47 | |
|---|
| 48 | # Minimum difference from matrix to localized point |
|---|
| 49 | maxdiff = 0.05 |
|---|
| 50 | |
|---|
| 51 | # Projection of the matrix |
|---|
| 52 | matProj = 'latlon' |
|---|
| 53 | |
|---|
| 54 | ####### ####### |
|---|
| 55 | ## MAIN |
|---|
| 56 | ####### |
|---|
| 57 | main = 'ORforcing_reconstruct.py' |
|---|
| 58 | fname = 'ORforcing_reconstruct.py' |
|---|
| 59 | |
|---|
| 60 | onc = NetCDFFile(filen, 'r') |
|---|
| 61 | |
|---|
| 62 | availProj = ['latlon'] |
|---|
| 63 | |
|---|
| 64 | ofilen = 'reconstruct_matrix_' + opts.year + '.nc' |
|---|
| 65 | |
|---|
| 66 | ncvars = onc.variables.keys() |
|---|
| 67 | |
|---|
| 68 | varstocheck = [indvar, lonvar, latvar] |
|---|
| 69 | for vn in varstocheck: |
|---|
| 70 | if not gen.searchInlist(ncvars, vn): |
|---|
| 71 | print gen.errormsg |
|---|
| 72 | print ' ' + main + ": file '" + filen + "' does not have variable '" + vn + \ |
|---|
| 73 | "' !!" |
|---|
| 74 | print ' available ones:', ncvars |
|---|
| 75 | quit(-1) |
|---|
| 76 | |
|---|
| 77 | oind = onc.variables[indvar] |
|---|
| 78 | olon = onc.variables[lonvar] |
|---|
| 79 | olat = onc.variables[latvar] |
|---|
| 80 | |
|---|
| 81 | vecdimn = oind.dimensions[0] |
|---|
| 82 | Xdimn = olon.dimensions[1] |
|---|
| 83 | Ydimn = olon.dimensions[0] |
|---|
| 84 | |
|---|
| 85 | # Fortran indices, first 1 |
|---|
| 86 | indv = oind[:] |
|---|
| 87 | lonv = olon[:] |
|---|
| 88 | latv = olat[:] |
|---|
| 89 | indv = indv - 1 |
|---|
| 90 | |
|---|
| 91 | dimx = lonv.shape[1] |
|---|
| 92 | dimy = lonv.shape[0] |
|---|
| 93 | |
|---|
| 94 | veclon = np.zeros((len(indv)), dtype=np.float) |
|---|
| 95 | veclat = np.zeros((len(indv)), dtype=np.float) |
|---|
| 96 | |
|---|
| 97 | # Construct veclon, veclat |
|---|
| 98 | for iid in range(len(indv)): |
|---|
| 99 | iy = int((indv[iid]-1)/dimx) |
|---|
| 100 | ix = indv[iid] - iy*dimx - 1 |
|---|
| 101 | veclon[iid] = lonv[iy,ix] |
|---|
| 102 | veclat[iid] = latv[iy,ix] |
|---|
| 103 | |
|---|
| 104 | if not gen.searchInlist(availProj, matProj): |
|---|
| 105 | print errormsg |
|---|
| 106 | print ' ' + fname + ": projection '" + matProj + "' not available !!" |
|---|
| 107 | print ' available ones:', availProj |
|---|
| 108 | quit(-1) |
|---|
| 109 | |
|---|
| 110 | ncdims = onc.dimensions |
|---|
| 111 | |
|---|
| 112 | if variable == 'all': |
|---|
| 113 | varns = ncvars |
|---|
| 114 | else: |
|---|
| 115 | varns = variable.split(',') |
|---|
| 116 | for vn in varns: |
|---|
| 117 | if not gen.searchInlist(ncvars, vn): |
|---|
| 118 | print gen.errormsg |
|---|
| 119 | print ' ' + fname + ": file '" + ncfile + "' does not have " + \ |
|---|
| 120 | " variable '" + vn + "' !!" |
|---|
| 121 | print ' available ones:', ncvars |
|---|
| 122 | quit(-1) |
|---|
| 123 | |
|---|
| 124 | if gen.searchInlist(olon.ncattrs(), 'units'): |
|---|
| 125 | xunits = olon.getncattr('units') |
|---|
| 126 | else: |
|---|
| 127 | xunits = '-' |
|---|
| 128 | if gen.searchInlist(olat.ncattrs(), 'units'): |
|---|
| 129 | yunits = olat.getncattr('units') |
|---|
| 130 | else: |
|---|
| 131 | yunits = '-' |
|---|
| 132 | |
|---|
| 133 | if type(Xmin) == type('2') and Xmin == 'all': |
|---|
| 134 | Xmin = np.min(lonv) |
|---|
| 135 | Xmax = np.max(lonv) |
|---|
| 136 | Ymin = np.min(latv) |
|---|
| 137 | Ymax = np.max(latv) |
|---|
| 138 | |
|---|
| 139 | # Matrix values |
|---|
| 140 | if matProj == 'latlon': |
|---|
| 141 | dimx = int((Xmax - Xmin+resX)/resX) |
|---|
| 142 | dimy = int((Ymax - Ymin+resY)/resY) |
|---|
| 143 | |
|---|
| 144 | matindt, matXt, matYt = Sci.module_scientific.reconstruct_matrix(vectorxpos=veclon, \ |
|---|
| 145 | vectorypos=veclat, dvec=veclon.shape[0], xmin=Xmin, xmax=Xmax, ymin=Ymin,ymax=Ymax,\ |
|---|
| 146 | dmatx=dimx, dmaty=dimy, matproj=matProj, maxdiff=maxdiff) |
|---|
| 147 | |
|---|
| 148 | matind = matindt.transpose() |
|---|
| 149 | # Fortran like, First 1 |
|---|
| 150 | matind = np.where(matind != -1, matind - 1, matind) |
|---|
| 151 | |
|---|
| 152 | # Creation of file |
|---|
| 153 | onewnc = NetCDFFile(ofilen, 'w') |
|---|
| 154 | |
|---|
| 155 | # Dimensions |
|---|
| 156 | newdim = onewnc.createDimension('x', dimx) |
|---|
| 157 | newdim = onewnc.createDimension('y', dimy) |
|---|
| 158 | |
|---|
| 159 | # Variable-dimension |
|---|
| 160 | newvar = onewnc.createVariable('lon', 'f8', ('y', 'x')) |
|---|
| 161 | newvar[:] = matXt.transpose() |
|---|
| 162 | ncvar.basicvardef(newvar, 'lon', 'Longitude', 'degrees_east') |
|---|
| 163 | |
|---|
| 164 | newvar = onewnc.createVariable('lat', 'f8', ('y', 'x')) |
|---|
| 165 | newvar[:] = matYt.transpose() |
|---|
| 166 | ncvar.basicvardef(newvar, 'lat', 'Latitude', 'degrees_north') |
|---|
| 167 | onewnc.sync() |
|---|
| 168 | |
|---|
| 169 | # Variable indices |
|---|
| 170 | newvar = onewnc.createVariable('vec1D_matind', 'i', ('y', 'x'), fill_value=-1) |
|---|
| 171 | newvar[:] = matind |
|---|
| 172 | ncvar.basicvardef(newvar, 'vec1D_matind', 'matrix with the equivalencies from 1D ' + \ |
|---|
| 173 | 'vector indices', '-') |
|---|
| 174 | ncvar.set_attribute(newvar, 'coordinates', 'lon lat') |
|---|
| 175 | |
|---|
| 176 | # Looking for equivalencies in the 1D vector |
|---|
| 177 | matlonlat = matind.copy() |
|---|
| 178 | for j in range(dimy): |
|---|
| 179 | for i in range(dimx): |
|---|
| 180 | if matind[j,i] != -1: |
|---|
| 181 | matlonlat[j,i] = indv[matind[j,i]] |
|---|
| 182 | |
|---|
| 183 | newvar = onewnc.createVariable('lonlat_matind', 'i', ('y', 'x'), fill_value=-1) |
|---|
| 184 | newvar[:] = matlonlat |
|---|
| 185 | ncvar.basicvardef(newvar, 'lonlat_matind', 'matrix with the equivalencies from ' + \ |
|---|
| 186 | '2D lon, lat matrices', '-') |
|---|
| 187 | ncvar.set_attribute(newvar, 'coordinates', 'lon lat') |
|---|
| 188 | onewnc.sync() |
|---|
| 189 | |
|---|
| 190 | # Getting variables |
|---|
| 191 | for vn in varns: |
|---|
| 192 | if not onewnc.variables.has_key(vn): |
|---|
| 193 | ovar = onc.variables[vn] |
|---|
| 194 | indn = ovar.dimensions |
|---|
| 195 | vardims = [] |
|---|
| 196 | shapevar = [] |
|---|
| 197 | for dn in indn: |
|---|
| 198 | if not gen.searchInlist(onewnc.dimensions, dn) and dn != Xdimn and \ |
|---|
| 199 | dn != Ydimn: |
|---|
| 200 | if onc.dimensions[dn].isunlimited(): |
|---|
| 201 | newdim = onewnc.createDimension(dn, None) |
|---|
| 202 | else: |
|---|
| 203 | newdim = onewnc.createDimension(dn, len(onc.dimensions[dn])) |
|---|
| 204 | |
|---|
| 205 | if dn == vecdimn: |
|---|
| 206 | vardims.append('y') |
|---|
| 207 | vardims.append('x') |
|---|
| 208 | shapevar.append(dimy) |
|---|
| 209 | shapevar.append(dimx) |
|---|
| 210 | else: |
|---|
| 211 | vardims.append(dn) |
|---|
| 212 | shapevar.append(len(onc.dimensions[dn])) |
|---|
| 213 | |
|---|
| 214 | if ovar.dtype == type(int(2)): |
|---|
| 215 | newvar= onewnc.createVariable(vn,ncvar.nctype(ovar.dtype),tuple(vardims),\ |
|---|
| 216 | fill_value=gen.fillValueI) |
|---|
| 217 | varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueI |
|---|
| 218 | elif ovar.dtype == type(np.int32(2)): |
|---|
| 219 | newvar= onewnc.createVariable(vn,ncvar.nctype(ovar.dtype),tuple(vardims),\ |
|---|
| 220 | fill_value=gen.fillValueI) |
|---|
| 221 | varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueI |
|---|
| 222 | elif ovar.dtype == type(np.int64(2)): |
|---|
| 223 | newvar= onewnc.createVariable(vn,ncvar.nctype(ovar.dtype),tuple(vardims),\ |
|---|
| 224 | fill_value=gen.fillValueI) |
|---|
| 225 | varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueI |
|---|
| 226 | elif ovar.dtype == type(np.float(2.)): |
|---|
| 227 | newvar= onewnc.createVariable(vn,ncvar.nctype(ovar.dtype),tuple(vardims),\ |
|---|
| 228 | fill_value=gen.fillValueF) |
|---|
| 229 | varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueF |
|---|
| 230 | elif ovar.dtype == type(np.float32(2.)): |
|---|
| 231 | newvar= onewnc.createVariable(vn,ncvar.nctype(ovar.dtype),tuple(vardims),\ |
|---|
| 232 | fill_value=gen.fillValueF) |
|---|
| 233 | varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueF |
|---|
| 234 | else: |
|---|
| 235 | print gen.errormsg |
|---|
| 236 | print ' ' + fname + ': variable type:', ovar.dtype, ' not ready !!' |
|---|
| 237 | quit(-1) |
|---|
| 238 | |
|---|
| 239 | print ' reconstructing:', vn, '...' |
|---|
| 240 | # Filling variable. It would be faster if we can avoid this loop... I'm feeling lazy! |
|---|
| 241 | if not gen.searchInlist(vardims,'x') and not gen.searchInlist(vardims,'y'): |
|---|
| 242 | newvar[:] = ovar[:] |
|---|
| 243 | else: |
|---|
| 244 | ovart = ovar[:].transpose() |
|---|
| 245 | if newvar.dtype == type(float(2.)) or newvar.dtype == type(np.float(2.)) \ |
|---|
| 246 | or newvar.dtype == type(np.float32(2)) or \ |
|---|
| 247 | newvar.dtype == type(np.float64(2)): |
|---|
| 248 | newvals = Sci.module_scientific.fill3dr_2dvec(matind=matindt, \ |
|---|
| 249 | inmat=ovart, id1=ovart.shape[0], id2=ovart.shape[1], \ |
|---|
| 250 | od1=newvar.shape[2], od2=newvar.shape[1], od3=newvar.shape[0]) |
|---|
| 251 | else: |
|---|
| 252 | newvals = Sci.module_scientific.fill3di_2dvec(matind=matindt, \ |
|---|
| 253 | inmat=ovart, id1=ovart.shape[0], id2=ovart.shape[1], \ |
|---|
| 254 | od1=newvar.shape[2], od2=newvar.shape[1], od3=newvar.shape[0]) |
|---|
| 255 | newvar[:] = newvals.transpose() |
|---|
| 256 | |
|---|
| 257 | # Attributes |
|---|
| 258 | for atn in ovar.ncattrs(): |
|---|
| 259 | if atn != '_FillValue' and atn != 'units': |
|---|
| 260 | atv = ovar.getncattr(atn) |
|---|
| 261 | ncvar.set_attribute(newvar, atn, atv) |
|---|
| 262 | ncvar.set_attribute(newvar, 'coordinates', 'lon lat') |
|---|
| 263 | onewnc.sync() |
|---|
| 264 | |
|---|
| 265 | # Global attributes |
|---|
| 266 | for atn in onc.ncattrs(): |
|---|
| 267 | atv = ovar.getncattr(atn) |
|---|
| 268 | ncvar.set_attribute(onewnc, atn, atv) |
|---|
| 269 | onewnc.sync() |
|---|
| 270 | ncvar.add_global_PyNCplot(onewnc, main, fname, '0.1') |
|---|
| 271 | onc.close() |
|---|
| 272 | |
|---|
| 273 | # Reconstructing times |
|---|
| 274 | otime = onewnc.variables['time'] |
|---|
| 275 | ncvar.set_attribute(otime, 'units', 'seconds since ' + opts.year + '-01-01 00:00:00') |
|---|
| 276 | |
|---|
| 277 | onewnc.sync() |
|---|
| 278 | onewnc.close() |
|---|
| 279 | |
|---|
| 280 | print fname + ": Successful writing of file '" + ofilen + ".nc' !!" |
|---|
![(please configure the [header_logo] section in trac.ini)](/LMDZ_WRF/chrome/site/your_project_logo.png){kind=logo}
