content: # -*- coding: utf-8 -*-
  """
  Created on Fri Jul 25 08:48:28 2014
  
  @author: david
  """
  
  #*************** IMPORT DEPENDANCIES*******************************************
  import numpy as np
  #import spec_gdal4 as spg
  from osgeo import gdal
  import os
  import csv
  #import h5py
  import datetime
  import numpy.ma as ma
  #from StringIO import StringIO
  #import shapely
  #import r2py
  
  from osgeo import gdal_array
  from osgeo import gdalconst
  from osgeo.gdalconst import *
  
  
  from osgeo import ogr
  from osgeo import osr
  
  
  from scipy.spatial import ConvexHull
  from scipy.signal import find_peaks_cwt
  from scipy.signal import savgol_filter
  from scipy import interpolate
  
  import matplotlib.pyplot as plt
  
  #from shapely.geometry import LineString
  
  ################# Functions ###################################################
  '''These here are functions that are not part of any specific class- these 
  are used by the data import classes for functions such as smoothing'''
  
  def smoothing(perc_out, block_start, block_end, kparam, weight, sparam):   
      #D
      sm_spline_block = perc_out[block_start:block_end,:]
      sm_x = sm_spline_block[:,0]
      sm_y = sm_spline_block[:,1]
      sm_len = sm_x.shape
        
      sm_weights = np.zeros(sm_len)+weight
     
      
      sm_spline = interpolate.UnivariateSpline(sm_x, 
                                               sm_y, 
                                               k=kparam,
                                               w=sm_weights,
                                               s=sparam)
                                               
      
      spline = sm_spline(sm_x)
      
      
      spline = np.column_stack((sm_x,spline))
      
      
      
      return spline
      
      
  def interpolate_gaps(array1, array2):          
     array_end = array1.shape[0]-1
     array1_endx = array1[array_end, 0]
     #get the start point of the second array
     array2_start = array2[0,0] 
     #get the length of the area to be interpolated
     x_len = array2_start-array1_endx+1
     #generate x values to use for the array
     xvals = np.linspace(array1_endx, array2_start, num=x_len)
     #y val for the start of  the interpolated area
     yval_array1 = array1[array_end,1]
     # y val for the end of interpolated area    
     yval_array2 = array2[0,1]
     #stack the values into a new array
     xin = np.append(array1_endx, array2_start)
     yin = np.append(yval_array1, yval_array2) 
     
     #numpy.interp(x, xp, fp)
     gap_filling = np.interp(xvals, xin, yin)
     filled_x = np.column_stack((xvals, gap_filling))
     print (filled_x.shape)
     return filled_x
  
  class absorption_feature():
      '''this class is used for the characterisation of spectral absortion features, 
      and their investigation using continuum removal'''
      def __init__(self, spectra, feat_start, feat_end, feat_centre): 
          self.wl = spectra[:,0]
          self.values = spectra[:,1]
          print ('CALL TO ABSORPTION FEATURE')
          # start of absorption feature
          self.feat_start = feat_start
          # end of absorption feature
          self.feat_end = feat_end
          # approximate 'centre' of feature
          self.feat_centre = feat_centre        
          
          #get the range of the data
          self.min_wl = self.wl[0]
          self.max_wl = self.wl[-1]
          
          print ('Absorption feature',self.feat_start,self.feat_end)
  
          #define feature name
          self.feat_name = str(self.feat_start)+'_'+str(self.feat_end)
           
          '''# if the feature is within the range of the sensor, do stuff
          if self.feat_start > self.min_wl and self.feat_end < self.max_wl:
              print 'can do stuff with this data'
              try:
                  self.abs_feature()
                  print ('Absorption feature analysis sussceful')
              except:
                  print ('ERROR analysing absorption feature', self.feat_name)
                  pass
          else:
              print ('Cannot define feature: Out of range')'''                  
      
      ########## Methods ##################################################
      def abs_feature(self):
          print ('Call to abs_feature made')
          # Meffod to calculate the end points of the absorption feature
          # Does this using the Qhull algorithim form scipy spatial
          
          #use the initial defintnion of the absorption feature as a staring point
          # get the indices for these
          
          cont_rem_stacked = None
          ft_def_stacked = None
          
          
          
          start_point = np.argmin(np.abs(self.wl-self.feat_start))
          end_point = np.argmin(np.abs(self.wl-self.feat_end))
          centre = np.argmin(np.abs(self.wl-self.feat_centre))
              
          #find the index minima of reflectance 
          minima = np.argmin(self.values[start_point:end_point])+start_point
          
          # if the minima = the start point then the start point is the minima 
          if minima == start_point:
              left = minima
          #if not then the left side of the feature is the maixima on the left of the minima 
          elif minima <= centre:
              left = start_point+np.argmax(self.values[start_point:centre])    
          else:
              left = start_point+np.argmax(self.values[start_point:minima])
              
          #right is the maxima on the right of the absorption feature
          if minima == end_point:
              right = minima     
          else:
              right = minima+np.argmax(self.values[minima:end_point])
          
          # use left and right to create a 2D array of points
          hull_in = np.column_stack((self.wl[left:right],self.values[left:right]))
          
          #determine the minima of the points
          hull_min = minima-left
          
          if hull_min <= 0:
              hull_min=0
          
          #find the wavelength at minima
          hull_min_wl = hull_in[hull_min,0]
          
          # define the wavelength ranges we'll use to select simplices
          ft_left_wl = hull_min_wl-((hull_min_wl-hull_in[0,0])/2)
          ft_right_wl = hull_min_wl+((hull_in[-1,0]-hull_min_wl)/2)
          
          #use scipy.spatial convex hull to determine the convex hull of the points
          hull = ConvexHull(hull_in)
          
          # get the simplex tuples from the convex hull
          simplexes = hull.simplices
          # create an empty list to store simplices potentially related to our feature
          feat_pos = []
          #iterate through the simplices
          for simplex in simplexes:
              #extract vertices from simplices
              vertex1 = simplex[0]
              vertex2 = simplex[1]
              
              #print 'VERT!',hull_in[vertex1,0],hull_in[vertex2,0]
              
              ''' We're only interested in the upper hull. Qhull moves counter-
              clockwise. Therefore we're only interested in those points where 
              vertex 1 is greater than vertex 2'''
              '''The above may be total bollocks'''
              
              if not vertex1 < vertex2:
                  '''We then use the wavelength ranges to determine which simplices
                  relate to our absorption feature'''
                  if hull_in[vertex2,0] <= ft_left_wl and \
                     hull_in[vertex2,0] >= self.wl[left] and \
                     hull_in[vertex1,0] >= ft_right_wl and \
                     hull_in[vertex1,0] <= self.wl[right]:
                      # append the vertices to the list
                      print (hull_in[vertex2,0])
                      print (hull_in[vertex1,0])
                      feat_pos.append((vertex2,vertex1))
                      print ('feat_pos length:',len(feat_pos), type(feat_pos))
                      #print feat_pos[0],feat_pos[1]
                  
                  else:
                      continue
          
          '''We only want one feature here. If there's more than one or less
          than one we're not interested as we're probably not dealing with
          vegetation'''
          # If there's less than one feature...
          if len(feat_pos) < 1:
              print ('Absorption feature cannot be defined:less than one feature')
              ft_def_stacked = None 
              ft_def_hdr = None 
              cont_rem_stacked = None
  
              
              
              
              
          elif len(feat_pos) == 1:
              feat_pos=feat_pos[0]
              print ('£££££',feat_pos, type(feat_pos))
          
          else:
              #if theres more than one fid the widest one. this is not optimal.
              if len(feat_pos) >1:
                  feat_width = []
                  
                  for pair in feat_pos:
                      feat_width.append(pair[1]-pair[0])
                      print ('feat width:', feat_width)
                  #feat_width = np.asarray(feat_width)
                  print (feat_width)
                  f_max = feat_width.index(max(feat_width))
                  print (f_max)
                  feat_pos = feat_pos[f_max]
                  print (type(feat_pos))
                  
                  
         
          
          if not feat_pos==None:
              feat_pos = feat_pos[0], feat_pos[1]
              print ('DOES MY FEAT_POS CONVERSION WORK?', feat_pos)
                             
              print ('Analysing absorption feature')
              #slice
              feature = hull_in[feat_pos[0]:feat_pos[1],:]
              print ('Feature shape',feature.shape,'start:',feature[0,0],'end:',feature[-1,0])
              #get the minima in the slice
              minima_pos = np.argmin(feature[:,1])
              #continuum removal
              contrem = self.continuum_removal(feature,minima_pos)
      
              # set up single value outputs
              # start of feature
              refined_start = feature[0,0]
              # end of feature
              refined_end = feature[-1,0]
              # wavelength at minima
              minima_WL = feature[minima_pos,0]
              # reflectance at minima
              minima_R = feature[minima_pos,1]
              # area of absorption feature
              feat_area = contrem[4]
              # two band normalised index of minima and start of feature
              left_tbvi = (refined_start-minima_R)/(refined_start+minima_R)
              # two band normalised index of minima and right of feature
              right_tbvi = (refined_end-minima_R)/(refined_end+minima_R)
              # gradient of the continuum line
              cont_gradient = np.mean(np.gradient(contrem[0]))
              # area of continuum removed absorption feature
              cont_rem_area = contrem[3]
              # maxima of continuum removed absorption feature
              cont_rem_maxima = np.max(contrem[1])
              # wavelength of maxima of continuum removed absorption feature
              cont_rem_maxima_wl = feature[np.argmax(contrem[1]),0]
              #area of left part of continuum removed feature
              cont_area_l = contrem[5]
              if cont_area_l == None:
                  cont_area_l=0
              #are aof right part of continuum removed feature
              cont_area_r = contrem[6]
  
              #stack these into a lovely array
              ft_def_stacked = np.column_stack((refined_start,
                                                refined_end,
                                                minima_WL,
                                                minima_R,
                                                feat_area,
                                                left_tbvi,
                                                right_tbvi,
                                                cont_gradient,
                                                cont_rem_area,
                                                cont_rem_maxima,
                                                cont_rem_maxima_wl,
                                                cont_area_l,
                                                cont_area_r))
                                                
              ft_def_hdr = str('"Refined start",'+
                                 '"Refined end",'+
                                 '"Minima Wavelenght",'+
                                 '"Minima Reflectance",'+
                                 '"Feature Area",'+
                                 '"Left TBVI",'+
                                 '"Right TBVI",'+
                                 '"Continuum Gradient",'+
                                 '"Continuum Removed Area",'+
                                 '"Continuum Removed Maxima",'+
                                 '"Continuum Removed Maxima WL",'+
                                 '"Continuum Removed Area Left",'+
                                 '"Continuum Removed Area Right",')
                                                
              #print ft_def_stacked.shape #save the stacked outputs as hdf
                  
                  
                                                    
              # stack the 2d continuum removed outputs
              cont_rem_stacked = np.column_stack((feature[:,0],
                                                  feature[:,1],
                                                  contrem[0],
                                                  contrem[1],
                                                  contrem[2]))
                                                  
              print ('CREM', cont_rem_stacked.shape)
                                              
          return ft_def_stacked, ft_def_hdr, cont_rem_stacked
              
              
              
      def continuum_removal(self,feature,minima):  
          #method to perform continuum r=<emoval
          #pull out endmenmbers
          end_memb = np.vstack((feature[0,:],feature[-1,:]))
          #interpolate between the endmembers using x intervals    
          continuum_line = np.interp(feature[:,0], end_memb[:,0], end_memb[:,1])
          #continuum removal
          continuum_removed = continuum_line/feature[:,1]
  
          #stack into coord pairs so we can measure the area of the feature
          ft_coords = np.vstack((feature,
                                 np.column_stack((feature[:,0],continuum_line))))
          #get the area
          area = self.area(ft_coords)
                  
          #get the area of the continuum removed feature
          cont_rem_2d = np.column_stack((feature[:,0],continuum_removed))        
          cont_r_area = self.area(cont_rem_2d)
  
          #band-normalised by area continuum removal
          cont_BNA = (1-(feature[:,1]/continuum_line))/area
          
          #continuum removed area on left of minima
          cont_area_left = self.area(cont_rem_2d[0:minima,:])
  
          #continuum removed area on right of minima
          cont_area_right = self.area(cont_rem_2d[minima:,:])
          
          return (continuum_line, 
                  continuum_removed, 
                  cont_BNA, 
                  cont_r_area,  
                  area,
                  cont_area_left,
                  cont_area_right)
           
      #define area of 2d polygon- using shoelace formula
      def area(self, coords2d):
          #setup counter
          total = 0.0
          #get the number of coorsinate pairs
          N = coords2d.shape[0]
          #iterate through these
          for i in range(N):
              #define the first coordinate pair
              vertex1 = coords2d[i]
              #do the second
              vertex2 = coords2d[(i+1) % N]
              #append the first & second distance to the toatal
              total += vertex1[0]*vertex2[1] - vertex1[1]*vertex2[0]
              #return area
              return abs(total/2)  
      
  class Indices():
      #class that does vegetation indices
      def __init__(self,spectra):
          self.wl = spectra[:,0]
          self.values = spectra[:,1]
          self.range = (np.min(self.wl),np.max(self.wl))
          '''So, the init method here checks the range of the sensor and runs
          the appropriate indices within that range, and saves them as hdf5. 
          The indices are all defined as methods of this class'''
      def visnir(self):
      # Sensor range VIS-NIR
          if self.range[0] >= 350 and \
             self.range[0] <= 500 and \
             self.range[1] >= 900:
                 vis_nir = np.column_stack((self.sr700_800(),
                                            self.ndvi694_760(),
                                            self.ndvi695_805(),
                                            self.ndvi700_800(),
                                            self.ndvi705_750(),
                                            self.rdvi(),
                                            self.savi(),
                                            self.msavi2(),
                                            self.msr(),
                                            self.msrvi(),
                                            self.mdvi(),
                                            self.tvi(),
                                            self.mtvi(),
                                            self.mtvi2(),
                                            self.vog1vi(),
                                            self.vog2(),
                                            self.prsi(),
                                            self.privi(),
                                            self.sipi(),
                                            self.mcari(),
                                            self.mcari1(),
                                            self.mcari2(),
                                            self.npci(),
                                            self.npqi(),
                                            self.cri1(),
                                            self.cri2(),
                                            self.ari1(),
                                            self.ari2(),
                                            self.wbi()))
                                            
                                            
                 vis_nir_hdr=str('"sr700_800",'+
                                 '"ndvi694_760",'+
                                 '"ndvi695_805",'+
                                 '"ndvi700_800",'+
                                 '"ndvi705_750",'+
                                 '"rdvi",'+
                                 '"savi",'+
                                 '"msavi2",'+
                                 '"msr",'+
                                 '"msrvi",'+
                                 '"mdvi",'+
                                 '"tvi",'+
                                 '"mtvi",'+
                                 '"mtvi2",'+
                                 '"vog1vi",'+
                                 '"vog2",'+
                                 '"prsi"'+
                                 '"privi",'+
                                 '"sipi",'+
                                 '"mcari",'+
                                 '"mcari1",'+
                                 '"mcari2",'+
                                 '"npci",'+
                                 '"npqi",'+
                                 '"cri1",'+
                                 '"cri2",'+
                                 '"ari1",'+
                                 '"ari2",'+
                                 '"wbi"')
          else:
              vis_nir = None
              vis_nir_hdr = None
          
          return vis_nir,vis_nir_hdr
              #Range NIR-SWIR
          
      def nir_swir(self):
          if self.range[0] <= 900 and self.range[1] >=2000:
              nir_swir = np.column_stack((self.ndwi(),
                                          self.msi(),
                                          self.ndii()))
              nir_swir_hdr = str('"ndwi",'+
                                 '"msi",'+
                                 '"ndii"')                                
              
          else:
              #continue
              print ('not nir-swir')
              nir_swir=None
              nir_swir_hdr=None            
              
              
          return nir_swir, nir_swir_hdr
          #range SWIR
      def swir(self):
          if self.range[1] >=2000:
              swir = np.column_stack((self.ndni(),
                                      self.ndli()))
              
              swir_hdr=str('"ndni",'+
                           '"ndli"')
          else:
              print ('swir-nir')
              swir = None
              swir_hdr = None
              #continue
          return swir,swir_hdr
          
          
      #||||||||||||||||||||| Methods |||||||||||||||||||||||||||||||||||||||||||||||        
      # function to run every permutation of the NDVI type index across the Red / IR
      # ...... VIS / NIR methods ....
      def multi_tbvi (self, red_start=650, red_end=750, ir_start=700, ir_end=850):
          # get the indicies of the regions we're going to use.
          # we've added default values here, but they can happily be overidden
          #start of red
          red_l =np.argmin(np.abs(self.wl-red_start))
          #end of red
          red_r = np.argmin(np.abs(self.wl-red_end))
          #start of ir
          ir_l = np.argmin(np.abs(self.wl-ir_start))
          #end of ir
          ir_r = np.argmin(np.abs(self.wl-ir_end))
          
          #slice
          left = self.values[red_l:red_r]
          right = self.values[ir_l:ir_r]
          
          #set up output
          values = np.empty(3)
          
          #set up counter
          l = 0    
          #loop throught the values in the red
          for lvalue in left:
              l_wl = self.wl[l+red_l]
              r = 0
              l = l+1
              #then calculate the index with each wl in the NIR
              for rvalue in right:
                  value = (rvalue-lvalue)/(rvalue+lvalue)
                  r_wl = self.wl[r+ir_l]
                  out = np.column_stack((l_wl,r_wl,value))
                  values = np.vstack((values, out))
                  out = None
                  r = r+1
          return values[1:,:]
          
      def sr700_800 (self, x=700, y=800):
          index = self.values[np.argmin(np.abs(self.wl-x))]/self.values[np.argmin(np.abs(self.wl-y))]
          return index
      
      def ndvi705_750 (self, x=705, y=750):
          index = (self.values[np.argmin(np.abs(self.wl-y))]-self.values[np.argmin(np.abs(self.wl-x))])/\
                  (self.values[np.argmin(np.abs(self.wl-y))]+self.values[np.argmin(np.abs(self.wl-x))])
          return index    
      
      def ndvi700_800 (self, x=700, y=800):
          index = (self.values[np.argmin(np.abs(self.wl-y))]-self.values[np.argmin(np.abs(self.wl-x))])/\
                  (self.values[np.argmin(np.abs(self.wl-y))]+self.values[np.argmin(np.abs(self.wl-x))])
          return index
          
      def ndvi694_760 (self, x=694, y=760):
          index = (self.values[np.argmin(np.abs(self.wl-y))]-self.values[np.argmin(np.abs(self.wl-x))])/\
                  (self.values[np.argmin(np.abs(self.wl-y))]+self.values[np.argmin(np.abs(self.wl-x))])
          return index
      
      def ndvi695_805 (self, x=695, y=805):
          index = (self.values[np.argmin(np.abs(self.wl-y))]-self.values[np.argmin(np.abs(self.wl-x))])/\
                  (self.values[np.argmin(np.abs(self.wl-y))]+self.values[np.argmin(np.abs(self.wl-x))])
          return index
          
      def npci (self, x=430, y=680):
          index = (self.values[np.argmin(np.abs(self.wl-y))]-self.values[np.argmin(np.abs(self.wl-x))])/\
                  (self.values[np.argmin(np.abs(self.wl-y))]+self.values[np.argmin(np.abs(self.wl-x))])
          return index
          
      def npqi (self, x=415, y=435):
          index = (self.values[np.argmin(np.abs(self.wl-y))]-self.values[np.argmin(np.abs(self.wl-x))])/\
                  (self.values[np.argmin(np.abs(self.wl-y))]+self.values[np.argmin(np.abs(self.wl-x))])
          return index
      
          #mSRvi
      #= (750-445)/(705+445)
      def msrvi (self):   
          x = 750
          y = 445
          z = 705    
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          z_val = self.values[np.argmin(np.abs(self.wl-z))]
          msrvi_val = (x_val-y_val)/(z_val+y_val)
          return msrvi_val
          
      #Vogelmann Red Edge 1
      #740/720
      def vog1vi (self):   
          x = 740
          y = 720    
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          vog1vi_val = (x_val/y_val)
          return vog1vi_val
          
      #Vogelmann Red Edge 2
      #= (734-747)/(715+726)
      def vog2 (self):   
          v = 734
          x = 747
          y = 715
          z = 726    
          v_val = self.values[np.argmin(np.abs(self.wl-v))]
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          z_val = self.values[np.argmin(np.abs(self.wl-z))]
          vog2_val = (v_val-x_val)/(y_val+z_val)
          return vog2_val
          
      #PRI
      # (531-570)/(531+570)
      def privi (self):   
          x = 531
          y = 570    
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          privi_val = (x_val-y_val)/(x_val+y_val)
          return privi_val
                  
      #SIPI
      #(800-445)/(800-680)
      def sipi (self):   
          x = 800
          y = 445
          z = 680    
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          z_val = self.values[np.argmin(np.abs(self.wl-z))]
          sipi_val = (x_val-y_val)/(x_val+z_val)
          return sipi_val
  
      #Water band index
      # WBI = 900/700
      def wbi (self):   
          x = 900
          y = 700    
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          wbi_val = (x_val/y_val)
          return wbi_val
          
      #mNDVI
      #= (750-705)/((750+705)-(445))
      def mdvi (self):
          x = 750
          y = 705
          z = 445    
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          z_val = self.values[np.argmin(np.abs(self.wl-z))]
          mdvi_val = (x_val-y_val)/((x_val+y_val)-z_val)
          return mdvi_val
          
      #Carotenid Reflectance Index
      #CRI1 = (1/510)-(1/550)
      def cri1 (self):   
          x = 510
          y = 550    
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          cri1_val = (1/x_val)-(1/y_val)
          return cri1_val
          
      #CRI2 = (1/510)-(1/700)
      def cri2 (self):   
          x = 510
          y = 700    
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          cri2_val = (1/x_val)-(1/y_val)
          return cri2_val
          
      #Anthocyanin
      #ARI1 = (1/550)-(1/700)
      def ari1 (self):   
          x = 550
          y = 700   
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          ari1_val = (1/x_val)-(1/y_val)
          return ari1_val
          
      #ARI2 = 800*((1/550)-(1/700)_))
      def ari2 (self):
          x = 510
          y = 700   
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          ari2_val = 800*((1/x_val)-(1/y_val))
          return ari2_val
          
      #MSR
      #=((800/670)-1)/SQRT(800+670)
      def msr (self):   
          x = 800
          y = 670   
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          msr_val = ((x_val/y_val)-1)/(np.sqrt(x_val+y_val))
          return msr_val
          
      #SAVI
      #= (1+l)(800-670)/(800+670+l)
      def savi (self, l=0.5):   
          x = 800
          y = 670
          l = 0.5
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          savi_val = ((1+l)*(x_val-y_val))/(x_val+y_val+l)
          return savi_val
          
      #MSAVI
      #=1/2(sqrt(2*800)+1)-SQRT(((2*800+1)sqr)-8*(800-670)
      def msavi2 (self):   
          x = 800
          y = 670
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          msavi2_top1 = (2*x_val+1)
          msavi2_top2 = (np.sqrt(np.square(2*x_val+1)-(8*(x_val-y_val))))
          msavi2_top = msavi2_top1-msavi2_top2
          msavi2_val = msavi2_top/2
          return msavi2_val
          
      #Modified clhoropyll absorption indec
      #MCARI = ((700-670)-0.2*(700-550))*(700/670)
      def mcari (self):   
          x = 700
          y = 670
          z = 550
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          z_val = self.values[np.argmin(np.abs(self.wl-z))]
          mcari_val = (x_val-y_val)-(0.2*(x_val-z_val)*(x_val/y_val))    
          return mcari_val
          
      #Triangular vegetation index
      #TVI 0.5*(120*(750-550))-(200*(670-550))
      def tvi (self):   
          x = 750
          y = 550
          z = 670    
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          z_val = self.values[np.argmin(np.abs(self.wl-z))]
          tvi_val = 0.5*((120*(x_val-y_val))-(200*(z_val+y_val)))
          return tvi_val
          
      #MCAsavRI1 = 1.2*(2.5*(800-67-)-(1.3*800-550)
      def mcari1 (self):   
          x = 800
          y = 670
          z = 550    
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          z_val = self.values[np.argmin(np.abs(self.wl-z))]
          mcari1_val = (1.2*((2.5*(x_val-y_val)))-(1.3*(x_val+z_val)))
          return mcari1_val
          
      #MTVI1
      #=1.2*((1.2*(800-550))-(2.5(670-550)))
      def mtvi (self):   
          x = 800
          y = 550
          z = 670    
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          z_val = self.values[np.argmin(np.abs(self.wl-z))]
          mtvi_val = (1.2*(12*(x_val-y_val)))-(2.5*(z_val-y_val))
          return mtvi_val
  
      def mcari2 (self):   
          x = 800
          y = 670
          z = 550
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          z_val = self.values[np.argmin(np.abs(self.wl-z))]
          mcari2_top = (1.5*(2.5*(x_val-y_val)))-(1.3*(x_val-z_val))
          mcari2_btm = np.sqrt((np.square(2*x_val)+1)-((6*x_val)-(5*(np.sqrt(y_val))))-0.5)
          mcari2_val = mcari2_top/mcari2_btm
          return mcari2_val
          
      #MTVI2=(1.5*(2.5(800-670)-2.5*(800-550))/sqrt((2*800+1s)sq)-((6*800)-(5*sqrt670))-0.5
      def mtvi2 (self):   
          x = 800
          y = 670
          z = 550    
      
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          z_val = self.values[np.argmin(np.abs(self.wl-z))]
          mtvi2_top = (1.5*(2.5*(x_val-z_val)))-(1.3*(x_val-z_val))
          mtvi2_btm = np.sqrt((np.square(2*x_val)+1)-((6*x_val)-(5*(np.sqrt(y_val))))-0.5)
          mtvi2_val = mtvi2_top/mtvi2_btm
          return mtvi2_val
          
      #Renormalised DVI
      #RDVI = (800-670)/sqrt(800+670)
      def rdvi (self):
          x = 800
          y = 670
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          rdvi_val = (x_val-y_val)/np.sqrt(x_val+y_val)
          return rdvi_val
          
      #Plant senescance reflectance index
      #PRSI = (680-500)/750
      def prsi (self):   
          x = 680
          y = 500
          z = 750
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          z_val = self.values[np.argmin(np.abs(self.wl-z))]
          prsi_val = (x_val-y_val)/z_val
          return prsi_val
          
      #||||||||||||||||||||||| SWIR methods ||||||||||||||||||||||||||||||||||||
  
      #Cellulose Absorption Index
      #CAI =0.5*(2000-2200)/2100
      def cai (self):   
          x = 2000
          y = 2200
          z = 2100
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          z_val = self.values[np.argmin(np.abs(self.wl-z))]
          cai_val = 0.5*(x_val-y_val)-z_val
          return cai_val
          
      #Normalized Lignin Difference
      #NDLI = (log(1/1754)-log(1/1680))/(log(1/1754)+log(1/1680))
      def ndli (self):   
          x = 1754
          y = 2680
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))] 
          ndli_val = (np.log(1/x_val)-np.log(1/y_val))/(np.log(1/x_val)+np.log(1/y_val))
          return ndli_val
          
      #Canopy N
      #NDNI =(log(1/1510)-log(1/1680))/(log(1/1510)+log(1/1680))
      def ndni (self):   
          x = 1510
          y = 1680
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          ndni_val = (np.log(1/x_val)-np.log(1/y_val))/(np.log(1/x_val)+np.log(1/y_val))
          return ndni_val
          
      #|||||||||||||||||||||| Full spectrum (VIS-SWIR)||||||||||||||||||||||||||||
      #Normalised Difference IR index
      #NDII = (819-1649)/(819+1649)#NDII = (819-1649)/(819+1649)
      def ndii (self):   
          x = 819
          y = 1649    
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          ndii_val = (x_val-y_val)/(x_val+y_val)
          return ndii_val
      
      #Moisture Stress Index
      #MSI = 1599/819http://askubuntu.com/questions/89826/what-is-tumblerd
      def msi (self):   
          x = 1599
          y = 810    
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          msi_val = (x_val/y_val)
          return msi_val
          
      #NDWI
      #(857-1241)/(857+1241)
      def ndwi (self):   
          x = 857
          y = 1241    
          x_val = self.values[np.argmin(np.abs(self.wl-x))]
          y_val = self.values[np.argmin(np.abs(self.wl-y))]
          ndwi_val = (x_val-y_val)/(x_val+y_val)
          return ndwi_val
  
  class red_edge():
      '''Class to derive red edge position using a number of different methods'''
      def __init__(self, spectra):
          self.wl = spectra[:,0]
          self.values = spectra[:,1]
          self.range = (np.min(self.wl),np.max(self.wl))
          '''Again, the mehtod that initialises this class uses the range of the 
          sensor to check to see if it falls within the red-edge reigion. If so,
          it will derive the red edge using the differnet methods and save these 
          as seprate hdf5 datasets in the appropriate group''' 
          if self.range[0] <= 670 and self.range[1] >=750:
              self.redge_vals = np.column_stack((self.redge_linear(),
                                                 self.redge_lagrange(),
                                                 self.redge_linear_extrapolation()))
              print (self.redge_vals)
              print (self.redge_linear,self.redge_lagrange,self.redge_linear_extrapolation)
              
              self.redge_hdr = str('"linear",'+
                                   '"lagrange",'+
                                   '"extrapolated"')
              
          else:
              print ('red_edge out of range')
              self.redge_vals = None
              self.redge_hdr = None
              
              
      ##################### METHODS #########################################
      #linear- defined by clevers et al 1994:
      def redge_linear(self):
          r670 = self.values[np.argmin(np.abs(self.wl-670))]
          r780 = self.values[np.argmin(np.abs(self.wl-780))]
          r700 = self.values[np.argmin(np.abs(self.wl-700))]
          r740 = self.values[np.argmin(np.abs(self.wl-740))]
          r_edge = (r670+r780)/2
          lin_rep =700+40*((r_edge-r700)/(r740-r700))
          print ('REDGE_LINEAR',lin_rep)
          return lin_rep
          
      #Lagrangian method, after Dawson & Curran 1998
      def redge_lagrange(self):        
          #select the red edge region of the first derviative and associate this
          #with wavelength
          x = 680
          y = 730   
          first_diff = np.diff(self.values, 1) 
          spec_in = np.column_stack((self.wl[1:], first_diff))
          l680 = np.argmin(np.abs(spec_in[:,0]-x))
          r680 = spec_in[l680,0]
          l730 = np.argmin(np.abs(spec_in[:,0]-y))
          r730 = spec_in[l730,0]
          redge_region_sel = np.where(np.logical_and(spec_in[:,0]>r680-1, 
                                                     spec_in[:,0]<r730+1))
          redge_region = spec_in[redge_region_sel]        
          #find the maximum first derivative, return index
          dif_max = np.argmax(redge_region[:,1], axis=0)
          #find band with the max derivative -1, return index     
          dif_max_less = (np.argmax(redge_region[:,1], axis=0))-1
          #find band with the max derivative +1, return index
          dif_max_more = (np.argmax(redge_region[:,1], axis=0))+1 
          
          
          if dif_max_more >= redge_region.shape[0]:
              dif_max_more = redge_region.shape[0]-1
          
          #use these indeces to slice the array    
          rmax = redge_region[dif_max]    
          rmax_less =redge_region[dif_max_less]
          rmax_more =redge_region[dif_max_more]
          #lagrangian interpolation with three points
          #this has been expanded to make the syntax easier    
          a = rmax_less[1]/(rmax_less[0]-rmax[0])*(rmax_less[0]-rmax_more[0])
          b = rmax[1]/(rmax[0]-rmax_less[0])*(rmax[0]-rmax_more[0])
          c = rmax_more[1]/(rmax_more[0]-rmax_less[0])*(rmax_more[0]-rmax[0])
          d = a*(rmax[0]+rmax_more[0])
          e = b*(rmax_less[0]+rmax_more[0])
          f = c*(rmax_less[0]+rmax[0])    
          lg_rep = (d+e+f)/(2*(a+b+c))
          print ('Lagrangian', lg_rep)
          return lg_rep
      
      
      
      
      
      #Linear extrapolation- after Cho & Skidmore 2006, Cho et al 2007
      def redge_linear_extrapolation(self):
          diff = np.diff(self.values)
          d680 = diff[np.argmin(np.abs(self.wl-680+1))]
          d694 = diff[np.argmin(np.abs(self.wl-694+1))]
          
          d724 = diff[np.argmin(np.abs(self.wl-724+1))]
          d760 = diff[np.argmin(np.abs(self.wl-760+1))]
          
          red_slope = ((d694-d680)/(694-680))
          ir_slope = ((d760-d724)/(760-724))
          
          red_inter = d680-(red_slope*680)
          
          ir_inter = d724-(ir_slope*724)
          
          wl = (ir_inter-red_inter)/(ir_slope-red_slope)
  
          print ('^!!!!!!!!! Linear:',wl)
          return np.abs(wl)
  
          
  
          
          
  
  class fluorescence():
      '''this class is inteded to look for evidence of photosynthetic flourescence
      currently this is limited to simple reflectance indices. This should be
      expanded to take in other more complex methods to invesitgae fluorescence'''
      
      def __init__(self, spectra):
          self.wl = spectra[:,0]
          self.values = spectra[:,1]
          self.range = (np.min(self.wl),np.max(self.wl))
          print ('call to fluor')
          '''The init method checks the range to establish if it overlaps with 
          region of chlorophyll flourescence. If so it will will perform the 
          analysis methods and output to hdf5'''
          
  
      def wl_selector(self, x):
          '''this method finds the index of the wavelength closest to that 
          specified for reflectance'''
          value = self.values[np.argmin(np.abs(self.wl-x))]
          return value
      
      def d_wl_selector(self, x):
          '''this method finds the index of the wavelength closest to that 
          specified for the first derivative'''
          diff = np.diff(self.values)
          value = diff[np.argmin(np.abs(self.wl-x))+1]
          return value
          
      def wl_max_d(self):
          '''method to extract wavelength of the maxima of the first derivative
          and return this'''
          start = np.argmin(np.abs(self.wl-650))
          end = np.argmin(np.abs(self.wl-760))
          diff = np.diff(self.values[start:end])
          maxdiff = np.argmax(diff)
          maxdiffwl = self.wl[maxdiff+start+1]
          return maxdiffwl, diff[maxdiff]
          
      def simple_ratios(self):
          ''' This method runs flourescence indices ratios and returns them as a
          stacked numpy array'''
          #r680/r630
          r680r630 = self.wl_selector(680)/self.wl_selector(630)
          print (r680r630)
          #r685/r630
          r685r630 = self.wl_selector(685)/self.wl_selector(630)
          print (r685r630)
          #r685/r655
          r685r655 = self.wl_selector(685)/self.wl_selector(655)
          print (r685r655)
          #r687/r630
          r687r630 = self.wl_selector(687)/self.wl_selector(630)
          print (r687r630)
          #r690/r630
          r690r630 = self.wl_selector(690)/self.wl_selector(630)
          print (r690r630)
          #r750/r800
          r750r800 = self.wl_selector(750)/self.wl_selector(800)
          print (r750r800)
          
          #sq(r685)/(r675-r690)
          sqr685 = np.square(self.wl_selector(685))/(self.wl_selector(675)-self.wl_selector(690))
          print (sqr685)
          
          #(r675-r690)/sq(r683) Zarco-Tejada 2000
          r675r690divsq683 = (self.wl_selector(675)-self.wl_selector(690))/np.square(self.wl_selector(683))
          print (r675r690divsq683)
           
          #d705/d722
          d705d722 = self.d_wl_selector(705)/self.d_wl_selector(722)
          print (d705d722)
          
          #d730/d706
          d730d706 = self.d_wl_selector(730)/self.d_wl_selector(706)
          print (d730d706)
          
          #(d688-d710)/sq(d697)
          d686d710sq697 = (self.d_wl_selector(688)-self.d_wl_selector(710))\
                          /np.square(self.d_wl_selector(697))
          print (d686d710sq697)
          
          #wl at max d / d720
          maxdd720 = self.wl_max_d()[1]/self.d_wl_selector(720)
          print (maxdd720)
          
          #wl at max d / d703
          maxdd703 = self.wl_max_d()[1]/self.d_wl_selector(703)
          print (maxdd703)
          
          #wl at max d / d(max d+12)
          print (self.wl_max_d()[0])
          maxd12 = self.wl_max_d()[1]/self.d_wl_selector(self.wl_max_d()[0]+12)
          print (maxd12)
          
          combined = np.vstack((r680r630,
                                r685r630,
                                r685r655,
                                r687r630,
                                r690r630,
                                r750r800,
                                sqr685,
                                r675r690divsq683,
                                d705d722,
                                d730d706,
                                d686d710sq697,
                                maxdd720,
                                maxdd703,
                                maxd12))
                                
                                
          fluo_hdr = str('"r680r630",'+
                         '"r685r630",'+
                         '"r685r655",'+
                         '"r687r630",'+
                         '"r690r630",'+
                         '"r750r800",'+
                         '"sqr685",'+
                         '"r675r690divsq683",'+
                         '"d705d722",'+
                         '"d730d706",'+
                         '"d686d710sq697",'+
                         '"maxdd720",'+
                         '"maxdd703",'+
                         '"maxd12"')
          
          return combined, fluo_hdr
          
      def dual_peak(self):
          '''This fuction loogs for a dual peak in the red-edge region. If it's
          there it measures the depth of the feature between the two peaks. 
          UNTESTED'''
          start = self.wl_selector(640)
          end = self.wl_selector(740)
          
          d1_region = np.diff(self.values[start:end])
          #d2_region = np.diff(self.values[start:end], n=2)
          
  
          
          peak_finder = find_peaks_cwt(d1_region, np.arange(3,10))
          
              
          peak_wl = wavelengths[peak_finder]
          
          fluor_peaks = []
  
          for peak in peak_finder:
              if peak_wl[peak] == self.wl[self.wl_selector(668)]:
                  print ('found flourescence peak at 668nm')
                  fluor_peaks.append(peak)
                  
              elif peak_wl[peak] == self.wl[self.wl_selector(735)]:
                  print ('found flourescence peak at 735nm')
                  fluor_peaks.append[peak]
                  
              else:
                  print ('unknown peak')
                  
          '''if len(fluor_peaks) == 2:
              something = 'something'''
              
  class load_asd():
                         
      def __init__(self, indir, output_dir):
          data_list = os.listdir(indir)
          print (data_list)
          
          #output_dir = os.path.join(indir,'output')
          if not os.path.exists(output_dir):
              os.mkdir(output_dirx)
          
          for directory in data_list:
              parent = os.path.join(indir, directory)
              spectra_dir = os.path.join(parent, 'raw_spectra')
              reading_info_dir = os.path.join(parent, 'reading_info')
              
              sensor_name = 'ASD FieldSpec Pro'
              sensor_type = 'SPR'
              sensor_units = 'nm'
              sensor_range = [350,2500]
              
              os.chdir(reading_info_dir)
              
              reading_info_file = open('reading_atributes.txt','rb')
              
              reading_info = csv.DictReader(reading_info_file)         
              
              reading_info_array = np.empty(12)
              
              readings_list = [row for row in reading_info]
              
              
              for reading in readings_list[:]:
                  
                  reading_filename = str(reading['reading_id']+'.txt')
                  reading_info_line = np.column_stack((reading['reading_id'],
                                                       reading['dartField'],
                                                       reading['transect'],
                                                       reading['transectPosition'],
                                                       reading['reading_type'],
                                                       reading['reading_coord_osgb_x'],
                                                       reading['reading_coord_osgb_y'],
                                                       reading['dateOfAcquisition'],
                                                       reading['timeOfAcquisition'],
                                                       reading['instrument_number'],
                                                       reading['dark_current'],
                                                       reading['white_ref']))
                  #print reading_info_line
                  
                  if reading['reading_type']== 'REF':
                      reading_info_array = np.vstack((reading_info_array,reading_info_line))
                      #print reading_info_array
                      
                      print ('*********** Loading File', reading_filename, '***********')
                      os.chdir(spectra_dir)
                      spec = np.genfromtxt(reading_filename, 
                                           delimiter=', ',
                                           skiprows=30)
                                           
                      spec = np.column_stack((spec[:,0],spec[:,1]*100))
                                
                      nir_start = 0
                      nir_end = 990
                      nir_weight = 3.5
                      nir_k = 4.9
                      nir_s =45
                      
                      swir1_start = 1080
                      swir1_end = 1438
                      swir1_weight = 8.5
                      swir1_k = 3.5
                      swir1_s = 35
                      
                      swir2_start = 1622
                      swir2_end = 2149
                      swir2_weight = 1.2
                      swir2_s = 92
                      swir2_k = 2.8
                      
                      #smoothing(perc_out, block_start, block_end, kparam, weight, sparam)
                      nir_smoothed = smoothing(spec, nir_start, nir_end, nir_k, nir_weight, nir_s)
                      swir1_smoothed = smoothing(spec, swir1_start, swir1_end, swir1_k, swir1_weight, swir1_s)
                      swir2_smoothed = smoothing(spec, swir2_start, swir2_end, swir2_k, swir2_weight, swir2_s)
                      
                      print ('Smoothed array shape', nir_smoothed.shape,swir1_smoothed.shape,swir2_smoothed.shape)
                      
                      
                      
                      nir_swir_gap = interpolate_gaps(nir_smoothed,swir1_smoothed)
                      swir2_gap = interpolate_gaps(swir1_smoothed,swir2_smoothed)
                      
                   
                      spec_smoothed = np.vstack((nir_smoothed,
                                                 nir_swir_gap,
                                                 swir1_smoothed,
                                                 swir2_gap,
                                                 swir2_smoothed))
                                        
                      print ('Spec SHAPE:', spec.shape)
                      
                      survey_dir = os.path.join(output_dir, directory)
                      if not os.path.exists(survey_dir):
                          os.mkdir(survey_dir)
                          
                      os.chdir(survey_dir)
                          
                      try:
                          abs470 = absorption_feature(spec_smoothed,400,518,484)
                          print (abs470.abs_feature()[0])
                          abs470_ftdef = abs470.abs_feature()[0]
                          print (abs470_ftdef)
                          abs470_crem = abs470.abs_feature()[2]
                          if not abs470_ftdef == None:
                              np.savetxt(reading_filename[0:-4]+'_abs470_ftdef.txt',
                                         abs470_ftdef,
                                         header=abs470.abs_feature()[1],
                                         delimiter=',')
                              np.savetxt(reading_filename[0:-4]+'_abs470_crem.txt',
                                         abs470_crem,
                                         delimiter=',')
                      except:
                          pass
                      try:
                          abs670 = absorption_feature(spec_smoothed,548,800,670)
                          abs670_ftdef = abs670.abs_feature()[0]
                          abs670_crem = abs670.abs_feature()[2]
                          if not abs670_ftdef == None:
                              np.savetxt(reading_filename[0:-4]+'_abs670_ftdef.txt',
                                         abs670_ftdef,
                                         header=abs670.abs_feature()[1],
                                         delimiter=',')
                              np.savetxt(reading_filename[0:-4]+'_abs670_crem.txt',
                                         abs670_crem,
                                         delimiter=',')
                      except:
                          pass
                      
                      try:
                          abs970 = absorption_feature(spec_smoothed,880,1115,970)
                          abs970_ftdef = abs970.abs_feature()[0]
                          abs970_crem = abs970.abs_feature()[2]
                          if not abs970_ftdef == None:
                              np.savetxt(reading_filename[0:-4]+'_abs970_ftdef.txt',
                                         abs970_ftdef,
                                         header=abs970.abs_feature()[1],
                                         delimiter=',')
                              np.savetxt(reading_filename[0:-4]+'_abs970_crem.txt',
                                         abs970_crem,
                                         delimiter=',')
                      except:
                          pass
                      
                      try:
                          abs1200 = absorption_feature(spec_smoothed,1080,1300,1190)
                          abs1200_ftdef = abs1200.abs_feature()[0]
                          abs1200_crem = abs1200.abs_feature()[2]
                          if not abs1200_ftdef == None:
                              np.savetxt(reading_filename[0:-4]+'_abs1200_ftdef.txt',
                                         abs1200_ftdef,
                                         header=abs1200.abs_feature()[1],
                                         delimiter=',')
                              np.savetxt(reading_filename[0:-4]+'_abs1200_crem.txt',
                                         abs1200_crem,
                                         delimiter=',')
                      except:
                          pass
                      try:               
                          abs1730 = absorption_feature(spec_smoothed,1630,1790,1708)
                          abs1730_ftdef = abs1730.abs_feature()[0]
                          abs1730_crem = abs1730.abs_feature()[2]
                          if not abs1730_ftdef == None:
                              np.savetxt(reading_filename[0:-4]+'_abs1730_ftdef.txt',
                                         abs1730_ftdef,
                                         header=abs1730.abs_feature()[1],
                                         delimiter=',')
                              np.savetxt(reading_filename[0:-4]+'_abs1730_crem.txt',
                                         abs1730_crem,
                                         delimiter=',')
                      except:
                          pass
                      
                      print (spec_smoothed.shape)
                      
                      try:
                          abs2100 = absorption_feature(spec_smoothed,2001,2196,2188)
                          abs2100_ftdef = abs2100.abs_feature()[0]
                          abs2100_crem = abs2100.abs_feature()[2]
                          if not abs2100_ftdef == None:
                              np.savetxt(reading_filename[0:-4]+'_abs2100_ftdef.txt',
                                         abs2100_ftdet,
                                         header=abs2100.abs_feature()[1],
                                         delimiter=',')
                              np.savetxt(reading_filename[0:-4]+'_abs2100_crem.txt',
                                         abs2100_crem,
                                         delimiter=',')
                      except:
                          pass
                      
                      veg_indices = Indices(spec_smoothed)
                      indices = np.column_stack((veg_indices.visnir()[0],
                                                 veg_indices.nir_swir()[0],
                                                 veg_indices.swir()[0]))
                      print (veg_indices.visnir()[1],veg_indices.nir_swir()[1],veg_indices.swir()[1])
                      hdr = str(veg_indices.visnir()[1]+','+veg_indices.nir_swir()[1]+','+veg_indices.swir()[1])
                      
                      np.savetxt(reading_filename[0:-4]+'_indices.txt',
                                 indices,
                                 header=hdr,
                                 delimiter=',')
                      
                      mtbvi = veg_indices.multi_tbvi()
                      np.savetxt(reading_filename[0:-4]+'_mtbvi.txt',
                                 mtbvi,
                                 delimiter=',')
                                 
                      redge = red_edge(spec_smoothed)
                      print (redge.redge_vals.shape)
                      print (redge.redge_vals)
                      np.savetxt(reading_filename[0:-4]+'_redge.txt',
                                 redge.redge_vals,
                                 delimiter=',')
                      
                      fluo = fluorescence(spec_smoothed)
                      np.savetxt(reading_filename[0:-4]+'_flou.txt',
                                 np.transpose(fluo.simple_ratios()[0]),
                                 header = fluo.simple_ratios()[1],
                                 delimiter=',')
                      
                       
                      np.savetxt(reading_filename[0:-4]+'_spec.txt',
                                 spec_smoothed,
                                 delimiter=',')
  
  
  class load_image():
      def __init__(self, wavlengths_dir,image_dir,out_dir):
          os.chdir(wavelengths_dir)
          wavelengths = np.genfromtxt('wavelengths.txt')
          print ('wavelengths array', wavelengths)
          os.chdir(image_dir)
          
          image_list = os.listdir(image_dir)
          
          for image in image_list:
              import_image = self.get_image(image)
              image_name = image[:-4]
              print ('IMAGE NAME:', image_name)
              row = 1
              
              img_array = import_image[0]
              print ('Image_array', img_array)
              projection = import_image[1]
              print ('Projection',projection)
              x_size = import_image[2]
              print ('Xdim',x_size)
              y_size = import_image[3]
              print ('Ydim', y_size)
              spatial = import_image[4]
              print (spatial)
              x_top_left = spatial[0]
              ew_pix_size = spatial[1]
              rotation_ew = spatial[2]
              y_top_left = spatial[3]
              rotation_y = spatial[4]
              ns_pixel_size = spatial[5]
              
              print ('Spatial', x_top_left,ew_pix_size,rotation_ew,y_top_left,rotation_y,ns_pixel_size)
              
              print ('IMAGE ARRAY SHAPE',img_array.shape)
              
              
              img_dims = img_array.shape
              
              print (img_dims[0],'/',img_dims[1])
              
              #indices+29
              indices_out = np.zeros((img_dims[0],img_dims[1],29), dtype=np.float32)
              #print indices_out
              #redge=3
              redge_out = np.zeros((img_dims[0],img_dims[1]),dtype=np.float32)
  
              #fluo=14
              fluo_out=np.zeros((img_dims[0],img_dims[1],14), dtype=np.float32)
              
              print ('fluo out', fluo_out.shape)
  
              
              ft470_out = np.zeros((img_dims[0],img_dims[1],13), dtype=np.float32)
              ft670_out = np.zeros((img_dims[0],img_dims[1],13), dtype=np.float32)
              ft970_out = np.zeros((img_dims[0],img_dims[1],13), dtype=np.float32)
              
              x470 = np.argmin(np.abs(wavelengths-400)) 
              y470 = np.argmin(np.abs(wavelengths-518))
              len470 = y470-x470
              cr470_out = np.zeros((img_dims[0],img_dims[1],len470), dtype=np.float32)
                             
              
              x670 = np.argmin(np.abs(wavelengths-548)) 
              y670 = np.argmin(np.abs(wavelengths-800))
              len670 = y670-x670
              cr670_out = np.zeros((img_dims[0],img_dims[1],len670), dtype=np.float32)
              print (cr670_out)
              
              x970 = np.argmin(np.abs(wavelengths-880)) 
              y970 = np.argmin(np.abs(wavelengths-1000))
              len970 = y970-x970
              cr970_out = np.zeros((img_dims[0],img_dims[1],len970), dtype=np.float32)
              
              #print cr970_out
              
              
              print (wavelengths)
              
              row = 0
              print ('***', row, img_dims[0])
              for i in range(0,img_dims[0]):
                  print (i)
                  column = 0
                  #print 'COL',column
                  
                  
                  
                  for j in range(0,img_dims[1]):
                      print ('COLUMN',column)
                      #print 'Pixel',pixel
                      name = '%s_pix-%s_%s' % (image_name,row,column)
                      print ('NAME',name)
                  
                      pixel = img_array[row,column,:]
                      
                      #smoothed = savgol_filter(pixel,5,2)
                      
                      #spec_smoothed = np.column_stack((wavelengths,smoothed))
                      
                      spec_smoothed = np.column_stack((wavelengths,pixel))
                      
                      
                      
                      print (spec_smoothed)
                      veg_indices = Indices(spec_smoothed)
                      indices = veg_indices.visnir()[0]
                      print ('(*&)(*)(*&&^)^)^)*&^)*^)*&', indices)
                      indices_out[row,column,:]=indices
                      
                      fluo = fluorescence(spec_smoothed)
                      fluo_out[row,column,:]=np.transpose(fluo.simple_ratios()[0])
                      
                      redge = red_edge(spec_smoothed)
                      print (redge.redge_vals.shape)
                      redge_out[row,column]= redge.redge_vals[0,2]
                      
                      
                      try:
                          abs470 = absorption_feature(spec_smoothed,400,518,484)
                          abs470_ftdef = abs470.abs_feature()[0]
                          
                          abs470_crem = abs470.abs_feature()[2]
                          abs470_crem = np.column_stack((abs470_crem[:,0],abs470_crem[:,4]))
                          
                          print ('!*!*!*!*!&!*!*', abs470_crem)
                          
                          crem470_fill = self.crem_fill(x470,y470,abs470_crem,wavelengths)
                          
                          ft470_out[row,column,:]=abs470_ftdef
                          cr470_out[row,column,:]=crem470_fill
                          
                      except:
                          pass
                      
                      try:
                          abs670 = absorption_feature(spec_smoothed,548,800,670)
                          abs670_ftdef = abs670.abs_feature()[0]
                          abs670_crem = abs670.abs_feature()[2]
                          abs670_crem = np.column_stack((abs670_crem[:,0],abs670_crem[:,4]))
                          ft670_out[row,column,:]=abs670_ftdef
                          
                          crem670_fill = self.crem_fill(x670,y670,abs670_crem,wavelengths)                            
                          
                          cr670_out[row,column,:]=crem670_fill
                      except:
                          pass
                      
                      try:
                          abs970 = absorption_feature(spec_smoothed,880,1000,970)
                          abs970_ftdef = abs970.abs_feature()[0]
                          abs970_crem = abs970.abs_feature()[2]
                          abs970_crem = np.column_stack((abs970_crem[:,0],abs970_crem[:,4]))
                          
                          crem970_fill = self.crem_fill(x970,y970,abs970_crem,wavelengths)
                          
                          ft970_out[row,column,:]=abs970_ftdef
                          cr970_out[row,column,:]=crem970_fill
                      except:
                          pass
                      
                      column = column+1
                      print (pixel.shape)
                  row = row+1
                      
              
        
              self.writeimage(out_dir,image+'_indices.tif',indices_out,spatial)
              self.writeimage(out_dir,image+'_fluo.tif',fluo_out,spatial)
              self.writeimage(out_dir,image+'_redge.tif',redge_out,spatial)
              self.writeimage(out_dir,image+'_ft470.tif',ft470_out,spatial)
              self.writeimage(out_dir,image+'_cr470.tif',cr470_out,spatial)
              self.writeimage(out_dir,image+'_ft670.tif',ft670_out,spatial)
              self.writeimage(out_dir,image+'_cr670.tif',cr670_out,spatial)
              self.writeimage(out_dir,image+'_ft970.tif',ft970_out,spatial)
              self.writeimage(out_dir,image+'_cr970.tif',cr970_out,spatial)
                  
      def crem_fill(self,xwl,ywl,bna,wavelengths):
          bna_out=np.zeros((ywl-xwl))
          bna_wvl = bna[:,0]
          bna_refl= bna[:,1]
          full_wl = wavelengths[xwl:ywl]
          index = np.argmin(np.abs(wavelengths-bna_wvl[0]))
          bna_out[index:]=bna_refl
  
          return bna_out
      
                  
      def get_image(self, image):
          print ('call to get_image')
          # open the dataset
          dataset = gdal.Open(image, GA_ReadOnly)
          print ('Dataset',dataset)
          # if there's nothign there print error
          if dataset is None: 
              print ('BORK: Could not load file: %s' %(image))
         # otherwise do stuff
          else:
              #get the format
              driver = dataset.GetDriver().ShortName
              #get the x dimension
              xsize = dataset.RasterXSize
              #get the y dimension
              ysize = dataset.RasterYSize
              #get the projection
              proj = dataset.GetProjection()
              #get the number of bands
              bands = dataset.RasterCount
              
              #get the geotransform Returns a list object. This is standard GDAL ordering:
                  #spatial[0] = top left x
                  #spatial[1] = w-e pixel size
                  #spatial[2] = rotation (should be 0)
                  #spatial[3] = top left y
                  #spatial[4] = rotation (should be 0)
                  #spatial[5] = n-s pixel size
              spatial = dataset.GetGeoTransform()
              
              
              
              #print some stuff to console to show  we're paying attention
              print ('Found raster in %s format. Raster has %s bands' %(driver,bands))
              print ('Projected as %s' %(proj))
              print ('Dimensions: %s x %s' %(xsize,ysize))
              
              #instantiate a counter
              count = 1
              
              #OK. This is the bit that catually loads the bands in in a while loop
              # Loop through bands as long as count is equal to or less than total
              while (count<=bands):
                  #show that your computer's fans are whining for a reason
                  print ('Loading band: %s of %s' %(count,bands))
                 
                  #get the band
                  band = dataset.GetRasterBand(count)
                  # load this as a numpy array
                  data_array = band.ReadAsArray()
                  '''data_array = ma.masked_where(data_array == 0, data_array)
                  data_array = data_array.filled(-999)'''
                  data_array = data_array.astype(np.float32, copy=False)
                  
                  # close the band object
                  band = None
               
                  #this bit stacks the bands into a combined numpy array
                  #if it's the first band copy the array directly to the combined one
                  if count == 1:
                      stacked = data_array
                  #else combine these 
                  else:
                      stacked = np.dstack((stacked,data_array))            
                  
                  #stacked = stacked.filled(-999)
                  
                  #just to check it's working
                  #print stacked.shape
                  
                  # increment the counter
                  count = count+1
                  
          #stacked = stacked.astype(np.float32, copy=False)
  
          return stacked,proj,xsize,ysize,spatial
          
      def writeimage(self,
                 outpath, 
                 outname,
                 image,
                 spatial):
  
          data_out = image
          print ('ROWS,COLS',image.shape)
          print ('Call to write image')
          os.chdir(outpath)
          print ('OUTPATH',outpath)
          print ('OUTNAME',outname)
          
          #load the driver for the format of choice
          driver = gdal.GetDriverByName("Gtiff") 
          #create an empty output file
          #get the number of bands we'll need
          try:        
              bands = image.shape[2]
          except:
              bands=1
          
          print ('BANDS OUT', bands)
          
          #file name, x columns, y columns, bands, dtype
          out = driver.Create(outname, image.shape[1], image.shape[0], bands, gdal.GDT_Float32)
          #define the location using coords of top-left corner
          # minimum x, e-w pixel size, rotation, maximum y, n-s pixel size, rotation
          out.SetGeoTransform(spatial)
  
          srs = osr.SpatialReference()
          #get the coodrinate system using the ESPG code
          srs.SetWellKnownGeogCS("EPSG:27700")
          #set pstackedstackedstackedtojection of output file 
          out.SetProjection(srs.ExportToWkt())
          
          band = 1
          
          if bands == 1:
              out.GetRasterBand(band).WriteArray(data_out)
              #set the no data value
              out.GetRasterBand(band).SetNoDataValue(-999)
              #apend the statistics to dataset
              out.GetRasterBand(band).GetStatistics(0,1)
              
              print ('Saving %s/%s' % (band,bands))
              
          else:
              while (band<=bands):
                  data = data_out[:,:,band-1]
                  #write values to empty array
                  out.GetRasterBand(band).WriteArray( data )    
                  #set the no data value
                  out.GetRasterBand(band).SetNoDataValue(-999)
                  #apend the statistics to dataset
                  out.GetRasterBand(band).GetStatistics(0,1)  
                  print ('Saving %s/%s' % (band,bands))
                  band = band+1  
          out = None     
          
          print ('Processing of %s complete' % (outname))       
      
          return outname
      
  
               
  if __name__ == "__main__": 
      #dir_path = os.path.dirname(os.path.abspath('...'))
      #data_root = os.path.join(dir_path, 'data')
  
      data_root = '/home/dav/data/temp/test/test_spec'
      
      for folder in os.listdir(data_root):
          input_dir = os.path.join(data_root,folder)
          print (input_dir)
          
          surveys_list = os.listdir(input_dir)
          print (surveys_list)
          
          for survey_dir in surveys_list:
              print (survey_dir)
              
              site_dir=os.path.join(input_dir,survey_dir)
              print (site_dir)
              image_path = os.path.join(site_dir, 'image')
              print (image_path)
              wavelengths_dir = os.path.join(site_dir, 'wavelengths')
              print (wavelengths_dir)
              out_dir = os.path.join(site_dir,'output')
              if not os.path.exists(out_dir):
                  os.mkdir(out_dir)        
              
              load_image(wavelengths_dir,image_path,out_dir)