source: anuga_work/development/ted_rigby/get_tuflow_data.py @ 6823

"""
--------------------------------------------------------------get_tuflow_data.py
ANUGA UTILITIES TO READ TUFLOW MID/MIF INPUT FILES
These functions assist Anuga to read in data from a previous tuflow model

            get_tuflow_bc            Reads in the inflow and dsbc timeseries
            get_wbnm_coords          Gets the outlet coords of each inflow subarea
            get_tuflow_bc_filetype   Gets the tuflow (ts1 or csv) bc filetype
            convert_asc_ts_to_tms    Converts an ascii timeseries to NetCDF format
            get_tuflow_tmf           Reads in the  surface classes and properties
            get_tuflow_mat           Reads in the surface class polygons
            get_tuflow_2d_po         Reads the tuflow po flowlines 
            
            

This script is imported by model_run.py

Code Version : 1.00 June 2008 Initial  release
Author       : E Rigby (0437 250 500) ted.rigby@rienco.com.au

-------------------------------------------------------------------------------
"""

import string

#START FUNCTION ================================================= get_tuflow_bc

def get_tuflow_bc(location,loc_filename,bc_filetype,verbose):
    'returns an ascii hydrograph for that location'
               
    if verbose :
        print "         ++++ get_tuflow_bc - Opening file %s to obtain temporal data for %s type bc at %s " %(loc_filename,bc_filetype,location) 

    if bc_filetype=="ts1":
    
        #   Is a ts1 florate file - open and read in the event and location specific flow hydrograph
        bctype = "inflow"
        
        if verbose :
            print "         ++++ get_tuflow_bc - bc type is  %s" %(bctype)
        ts1_fid=open(loc_filename, 'r') # open the associated hydrograph file for that loc
        ts1_lines=ts1_fid.readlines()   # read all file into string list
        ts1_fid.close()                 # close the loc specific hydrograph file
        subsline=ts1_lines[4]           # select the line with subarea headings

        if verbose :
            print "         ++++ get_tuflow_bc - Subarea header line is ", subsline [:50]
        index=subsline.index(location)+len(location) -12  # locate start pos of selected subarea header col
        tsecs=[]
        qmecs=[]
        for ts1_line in ts1_lines[5:]:
            t=float(ts1_line[:12])*60.0          # time convert time to secs - was minutes
            if (ts1_line[index:index+12]) == "            ":
                q=0                              # end of this hydrograph not all full duration
            else:
                q=float(ts1_line[index:index+12])# flowrate in m3/sec
            tsecs.append(t)
            qmecs.append(q)
        if verbose :
            print "         ++++ get_tuflow_bc - ts1 flow hydrograph at %s (first 5 values)" %(location)
            for i in range (5):
                print "         ++++ ", tsecs[i],qmecs[i]
                
        return tsecs,qmecs

    else:

        #    Is a dsbc csv file - open and read in the event and location specific elevation hydrograph
        bctype = "dsbc"

        if verbose :
            print "         ++++ bc type is  %s" %(bctype)
        csv_fid=open(loc_filename, 'r')  # open the associated elev hydrograph file for that loc
        csv_lines=csv_fid.readlines()    # read all file into string list
        csv_fid.close()                  # close the loc specific hydrograph file
        tsecs=[]
        elev =[]
        for csv_line in csv_lines[1:]:   # read the dsbc data
            sep_col=csv_line.find(',')   # locate seperator - not fixed field widths
            t=float(csv_line[:sep_col])* 3600.0 # time convert to seconds for Anuga - was hours
            e=float(csv_line[sep_col+1:])# elevation at bdry in (m)
            tsecs.append(t)
            elev.append(e)
        if verbose :
            for i in range (5):
                print "         ++++ csv elev hydrograph at ", location,tsecs[i],elev[i]

        return tsecs,elev

#END FUNCTION =================================================== get_tuflow_bc       



#START FUNCTION =============================================== get_wbnm_coords

