source: lmdz_wrf/trunk/tools/compute_accum.py @ 1333

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()