Changeset 2168 in lmdz_wrf for trunk/tools

Timestamp:

Oct 5, 2018, 10:19:51 PM (7 years ago)

Author:

lfita

Message:

Adding:

• `retrieve_stations': Function to retrieve temporal values at given stations provided by a secondary netcdf

Location:

trunk/tools

Files:

• nc_var.py (modified) (4 diffs)
• nc_var_tools.py (modified) (9 diffs)

trunk/tools/nc_var.py

@@ -61 +61 @@
 ## e.g. # nc_var.py -o merge_files -S 'plev|plev,time|time:merged.nc' -f 'ncobs/AliceSprings/snd_94326_198201-198201.nc,ncobs/AliceSprings/snd_94326_198202-198202.nc' -v 'plev,time'
 ## e.g. # nc_var.py -o temporal_stats -S 'Time:WRFtime:day@1@min,LTday@-3@1@min:bottom_top@ZNU,south_north@XLAT,west_east@XLONG' -f wrfout_d01_1995-01-01_00:00:00 -v T2,Q2
+## e.g. # nc_var.py -o retrieve_stations -f wrfout_d01_1995-01-01_00:00:00 -S 'tmin_percentiles.nc:stname:None:stlon:stlat:None:nearest:west_east:XLONG:south_north:XLAT:HGT:Time:WRFtime' -v T2,QVAPOR
 
 from optparse import OptionParser
@@ -158 +159 @@
 # pinterp: Function to vertically interpolate using subroutines from the p_interp.F90 NCAR program
 # remapnn: Function to remap to the nearest neightbor a variable using projection from another file
+# retrieve_stations: Function to retrieve temporal values at given stations provided by a secondary netcdf
 # seasmean: Function to compute the seasonal mean of a variable
 # sellonlatbox: Function to select a lotlan box from a data-set
@@ -219 +221 @@
   'netcdf_fold_concatenation_HMT', 'reproject', 'Partialmap_Entiremap',              \
   'Partialmap_EntiremapFor', 'Partialmap_EntiremapForExact',                         \
-  'pinterp', 'reconstruct_matrix_from_vector', 'remapnn', 'rm_FillValue',            \
+  'pinterp', 'reconstruct_matrix_from_vector', 'remapnn', 'retrieve_stations',       \
+  'rm_FillValue',                                                                    \
   'seasmean', 'sellonlatbox', 'sellonlatlevbox', 'selvar', 'setvar_asciivalues',     \
   'setvar_nc', 'sorttimesmat', 'spacemean', 'SpatialWeightedMean',                   \
@@ -470 +473 @@
 elif oper == 'reproject':
     ncvar.reproject(opts.values, opts.ncfile, opts.varname)
+elif oper == 'retrieve_stations':
+    ncvar.retrieve_stations(opts.values, opts.ncfile, opts.varname)
 elif oper == 'rm_FillValue':
     ncvar.rm_FillValue(opts.values, opts.ncfile, opts.varname)

trunk/tools/nc_var_tools.py

@@ -112 +112 @@
 # maskvar: Function to mask a variable using a mask. It is assumed that mask[...,dimM,dimJ,dimK] and var[...,dimM,dimJ,dimK] share the last dimensions
 # merge_files: Function to merge variables from two files
-# mindist: Function to proivde the minimum distance toa pair of lon,lat
+# mindist: Function to provide the minimum distance toa pair of lon,lat
 # ModelChar: Class object to determine major characterisitcs of a given model output
 # model_characteristics: Functino to provide major characterisitcs of a given model output
@@ -135 +135 @@
 # remapnn_old: Function to remap to the nearest neightbor a variable using projection from another fil
 # reproject: Function to reproject values to another one
+# retrieve_stations: Function to retrieve temporal values at given stations provided by a secondary netcdf
 # rm_FillValue: Operation to remove the _FillValue from a variable inside a netCDF file
 # seasmean: Function to compute the seasonal mean of a variable
@@ -25381 +25382 @@
     stats = []
     for timest in timestatsS:
+        print 'Lluis timest:', timest, 'split:', timest.split('@')
         timestatsv = []
         period = timest.split('@')[0]
@@ -25512 +25514 @@
 
         # Adding period dimension and variables
-        newdim = onewnc.createDimension(periodn+'_time', Ntslc)
-        newvart = onewnc.createVariable(periodn+'_time', 'f8',                       \
-          (periodn+'_time'))
-        basicvardef(newvart, periodn+'_time', 'time for ' + periodn, timeu)
-        newvart.setncattr('calendar', timec)
-        newvart.setncattr('bounds', periodn+'_time_bounds')
-        newvarb = onewnc.createVariable(periodn+'_time_bounds', 'f8',                \
-          (periodn+'_time', 'bounds'))
-        basicvardef(newvarb, periodn+'_time_bounds', 'time boundaries for '+ periodn,\
-          timeu)
-        newvarb.setncattr('calendar', timec)
-        for it in range(Ntslc):
-            timeslcev = tslcs[it]
-            islc = timevc[timeslcev[0]]
-            eslc = timevc[timeslcev[1]]
-
-            newvart[it] = (eslc + islc) / 2.
-            newvarb[it,:] = [islc, eslc]
-        onewnc.sync()
+        if not gen.searchInlist(onewnc.dimensions, periodn+'_time'):
+            newdim = onewnc.createDimension(periodn+'_time', Ntslc)
+            newvart = onewnc.createVariable(periodn+'_time', 'f8',                   \
+              (periodn+'_time'))
+            basicvardef(newvart, periodn+'_time', 'time for ' + periodn, timeu)
+            newvart.setncattr('calendar', timec)
+            newvart.setncattr('bounds', periodn+'_time_bounds')
+            newvarb = onewnc.createVariable(periodn+'_time_bounds', 'f8',            \
+              (periodn+'_time', 'bounds'))
+            basicvardef(newvarb, periodn+'_time_bounds', 'time boundaries for '+     \
+              periodn, timeu)
+            newvarb.setncattr('calendar', timec)
+            for it in range(Ntslc):
+                timeslcev = tslcs[it]
+                islc = timevc[timeslcev[0]]
+                eslc = timevc[timeslcev[1]]
+    
+                newvart[it] = (eslc + islc) / 2.
+                newvarb[it,:] = [islc, eslc]
+            onewnc.sync()
 
         for varn in varnames:
@@ -25576 +25579 @@
                 newvarn = varn+str(amount)+periodn+statn
 
-            newvar = onewnc.createVariable(newvarn, 'f', tuple(dimntstat),          \
+            newvarv = onewnc.createVariable(newvarn, 'f', tuple(dimntstat),         \
               fill_value=gen.fillValueR)
             varunits = get_varunits(ovar)
@@ -25584 +25587 @@
               CFvarattr[4].replace('|',' ')
             add_varattrs(onc,onewnc,[varn],[newvarn])
-            basicvardef(newvar, newvarn, Lname, varunits)
-            newvar.setncattr('cell_method', 'time: ' + Lstatn[statn])
-            newvar.setncattr('bounds', periodn+'_time_bounds')
-            if periodn == 'LTday': set_attributek(newvar, 'UTCshift', UTCshift, 'I')
-
-            tstat = np.ones(tuple(shapetstat), dtype=np.float)*gen.fillValueR
+            basicvardef(newvarv, newvarn, Lname, varunits)
+            newvarv.setncattr('cell_method', 'time: ' + Lstatn[statn])
+            newvarv.setncattr('bounds', periodn+'_time_bounds')
+            if periodn == 'LTday': set_attributek(newvarv, 'UTCshift', UTCshift, 'I')
+
+            newvarN = onewnc.createVariable(newvarn+'_Nsteps', 'i', (periodn+'_time'))
+            basicvardef(newvarN, newvarn+'_Nsteps', 'Number of time steps used to ' +\
+              'compute ' + Lname, '1')
+
+            #tstat = np.ones(tuple(shapetstat), dtype=np.float)*gen.fillValueR
 
             # Looping into the period
             for islc in range(Ntslc):
                 timeslcv = tslcs[islc]
+                newvarN[islc] = timeslcv[1]-timeslcv[0]+1
                 
                 timeslc = []
@@ -25611 +25619 @@
                 tvals = ovar[tuple(timeslc)]
                 if statn == 'min':
-                    tstat[tuple(tstatslc)] = np.min(tvals, axis=itdim)
+                    newvarv[tuple(tstatslc)] = np.min(tvals, axis=itdim)
                 elif statn == 'max':
-                    tstat[tuple(tstatslc)] = np.max(tvals, axis=itdim)
+                    newvarv[tuple(tstatslc)] = np.max(tvals, axis=itdim)
                 elif statn == 'mean':
-                    tstat[tuple(tstatslc)] = np.mean(tvals, axis=itdim)
+                    newvarv[tuple(tstatslc)] = np.mean(tvals, axis=itdim)
                 elif statn == 'mean2':
-                    tstat[tuple(tstatslc)] = np.mean(tvals**2, axis=itdim)
+                    newvarv[tuple(tstatslc)] = np.mean(tvals**2, axis=itdim)
                 elif statn == 'std':
-                    tstat[tuple(tstatslc)] = np.std(tvals, axis=itdim)
+                    newvarv[tuple(tstatslc)] = np.std(tvals, axis=itdim)
                 else:
                     print errormsg
@@ -25627 +25635 @@
                     print '    available ones:', statnor
                     quit(-1)
-
-            newvar[:] = tstat
+                onewnc.sync()
+
+            #newvar[:] = tstat
             onewnc.sync()
 
@@ -25646 +25655 @@
 #temporal_stats(vals, '/home/lluis/PY/wrfout_d01_1995-01-01_00:00:00', 'T2')
 
-#quit()
+def retrieve_stations(values, ncfile, variable):
+    """ Function to retrieve temporal values at given stations provided by a 
+          secondary netcdf
+      values=[stnetcdf]:[stvnamen]:[stvidn]:[stvlonn]:[stvlatn]:
+          [stvheight]:[method]:[datadlonn]:[datavlonn]:[datadlatn]:[datavlatn]:
+          [datavheightn]:[datadtimen]:[datavtimen]
+        [stnetcdf]: name of the file with the stations
+        [stvnamen]: name of the variable in [stnetcdf] with the names of the stations
+          (None for inexstent)
+        [stvidn]: name of the variable in [stnetcdf] with the ids of the stations 
+           (None for inexstent)
+        [stvlonn]: name of the variable in [stnetcdf] with the longitudes of the 
+           stations (None for inexstent)
+        [stvlatn]: name of the variable in [stnetcdf] with the latitudes of the 
+          stations (None for inexstent)
+        [stvheightn]: name of the variable in [stnetcdf] with the heights of the 
+          stations (None for inexstent)
+        [method]: methodology to use to retrieve the values
+          'nearest': nearest grid point
+          'mean': mean of the closest 4 grid points
+          'dist': inverse distance weighted mean of the closest 4 grid points
+        [datadlonn]: name of dimension in datafile for the lon axis
+        [datavlonn]: name of variable in datafile with the longitudes
+        [datadlatn]: name of dimension in datafile for the lat axis
+        [datavlatn]: name of variable in datafile with the latitudes
+        [datavheightn]: name of variable in datafile with the heighs (None for 
+          inexistent)
+        [datadtimen]: name of dimension in datafile for the time axis
+        [datavtimen]: name of variable in datafile for the time values ('WRFtime' 
+          for WRF times)
+      ncfile= file from which retrieve values of stations
+      variable= ',' list of values to retrieve values ('all' for all variables)
+    """
+    fname = 'retrieve_stations'
+    availmethod = ['nearest', 'mean', 'dist']
+    methoddesc = {'nearest': 'nearest grid point', 'mean': 'mean of the closest 4 '+ \
+      'grid points', 'dist': 'inverse distance weighted mean of the closest 4 ' +    \
+      'grid points'}
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print retrieve_stations.__doc__
+        quit()
+
+    expectargs = '[stnetcdf]:[stvnamen]:[stvidn]:[stvlonn]:[stvlatn]:[stvheightn]:'+ \
+      '[method]:[datadlonn]:[datavlonn]:[datadlatn]:[datavlatn]:[datavheightn]:' +   \
+      '[datadtimen]:[datavtimen]'
+    gen.check_arguments(fname, values, expectargs, ':')
+
+    stnetcdf = values.split(':')[0]
+    stvnamen = values.split(':')[1]
+    stvidn = values.split(':')[2]
+    stvlonn = values.split(':')[3]
+    stvlatn = values.split(':')[4]
+    stvheightn = values.split(':')[5]
+    method = values.split(':')[6]
+    datadlonn = values.split(':')[7]
+    datavlonn = values.split(':')[8]
+    datadlatn = values.split(':')[9]
+    datavlatn = values.split(':')[10]
+    datavheightn = values.split(':')[11]
+    datadtimen = values.split(':')[12]
+    datavtimen = values.split(':')[13]
+
+    # Retrieving station information
+
+    # Getting variables
+    filestinf = {}
+    stinf = []
+    if stvnamen != 'None':
+        filestinf['station_name'] = stvnamen
+        stinf.append('station_name')
+    else:
+        filestinf['station_name'] = None
+    if stvidn != 'None':
+        filestinf['station_id'] = stvidn
+        stinf.append('station_id')
+    else:
+        filestinf['station_id'] = None
+    if stvlonn != 'None':
+        filestinf['station_lon'] = stvlonn
+        stinf.append('station_lon')
+    else:
+        print errormsg
+        print '  ' + fname + ": at least a variable with longitudes of the " +       \
+          "stations is required !!"
+        quit(-1)
+    if stvlatn != 'None':
+        filestinf['station_lat'] = stvlatn
+        stinf.append('station_lat')
+    else:
+        print errormsg
+        print '  ' + fname + ": at least a variable with latitudes of the " +        \
+          "stations is required !!"
+        quit(-1)
+    if stvheightn != 'None':
+        filestinf['station_height'] = stvheightn
+        stinf.append('station_height')
+    else:
+        filestinf['station_height'] = None
+
+    # Getting values
+    if not os.path.isfile(stnetcdf):
+        print errormsg
+        print '  ' + fname + ": file with stations information  '" + stnetcdf +      \
+          "' does not exist !!"
+        quit(-1)
+    ostnc = NetCDFFile(stnetcdf, 'r')
+
+    stationsinf = {}
+    for sti in stinf:
+        varn = filestinf[sti]
+        if not ostnc.variables.has_key(varn):
+            print errormsg
+            print '  ' + fname + ": stations does not have variable '" + varn +      \
+              "' for '" + sti + "' !!"
+            varns = ostnc.variables.keys()
+            varns.sort()
+            print '    available ones:', varns
+            quit(-1)
+
+        ovar = ostnc.variables[varn]
+        if sti == 'station_name':
+            stvals = get_str_nc(ovar, ovar.shape[1])
+            stnames = []
+            for ist in range(len(stvals)):
+                stnames.append(''.join(stvals[ist]))
+            stnames.remove('')
+            stationsinf[sti] = stnames
+        else:
+            stationsinf[sti] = ovar[:]
+
+    #print infmsg
+    #print '    ' + fname  + ': stations _______'
+    #gen.printing_dictionary(stationsinf)
+
+    # Getting stations points
+    donc = NetCDFFile(ncfile, 'r')
+
+    if not donc.variables.has_key(datavlonn):
+        print errormsg
+        print '  ' + fname + ": data file does not have '" + datavlonn + "' for " +  \
+          "longitudes !!"
+        varns = donc.variables.keys()
+        varns.sort()
+        print '    available ones:', varns
+        quit(-1)
+
+    if not donc.variables.has_key(datavlatn):
+        print errormsg
+        print '  ' + fname + ": data file does not have '" + datavlatn + "' for " +  \
+          "latitudes !!"
+        varns = donc.variables.keys()
+        varns.sort()
+        print '    available ones:', varns
+        quit(-1)
+
+    if datavheightn != 'None' and not donc.variables.has_key(datavheightn):
+        print errormsg
+        print '  ' + fname + ": data file does not have '" + datavheightn + "' for "+\
+          "Heights !!"
+        varns = donc.variables.keys()
+        varns.sort()
+        print '    available ones:', varns
+        quit(-1)  
+
+    if datavtimen != 'WRFtime' and not donc.variables.has_key(datavtimen):
+        print errormsg
+        print '  ' + fname + ": data file does not have '" + datavtimen + "' for " + \
+          "times !!"
+        varns = donc.variables.keys()
+        varns.sort()
+        print '    available ones:', varns
+        quit(-1)  
+
+    olon = donc.variables[datavlonn]
+    lon0 = olon[:]
+    olat = donc.variables[datavlatn]
+    lat0 = olat[:]
+    oheight = donc.variables[datavheightn]
+    height0 = oheight[:]
+
+    lon, lat = gen.lonlat2D(lon0, lat0)
+    ddimx = lon.shape[1]
+    ddimy = lon.shape[0]
+
+    if datavheightn != 'None':
+        if len(height0.shape) == 3: height = height0[0,:,:]
+        elif len(height0.shape) == 2: height = height0[:]
+        else:
+            print errormsg
+            print '  ' + fname + ": height variable has not 2D, 3D shape:", height0.shape
+            quit(-1)
+    else:
+        height = None
+
+    if datavtimen != 'WRFtime':
+        otime = donc.variables[datavtimen]
+        timev = otime[:]
+        timeu = time.units
+        if gen.searchInlist(otime.ncattrs(), 'calendar'):
+            timec = time.calendar
+        else:
+            print warnsmg
+            print '    ' + fname + ": variable time without calendar attribute !"
+            print "    assuming 'standard/greogorian'"
+            timec = 'standard'
+    else:
+        otimev = donc.variables['Times']
+        timewrfv = otimev[:]
+        timev, timeu = compute_WRFtime(timewrfv, refdate='19491201000000',           \
+          tunitsval='minutes')
+        timec = 'standard'
+
+    if variable == 'all':
+        varns = donc.variables.keys()
+        # Removing that variables without lon,lat dimensoins
+        for vn in varns:
+            ovn = donc.variables[vn]
+            vdims = ovn.dimensions
+            if not gen.searchInlist(vdims, datadlonn) or not                         \
+              gen.searchInlist(vdims, datadlonn): varns.remove(vn)
+    else:
+        varns = gen.str_list(variable, ',')
+
+    stns = stationsinf['station_name']
+    stlons = stationsinf['station_lon']
+    stlats = stationsinf['station_lat']
+
+    Nstations = len(stns)
+    print infmsg
+    print '    ' + fname + ': to process', Nstations, 'stations !!'
+
+    # File creation
+    ## 
+    ofile = fname + '.nc'
+    onewnc = NetCDFFile(ofile, 'w')
+
+    # Creation of dimensions
+    newdim = onewnc.createDimension('Lstring',256)
+    newdim = onewnc.createDimension(datadtimen,None)
+    newdim = onewnc.createDimension('station',Nstations)
+    newdim = onewnc.createDimension('time',None)
+
+    # Create variable-dimensions
+    newvar = onewnc.createVariable('station_name','c', ('station', 'Lstring'))
+    writing_str_nc(newvar, stns, 256)
+    basicvardef(newvar, 'station_name', 'Name of the station', '-')    
+
+    newvar = onewnc.createVariable('time','f8', ('time'))
+    newvar[:] = timev
+    basicvardef(newvar, 'time', 'Time', timeu)
+    newvar.setncattr('calendar', timec)
+
+    newvari = onewnc.createVariable('station_i','i', ('station'))
+    basicvardef(newvari, 'station_i', 'x-axis index of the station in the grid', '-')
+
+    newvarj = onewnc.createVariable('station_j','i', ('station'))
+    basicvardef(newvarj, 'station_j', 'y-axis index of the station in the grid', '-')
+
+    newvarlon = onewnc.createVariable('station_lon','f', ('station'))
+    basicvardef(newvarlon, 'station_lon', 'Longitude of the station in the grid',    \
+      'degrees_east')
+
+    newvarlat = onewnc.createVariable('station_lat','f', ('station'))
+    basicvardef(newvarlat, 'station_lat', 'Latitude of the station in the grid',     \
+      'degrees_north')
+
+    newvardist = onewnc.createVariable('station_dist','f', ('station'))
+    basicvardef(newvardist, 'station_dist', 'Distance of the station in the grid',   \
+      'degrees')
+
+    if datavheightn != 'None':
+        newvarhgt = onewnc.createVariable('station_hight','f', ('station'))
+        basicvardef(newvarhgt, 'station_height', 'Height of the station in the grid',\
+          'm')
+
+    if method == 'nearest':
+        newdim = onewnc.createDimension('Nweight', 1)
+    else:
+        newdim = onewnc.createDimension('Nweight', 4)
+    newvarwght = onewnc.createVariable('station_weights','f',('station','Nweight'))
+    basicvardef(newvarwght, 'station_weights', 'Weights to retrieve station ' +      \
+      'values from grid', '1')
+    newvarwght.setncattr('method',method)
+
+    stdataij = {}
+    stslice = {}
+    stweight = {}
+
+    if method == 'dist': stdist = {}
+
+    pointdatavals = datavlonn+':'+datavlatn+':'+datadlonn+'|-1,'+datadlatn+'|-1,' +  \
+      datadtimen+'|0'
+    for ist in range(Nstations):
+        
+        ijst, dst = get_point(pointdatavals, ncfile, str(stlons[ist])+','+           \
+          str(stlats[ist]))
+        stdataij[stns[ist]] = [ijst, dst]
+        dlon = lon[ijst[0],ijst[1]]
+        dlat = lat[ijst[0],ijst[1]]
+
+        # Writing in file station values from grid
+        newvari[ist] = ijst[1]
+        newvarj[ist] = ijst[0]
+        newvarlon[ist] = dlon
+        newvarlat[ist] = dlat
+        newvardist[ist] = dst
+        if height is not None: newvarhgt[ist] = height[ijst[0],ijst[1]]
+        onewnc.sync()
+
+        ptvslice = []
+        if method == 'nearest':
+            ptvslice = ijst + []
+            weights = np.ones((1), dtype=np.float)
+        elif method == 'mean' or method == 'dist':
+            ddx = lon[ijst[0],ijst[1]+1] - dlon
+            ddy = lat[ijst[0]+1,ijst[1]] - dlat
+            signx = np.abs(ddx)/ddx
+            signy = np.abs(ddy)/ddy
+
+            if stlons[ist] >= dlon: incx = 1*signx
+            else: incx = -1*signx
+            if stlats[ist] >= dlat: incy = 1*signy
+            else: incy = -1*signy
+ 
+            # Avoiding outside shape
+            xij = np.min([ijst[1]+incx, ddimx-1])
+            yij = np.min([ijst[0]+incy, ddimy-1])
+            xij = int(np.max([0, xij]))
+            yij = int(np.max([0, yij]))
+
+            ptvslice.append(ijst)
+            ptvslice.append([ijst[0],xij])
+            ptvslice.append([yij,ijst[1]])
+            ptvslice.append([yij,xij])
+
+            weights = np.ones((4), dtype=np.float)
+            if method == 'dist':
+                dists = np.zeros((4), dtype=np.float)
+                dists[0] = dst
+                for ineig in range(1,4):
+                    ij = ptvslice[ineig]
+                    distv = (lon[ij[1],ij[0]]-stlons[ist])**2 +                      \
+                      (lat[ij[1],ij[0]]-stlats[ist])**2
+                    dists[ineig] = np.sqrt(distv)
+                weights = 1./dists
+                stdist[stns[ist]] = dists
+                weights = weights/np.sum(weights)
+        else:
+            print errormsg
+            print '  ' + fname + ": method '" + + "' not ready !!"
+            print '    available ones:', availmethod
+            quit(-1)
+        stslice[stns[ist]] = ptvslice
+        stweight[stns[ist]] = weights
+        newvarwght[ist,:] = weights
+
+    for ist in range(Nstations):
+        ptvslice = stslice[stns[ist]]
+
+        if type(ptvslice[0]) == type(np.int64(1)): Nstpts = 1
+        else: Nstpts = len(ptvslice)
+        wghts = stweight[stns[ist]]
+
+        for varn in varns:
+            if not donc.variables.has_key(varn):
+                print errormsg
+                print '  ' + fname + ": data file does not have variable '" + varn + \
+                  "' !!"
+                fvarns = donc.variables.keys()
+                fvarns.sort()
+                print '    available ones:', varns
+                quit(-1)
+
+            if ist == 0:
+                print varn + " ..."
+
+            ovarn = donc.variables[varn]
+
+            vslice = []
+            idd = 0
+            ilon = -1
+            ilat = -1
+            vstshape = [Nstpts]
+            vstslice = [Nstpts]
+            vardimns = ['station']
+            stgridvals = [Nstations]
+            for dn in ovarn.dimensions:
+                if dn == datadlonn: 
+                    vslice.append(0)
+                    ilon = idd
+                elif dn == datadlatn:
+                    vslice.append(0)
+                    ilat = idd
+                else: 
+                    dshape = ovarn.shape[idd]
+                    vslice.append(slice(0,dshape))
+                    vstshape.append(dshape)
+                    vstslice.append(slice(0,dshape))
+                    vardimns.append(ovarn.dimensions[idd])
+                    if ist == 0 and not gen.searchInlist(onewnc.dimensions, dn):
+                        newdim=onewnc.createDimension(dn, dshape)
+                    stgridvals.append(slice(0,dshape))
+                idd = idd + 1
+
+            if ist == 0: 
+                varattrs = ovarn.ncattrs()
+                if gen.searchInlist(varattrs, '_fillValue'):
+                    fillval = ovarn.getncattr('_fillValue')
+                    newvar= onewnc.createVariable(varn, ovarn.dtype, tuple(vardimns),\
+                      fill_value=fillval)
+                    varatttrs.remove('_fillValue')
+                else:
+                    newvar= onewnc.createVariable(varn, ovarn.dtype, tuple(vardimns))
+                for attrn in varattrs:
+                    attrv = ovarn.getncattr(attrn)
+                    newvar.setncattr(attrn,attrv)
+
+                onewnc.sync()
+            else:
+                newvar = onewnc.variables[varn]
+
+            varvalues = np.zeros((vstshape), dtype=ovar.dtype)
+            vwghts = np.ones((vstshape), dtype=np.float)
+            if Nstpts > 2:
+                for istpts in range(Nstpts):
+                    vstslice[0] = istpts
+                    istp = ptvslice[istpts]
+                    vslice[ilat] = istp[0]
+                    vslice[ilon] = istp[1]
+
+                    varvalues[tuple(vstslice)] = ovarn[tuple(vslice)]
+                    vwghts[tuple(vstslice)] = wghts[istpts]
+            else:
+                vslice[ilat] = ptvslice[0]
+                vslice[ilon] = ptvslice[1]
+                varvalues[0] = ovarn[tuple(vslice)]
+
+                # Shapping weights
+                vwghts = np.ones(varvalues.shape, dtype=np.float)
+
+            stvarv = np.sum(varvalues*vwghts, axis=0)
+
+            stgridvals[0] = ist
+            if ist == 1:
+                print 'Lluis shapes ____'
+                print 'newvar:', newvar.shape
+                print 'vslice:', vslice
+                print 'vstslice:', vstslice
+                print 'stgridvals:', stgridvals
+                print 'varvalues:', varvalues.shape
+                print 'weights:', vwghts.shape
+                print 'stvarv:', stvarv.shape
+                vvv =  varvalues*vwghts
+                print ' Lluis shape prod:', vvv.shape
+
+            newvar[tuple(stgridvals)] = stvarv[:]
+            onewnc.sync()
+
+            #quit()
+    
+    # Adding PyNCplot
+    add_global_PyNCplot(onewnc, 'nc_var_tools', fname, '1.0')
+    onewnc.setncattr('description', 'Values at different stations from a grid file')
+    onewnc.setncattr('stations_file', stnetcdf)
+    onewnc.setncattr('grid_file', ncfile)
+    onewnc.setncattr('method', method)
+    onewnc.setncattr('method_description', methoddesc[method])
+    onewnc.sync()
+    
+    add_globattrs(ostnc, onewnc, 'all')
+    onewnc.sync()
+
+    ostnc.close()
+    onewnc.close()
+    print fname + ": successfull creation of file '" + ofile + "' !!"
+
+    return
+
+#values = '/home/lluis/estudios/WRFsensSFC/dyn-stat/stat/tmin_percentiles.nc:stname:None:stlon:stlat:None:dist:west_east:XLONG:south_north:XLAT:HGT:Time:WRFtime'
+
+#retrieve_stations(values, '/home/lluis/PY/wrfout_d01_1995-01-01_00:00:00', 'T2')
+
+#quit()