Changeset 738 in lmdz_wrf

Timestamp:

May 2, 2016, 7:35:35 PM (9 years ago)

Author:

lfita

Message:

Adding sellonlatbox', nwew version proper done of sellonlatboxold'
Adding `|S1' to 'set_attributek'

File:

• trunk/tools/nc_var_tools.py (modified) (2 diffs)

trunk/tools/nc_var_tools.py

@@ -9298 +9298 @@
       varn= ',' list of names of the variables ('all', for all variables)
     """
+    import numpy.ma as ma
     fname = 'sellonlatbox'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print sellonlatbox.__doc__
+        quit()
+
+    arguments = '[lonName],[latName],[lonSW],[latSW],[lonNE],[latNE]'
+    check_arguments(fname,values,arguments,',')
+
+    ofile = 'sellonlatbox_' + varn + '.nc'
+
+    lonn = values.split(',')[0]
+    latn = values.split(',')[1]
+
+    lonSW = np.float(values.split(',')[2])
+    latSW = np.float(values.split(',')[3])
+    lonNE = np.float(values.split(',')[4])
+    latNE = np.float(values.split(',')[5])
+
+    objfile = NetCDFFile(ncfile, 'r')
+    if  not objfile.variables.has_key(lonn):
+        print errormsg
+        print '  ' + fname + ": file '" + ncfile + "' des not have longitude " +     \
+          "variable '" + lonn + "' !!"
+        quit(-1)
+    if  not objfile.variables.has_key(latn):
+        print errormsg
+        print '  ' + fname + ": file '" + ncfile + "' des not have latitude " +      \
+          "variable '" + lonn + "' !!"
+        quit(-1)
+
+
+    lonobj = objfile.variables[lonn]
+    latobj = objfile.variables[latn]
+
+    loninf = variable_inf(lonobj)
+    latinf = variable_inf(latobj)
+
+    lonv = lonobj[:]
+    latv = latobj[:]
+
+    if varn == 'all':
+        varns = objfile.variables
+    elif varn.find(',') != -1:
+        varns = varn.split(',')
+    else:
+        varns = [varn]
+
+# Lon/lat box
+    nalon = lonv.min()
+    nalat = latv.min()
+    xalon = lonv.max()
+    xalat = latv.max()
+
+    if lonSW < nalon:
+        print errormsg
+        print '  ' + fname + ': longitude SW corner:', lonSW,' too small!!'
+        print '    min lon data:', nalon, 'increase SW corner longitude'
+        quit(-1)
+    if latSW < nalat:
+        print errormsg
+        print '  ' + fname + ': latitude SW corner:', latSW,' too small!!'
+        print '    min lat data:', nalat, 'increase SW corner latitude'
+        quit(-1)
+    if lonNE > xalon:
+        print errormsg
+        print '  ' + fname + ': longitude NE corner:', lonNE,' too large!!'
+        print '    max lon data:', xalon, 'decrease NE corner longitude'
+        quit(-1)
+    if latNE > xalat:
+        print errormsg
+        print '  ' + fname + ': latitude NE corner:', latNE,' too large!!'
+        print '    max lat data:', xalat, 'decrease NE corner latitude'
+        quit(-1)
+
+    malon = ma.masked_outside(lonv,lonSW,lonNE)
+    malat = ma.masked_outside(latv,latSW,latNE)
+
+    if len(lonv.shape) == 1:
+        dimx = lonv.shape[0]
+        dimy = latv.shape[0]
+        malonlat = np.zeros((dimy,dimx), dtype=lonv.dtype)
+        for i in range(dimx):
+            malonlat[:,i] = malon[i] + malat
+    elif len(lonv.shape) == 2:
+        dimx = lonv.shape[1]
+        dimy = lonv.shape[0]
+        malonlat = malon.mask + malat.mask
+        malonv = ma.array(lonv, mask=malonlat)
+        malatv = ma.array(latv, mask=malonlat)
+    elif len(lonv.shape) == 3:
+        print warnmsg
+        print '  ' + fname + ': assuming that lon/lat is constant to the first dimension!!'
+        print '    shape lon:', lonv.shape, 'lat:', latv.shape
+        dimx = lonv.shape[2]
+        dimy = lonv.shape[1]
+        malonlat = np.zeros((dimy,dimx), dtype=bool)
+        malonlat = malon.mask[0,:,:] + malat.mask[0,:,:]
+        malonv = ma.array(lonv, mask=malon)
+        malatv = ma.array(latv, mask=malon)
+    else:
+        print errormsg
+        print '  ' + fname + ': matrix shapes not ready!!'
+        print '    shape lon:',lonv.shape,'lat:',latv.shape
+        quit(-1)
+
+# Size of the new domain
+    nlon = malon.min()
+    nlat = malat.min()
+    xlon = malon.max()
+    xlat = malat.max()
+
+    print '  ' + fname + ': data box: SW', nlon, nlat,'NE:',xlon,xlat 
+
+    ilon = dimx
+    for j in range(dimy):
+        for i in range(dimx):
+            if malonlat[j,i] == False and i < ilon: 
+              ilon = i
+    elon = 0
+    for j in range(dimy):
+        for i in range(dimx):
+            if malonlat[j,i] == False and i > elon: 
+              elon = i
+
+    ilat = dimy
+    for j in range(dimy):
+        for i in range(dimx):
+            if malonlat[j,i] == False and j < ilat: 
+              ilat = j
+    elat = 0
+    for j in range(dimy):
+        for i in range(dimx):
+            if malonlat[j,i] == False and j > elat: 
+              elat = j
+
+    newdimx = elon - ilon + 1
+    newdimy = elat - ilat + 1
+    print '  ' + fname + ': new lon size:', ilon, ',', elon, 'new dimx:', newdimx
+    print '  ' + fname + ': new lat size:', ilat, ',', elat, 'new dimy:', newdimy
+
+    newncobj = NetCDFFile(ofile, 'w')
+
+# Creation of lon,lat related dims
+    newdim = newncobj.createDimension(lonn, newdimx)
+    newdim = newncobj.createDimension(latn, newdimy)
+
+    for vn in varns:
+        print 'Adding "' + vn + "' ..."
+        if  not objfile.variables.has_key(lonn):
+            print errormsg
+            print '  ' + fname + ": file '" + ncfile + "' des not have variable '" + \
+              vn + "' !!"
+            quit(-1)
+        if not searchInlist(objfile.variables, vn):
+            print errormsg
+            print '  ' + fname + ': variable name "' + vn + '" is not in file "' +   \
+              ncfile + '" !!!!!'
+            quit(-1)
+
+        varobj = objfile.variables[vn]
+        varinf = variable_inf(varobj)
+
+# Adding variable dimensions
+        varslice = []
+        for dimn in varinf.dimns:
+            if not searchInlist(newncobj.dimensions, dimn):
+                newncobj.createDimension(dimn, len(objfile.dimensions[dimn]))
+            if dimn == lonn:
+                varslice.append(slice(ilon,elon+1))
+            elif dimn == latn:
+                varslice.append(slice(ilat,elat+1))
+            else:
+                varslice.append(slice(0,len(objfile.dimensions[dimn])))
+
+        if varinf.FillValue is not None:
+            newvar = newncobj.createVariable(vn, nctype(varinf.dtype), varinf.dimns)
+        else:
+            newvar = newncobj.createVariable(vn, nctype(varinf.dtype), varinf.dimns, \
+              fill_value = varinf.FillValue)
+        newvar[:] = varobj[tuple(varslice)]
+
+        for atrn in varinf.attributes:
+            if atrn != '_FillValue':
+                attrv = varobj.getncattr(atrn)
+                newattr = set_attributek(newvar, atrn, attrv, type(attrv))
+
+        newncobj.sync()
+
+    for atrn in objfile.ncattrs():
+        attrv = objfile.getncattr(atrn)
+        newattr = set_attributek(newncobj, atrn, attrv, type(attrv))
+
+# global attributes
+    newncobj.setncattr('author', 'L. Fita')
+    newncobj.setncattr('institution', 'Laboratire de Meteorologie Dynamique')
+    newncobj.setncattr('university', 'Pierre Marie Curie - Jussieu')
+    newncobj.setncattr('center', 'Centre National de Recherches Scientifiques')
+    newncobj.setncattr('city', 'Paris')
+    newncobj.setncattr('country', 'France')
+    newncobj.setncattr('script', 'nc_var_tools.py')
+    newncobj.setncattr('function', 'sellonlatbox')
+    newncobj.setncattr('version', '1.0')
+
+    objfile.close()
+
+    newncobj.sync()
+    newncobj.close()
+
+    print '  ' + fname + ': successful creation of file "' + ofile + '" !!!'
+
+    return
+
+
+def sellonlatboxold(values, ncfile, varn):
+    """ Function to select a lotlan box from a data-set
+      values= [lonName],[latName],[lonSW],[latSW],[lonNE],[latNE]
+        [lonName]: name of the variable with the longitudes
+        [latName]: name of the variable with the latitudes 
+        [lonSW],[latSW]: longitude and latitude of the SW vertex of the box
+        [lonNE],[latNE]: longitude and latitude of the NE vertex of the box
+      ncfile= netCDF file
+      varn= ',' list of names of the variables ('all', for all variables)
+    """
+    fname = 'sellonlatboxold'
 
     if values == 'h':
@@ -18293 +18519 @@
 
     if vartype == type('c'):
+         ncs = 'c'
+    elif vartype == '|S1':
          ncs = 'c'
     elif vartype == type(int(1)):