def get_wbnm_coords(wbnm_filename,subarea,verbose):
    'reurns the x and y outlet coordinates of the subarea'
    
    wbn_fid=open(wbnm_filename, 'r')           # open the wbnm run file
    for wbn_line in  wbn_fid:                  # read file a line at a time
        
        if wbn_line[:12].find(subarea.strip())>=0:     # found the subarea line
            xcoord = float(wbn_line[36:47])    # Subareas outlet xcoord
            ycoord = float(wbn_line[48:59])    # Subareas outlet xcoord
            wbn_fid.close()                    # close the wbn run file
            if verbose :
                print '         ++++ get_wbnm_coords - Subareal outlet for %s found at %7.2f %7.2f ' %(subarea,xcoord,ycoord)   
            return xcoord, ycoord              # return the outlet coords   
                        
        if wbn_line.find("#####END_TOPOLOGY")>0:   # did not find the subarea in topology
            print "        >>>> ERROR - get_wbnm_coords - Could not find %s in %s" %(subarea,wbnm_filename),
            sys.exit(0)                        # terminate the program
    


#END FUNCTION =================================================== get_wbnm_coords


#START FUNCTION =========================================== get_tuflow_bc_filetype

def get_tuflow_bc_filetype(loc_filename,verbose):
    'returns type ("ts1") ("csv") of bc_file'

    if loc_filename.find(".ts1")>=0:    # Is file containing hydrographs for each subarea
        bc_filetype="ts1"
        return bc_filetype
        
    elif loc_filename.find(".csv")>=0:  # Is time(hrs),elev(m) dsbc.csv file
        bc_filetype="csv"
        return bc_filetype
        
    else:
        print "         >>>> ERROR get_tuflow_bc_filetype - Unrecognised bc filetype %s - aborting " %(bc_filetype)
        sys.exit(0)                    # terminate the program



#END FUNCTION =========================================== get_tuflow_bc_filetype

#START FUNCTION ========================================== convert_asc_ts_to_tms
def convert_asc_ts_to_tms(tms_dir,ts_name, time_units, time, value,verbose):
    'saves the [time][value] timeseries into a netcdf tms object '

    import os
    from os import sep
    from Numeric import array, Float
    from Scientific.IO.NetCDF import NetCDFFile
      
    if verbose :
        print "         ++++ asc_ts_to_tms - Converting %s time series in %s to tms " %(ts_name,time_units)
        for i in range (5):     # len(time) would be full list
            print "         ++++ ", time[i],value[i]
    
    tms_filename = tms_dir + sep + ts_name.strip() + '.tms' # buildname of tms file to be created
    if verbose :
        print "         ++++ asc_ts_to_tms - Created tms file named %s to store time series in" %(tms_filename)
 
    
# compute time conversion factor to bring to secs

    if time_units=='days':
        factor=60.0*60.0*24.0
    elif time_units=='hours':
        factor=60.0*60.0
    elif time_units=='minutes':
        factor=60.0
    else:
        factor=1.0
    for i in range (len(time)):
        time[i]=time[i]*factor
       
    # convert to and create NetCDF tms file of time series
    

    T=array(time,Float)          # secs
    V=array(value,Float)         # whatever

    tms_fid=NetCDFFile(tms_filename,'w')
    tms_fid.institution='Rienco Consulting'
    tms_fid.description='Utility to create NetCDF (tms) file of ascii time series'
    tms_fid.starttime=0.0
    tms_fid.createDimension('number_of_timesteps', len(T))
    tms_fid.createVariable('time',Float,('number_of_timesteps',) )
    tms_fid.variables['time'][:] = T
    tms_fid.createVariable('value',Float,('number_of_timesteps',) )
    tms_fid.variables['value'][:]= V
    tms_fid.close()
    
    if verbose :
            print "         ++++ asc_ts_to_tms - NetCDF time series created and stored in %s - time units are %s " %(tms_filename,time_units)
        
    return tms_filename
   

#END FUNCTION ========================================== convert_asc_ts_to_tms

