Changeset 574 in lmdz_wrf for trunk/tools

Timestamp:

Jul 6, 2015, 5:32:42 PM (10 years ago)

Author:

lfita

Message:

Adding `getvalues_lonlat' to get values to the closest longitude, latitude point

File:

• trunk/tools/nc_var_tools.py (modified) (1 diff)

trunk/tools/nc_var_tools.py

@@ -16402 +16402 @@
     return newstr
 
-#quit()
+def getvalues_lonlat(values, ncfile):
+    """ Function to retrieve the values from the closest grid point to a set of longitude, latitude values
+      values=[longitude]:[latitude]:[dnames]:[vdnames]
+        [longitude]: value for the longitude
+        [latitude]: value for the latitude
+        [dnames]: [xdimname],[ydimname] names of the dimensions on the 'x' and 'y' axis
+        [dvnames]: [xvdimname],[yvdimname] names of the variables with the values for
+          the dimensions on the 'x' and 'y' axis
+      ncfile= netcdf file to use
+    """
+    import subprocess as sub
+    fname = 'getvalues_lonlat'
+
+    ofile = fname + '.nc'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print getvalues_lonlat.__doc__
+        quit()
+
+    arguments = '[longitude]:[latitude]:[dnames]:[vdnames]'
+    check_arguments(fname, len(values.split(':')), arguments, ':',                   \
+      len(arguments.split(',')))
+
+    longitude = np.float(values.split(':')[0])
+    latitude = np.float(values.split(':')[1])
+    dnames = values.split(':')[2].split(',')
+    dvnames = values.split(':')[3].split(',')
+
+    onc = NetCDFFile(ncfile, 'r')
+
+    for dn in dnames:
+        if not searchInlist(onc.dimensions, dn):
+            print errormsg
+            print '  ' + fname + ": file '" + ncfile + "' does not have dimension '"+\
+              dn + "' !!"
+            quit(-1) 
+
+    for vdn in dvnames:
+        if not searchInlist(onc.variables, vdn):
+            print errormsg
+            print '  ' + fname + ": file '" + ncfile + "' does not have variable '"+ \
+              vdn + "' !!"
+            quit(-1) 
+
+# Getting grid point localization
+##
+    lon2d, lat2d = lonlat2D(onc.variables[dvnames[0]], onc.variables[dvnames[1]])
+    
+    dimx = lon2d.shape[1]
+    dimy = lon2d.shape[0]
+
+    mappairs = np.ones((2), dtype=int)
+    dist=np.sqrt((longitude - lon2d)**2. + (latitude - lat2d)**2.)
+    mindist = np.min(dist)
+    minydist, minxdist = index_mat(dist,mindist)
+    mappairs[:] = [minydist, minxdist]
+
+    print '  ' + fname + ': nearest point: ',mappairs, 'distance:', mindist
+
+    onc.close()
+
+    vals = dnames[1] + ',' + str(mappairs[1]) + ',' + str(mappairs[1]) + ',1@'
+    vals = vals + dnames[0] + ',' + str(mappairs[0]) + ',' + str(mappairs[0]) + ',1'
+
+    DataSetSection_multidims(vals, ncfile)
+    of = ncfile.split('.')[0] + '_' +  dnames[1] + '_B' + str(mappairs[1]) + '-E' +  \
+      str(mappairs[1]) + '-I1_' + dnames[0] + '_B' + str(mappairs[0]) + '-E' +       \
+      str(mappairs[0]) + '-I1.nc'
+
+    sub.call(['mv',of,ofile])
+
+    onc = NetCDFFile(ofile, 'a')
+
+    newdim = onc.createDimension('lLd',1)
+    newvar = onc.createVariable('lon_point','f',('lLd'))
+    newattr = basicvardef(newvar,'lon_point','longitude of selection of the values', \
+      'degrees_East')
+    newvar[:] = longitude
+
+    newvar = onc.createVariable('lat_point','f',('lLd'))
+    newattr = basicvardef(newvar,'lat_point','latitude of selection of the values',  \
+      'degrees_North')
+    newvar[:] = latitude
+
+    newvar = onc.createVariable('flon_point','f',('lLd'))
+    newattr = basicvardef(newvar,'flon_point','closest longitude of selection of ' + \
+      'the values', 'degrees_East')
+    newvar[:] = lon2d[tuple(mappairs)]
+
+    newvar = onc.createVariable('flat_point','f',('lLd'))
+    newattr = basicvardef(newvar,'flat_point','closest latitude of selection of ' +  \
+      'the values', 'degrees_North')
+    newvar[:] = lat2d[tuple(mappairs)]
+
+    newvar = onc.createVariable('fdist_point','f',('lLd'))
+    newattr = basicvardef(newvar,'fdist_point','distance to the point of selection '+\
+      'of the values', 'degrees')
+    newvar[:] = mindist
+
+    onc.sync()
+    onc.close()
+
+    print fname + ": successful writing of file '" + ofile + "' !!"
+
+    return 
+
+vals='20.5:36.0:west_east,south_north:XLONG_M,XLAT_M'
+ncf='/home/lluis/etudes/domains/medic051215_2km/geo_em.d01.nc'
+
+getvalues_lonlat(vals, ncf)
+quit()
 
 """operation to make: