def readheader(fid,nheader) :
#Reading the header
        for j in range(nheader):
	        line = fid.readline()


import numpy as nm
import matplotlib.pyplot as plt # for plotting
#import os
import time as tm
import os

nplot = input("Choose a variables to plot \n \
 1  - -w*/u* versus mixed-layer depth \n \
 2  -  None \n ")

#path = '../results/Seida1_2007-2008/time_series/' # path to files
# Number of experiments is limited to 10
#path = [os.path.expanduser('~/Files/main/observdata/Kuivajarvi/')]
path = ['../results/Kuivajarvi2013/time_series/'] #, \
#        '../results/Kuivajarvi2013/noseiches/'] # path to files
npaths = len(path)
#nameexp = ['obs']
nameexp = ['base'] #,'noseiches']

nlevs_water = 16

# User controls
if nplot == 1: # Temperature of surface and soil
        basename = 'conv_series  1  1'
        fileout = 'w*_u*_vs_ML_ser'
        nheader = 18 #Number of lines in the header
        ncols = 18 #Total number of columns
        ncol = [15,12] # ,20,21,22,23] #The number of columns in the file to display
        legend = [''] #
        labelx = '$-w_*/u_*$, n/d'
        labely = 'Mixed-layer depth, m'
        ctime = 5
        timescale = 365/12*24 # 365/12*24 - number of hours in a month
        time0 = 4 #Initial time
        #revert_yaxis = 1
        missingval = [-1e+3,1e+3]
#if nplot == 2: # Temperature at different depths
#        basename = 'water_temp  1  1'
#        fileout  = 'temp_wat'
#        nheader  = 6
#        ncols    = 37 #Total number of columns
#        ncol     = [7,9,11,13,15,17,19,21,23,25,27,29,31,33,35,37] #The number of columns in the file to display
#        legend   = ['T_0.2m','T_0.5m','T_1m','T_1.5m','T_2m','T_2.5m', \
#        'T_3m','T_3.5m','T_4m','T_4.5m','T_5m','T_6m','T_7m','T_8m','T_10m','T_12m']
#        labelx   = 'Time, months'
#        labely   = 'Depth, m'
#        timescale = 365/12 # 365/12*24 - number of hours in a month
#        ctime = 5
#        time0 = 4 #Initial time
#if nplot == 3: # Mixed-layer depth
#        basename = 'conv_series  1  1'
#        fileout = 'Mixed_layer_ser'
#        nheader = 15 #Number of lines in the header
#        ncols = 15 #Total number of columns
#        ncol = [11,12] # ,20,21,22,23] #The number of columns in the file to display
#        legend = ['MLD from $N^2$ profile','MLD from TKE profile'] #
#        labelx = 'Time, months'
#        labely = 'Depth, m'
#        ctime = 5
#        timescale = 365/12*24 # 365/12*24 - number of hours in a month
#        time0 = 4 #Initial time
#        revert_yaxis = 1
#if nplot == 4: # u_*/w_*
#        basename = 'conv_series  1  1'
#        fileout = 'w*_u*_ser'
#        nheader = 15 #Number of lines in the header
#        ncols = 15 #Total number of columns
#        ncol = [15] # ,20,21,22,23] #The number of columns in the file to display
#        legend = ['-w*/u*'] #
#        labelx = 'Time, months'
#        labely = 'Depth, m'
#        ctime = 5
#        timescale = 365/12*24 # 365/12*24 - number of hours in a month
#        time0 = 4 #Initial time
#        #revert_yaxis = 1
#elif nplot == 5: # Methane
#        basename = 'methane_water  1  1'
#        fileout = 'Methane_water_ser'
#        nheader = 6 #Number of lines in the header
#        ncols = 39 #Total number of columns
#        ncol = [7,11,19,27,31,37] #The number of columns in the file to display
#        legend = ['CH4_0.2m','CH4_1m','CH4_3m','CH4_5m','CH4_7m','CH4_12m']
#        labelx = 'Time, months'
#        labely = 'Methane concentration, $\mu g/l$'
#        ctime = 5
#        timescale = 365/12 # 365/12*24 - number of hours in a month
#        time0 = 4 #Initial time
#elif nplot == 6: # Oxygen
#        basename = 'oxygen_water  1  1'
#        fileout = 'Oxygen_water_ser'
#        nheader = 6 #Number of lines in the header
#        ncols = 39 #Total number of columns
#        ncol = [7,11,19,27,31,37] #The number of columns in the file to display
#        legend = ['O2_0.2m','O2_1m','O2_3m','O2_5m','O2_7m','O2_12m']
#        labelx = 'Time, months'
#        labely = 'Oxygen concentration, mg/l'
#        ctime = 5
#        timescale = 365/12 # 365/12*24 - number of hours in a month
#        time0 = 4 #Initial time
#elif nplot == 7: # Carbon dioxide
#        basename = 'co2_water  1  1'
#        fileout = 'co2_water_ser'
#        nheader = 6 #Number of lines in the header
#        ncols = 39 #Total number of columns
#        ncol = [7,11,19,27,31,37] #The number of columns in the file to display
#        legend = ['CO2_0.2m','CO2_1m','CO2_3m','CO2_5m','CO2_7m','CO2_12m']
#        labelx = 'Time, months'
#        labely = 'CO2 concentration, mg/l'
#        ctime = 5
#        timescale = 365/12 # 365/12*24 - number of hours in a month
#        time0 = 4 #Initial time
#elif nplot == 8: # water surface level at souther/northern parts
#        basename = 'layers  1  1'
#        fileout = 'hS-hN_ser'
#        nheader = 15 #Number of lines in the header
#        ncols = 15 #Total number of columns
#        ncol = [7,8] #The number of columns in the file to display
#        legend = ['Mean level S','Mean level N']
#        labelx = 'Time, minutes'
#        labely = 'Deviation above mean level, m'
#        ctime = 5
#        timescale = 1./60 # 365/12*24 - number of hours in a month
#        time0 = 0 #Initial time
#elif nplot == 9: # Ebullition at the surface
#        basename = 'methane_series  1  1'
#        fileout = 'ch4_ebul_soilcols_ser'
#        nheader = 27 #Number of lines in the header
#        ncols = 30 #Total number of columns
#        ncol = [26,27,28,29,30] #The number of columns in the file to display
#        legend = ['Fb_1','Fb_2','Fb_3','Fb_4','Fb_5']
#        labelx = 'Time, months'
#        labely = 'Surface bubble flux, $mg/(m^2day)$'
#        ctime = 5
#        timescale = 365/12*24 # 365/12*24 - number of hours in a month
#        time0 = 4 #Initial time
#        #valscale = 1e+6 # converting from mol/(m^2*s) to mcmol/(m^2*s)
#elif nplot == 10: # Methane fluxes
#        basename = 'methane_series  1  1'
#        fileout = 'methane_flux_ser'
#        nheader = 27 #Number of lines in the header
#        ncols = 30 #Total number of columns
#        ncol = [14,16,17] #The number of columns in the file to display
#        legend = ['Methane ebullition flux at the surface', \
#        'Methane diffusion flux at the lake bottom','Methane diffusion flux at the lake surface']
#        labelx = 'Time, months'
#        labely = 'Methane flux, $\mu mol/(m^2*s)$'
#        ctime = 5
#        timescale = 365/12*24 # 365/12*24 - number of hours in a month
#        time0 = 4 #Initial time
#        valscale = 1e+6 # converting from mol/(m^2*s) to mcmol/(m^2*s)
#elif nplot == 11: # Methane production in sediments
#        basename = 'rprod_soil  1  1'
#        fileout = 'methane_prod_new_ser'
#        nheader = 6 #Number of lines in the header
#        ncols = 15 #Total number of columns
#        ncol = [7,9,11,13,15] #The number of columns in the file to display
#        legend = ['Methane production in soil col 1', \
#                  'Methane production in soil col 2', \
#                  'Methane production in soil col 3', \
#                  'Methane production in soil col 4', \
#                  'Methane production in soil col 5' ]
#        labelx = 'Time, months'
#        labely = 'Methane production rate, $mol/(m^2s)$'
#        ctime = 5
#        timescale = 365/12 # 365/12*24 - number of hours in a month
#        time0 = 4 #Initial time
#if nplot == 12: # Temperature measured
#        basename = 'KJdata2013'  #'obs'
#        fileout = 'Temp_obs'
#        nheader = 0 #Number of lines in the header
#        ncols = 21 #Total number of columns
#        ncol = [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17] # ,20,21,22,23] #The number of columns in the file to display
#        legend = ['T_0.2m','T_0.5m','T_1m','T_1.5m','T_2m','T_2.5m', \
#         'T_3m','T_3.5m','T_4m','T_4.5m','T_5m','T_6m','T_7m','T_8m', \
#         'T_10m','T_12m'] 
#        labelx = 'Time, months'
#        labely = 'Temperature, Celsius'
#        ctime = 1
#        timescale = 365/12 # 365/12*24 - number of hours in a month
#        time0 = 0 #Initial time
#if nplot == 13: # Carbon dioxide measured
#        basename = 'KJdata2013'
#        fileout = 'co2_obs'
#        nheader = 1 #Number of lines in the header
#        ncols = 21 #Total number of columns
#        ncol = [18,19,20,21]  #The number of columns in the file to display
#        legend = ['CO2_0.2m','CO2_1.5m','CO2_2.5m','CO2_7m'] 
#        labelx = 'Time, months'
#        labely = 'Carbon dioxide, ppm'
#        ctime = 1
#        timescale = 365/12 # 365/12*24 - number of hours in a month
#        time0 = 0 #Initial time
#if nplot == 14: # CO_2 measured manually
#        basename = 'ch4&co2man2013'
#        fileout = 'co2_obs_man'
#        nheader = 1 #Number of lines in the header
#        ncols = 28 #Total number of columns
#        ncol = [6,9,12,15,18,21,24,27]  #The number of columns in the file to display
#        legend = ['CO2_0.2m','CO2_1m','CO2_3m','CO2_5m','CO2_7m', \
#                  'CO2_9m','CO2_11m','CO2_12m'] 
#        labelx = 'Time, months'
#        labely = '$CO_2$ concentration, $mg/l$'
#        ctime = 3
#        timescale = 365/12 # 365/12*24 - number of hours in a month
#        time0 = 0 #Initial time
#if nplot == 15: # Methane measured manually
#        basename = 'ch4&co2man2013'
#        fileout = 'ch4_obs_man'
#        nheader = 1 #Number of lines in the header
#        ncols = 28 #Total number of columns
#        ncol = [5,8,11,14,17,20,23,26]  #The number of columns in the file to display
#        legend = ['CH4_0.2m','CH4_1m','CH4_3m','CH4_5m','CH4_7m', \
#                  'CH4_9m','CH4_11m','CH4_12m'] 
#        labelx = 'Time, months'
#        labely = 'Methane concentration, $\mu g/l$'
#        ctime = 3
#        timescale = 365/12 # 365/12*24 - number of hours in a month
#        time0 = 0 #Initial time
#if nplot == 16: # Oxygen measured manually
#        basename = 'o2man2013'
#        fileout = 'o2_obs'
#        nheader = 1 #Number of lines in the header
#        ncols = 24 #Total number of columns
#        ncol = [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24]  #The number of columns in the file to display
#        legend = ['CH4_0m','CH4_0.5m','CH4_1m','CH4_1.5m','CH4_2m', \
#                  'CH4_2.5m','CH4_3m','CH4_3.5m','CH4_4m','CH4_4.5m','CH4_5m', \
#                  'CH4_5.5m','CH4_6m','CH4_6.5m','CH4_7m','CH4_7.5m','CH4_8m', \
#                  'CH4_8.5m','CH4_9m','CH4_10m','CH4_11m','CH4_12m','CH4_12.5m'] 
#        labelx = 'Time, months'
#        labely = 'Oxygen concentration, $mg/l$'
#        ctime = 1
#        timescale = 365/12 # 365/12*24 - number of hours in a month
#        time0 = 0 #Initial time
#if nplot == 17: # Methane flux measured
#        basename = 's13_Kuic_all'
#        fileout =  'ch4flux_obs'
#        nheader = 1 #Number of lines in the header
#        ncols = 41 #Total number of columns
#        ncol = [24]  #The number of columns in the file to display
#        legend = ['CH4_flux_surf'] 
#        labelx = 'Time, months'
#        labely = 'Surface methane flux, $\mu mol/(m^2*s)$'
#        ctime = 1
#        timescale = 365/12 # 365/12*24 - number of hours in a month
#        time0 = 0 #Initial time
#        missingval = -999.
#if nplot == 18: # Carbon dioxide flux measured
#        basename = 's13_Kuic_all'
#        fileout =  'co2flux_obs'
#        nheader = 1 #Number of lines in the header
#        ncols = 41 #Total number of columns
#        ncol = [5]  #The number of columns in the file to display
#        legend = ['CO2_flux_surf'] 
#        labelx = 'Time, months'
#        labely = 'Surface CO_2 flux, mol/(m^2*s)$'
#        ctime = 1
#        timescale = 365/12 # 365/12*24 - number of hours in a month
#        time0 = 0 #Initial time
#        missingval = -999.
#if nplot == 19: # Richardson number
#        basename = 'Ri  1  1'
#        fileout = 'Ri_timser'
#        nheader = 6 #Number of lines in the header
#        ncols = 39
#        ncol = [11,19,27,31,37]
#        legend = ['Ri_1m','Ri_3m','Ri_5m','Ri_7m','Ri_12m']
#        labelx = 'Time, months'
#        labely = 'Ri'
#        timescale = 365/12 # 365/12*24 - number of hours in a month
#        ctime = 5
#        time0 = 4 #Initial time
#if nplot == 20: # Heat conductance
#        basename = 'lamw  1  1'
#        fileout = 'lamw_timser'
#        nheader = 6 #Number of lines in the header
#        ncols = 39
#        ncol = [11,19,27,31,37]
#        legend = ['HC_1m','HC_3m','HC_5m','HC_7m','HC_12m']
#        labelx = 'Time, months'
#        labely = 'Heat conductance, $m^2/s$'
#        timescale = 365/12 # 365/12*24 - number of hours in a month
#        ctime = 5
#        time0 = 4 #Initial time

