Index: trunk/tools/compute_accum.py
===================================================================
--- trunk/tools/compute_accum.py	(revision 653)
+++ trunk/tools/compute_accum.py	(revision 653)
@@ -0,0 +1,169 @@
+import numpy as np
+from netCDF4 import Dataset as NetCDFFile
+import os
+import re
+from optparse import OptionParser
+import numpy.ma as ma
+
+# version
+version=1.1
+
+# Filling values for floats, integer and string
+fillValueF = 1.e20
+fillValueI = -99999
+fillValueS = '---'
+
+# Length of the string variables
+StringLength = 200
+
+# Size of the map for the complementary variables/maps
+Ndim2D = 100
+
+main = 'create_OBSnetcdf.py'
+errormsg = 'ERROR -- error -- ERROR -- error'
+warnmsg = 'WARNING -- warning -- WARNING -- warning'
+
+fillValue = 1.e20
+
+def searchInlist(listname, nameFind):
+    """ Function to search a value within a list
+    listname = list
+    nameFind = value to find
+    >>> searInlist(['1', '2', '3', '5'], '5')
+    True
+    """
+    for x in listname:
+      if x == nameFind:
+        return True
+    return False
+
+
+def set_attribute(ncvar, attrname, attrvalue):
+    """ Sets a value of an attribute of a netCDF variable. Removes previous attribute value if exists
+    ncvar = object netcdf variable
+    attrname = name of the attribute
+    attrvalue = value of the attribute
+    """
+    import numpy as np
+    from netCDF4 import Dataset as NetCDFFile
+
+    attvar = ncvar.ncattrs()
+    if searchInlist(attvar, attrname):
+        attr = ncvar.delncattr(attrname)
+
+    attr = ncvar.setncattr(attrname, attrvalue)
+
+    return ncvar
+
+def basicvardef(varobj, vstname, vlname, vunits):
+    """ Function to give the basic attributes to a variable
+    varobj= netCDF variable object
+    vstname= standard name of the variable
+    vlname= long name of the variable
+    vunits= units of the variable
+    """
+    attr = varobj.setncattr('standard_name', vstname)
+    attr = varobj.setncattr('long_name', vlname)
+    attr = varobj.setncattr('units', vunits)
+
+    return
+
+####### ###### ##### #### ### ## #
+
+fold='/media/ExtDiskD/DATA/obs/HyMeX/IOP15/HyMeX_rainfal/V2'
+
+lonvals = []
+latvals = []
+
+for day in range(19,23):
+    if day == 22:
+        Nhour = 1
+    else:
+        Nhour = 24
+    for hour in range(0,Nhour):
+        print day, hour
+        filen = fold + '/geo_v2_201210' + str(day) + '{:02d}'.format(hour) + '_RR1.nc'
+        onc = NetCDFFile(filen, 'r')
+        opr = onc.variables['pracc']
+        olon = onc.variables['lon']
+        olat = onc.variables['lat']
+
+        lonv = olon[:]
+        latv = olat[:]
+       
+        Npts = len(lonv)
+        if day == 19 and hour == 0:
+            Nptsorig = Npts
+            lonvals = list(lonv)
+            latvals = list(latv)
+            heights = onc.variables['height'][:]
+            praccsum = np.zeros((Npts), dtype=np.float)
+            NOTfound = np.zeros((Npts), dtype=bool)
+            nNOTfound = 0
+
+        found = np.zeros((Nptsorig), dtype=bool)
+
+        nNotfounddate = 0
+        for ip in range(Npts):
+            if searchInlist(lonvals, lonv[ip]) and searchInlist(latvals, latv[ip]):
+                il = lonvals.index(lonv[ip])
+                if lonvals[il] == lonv[ip] and latvals[il] == latv[ip]:
+                    praccsum[il] = praccsum[il] + opr[ip]
+                    found[il] = True
+#                    print lonv[ip], latv[ip], ':', lonvals[il], latvals[il]
+            else:
+#                print 'Not in the first file!'
+                nNotfounddate = nNotfounddate + 1
+                    
+        print 'On',day,'at',hour,'not found;',nNotfounddate,'stations from',Npts
+        for ip in range(Nptsorig):
+            if not found[ip] and not NOTfound[ip]: 
+                NOTfound[ip] = True
+                nNOTfound = nNOTfound + 1
+
+        onc.close()
+
+print 'From originally',Nptsorig,'there are not continuosly records from:',nNOTfound
+#finalaccum = ma.masked_array(praccsum)
+for ip in range(Nptsorig):
+    if NOTfound[ip]:
+        praccsum[ip] = fillValueF
+
+# File creation
+newnc = NetCDFFile('accum.nc', 'w')
+
+# Dimensions definition
+newdim = newnc.createDimension('Npts',Nptsorig)
+
+# Variable creation
+newvar = newnc.createVariable('lon','f8',('Npts'))
+basicvardef(newvar, 'lon', 'Longitudes', 'degrees_East')
+newvar[:] = lonvals
+set_attribute(newvar,'axis','X')
+
+newvar = newnc.createVariable('lat','f8',('Npts'))
+basicvardef(newvar, 'lat', 'Latitudes', 'degrees_North')
+newvar[:] = latvals
+set_attribute(newvar,'axis','Y')
+
+# Total accumulation
+newvar = newnc.createVariable('pracc','f',('Npts'),fill_value=fillValueF)
+basicvardef(newvar, 'pracc', 'accumulated precipitation', 'kgm-2')
+newvar[:] = praccsum
+set_attribute(newvar,'coordinates','lon lat')
+
+# Global attributes
+##
+filen = fold + '/geo_v2_2012101900_RR1.nc'
+onc = NetCDFFile(filen, 'r')
+
+attrs = onc.ncattrs()
+
+for attr in attrs:
+    attrv = onc.getncattr(attr)
+    set_attribute(newnc, attr, attrv)
+
+onc.close()
+
+newnc.sync()
+newnc.close()
