source: lmdz_wrf/trunk/tools/create_stations-list_simulations.py @ 2831

# -*- coding: iso-8859-15 -*-
# L. Fita, LMD. Feburary 2015. Python script to generate two lists for simulation purposes
#  for WRF: tslits, lidarlist
#  for LMDZ: pointlocations.txt
## e.g. # create_stations-list_simulations.py -a -2,1.3,37.5,40.0:3,6.7,41.7,45.4 -m WRF -n 0.30 -s AEMET_sfcst.dat,sfcstations.dat,MeteoFrance_sfcst.dat -z lidarstations.dat,dropsondes.dat,soundingsstations.dat
## e.g. # create_stations-list_simulations.py -a -2,1.3,37.5,40.0:3,6.7,41.7,45.4 -m LMDZ -n 0.30 -s AEMET_sfcst.dat,sfcstations.dat,MeteoFrance_sfcst.dat -z lidarstations.dat,dropsondes.dat,soundingsstations.dat
import numpy as np
import os
import re
from optparse import OptionParser

main = 'create_stations-list_simulations.py'
errormsg = 'ERROR -- error -- ERROR -- error'
warnmsg = 'WARNING -- warning -- WARNING -- warning'

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 remove_NONascii(string):
    """ Function to remove that characters which are not in the standard 127 ASCII
      string= string to transform
    >>> remove_NONascii('Lluís')
    Lluis
    """
    fname = 'remove_NONascii'

    newstring = string

    RTFchar= ['á', 'é', 'í', 'ó', 'ú', 'à ', 'Ú', 'ì', 'ò', 'ù', 'â', 'ê', 'î', 'Î',  \
      'û', 'À', 'ë', 'ï', 'ö', 'ÃŒ', 'ç', 'ñ','Ê', 'œ', 'Á', 'É', 'Í', 'Ó', 'Ú', 'À', \
      'È', 'Ì', 'Ò', 'Ù', 'Â', 'Ê', 'Î', 'Ô', 'Û', 'Ä', 'Ë', 'Ï', 'Ö', 'Ü', 'Ç', 'Ñ',\
      'Æ', 'Œ', '\n', '\t']
    ASCchar= ['a', 'e', 'i', 'o', 'u', 'a', 'e', 'i', 'o', 'u', 'a', 'e', 'i', 'o',  \
      'u', 'a', 'e', 'i', 'o', 'u', 'c', 'n','ae', 'oe', 'A', 'E', 'I', 'O', 'U', 'A', \
      'E', 'I', 'O', 'U', 'A', 'E', 'I', 'O', 'U', 'A', 'E', 'I', 'O', 'U', 'C', 'N',\
      'AE', 'OE', '', ' ']

    Nchars = len(RTFchar)
    for ichar in range(Nchars):
        foundchar = string.find(RTFchar[ichar])
        if foundchar != 0:
            newstring = newstring.replace(RTFchar[ichar], ASCchar[ichar])

    return newstring

def within(pos, pt, mdist):
    """ Function to determine if a given point is closer to a given distance from a 
      series of points
      pos= series of positions
      pt= point
      mdist= distance
    >>> within([[4.0,42.0] ,[5.0,43.0]], [4.5,42,5], 0.5)
    True
    """
    fname = 'within'

    if mindist != 0.:
        if type(pos) == type([]):
            posv = np.array(pos)
        else:
            posv = pos

        for ip in range(posv.shape[0]):
            dist = np.sqrt( (posv[ip,0]-pt[0])**2. + (posv[ip,1]-pt[1])**2.)
            if (dist <= mdist):
                return True
    elif mindist == 360.:
        return True

    return False

def writting_stations(oout, fobs, ars, model, mindist):
    """ Function to write the given stations in a file
      oout= object of the file to fill up
      fobs= list of files with the observations
      ars= list of areas of interest
      model= for which model are the stations
      mindist= minimal distance between stations (in degrees)
    """
    fname = 'writting_stations'

    Nars = ars.shape[1]

    iobs = 1
    obsstn = []
    posst = []

    print fname + 'model:',model

    for filen in fobs:
        print '  ' + fname + ": opening '" + filen + "'"
        ofo = open(filen, 'r')
        for line in ofo:
            if line[0:1] != '#' and len(line) > 3:
                vals = line.split(' ')
                stlon = np.float(vals[1])
                stlat = np.float(vals[2].replace('\n',''))
                for ia in range(Nars):
                    if stlon>=ars[0,ia] and stlon<=ars[1,ia] and stlat>=ars[2,ia]    \
                      and stlat<=ars[3,ia] and not searchInlist(obsstn, vals[0]) and \
                      not within(posst, [stlon, stlat], mindist):
                          if len(vals[0]) > 24:
                              stname = vals[0][0:24]
                          else:
                              stname = vals[0]
                          if model == 'WRF':
                              oout.write('{0:24s}'.format(stname) + '  h' +          \
                                str(iobs).zfill(4) + ' ' + '{0:7.5g}'.format(stlat) +\
                                '  ' + '{0:6.4g}'.format(stlon) + '\n')
                          elif model == 'LMDZ':
                              oout.write('{0:5g}, '.format(iobs) +                   \
                                '{0:7.4g}, '.format(stlon) +                         \
                                '{0:7.4g}, '.format(stlat) + "'" + stname + "'\n")

                          obsstn.append(vals[0])
                          posst.append([stlon, stlat])
                          iobs = iobs + 1