if 'legend' in locals():
        nleg = len(legend)
        legendbase = []*nleg
        legendbase[:] = legend[:]

colorline = ['b','g','r','c','m','y','k']
typeline = ['-','--','-.',':','.',',','o','v','^','1']
ncolors = len(colorline)

# Changing default matplotlib settings
plt.rc('xtick',labelsize=16)
plt.rc('ytick',labelsize=16)
plt.rc('lines',linewidth=2)

ncol = nm.array(ncol,dtype='int32')
ncol[:] = ncol[:] - 1
nvardisp = len(ncol)

#nfiles = 10 #number of files
#time = [0]*nfiles
##Files names suffixes (timestamps)
#time[0] = '2003010100'
#time[1] = '2003010200'
#time[2] = '2003010300'
#time[3] = '2003010400'
#time[4] = '2003010500'
#time[5] = '2003010600'
#time[6] = '2003010700'
#time[7] = '2003010800'
#time[8] = '2003010900'
#time[9] = '2003011000'

#for i in range(nfiles):
for ip in range(npaths):
        fid = open(path[ip] + basename + '.dat','r')
        readheader(fid,nheader)

        #Counting the number of raws in a file (assuming the same in all files)
        array = fid.read()
        nraws = array.count('\n')
        print('The number of values is '+str(nraws))
        time = nm.zeros(nraws)
        vals = nm.zeros((nraws,nvardisp)) #[[0]*nvardisp]*nraws
        fid.seek(0) #rewinding file
        readheader(fid,nheader)

        #Reading array of depths and values to be plotted
        for j in range(nraws) :
	        line = fid.readline()
                line = line.split()
	        time[j] = float(line[ctime-1])
                for i in range(nvardisp):
                        vals[j,i] = float(line[ncol[i]])

        time[:] = time[:]/timescale + time0
        for j in range(nvardisp) :
               #vard = [0]*nraws
               #for i in range(nraws):
               #        vard[i] = vals[i,j]
               # Assigning missing values to NaN
               if 'missingval' in locals():
                        for i in range(nraws):
                                for k in range(len(missingval)):
                                        if vals[i,j] == float(missingval[k]):
                                                vals[i,j] = 'NaN'
                                                break
               if 'valscale' in locals():
                        vals[:,j] = vals[:,j]*valscale
        jcolor = ip-(ip//ncolors)*ncolors # periodic use of colors in one experiment
        #if 'plotid' in locals():
                # Adding to existent plot new curves
                #plt.plot(vals[:,1],vals[:,2],colorline[jcolor]+typeline[ip],figure=plotid)
        #else:
                # Creating plot 
        plotid = plt.scatter(vals[:,0],vals[:,1]) #,colorline[jcolor]+typeline[ip])
        
        # Adding to legend new entries from other experiments
        if ip > 0 and 'legend' in locals():
               for i in range(nleg):
                       legend.append(legendbase[i]+'_'+nameexp[ip])
        else:
               for i in range(nleg):
                       legend[i] = legend[i]+'_'+nameexp[ip]



plt.xlabel(labelx,fontsize=20)
plt.ylabel(labely,fontsize=20)
if 'legend' in locals():
        plt.legend(legend,fontsize=20,loc='best')

#Revert y axis if needed
if 'revert_yaxis' in locals():
        plt.gca().invert_yaxis()

plt.savefig(path[0] + fileout + '_' + nameexp[0] + '.png')
plt.show()