#START FUNCTION =============================================== get_tuflow_tmf
def get_tuflow_tmf(tmf_filename,verbose):
    'read in the tuflow roughness values and return as python list of values'

    
    tmf_fid=open(tmf_filename, 'r')           # open the tuflow material (roughness) file
    tmf_lines=tmf_fid.readlines()             # read all file into string list
    tmf_fid.close()                           # close the tmf file

    mat_roughness_data_list=[]

    if verbose :
        print "         ++++ get_tuflow_tmf - Surface roughness data read in from tmf and stored in mat_roughness_data_list " 

    for tmf_line in tmf_lines:
        #if tmf_line[0]=='!' or tmf_line[0]==' ' or tmf_line[0]=='\n':
        if tmf_line[0]=='!' or tmf_line[0] in string.whitespace:
            pass
        else:                                                 # decode and store the list values
            tmf_line=tmf_line.replace('!',',')                # get rid of the ! so can use split
            tmf_line_fields = tmf_line.split(',')             # break the line into fields
            tmf_line_fields[0] = int(tmf_line_fields[0])      # The tuflow surface RID
            tmf_line_fields[1] = float(tmf_line_fields[1])    # The tuflow fixed n value
            tmf_line_fields[2] = int(tmf_line_fields[2])      # Not used
            tmf_line_fields[3] = int(tmf_line_fields[3])      # Not used
            tmf_line_fields[4] = float(tmf_line_fields[4])    # The tuflow d1 value
            tmf_line_fields[5] = float(tmf_line_fields[5])    # The tuflow n1 value
            tmf_line_fields[6] = float(tmf_line_fields[6])    # The tuflow d2 value
            tmf_line_fields[7] = float(tmf_line_fields[7])    # The tuflow n2 value
            tmf_line_fields[8] = str(tmf_line_fields[8])      # The tuflow surface description
                                        
            if verbose :
                print "         ++++ get_tuflow_tmf  appending  "  ,tmf_line_fields

            mat_roughness_data_list.append(tmf_line_fields)   # there are nine data fields including the description
                        
    
                
    return mat_roughness_data_list            # a list of surface data fields as in the tuflow tmf file

#END FUNCTION ================================================ get_tuflow_tmf


#START FUNCTION =============================================== get_tuflow_mat
def get_tuflow_mat(mat_filename,verbose):
    'reads in the roughness polys from the tuflow 2d_mat.mif file and IDs from the mid'

    # Read in the polygons from the 2d_mat.mif file
    matmif_fid=open(mat_filename, 'r')       # open the tuflow material (roughness polys) file
    matmif_lines=matmif_fid.readlines()      # read all file into string list
    matmif_fid.close()                       # close the 2d_mat.mif file

    
    poly_coords_list=[]                               # float list of roughness polygon coord lists
    point_coords=[]                                   # float list of single coord x y pair
    
    for i in range (len(matmif_lines)):
        
        if matmif_lines[i][:12].find('Region')>=0: 
            nlines=int(matmif_lines[i+1] )           # read next line get no of coords (nlines) in poly
            poly=[]                                  # float list of coords x,y
            
            for j in range (nlines):                 # read next nlines of xy coords into poly[]
                x,y = matmif_lines[i+2+j].split()    # split into two strings
                x = float(x)                         # convert to numeric coords
                y = float(y)
                point_coords = [x,y]
                poly.append(point_coords)            # add to poly as float xy pairs

            poly_coords_list.append(poly)            # append this poly coord list to list of polys
            i=i+1+j                                  # reset counter to previous to last line read to avoid problem when at end of file

        else:
            # unwanted - skip line
            pass
     

    # Read in the polygon material (roughness) IDs from the 2d_mat.mid file
    matmid_fid=open(mat_filename[:-4] + '.mid', 'r') # open the tuflow material (roughness polys) file
    matmid_lines=matmid_fid.readlines()              # read all file into string list
    matmid_fid.close()                               # close the 2d_mat.mid file

    poly_RID_list=[]                                 # integer list identifying roughness polygons ID
    for matmid_line in matmid_lines:
        RID=int(matmid_line[:(matmid_line.index(','))])  # get the RID of each poly
        poly_RID_list.append(RID)

    if len(poly_RID_list)== len(poly_coords_list):   # all as it should be
        if verbose :
            print "         ++++ get_tuflow_mat - Read in %i surface roughness polygons from %s " %(len(poly_RID_list), mat_filename)
        return poly_RID_list,poly_coords_list        # [intIDs],[list of poly floatcoords list]
    else:
        print "         >>>> ERROR get_tuflow_mat - Poly regions list and poly ID list not of same length in get_tuflow_mat??"
        print "         >>>> ERROR get_tuflow_mat - Poly ID list %d and poly list %d" %(len(mat_ID_list),len(mat_poly_list))
        sys.exit(0)                                  # terminate the program

#END FUNCTION ================================================ get_tuflow_mat