#                          quit()
                #quit()
        ofo.close()

    print '  ' + fname + ':', iobs, 'stations added'

    return obsstn

####### ###### ##### #### ### ## #
mods = ['LMDZ', 'WRF']

parser = OptionParser()
parser.add_option("-a", "--areas", dest="areas", 
  help="':' list of areas (minLon,maxLon,minLat,maxLat)", metavar="VALUES")
parser.add_option("-m", "--model", dest="model", type='choice', choices=mods,
  help="name of the model for the files", metavar="FILES")
parser.add_option("-n", "--mindist", dest="mindist", 
  help="accepted minimal distance betwen stations (degrees)", metavar="VALUES")
parser.add_option("-s", "--sfcobs", dest="sobs", 
  help="',' list of files with surface observations (files as obsname lon lat)", 
  metavar="FILE")
parser.add_option("-z", "--heightobs", dest="hobs", 
  help="',' list of files with vertical-observations (files as obsname lon lat)", 
  metavar="FILES")
(opts, args) = parser.parse_args()

########    ####### 
## MAIN
    #######

ofile = 'newlist.dat'

if opts.areas is None:
    print errormsg
    print '  ' + main + ': no areas are provided!!'
    quit(-1)
else:
    ars = opts.areas.split(':')
    Nareas = len(ars)
    areas = np.zeros((4,Nareas), dtype = np.float)
    iar = 0
    for ar in ars:
        arvals = ar.split(',')
        for i in range(4):
            areas[i,iar] = np.float(arvals[i])
        iar = iar + 1

if opts.hobs is None:
    print errormsg
    print '  ' + main + ': no height observational files are provided!!'
    quit(-1)
else:
    hobs = opts.hobs.split(',')
    for hf in hobs:
        if not os.path.isfile(hf):
            print errormsg
            print '   ' + main + ": hieght obsevational file '" + hf +               \
              "' does not exist !!"
            quit(-1)

if opts.sobs is None:
    print errormsg
    print '  ' + main + ': no surface observational files are provided!!'
    quit(-1)
else:
    sobs = opts.sobs.split(',')
    for sf in sobs:
        if not os.path.isfile(sf):
            print errormsg
            print '   ' + main + ": surface obsevational file '" + sf +              \
              "' does not exist !!"
            quit(-1)

if opts.mindist is None:
    print warnmsg
    print '  ' + main + ': no minimal distance is provided!!'
    print '    asuming mindist=0.'
    mindist = 0.
else:
    mindist = np.float(opts.mindist)

if opts.model is None:
    print errormsg
    print '  ' + main + ': no model is provided!!'
    print '    available options:',mods
    quit(-1)
else:
    print main + ": '" + opts.model + "' model"
    if opts.model == 'WRF':
        sfilename = 'tslist'
        hfilename = 'lidarlist'
        for filen in [sfilename, hfilename]:
            oscff = open(filen, 'w')
            oscff.write('#-----------------------------------------------#\n')
            oscff.write('# 24 characters for name | pfx |  LAT  |   LON  |\n')
            oscff.write('#-----------------------------------------------#\n')
            oscff.close()

# Height observations
##
        objout = open(hfilename, 'a')
#        print main + ": creation of height observational file '" + hfilename + "' ..."
        hstations = writting_stations(objout, hobs, areas, opts.model, mindist)
        print main + "' succesfull written of height location observations '" +      \
          hfilename + "' !!"
        objout.close()

# surface observations
##
        objout = open(sfilename, 'a')
        sstations = writting_stations(objout, sobs, areas, opts.model, mindist)
        print main + "' succesfull written of surface location observations '" +     \
          sfilename + "' !!"
        objout.close()

    elif opts.model == 'LMDZ':
# Only one file as output
        shfilename = 'pointlocations.txt'
        shobs = hobs + sobs

        objout = open(shfilename, 'w')
        sstations = writting_stations(objout, shobs, areas, opts.model, mindist)
        print main + "' succesfull written of location observations '" +             \
          shfilename + "' !!"
        objout.close()
    else:
        print errormsg
        print '  ' + main + ": model '" + opts.model + "' not ready !!"
        quit(-1)