text: #=============================================================================
  # These functions create figures for visualising the m/n selection using the chi
  # mapping tool.
  #
  # It creates separate plots for each basin in the DEM.
  #
  # Authors:
  #     Simon M. Mudd
  #     Fiona J. Clubb
  #=============================================================================
  #=============================================================================
  # IMPORT MODULES
  #=============================================================================
  # set backend to run on server
  import matplotlib
  matplotlib.use('Agg')
  
  import numpy as np
  import LSDPlottingTools as LSDP
  import matplotlib.pyplot as plt
  from matplotlib import rcParams
  import matplotlib.ticker as ticker
  import pandas as pd
  from matplotlib import colors
  import math
  import os
  import subprocess
  from matplotlib import cm
  #from shapely.geometry import Polygon
  from LSDMapFigure import PlottingHelpers as Helper
  from LSDMapFigure.PlottingRaster import MapFigure
  from LSDMapFigure.PlottingRaster import BaseRaster
  from LSDPlottingTools import LSDMap_SAPlotting as SA
  from LSDPlottingTools import joyplot
  
  #===========================================
  # Function to make a figure object
  #===========================================
  def makefigure(size_format = "EPSL", aspect_ratio=16./9.):
      """
      This function makes a figure object based on the specified size_format
      of specific journals. This could be implemented at a global level.
      Sets RC params for figures
  
      MDH
  
      """
  
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
      rcParams['text.usetex'] = False
  
      # make the figure
      if type(size_format) is float:
          fig = plt.figure(1, facecolor='white',figsize=(size_format,size_format/aspect_ratio))
      elif size_format == "geomorphology":
          fig = plt.figure(1, facecolor='white',figsize=(6.25,6.25/aspect_ratio))
      elif size_format == "big":
          fig = plt.figure(1, facecolor='white',figsize=(16,16./aspect_ratio))
      elif size_format == "EPSL":
          fig = plt.figure(1, facecolor='white',figsize=(7.48,7.48/aspect_ratio))
      else:
          fig = plt.figure(1, facecolor='white',figsize=(4.92126,4.92126/aspect_ratio))
  
      return fig
  
  #=============================================================================
  #=============================================================================
  # ANALYSIS FUNCTIONS
  # These functions analyse the data (for example, doing the outlier checking)
  #=============================================================================
  def SimpleMaxMLECheck(BasinDF):
      """
      This function checks through the basin dataframe and returns a dict with the basin key and the best fit
      m/n
  
      Args:
          BasinDF: pandas dataframe from the basin csv file.
  
      Returns:
          m_over_n_dict: dictionary with the best fit m/n for each basin, the key is the basin key
          and the value is the best fit m/n
  
      Author: FJC
      """
      # read in the basin csv
      basin_keys = list(BasinDF['basin_key'])
      #print BasinDF
  
      # remove the first 2 columns from DF
      del BasinDF['basin_key']
      del BasinDF['outlet_jn']
  
      # now find the index and value of the max MLE in each row
      MOverNs = list(BasinDF.idxmax(axis=1))
      MOverNs = [float(x.split()[-1]) for x in MOverNs]
      print ("MAX MOVERNS")
      print (MOverNs)
  
      # zip into a dictionary
      MOverNDict = dict(zip(basin_keys, MOverNs))
      return MOverNDict
  
  def GetMOverNRangeMCPoints(BasinDF, start_movern=0.2, d_movern=0.1, n_movern=7):
      """
      This function checks through the MC points basin dataframe and returns the best fit and
      range of m/n values. The range is given as a list in the column (all values of m/n where
      the 3rd quartile is above the first quartile value for the best fit m/n)
  
      Args:
          BasinDF: pandas dataframe from the basin MC points csv file.
          start_movern (float): starting m/n value
          d_movern (float): step in m/n value
          n_movern (float): number of m/n values
  
      Returns:
          dataframe with the basin key, best fit m/n, and range of m/n for each basin
  
      Author: FJC
      """
      from scipy import stats
  
      # get the medians from the dataframe
      MedianDF = BasinDF.filter(regex='median')
  
      # find the median with the highest MLE for each basin
      MaxMedians = list(MedianDF.idxmax(axis=1))
      MaxMedians = [str(x) for x in MaxMedians]
      Median_MOverNs = []
      for median in MaxMedians:
          if "=" in median:
              Median_MOverNs.append(float(median.split("=")[-1]))
          else:
              Median_MOverNs.append(np.nan)
      #Median_MOverNs = [float(x.split("=")[-1]) for x in MaxMedians]
  
      # now find the first quartile that corresponds to this median
      FirstQDF = BasinDF.filter(regex='FQ')
      FirstQF_MLEs = []
      for i, median in enumerate(Median_MOverNs):
          if not np.isnan(median):
              # find the right column
              this_col = "FQ_MLE_m_over_n="+str(median)
              FirstQDF_mask = FirstQDF[this_col]
              FirstQF_MLEs.append(float(FirstQDF_mask.iloc[i]))
          else:
              FirstQF_MLEs.append(np.nan)
  
      # now, for each basin, find the columns in the 3rd quartile which are higher than the first Q MLE
      ThirdQDF = BasinDF.filter(regex='TQ')
  
      # change the column names to just have the m/n values
      column_names = list(ThirdQDF)
      column_names = [x.split("=")[-1] for x in column_names]
      ThirdQDF.columns = column_names
  
      # add the threshold first Q MLEs to the dataframe
      ThirdQDF['threshold'] = pd.Series(FirstQF_MLEs, index=ThirdQDF.index)
      print (ThirdQDF)
      # change DF to a boolean where values are greater than the threshold
      TempDF = ThirdQDF.drop('threshold', 1).gt(ThirdQDF['threshold'], 0)
      # get the column names where the values are greater than the threshold for each basin
      TempDF['Range_MOverNs'] = TempDF.apply(lambda x: ','.join(x.index[x]),axis=1)
  
      end_movern = float(start_movern)+float(d_movern)*(float(n_movern)-1)
  
      # get the m/n values greater than the threshold to a list, then get the min and max
      Min_MOverNs = []
      Max_MOverNs = []
      Range_MOverNs = list(TempDF['Range_MOverNs'])
      for i in range (len(Range_MOverNs)):
          movern_str = Range_MOverNs[i].split(",")
          if '' in movern_str:
              Min_MOverNs.append(np.nan)
              Max_MOverNs.append(np.nan)
          else:
              movern_floats = [float(x) for x in movern_str]
              # we need to linearly interpolate between the nearest m/ns for the min
              # and the max.
              Min_MOverN = min(movern_floats)
              Max_MOverN = max(movern_floats)
              
              
  
              # ok, for the minimum, get the previous m/n value
              if float(Min_MOverN) - float(d_movern) < float(start_movern):
                  new_min_movern = Min_MOverN
              else:
                  print("The minimum theta is: "+str(Min_MOverN))
                  this_key = round((float(Min_MOverN)-float(d_movern)),4)
                  print("The other theta is: "+str(this_key))
                  movern_list = [float(Min_MOverN)-float(d_movern), Min_MOverN]
                  mle_list = [ThirdQDF[str(this_key)][i],ThirdQDF[str(float(Min_MOverN))][i]]
                  slope, intercept, r_value, p_value, std_err = stats.linregress(movern_list, mle_list)
                  new_min_movern = (ThirdQDF['threshold'][i] - intercept)/slope
  
              # for the maximum get the next m/n value
              if float(Max_MOverN) + float(d_movern) > float(end_movern):
                  new_max_movern = Max_MOverN
              else:
                  movern_list = [Max_MOverN, float(Max_MOverN)+float(d_movern)]
  
                  # Alright, In some apparently random cases depending on computer architecture, float to string conversion can be F***** up.
                  # Here is a way to deal with it: it first try the regular code that works for many cases
                  # And only correct the string if needed.
                  # Slightly hacky, it might not work in the extremely rare cases where you try to constrain you concavity at 0.001 precision AND you have the wrong cpu or motherboard or whatever decides to convert 0.75 to 0.7500000000000004.
                  try:
                      # Regular case
                      mle_list = [ThirdQDF[str(Max_MOverN)][i], ThirdQDF[(str(float(Max_MOverN)+float(d_movern)))][i]]
                  except KeyError:
                      # DID NOT WORK, fixing
                      # DEaling with annoying string issue. No time to look for clean solution so here is a hacky way:
                      str_max_movern = "%0.2f" % Max_MOverN #  str(Max_MOverN)# converting
                      dtrdiff ="%0.2f" % (float(Max_MOverN)+float(d_movern)) # str(float(Max_MOverN)+float(d_movern))
                      # # if Conversion failed it screws the string into domething like 0.300000000000000004 instead of 0.3
                      # if(len(str_max_movern)>4):
                      #     #Fixing the extra 0
                      #     str_max_movern = str_max_movern[0:4]
                      # # Code won't work if 0.70000000000004 is converted to 0.70 and need in these case to be reconverted to 0.7 AAAAAAAAAAAAAAAAAAAAAAAAAAAH
                      if(str_max_movern[-1] == "0"):
                          str_max_movern = str_max_movern[:-1]
  
                      # # Same procedure for string 2, like EXACTLY the same
                      # dtrdiff = str(float(Max_MOverN)+float(d_movern))
                      # if(len(dtrdiff)>4):
                      #     #Fixing the extra 0
                      #     dtrdiff = dtrdiff[0:4]
                      if(dtrdiff[-1] == "0"):
                          dtrdiff = dtrdiff[:-1]
  
                      # # Done. Correcting with the correct concavity.
                      # try:
                      mle_list = [ThirdQDF[str_max_movern][i], ThirdQDF[(dtrdiff)][i]]
                      # except KeyError:
  
                  # Out of the bug area
                  slope, intercept, r_value, p_value, std_err = stats.linregress(movern_list, mle_list)
                  new_max_movern = (ThirdQDF['threshold'][i] - intercept)/slope
  
              Min_MOverNs.append(new_min_movern)
              Max_MOverNs.append(new_max_movern)
  
  
      # write the output dataframe
      OutputDF = pd.DataFrame()
      OutputDF['basin_key'] = BasinDF['basin_key']
      OutputDF['Median_MOverNs'] = pd.Series(Median_MOverNs)
      OutputDF['FirstQ_threshold'] = pd.Series(FirstQF_MLEs)
      OutputDF['Min_MOverNs'] = pd.Series(Min_MOverNs)
      OutputDF['Max_MOverNs'] = pd.Series(Max_MOverNs)
  
      return OutputDF
  
  def GetRangeMOverNChiResiduals(DataDirectory, fname_prefix, basin_list=[0], parallel=False):
      """
      This function reads in the CSV files with the chi residuals data and
      calculates the median best fit m/n along with the min and max from the
      first and third quartiles. These are returned as a pandas dataframe because
      I love pandas <3 <3
      SMM Note: this isn't really used since the residual analysis doesn't work. 
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_list: a list of the basins to make the plots for. If an empty list is passed then
          all the basins will be analysed. Default = basin 0.
  
      Returns:
          pandas dataframe with m/n data from the chi residuals analysis
  
      Author: FJC
      """
      # first let's read in the csv files with the residuals data
      if not parallel:
          dfs = Helper.ReadChiResidualsCSVs(DataDirectory,fname_prefix)
      else:
          dfs = Helper.AppendChiResidualsCSVs(DataDirectory,fname_prefix)
  
      # get the best fit m/n from the dataframes
      movern_data = []
      for i in range(len(dfs)):
          dfs[i] = dfs[i][dfs[i]['basin_key'].isin(basin_list)]
          ThisDF = dfs[i][dfs[i] < 0]
  
          indices = []
          for i, row in ThisDF.iterrows():
              # find the first index where it is not a nan
              index = row.first_valid_index()
              if index == None:
                  index = row.index[-1]
              indices.append(float(index.split("=")[-1]))
          movern_data.append(indices)
  
      # the movern_data is a list of lists containing the information.
      # movern_data[0] = median
      # movern_data[1] = first quartile
      # movern_data[2] = third quartile
      # the second index is the basin key.
      # so to get the median of basin 3 you would use movern_data[0][3]
  
      # write the output dataframe
      OutputDF = pd.DataFrame()
      OutputDF['basin_key'] = dfs[0]['basin_key']
      OutputDF['Median_MOverNs'] = pd.Series(movern_data[0])
      OutputDF['FirstQ_MOverNs'] = pd.Series(movern_data[1])
      OutputDF['ThirdQ_MOverNs'] = pd.Series(movern_data[2])
  
      return OutputDF
  
  def GetBestFitMOverNFromDisorder(BasinDF):
      """
      This function looks in the basin DF from the disorder CSV files
      and returns a dict with the basin ID and the best-fit m/n for each one.
      The best fit m/n is the one with the minimum disorder statistic.
  
      Args:
          BasinDF: basin disorder dataframe
  
      Returns: dict with the best fit m/n for each basin, key is the basin ID and value is the best fit m/n
  
      Author: FJC
      """
      # read in the basin csv
      basin_keys = list(BasinDF['basin_key'])
      #print BasinDF
  
      # remove the first 2 columns from DF
      del BasinDF['basin_key']
      del BasinDF['outlet_jn']
  
      # now find the index and value of the max MLE in each row
      MOverNs = list(BasinDF.idxmin(axis=1))
      MOverNs = [float(x.split()[-1]) for x in MOverNs]
      print ("Best-fit m/ns disorder")
      print (MOverNs)
  
      # zip into a dictionary
      MOverNDict = dict(zip(basin_keys, MOverNs))
      return MOverNDict
  
  
  def CompareChiAndSAMOverN(DataDirectory, fname_prefix, basin_list=[0], start_movern=0.2, d_movern=0.1, n_movern=7):
      """
      This function compiles
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_list: a list of the basins to make the plots for. If an empty list is passed then
          all the basins will be analysed. Default = basin 0.
          start_movern (float): the starting m/n value. Default is 0.2
          d_movern (float): the increment between the m/n values. Default is 0.1
          n_movern (float): the number of m/n values analysed. Default is 7.
  
      Author: SMM
      """
  
  
      from LSDPlottingTools import LSDMap_PointTools as PointTools
      from LSDPlottingTools import LSDMap_SAPlotting as SAPlot
  
      # read in binned SA data
      binned_csv_fname = DataDirectory+fname_prefix+'_SAbinned.csv'
      print("I'm reading in the csv file "+binned_csv_fname)
      binnedPointData = PointTools.LSDMap_PointData(binned_csv_fname)
  
      # get the basin keys and check if the basins in the basin list exist
      basin = binnedPointData.QueryData('basin_key')
      basin = [int(x) for x in basin]
      Basin = np.asarray(basin)
      these_basin_keys = np.unique(Basin)
  
      print("The unique basin keys are: ")
      print(these_basin_keys)
  
      final_basin_keys = []
      # A bit of logic for checking keys
      if (len(basin_list) == 0):
          final_basin_keys = these_basin_keys
      else:
          for basin in basin_list:
              if basin not in these_basin_keys:
                  print("You were looking for basin "+str(basin)+ " but it isn't in the basin keys.")
              else:
                  final_basin_keys.append(basin)
  
      print("The final basin keys are:")
      print(final_basin_keys)
  
      # Now get all the m/n values from the basin list
      SA_movern_dict = SAPlot.BinnedRegressionDriver(DataDirectory, fname_prefix,
                                                     basin_keys = final_basin_keys)
  
      # Now do the same with the chi-derived m/n data
      (Outlier_counter, removed_sources_dict, best_fit_movern_dict, MLEs_dict) = CheckMLEOutliers(DataDirectory, fname_prefix, basin_list=final_basin_keys,
                                   start_movern=start_movern, d_movern=d_movern,
                                   n_movern=n_movern)
  
      # Now print to files
      print("The best fit m.n for chi analysis are: ")
      print(best_fit_movern_dict)
  
      # Now plot the figure
      PlotMOverNDicts(DataDirectory, fname_prefix, SA_movern_dict,best_fit_movern_dict, FigFormat = "png", size_format = "ESURF")
  
  def CompareMOverNEstimatesAllMethods(DataDirectory, fname_prefix, basin_list=[0], start_movern=0.2, d_movern=0.1, n_movern=7, parallel=False, Chi_disorder=False):
      """
      This function reads in all the files with the data for the various methods of estimating
      the best fit m/n and produces a summary csv file which has the best fit m/n and uncertainty
      for each basin. This csv can then be used to plot the data.
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_list: a list of the basins to make the plots for. If an empty list is passed then
          all the basins will be analysed. Default = basin 0.
          start_movern (float): the starting m/n value. Default is 0.2
          d_movern (float): the increment between the m/n values. Default is 0.1
          n_movern (float): the number of m/n values analysed. Default is 7.
          Chi_disorder (bool): If true, will include the chi disorder stats
  
      Returns:
          writes a csv with the best fit m/n info for each basin
  
      Author: FJC
      """
      # check if a directory exists for the summary plots. If not then make it.
      summary_directory = DataDirectory+'summary_plots/'
      if not os.path.isdir(summary_directory):
          os.makedirs(summary_directory)
          
      # First, we need to check if the full chi actually is there
      Fname_bootstrap_stats = DataDirectory+fname_prefix+'_movernstats_basinstats.csv'
      from pathlib import Path
  
      bootstrap_exists = False
      my_file = Path(Fname_bootstrap_stats)
      if my_file.is_file():
          bootstrap_exists = True
          print("The bootstrap data exists. I hope it was worth the wait. ")
      else:
          print("I did not find bootstrap data, so I am assuming only S-A data and disorder data.")
          
      
  
      # read in the full chi dataframe
      if not parallel:
          FullChiBasinDF = Helper.ReadBasinStatsCSV(DataDirectory,fname_prefix)
      else:
          FullChiBasinDF = Helper.AppendBasinCSVs(DataDirectory,fname_prefix)
  
      # Let's get the list of basins
      if basin_list == []:
          print("You didn't supply a list of basins, so I'm going to analyse all of them.")
          basin_list = list(FullChiBasinDF['basin_key'])
  
      # We're going to write our summary csv using pandas, woooo
      # First of all we need to set up the dataframe
      OutDF = pd.DataFrame()
      OutDF['basin_key'] = pd.Series(basin_list)
  
      # get the best fit m/n for each basin in the list from the full chi method
      if bootstrap_exists:
          FullChiMOverNDict = SimpleMaxMLECheck(FullChiBasinDF)
          OutDF['Chi_MLE_full'] = OutDF['basin_key'].map(FullChiMOverNDict)
  
      # get the best fit m/n from the points method
      if bootstrap_exists:
          if not parallel:
              PointsChiBasinDF = Helper.ReadMCPointsCSV(DataDirectory,fname_prefix)
          else:
              PointsChiBasinDF = Helper.AppendBasinPointCSVs(DataDirectory,fname_prefix)
          PointsChiBasinDF = PointsChiBasinDF[PointsChiBasinDF['basin_key'].isin(basin_list)]
  
          UncertaintyDF = GetMOverNRangeMCPoints(PointsChiBasinDF, start_movern, d_movern, n_movern)
          OutDF['Chi_MLE_points'] = UncertaintyDF['Median_MOverNs']
          OutDF['Chi_MLE_points_min'] = UncertaintyDF['Min_MOverNs']
          OutDF['Chi_MLE_points_max'] = UncertaintyDF['Max_MOverNs']
  
      print ("Getting the m/n from the SA data")
  
      # Updated on the 21/01/2019 by Boris: I added error management on that part. 
      # Cedric Roerig from Giessen is having issue with plotting routines and the code complains data is empty raising ValueError
      # This should sort it
      # Data is replaced by zeros for SA
      # Let me know if any issue appears because of that
      # B.G.
  
      # get the best fit m/n from the raw SA data
      try:
          RawSADF = SA.LinearRegressionRawData(DataDirectory,fname_prefix,basin_list)
          OutDF['SA_raw'] = RawSADF['regression_slope']
          OutDF['SA_raw_sterr'] = RawSADF['std_err']
          OutDF['SA_raw_R2'] = RawSADF['R2']
          OutDF['SA_raw_p'] = RawSADF['p_value']
      except ValueError:
          print("Your Slope-area data is somehow empty or not conform... I am skipping it and data will be 0")
          OutDF['SA_raw'] = pd.Series(data = np.zeros(OutDF.shape[0]), index = OutDF.index)
          OutDF['SA_raw_sterr'] = pd.Series(data = np.zeros(OutDF.shape[0]), index = OutDF.index)
          OutDF['SA_raw_R2'] = pd.Series(data = np.zeros(OutDF.shape[0]), index = OutDF.index)
          OutDF['SA_raw_p'] = pd.Series(data = np.zeros(OutDF.shape[0]), index = OutDF.index)
  
      # get the SA tributary information
      try:
          SATribsDF = SA.GetRangeMOverNRawDataByChannel(DataDirectory,fname_prefix,basin_list)
          OutDF['SA_tribs'] = SATribsDF['median_movern']
          OutDF['SA_tribs_min'] = SATribsDF['FirstQ_movern']
          OutDF['SA_tribs_max'] = SATribsDF['ThirdQ_movern']
      except ValueError:
          print("Your Slope-area data is somehow empty or not conform... I am skipping it and data will be 0")
          OutDF['SA_tribs'] = pd.Series(data = np.zeros(OutDF.shape[0]), index = OutDF.index)
          OutDF['SA_tribs_min'] =pd.Series(data = np.zeros(OutDF.shape[0]), index = OutDF.index)
          OutDF['SA_tribs_max'] = pd.Series(data = np.zeros(OutDF.shape[0]), index = OutDF.index)
  
  
      # get the best fit m/n from the segmented SA data
      try:
          SASegmentedDF = SA.GetRangeMOverNSegmentedData(DataDirectory,fname_prefix,basin_list)
          OutDF['SA_segments'] = SASegmentedDF['median_movern']
          OutDF['SA_segments_min'] = SASegmentedDF['FirstQ_movern']
          OutDF['SA_segments_max'] = SASegmentedDF['ThirdQ_movern']
      except ValueError:
          print("Your Slope-area data is somehow empty or not conform... I am skipping it and data will be 0")
          OutDF['SA_segments'] = pd.Series(data = np.zeros(OutDF.shape[0]), index = OutDF.index)
          OutDF['SA_segments_min'] = pd.Series(data = np.zeros(OutDF.shape[0]), index = OutDF.index)
          OutDF['SA_segments_max'] = pd.Series(data = np.zeros(OutDF.shape[0]), index = OutDF.index)
  
      if Chi_disorder:
          if not parallel:
              DisorderDF = Helper.ReadDisorderUncertCSV(DataDirectory,fname_prefix)
          else:
              DisorderDF = Helper.AppendDisorderCSV(DataDirectory,fname_prefix)
          OutDF['Chi_disorder'] = DisorderDF['median']
          OutDF['Chi_disorder_min'] = DisorderDF['first_quartile']
          OutDF['Chi_disorder_max'] = DisorderDF['third_quartile']
  
      # print the SA segment data
      OutSAname = "_SA_segment_summary.csv"
      out_sa_name = summary_directory+fname_prefix+OutSAname
      SASegmentedDF.to_csv(out_sa_name,index=False)
  
      # now write the output dataframe to a csv file
      OutCSVname = "_movern_summary.csv"
      outname = summary_directory+fname_prefix+OutCSVname
      OutDF.to_csv(outname,index=False)
  
  
  def CheckMLEOutliers(DataDirectory, fname_prefix, basin_list=[0], start_movern=0.2, d_movern=0.1, n_movern=7, parallel=False):
      """
      This function uses the fullstats files to search for outliers in the
      channels. It loops through m/n values and for each m/n value calculates which
      tributaries have outlying MLE and RMSE values. Each time a tributary has an outlying value,
      an outlier counter is incremented so that for each basin we end up with the number
      of times each of its tributaries has been selected as an outlier. This outlier counter is
      returned. From this outlier counter, the tributaries with outliers are ranked in
      descending order of outlier counts. These are iteratively removed. So for example,
      if one tributary has 4 outlier counts and another has 2, the first trib is removed,
      and MLE is recalculated, and then the second is also removed and MLE is recalculated again.
      The removed tributaries and the order in which they are removed is contained within
      the dictionary removed_sources_dict. Finally, a dictionary with the m/n value
      for each basin with the lowest MLE is returned. Each basin has a list where the first
      element is the MLE with no tributaries removed, the second with the first outlier
      count removed, etc. Note that if two tributaries ahve the same outlier counts
      then these are removed at the same time. The key into the dictionaries are the basin numbers.
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_list: a list of the basins to make the plots for. If an empty list is passed then
          all the basins will be analysed. Default = basin 0.
          start_movern (float): the starting m/n value. Default is 0.2
          d_movern (float): the increment between the m/n values. Default is 0.1
          n_movern (float): the number of m/n values analysed. Default is 7.
  
      Returns:
          Outlier_counter (dict): This is a dictionary where the key is the basin
          number and the values are lists that have the cumulative number of times
          each tributary in the basin is detected to be an outlier.
          removed_sources_dict (dict): This is a dictionary where the key is the basin
          number and the values are lists of lists where the top level list is the steps of removing channels
          and the second level are the channels themselves. The reason there are nested
          lists and not individual channels being removed is because sometimes two
          channels have the same outlier counts and so they are removed at the same time.
          best_fit_movern_dict (dict): This is a dictionary where the key is the basin
          number and the values lists of the m/n ratio that has the lowest MLE. These are lists beacuse
          each element in the list represents the number of outlying tributaries removed.
          The first element is where no tributaries are removed.
          MLEs_dict (dictionary of arrays): The key is the basin number and the value is an
          array containing all the MLE values for each m/n and
          each iterated removed tributary.
          To get the MLEs of, say, the first removal for the 3rd basin you would use MLEs_dict[2][:,1]
  
      Author: SMM
      """
  
      # Get a vector of the m over n values
      print("start theta is: "+str(start_movern))
      print("d theta is: "+str(d_movern))
      print("n theta is: "+str(n_movern))
      end_movern = float(start_movern)+float(d_movern)*(float(n_movern)-1)
      print("end theta is: "+str(end_movern))                                                  
      m_over_n_values = np.linspace(start_movern,end_movern,n_movern)
      movern_strs = []
      for movern in m_over_n_values:
          movern_str = '%.2f'%movern
          if movern_str.endswith('0'):
              movern_str = movern_str[:-1]
              movern_strs.append(movern_str)
  
      # we open the first file just so that we can get a counter list
      print ("PARALLEL = ", parallel)
      if not parallel:
          FirstDF = Helper.ReadFullStatsCSV(DataDirectory,fname_prefix,movern_strs[0])
      else:
          FirstDF = Helper.AppendFullStatsCSVs(DataDirectory,movern_strs[0],fname_prefix)
  
      # get the number of basins
      basin_keys = list(FirstDF['basin_key'])
      basin_keys = [float(x) for x in basin_keys]
  
      # get the list of basins
      if basin_list == []:
          print("You didn't give me a list of basins, so I'll just run the analysis on all of them!")
          basin_list = basin_keys
          basin_set = set(basin_list)
          basin_list = list(basin_set)
          basin_list = [int(i) for i in basin_list]
  
      # make a data object that will hold the counters
      Outlier_counter = {}
      # loop through the basins
      for basin in basin_list:
          # mask the data so you only get the correct basin
          FirstDF_basin = FirstDF[FirstDF['basin_key'] == basin]
  
          trib_values = list(FirstDF_basin['test_source_key'])
          n_nodes = len(trib_values)
          # make the counter with zeros
          this_counter = np.zeros(n_nodes)
          Outlier_counter[basin] = this_counter
  
      # Now we loop through all the files, calculating the outliers
      for m_over_n in movern_strs:
  
          print("This_m_over_n is: "+m_over_n)
  
          #load the file
          if not parallel:
              FullStatsDF = Helper.ReadFullStatsCSV(DataDirectory,fname_prefix,m_over_n)
          else:
              FullStatsDF = Helper.AppendFullStatsCSVs(DataDirectory,m_over_n,fname_prefix)
  
  
          # loop through the basins
          for basin in basin_list:
  
              # mask the data so you only get the correct basin
              FullStatsDF_basin = FullStatsDF[FullStatsDF['basin_key'] == basin]
  
              # extract the relevant data
              MLE_values = list(FullStatsDF_basin['MLE'])
              RMSE_values = list(FullStatsDF_basin['RMSE'])
              trib_values = list(FullStatsDF_basin['test_source_key'])
              print ('N values: '+str(len(MLE_values))+' '+str(len(RMSE_values))+' '+str(len(trib_values)))
  
              # now get the outliers
              MLE_array = np.asarray(MLE_values)
              RMSE_array = np.asarray(RMSE_values)
  
              # Get the outliers using the MAD-based outlier function
              RMSE_outliers = LSDP.lsdstatsutilities.is_outlier(RMSE_array)
              MLE_outliers = LSDP.lsdstatsutilities.is_outlier(RMSE_array)
  
              # now check each of the outlier arrays to see if e need to flip the array
              RMSE_index_min = np.argmin(RMSE_array)
              #RMSE_min = np.min(RMSE_array)
  
              # if the max MLE is an outlier, flip the outlier vector
              if (RMSE_outliers[RMSE_index_min]):
                  RMSE_outliers = [not i for i in RMSE_outliers]
  
              MLE_index_max = np.argmax(MLE_array)
  
              # if the max MLE is an outlier, flip the outlier vector
              if (MLE_outliers[MLE_index_max]):
                  MLE_outliers = [not i for i in MLE_outliers]
  
              # turn the outliers vector into an integer
              int_Outlier = [int(i) for i in MLE_outliers]
              #print("Integer outliers are: ")
              #print(int_Outlier)
  
              # add this outlier counter to the outlier dict
              Outlier_counter[basin] = Outlier_counter[basin]+int_Outlier
  
      # Now try to calculate MLE by removing outliers
  
      # Set up dicts where the keys are the basin numbers and the
      # elements tell the m/n values of the best fit MLE, with each element
      # representing incrementally removed tributaries
      # The removed_sources_dict refers to the sources being removed each step
      best_fit_movern_dict = {}
      removed_sources_dict = {}
      MLEs_dict = {}
      for basin_number in basin_list:
          remove_list_index,movern_of_max_MLE,MLEs = Iteratively_recalculate_MLE_removing_outliers_for_basin(Outlier_counter,
                                                                                  DataDirectory,
                                                                                  fname_prefix,
                                                                                  basin_number,
                                                                                  start_movern,
                                                                                  d_movern,
                                                                                  n_movern, parallel)
          best_fit_movern_dict[basin_number] = movern_of_max_MLE
          removed_sources_dict[basin_number] = remove_list_index
          MLEs_dict[basin_number] = MLEs
  
      print("Here are the vitalstatisix, chief: ")
      print(best_fit_movern_dict)
  
      return Outlier_counter, removed_sources_dict, best_fit_movern_dict, MLEs_dict
  
  def Iteratively_recalculate_MLE_removing_outliers_for_basin(Outlier_counter, DataDirectory, fname_prefix, basin_number, start_movern=0.2, d_movern=0.1, n_movern=7, parallel=False):
      """
      This function drives the calculations for removing outliers incrementally
      from the MLE calculations. This is specific to a basin.
      It calls functions for specific m/n values
  
      Args:
          Outlier_counter (dict): The dictionary containing the outlier lists for each basin
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_number (list): a list of the basins to make the plots for. If an empty list is passed then
          all the basins will be analysed. Default = basin 0.
          start_movern (float): the starting m/n value. Default is 0.2
          d_movern (float): the increment between the m/n values. Default is 0.1
          n_movern (float): the number of m/n values analysed. Default is 7.
  
      Returns:
          remove_list_index (list of list): This is the sequence of tributaries that will be removed
          movern_of_max_MLE (list): The m over n values of the maximum MLE sequentially for removed tributaries
          MLEs (array): This is an array containing all the MLE values for each m/n and
          each iterated removed tributary. To get the MLEs of, say, the first removal you would use MLEs[:,1]
  
      Author: SMM
      """
  
      # get the outlier counter for this basin
      thisBasinOutlierCounter = Outlier_counter[basin_number]
  
      print("I am going to recalculate MLE for basin: "+str(basin_number))
      #print("The counter list is:")
      #print(thisBasinOutlierCounter)
  
      # Get the sroted version and the indices into the sorted version
      sort_index = np.argsort(thisBasinOutlierCounter)
      sorted_outliers = np.sort(thisBasinOutlierCounter)
  
      # make sure the sourted outliers are ints
      int_sorted_outliers = [int(i) for i in sorted_outliers]
  
      # we need to reverse these lists so that the biggest outlier counts come first
      sort_index = sort_index[::-1]
      sorted_outliers = sorted_outliers[::-1]
  
      # get all the duplicates, with the total number of duplicates for each counter
      # This uses the unbelievably handy collections.Counter tool
      from collections import Counter
      all_counter_dict=Counter(int_sorted_outliers)
  
      # pop out the zero duplicates: we don't exclude non-outlier data
      all_counter_dict.pop(0, None)
  
      # now we need to iteratively remove the offending counters.
      remove_list = []
      remove_list_index = []
      last_count = -1
  
      # Only enter the loop if there are positive counts
      if all_counter_dict != 0:
  
          # Now loop through the sorted outlier count
          for idx,sorted_outlier_count in enumerate(sorted_outliers):
  
              # first check if this has more than one count
              if sorted_outlier_count == 0:
                  all_removed = True
              else:
                  # if this is the first element, add it to the remove list
                  #if len(remove_list) == 0:
                  #    remove_list.append([])
                  #    remove_list_index.append([])
  
                  # get the number of counts
                  #n_counts = all_counter_dict[sorted_outlier_count]
  
                  # either append the count and index to the current count
                  # or make a new list for the next count
                  if sorted_outlier_count != last_count:
                      remove_list.append([])
                      remove_list_index.append([])
  
                  remove_list[-1].append(int(sorted_outlier_count))
                  remove_list_index[-1].append(sort_index[idx])
  
                  last_count = sorted_outlier_count
  
      # now print out the lists
      #print("remove list is: ")
      #print(remove_list)
  
      #print("and index of this list is: ")
      #print(remove_list_index)
  
      # now we loop through m over n, incrementally removing the tributaries
      # this is done by copying the MLE vector and then replacing
      # the offending channels with an MLE of 1
      # Get a vector of the m over n values
      end_movern = float(start_movern)+float(d_movern)*(float(n_movern)-1)
      m_over_n_values = np.linspace(start_movern,end_movern,n_movern)
  
      # Now get the movern values
      movern_of_max_MLE, MLEs = Calculate_movern_after_iteratively_removing_outliers(m_over_n_values,
                                                                               DataDirectory,
                                                                               fname_prefix,
                                                                               basin_number,
                                                                               remove_list_index, parallel)
  
      # Returns the remove_list_index, which is a list where each element
      # is a list of tributaries removed in an iteration,
      # And the maximum MLE values in each iteration.
      return remove_list_index,movern_of_max_MLE, MLEs
  
  
  
  def Calculate_movern_after_iteratively_removing_outliers(movern_list, DataDirectory,
                                                           fname_prefix, basin_number,
                                                           remove_list_index, parallel=False):
      """
      This function takes the remove list index, which contains information about
      the sequence of tributaries to be removed, and then recalculates MLE by incrementally
      removing tributaries.
  
      Args:
          movern_list (float): m/n value.
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_number (int): The basin you want
          remove_list_index (list of lists): This contains information about what tributaries to remove
  
      Returns:
          movern_of_max_MLE (list): A list containing the m/n values of the basin after outlying
          tributaries have been removed. Each element in the list represents an increment of
          tributary removed. The first element is with no tributaries removed.
          MLEs (array): This is an array containing all the MLE values for each m/n and
          each iterated removed tributary. To get the MLEs of, say, the first removal you would use MLEs[:,1]
  
      Author: SMM
      """
      # Loop through m over n values and recalculate MLE values after removing
      # the outlying data
      All_MLE = []
      for m_over_n in movern_list:
          movern_str = '%.2f'%m_over_n     # need to convert to a string here so files are read correctly. FJC 16/02/18
          if movern_str.endswith('0'):
              movern_str = movern_str[:-1]
      #    print movern_str
  
          MLE_vals = RecalculateTotalMLEWithRemoveList(DataDirectory, fname_prefix,
                                                       movern_str,basin_number, remove_list_index, parallel)
          All_MLE.append(MLE_vals)
  
  
      MLEs = np.asarray(All_MLE)
      index_of_maximums = np.argmax(MLEs,0)
  
      movern_of_max_MLE = []
      for index in index_of_maximums:
          movern_of_max_MLE.append(movern_list[index])
  
      # This is required because linspace gives floating point errors
      movern_of_max_MLE = np.around(movern_of_max_MLE,4)
  
      print("The MLEs for no removal are: ")
      print(MLEs[:,0])
  
      # Return the m/n ratio with the biggest MLE, but also the array
      # with all the MLE values
      return movern_of_max_MLE, MLEs
  
  def RecalculateTotalMLEWithRemoveList(DataDirectory, fname_prefix,
                                        movern,basin_number, remove_list_index, parallel=False):
      """
      This function takes the remove list index and then recalculates MLE by incrementally
      removing tributaries
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          movern (float): m/n value.
          basin_number (int): The basin you want
          remove_list_index (list of lists): This contains information about what tributaries to remove
  
      Returns:
          MLE_vals (list): The MLE data with incrementally removed tributaries
  
      Author: SMM
      """
      #load the file
      if not parallel:
          FullStatsDF = Helper.ReadFullStatsCSV(DataDirectory,fname_prefix,movern)
      else:
          FullStatsDF = Helper.AppendFullStatsCSVs(DataDirectory,movern,fname_prefix)
  
      # mask the data so you only get the correct basin
      FullStatsDF_basin = FullStatsDF[FullStatsDF['basin_key'] == basin_number]
  
      # extract the relevant data
      MLE_values = list(FullStatsDF_basin['MLE'])
      #trib_values = list(FullStatsDF_basin['test_source_key'])
  
      # get the MLE_values as an array
      MLE_values = np.asarray(MLE_values)
  
      #print("The total MLE is: ")
      #print(np.prod(MLE_values))
  
      MLE_vals = []
      MLE_vals.append(np.prod(MLE_values))
  
      # now loop through the remove list
      remove_list = []
      for stuff_to_remove in remove_list_index:
          this_MLE = MLE_values
  
          # need to extend the remove list with these tribs
          remove_list.extend(stuff_to_remove)
  
          # now set the MLE of the removed tribs to 1
          for idx in remove_list:
              this_MLE[idx] = 1
  
          #print("\n REMOVING; The total MLE is: ")
          #print(np.prod(this_MLE))
          MLE_vals.append(np.prod(this_MLE))
  
      return MLE_vals
  
  #=============================================================================
  # PLOTTING FUNCTIONS
  # Make plots of the m/n analysis
  #=============================================================================
  
  def MakePlotsWithMLEStats(DataDirectory, fname_prefix, basin_list = [0],
                    start_movern = 0.2, d_movern = 0.1, n_movern = 7, parallel=False):
      """
      This function makes a chi-elevation plot for each basin and each value of
      m/n and prints the MLE value between the tributaries and the main stem.
      The plot with the maximum value of MLE is highlighted in red (suggesting
      that this should be the appropriate m/n value for this basin). Channels
      are coloured by elevation.
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_list: a list of the basins to make the plots for. If an empty list is passed then
          all the basins will be analysed. Default = basin 0.
          start_movern (float): the starting m/n value. Default is 0.2
          d_movern (float): the increment between the m/n values. Default is 0.1
          n_movern (float): the number of m/n values analysed. Default is 7.
  
      Returns:
          Plot of each m/n value for each basin.
  
      Author: SMM, modified by FJC
      """
  
      # check if a directory exists for the chi plots. If not then make it.
      MLE_directory = DataDirectory+'basic_chi_plots/'
      if not os.path.isdir(MLE_directory):
          os.makedirs(MLE_directory)
  
      profile_suffix = "_movern.csv"
      basin_stats_suffix = "_movernstats_basinstats.csv"
  
      movern_profile_file = fname_prefix+profile_suffix
      movern_basin_stats_file = fname_prefix+basin_stats_suffix
  
      end_movern = start_movern+d_movern*(n_movern-1)
      m_over_n_values = np.linspace(start_movern,end_movern,n_movern)
  
      # get the maximum MLE of each basin
      if not parallel:
          pd_DF = Helper.ReadBasinStatsCSV(DataDirectory, fname_prefix)
      else:
          pd_DF = Helper.AppendBasinStatsCSVs(DataDirectory,fname_prefix)
  
      shp = pd_DF.shape
      max_MLEs = []
      max_MLEs_index = []
      for i in range(0,shp[0]):
          #print("I is: "+str(i))
          a = pd_DF.loc[[i]]
          b = np.asarray(a)
          c = b[0,2:]
          max_MLEs.append(max(c))
          max_MLEs_index.append(np.argmax(c))
      m_over_n_of_max = []
  
      for idx in max_MLEs_index:
          m_over_n_of_max.append(m_over_n_values[idx])
  
      print("The m over n of these max are: ")
      print(m_over_n_of_max)
  
      n_basins = len(max_MLEs)
  
      label_size = 12
  
      # Set up fonts for plots
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
      size_format = "default"
  
      # make a figure
      fig = makefigure(size_format)
  
      gs = plt.GridSpec(100,100,bottom=0.15,left=0.1,right=1.0,top=1.0)
      ax = fig.add_subplot(gs[5:100,10:95])
  
      # load the m_over_n data file
      thisPointData = LSDP.LSDMap_PointData(DataDirectory,movern_profile_file,parallel=parallel)
      allBasinStatsData = LSDP.LSDMap_PointData(DataDirectory+movern_basin_stats_file)
  
      print("m over n values are: ")
      print(m_over_n_values)
  
      mn_legends = []
      for mn in m_over_n_values:
          mn_legends.append("m_over_n = "+str(round(mn,4)))
  
      print("The mn labels are: ")
      print(mn_legends)
  
      # get the data form the profiles
      elevation = thisPointData.QueryData('elevation')
      elevation = [float(x) for x in elevation]
      basin = thisPointData.QueryData('basin_key')
      basin = [int(x) for x in basin]
      source = thisPointData.QueryData('source_key')
      source = [int(x) for x in source]
  
      # get the basin keys
      allstats_basinkeys = allBasinStatsData.QueryData("basin_key")
      allstats_basinkeys = [int(x) for x in allstats_basinkeys]
  
      # need to convert everything into arrays so we can mask different basins
      Elevation = np.asarray(elevation)
      Basin = np.asarray(basin)
      #Source = np.asarray(source)
  
      # Loop through m/n values aggregating data
      for idx,mn in enumerate(m_over_n_values):
  
          counter = str(idx).zfill(3)
          print("Counter is: "+counter)
          # first get the chi values for this m_over_n
          #print ("mn is: " + str(mn))
          #print("index is: "+str(idx))
          mn_legend = "m_over_n = "+str(mn)
          #print("I am looking for the data element: "+mn_legend)
          this_chi = thisPointData.QueryData(mn_legend)
  
          # get the MLE value for this m/n
          this_MLE = allBasinStatsData.QueryData(mn_legend)
  
          # convert to a numpy array for masking
          Chi = np.asarray(this_chi)
  
          # some info about the chi and elevation values
          #max_chi = np.amax(Chi)
          #max_Elevation = np.amax(Elevation)
          #min_Elevation = np.amin(Elevation)
  
          #z_axis_min = int(min_Elevation/10)*10
          #z_axis_max = int(max_Elevation/10)*10+10
          #chi_axis_max = int(max_chi/5)*5+5
  
          # Now mask the data. Initially we will do only basin 0
          if basin_list == []:
              print("You didn't give me any basins so I assume you want all of them.")
              basin_list = range(0,n_basins-1)
  
          for basin_key in basin_list:
              # now we need to find out if this basin is in the allstats file,
              # and if so what index it is
              #this_basin_index = -99
              #for ii,bk in enumerate(allstats_basinkeys):
              #    #print("index: "+str(ii)+" and basin_key is: "+str(bk))
              #    if (bk == basin_key):
              #        this_basin_index = idx
  
              #if(this_basin_index != -99):
              #    MLE = this_MLE[this_basin_index]
              #else:
              #    MLE = "NaN"
              MLE = this_MLE[basin_key]
              #MLE_str = str(MLE)
              #short_MLE = str("%03.02e" % round(MLE,2))
              short_MLE = str(round(MLE,3))
              print("The short MLE is: "+short_MLE)
  
              #print("The MLE of this basin for this m over n is: "+short_MLE)
  
              # this gets the mask (for the chosen basin)
              m = np.ma.masked_where(Basin!=basin_key, Basin)
  
              # this is the masked chi value
              maskX = np.ma.masked_where(np.ma.getmask(m), Chi)
              maskElevation = np.ma.masked_where(np.ma.getmask(m), Elevation)
  
              # now plot the data with a colourmap
              ax.scatter(maskX,maskElevation,s=2.5, c=maskElevation,cmap="terrain",edgecolors='none')
  
              # some formatting of the figure
              ax.spines['top'].set_linewidth(1)
              ax.spines['left'].set_linewidth(1)
              ax.spines['right'].set_linewidth(1)
              ax.spines['bottom'].set_linewidth(1)
  
              # make the lables
              ax.set_xlabel("$\chi$ (m)")
              ax.set_ylabel("Elevation (m)")
  
              # This affects all axes because we set share_all = True.
              #ax.set_ylim(z_axis_min,z_axis_max)
              #ax.set_ylim(z_axis_min,z_axis_max)
              #ax.set_xlim(0,chi_axis_max)
              #ax.text("Basin = " +mn_legend+", $m/n$ = "+str(mn))
  
              #newline = "\n"
              title_string = "Basin "+str(basin_key)+", $m/n$ = "+str(mn)
              title_string2 = "MLE = "+short_MLE
              ax.text(0.05, 0.95, title_string,
                      verticalalignment='top', horizontalalignment='left',
                      transform=ax.transAxes,
                      color='black', fontsize=10)
              # print("The basin index is: "+str(basin_key)+" and the max index is: "+str(max_MLEs_index[basin_key]))
              # if( idx == max_MLEs_index[basin_key]):
              #     print("This m/n is: "+str(mn)+" and it is the maximum MLE")
              #     ax.text(0.05, 0.88, title_string2+", maximum MLE in basin.",
              #         verticalalignment='top', horizontalalignment='left',
              #         transform=ax.transAxes,
              #         color='red', fontsize=10)
              # else:
              #     ax.text(0.05, 0.88, title_string2,
              #         verticalalignment='top', horizontalalignment='left',
              #         transform=ax.transAxes,
              #         color='black', fontsize=10)
  
              #save the plot
              newFilename = MLE_directory+"Chi_profiles_basin_"+str(basin_key)+"_"+counter+".png"
  
              # This gets all the ticks, and pads them away from the axis so that the corners don't overlap
              ax.tick_params(axis='both', width=1, pad = 2)
              for tick in ax.xaxis.get_major_ticks():
                  tick.set_pad(2)
  
              FigFormat = "png"
              plt.savefig(newFilename,format=FigFormat,dpi=300)
              ax.cla()
  
  def MakeChiPlotsMLE(DataDirectory, fname_prefix, basin_list=[0], start_movern=0.2, d_movern=0.1, n_movern=7,
                      size_format='ESURF', FigFormat='png', animate=False, keep_pngs=False, parallel=False):
      """
      This function makes chi-elevation plots for each basin and each value of m/n
      where the channels are coloured by the MLE value compared to the main stem.
      The main stem is plotted in black.
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_list: a list of the basins to make the plots for. If an empty list is passed then
          all the basins will be analysed. Default = basin 0.
          start_movern (float): the starting m/n value. Default is 0.2
          d_movern (float): the increment between the m/n values. Default is 0.1
          n_movern (float): the number of m/n values analysed. Default is 7.
          size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
          FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
          animate (bool): If this is true then it creates a movie of the chi-elevation plots coloured by MLE.
          keep_pngs (bool): If this is false and the animation flag is true, then the pngs are deleted and just the video is kept.
  
      Returns:
          Plot of each m/n value for each basin.
  
      Author: FJC
      """
      from matplotlib.ticker import FormatStrFormatter
  
      # check if a directory exists for the chi plots. If not then make it.
      MLE_directory = DataDirectory+'chi_plots/'
      if not os.path.isdir(MLE_directory):
          os.makedirs(MLE_directory)
  
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # make a figure
      fig = makefigure(size_format)
  
      gs = plt.GridSpec(100,100,bottom=0.15,left=0.1,right=1.0,top=1.0)
      ax = fig.add_subplot(gs[10:95,5:80])
      #colorbar axis
      ax2 = fig.add_subplot(gs[10:95,82:85])
  
      # read in the csv files
      if not parallel:
          ProfileDF = Helper.ReadChiProfileCSV(DataDirectory, fname_prefix)
          BasinStatsDF = Helper.ReadBasinStatsCSV(DataDirectory, fname_prefix)
      else:
          ProfileDF = Helper.AppendMovernCSV(DataDirectory, fname_prefix)
          BasinStatsDF = Helper.AppendBasinCSVs(DataDirectory, fname_prefix)
  
      # get the number of basins
      basin_keys = list(BasinStatsDF['basin_key'])
      basin_keys = [int(x) for x in basin_keys]
  
      # get the list of basins
      if basin_list == []:
          print("You didn't give me a list of basins, so I'll just run the analysis on all of them!")
          basin_list = basin_keys
  
      # loop through each m over n value
      end_movern = float(start_movern)+float(d_movern)*(float(n_movern)-1)
      m_over_n_values = np.linspace(start_movern,end_movern,n_movern)
  
      # best fit moverns
      best_fit_moverns = SimpleMaxMLECheck(BasinStatsDF)
  
      for m_over_n in m_over_n_values:
          # read in the full stats file
  
          #Stupid floating point representation issues
          movern_str = "%.2f" % round(m_over_n,2)
          if movern_str.endswith('0'):
              movern_str = movern_str[:-1]
  
          print("This m/n is: "+movern_str)
          if not parallel:
              FullStatsDF = Helper.ReadFullStatsCSV(DataDirectory,fname_prefix,movern_str)
          else:
              FullStatsDF = Helper.AppendFullStatsCSVs(DataDirectory,movern_str,fname_prefix)
  
          # loop through all the basins in the basin list
          for basin_key in basin_list:
              print("This basin key is: "+str(basin_key))
  
              # mask the data frames for this basin
              ProfileDF_basin = ProfileDF[ProfileDF['basin_key'] == basin_key]
              FullStatsDF_basin = FullStatsDF[FullStatsDF['basin_key'] == basin_key]
  
              #print FullStatsDF_basin
              print ("Getting the reference_source_key")
  
              print(FullStatsDF_basin.iloc[0]['reference_source_key'])
  
              # get the data frame for the main stem
              ProfileDF_MS = ProfileDF_basin[ProfileDF_basin['source_key'] == FullStatsDF_basin.iloc[0]['reference_source_key']]
  
              # get the data frame for the tributaries
              ProfileDF_basin = ProfileDF_basin[ProfileDF_basin['source_key'] != FullStatsDF_basin.iloc[0]['reference_source_key']]
              # merge with the full data to get the MLE for the tributaries
              ProfileDF_tribs = ProfileDF_basin.merge(FullStatsDF_basin, left_on = "source_key", right_on = "test_source_key")
  
              # get the chi and elevation data for the main stem
              movern_key = 'm_over_n = %s' % movern_str
              MainStemX = list(ProfileDF_MS[movern_key])
              MainStemElevation = list(ProfileDF_MS['elevation'])
  
              # get the chi, elevation, and MLE for the tributaries
              TributariesX = list(ProfileDF_tribs[movern_key])
              TributariesElevation = list(ProfileDF_tribs['elevation'])
              TributariesMLE = list(ProfileDF_tribs['MLE'])
  
              # get the colourmap to colour channels by the MLE value
              #NUM_COLORS = len(MLE)
              MLE_array = np.asarray(TributariesMLE)
              this_cmap = plt.cm.coolwarm
              cNorm  = colors.Normalize(vmin=np.min(MLE_array), vmax=np.max(MLE_array))
              plt.cm.ScalarMappable(norm=cNorm, cmap=this_cmap)
  
              # now plot the data with a colourmap
              sc = ax.scatter(TributariesX,TributariesElevation,c=TributariesMLE,cmap=this_cmap, norm=cNorm, s=2.5, edgecolors='none')
              ax.plot(MainStemX,MainStemElevation,lw=2, c='k')
  
              # some formatting of the figure
              ax.spines['top'].set_linewidth(1)
              ax.spines['left'].set_linewidth(1)
              ax.spines['right'].set_linewidth(1)
              ax.spines['bottom'].set_linewidth(1)
  
              # make the lables
              ax.set_xlabel("$\chi$ (m)")
              ax.set_ylabel("Elevation (m)")
  
              # the best fit m/n
              best_fit_movern = best_fit_moverns[basin_key]
  
              # label with the basin and m/n
              title_string = "Basin "+str(basin_key)+", "+ r"$\theta$ = "+movern_str
              if best_fit_movern == m_over_n:
                  ax.text(0.05, 0.95, title_string,
                          verticalalignment='top', horizontalalignment='left',
                          transform=ax.transAxes,
                          color='red', fontsize=10)
              else:
                  ax.text(0.05, 0.95, title_string,
                          verticalalignment='top', horizontalalignment='left',
                          transform=ax.transAxes,
                          color='black', fontsize=10)
  
              # add the colorbar
              colorbarlabel = "$MLE$"
              cbar = plt.colorbar(sc,cmap=this_cmap,spacing='uniform', orientation='vertical',cax=ax2)
              cbar.set_label(colorbarlabel, fontsize=10)
              ax2.set_ylabel(colorbarlabel, fontname='Liberation Sans', fontsize=10)
              ax2.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
  
              #save the plot
              newFilename = MLE_directory+"MLE_profiles"+str(basin_key)+"_"+movern_str+"."+str(FigFormat)
  
              # This gets all the ticks, and pads them away from the axis so that the corners don't overlap
              ax.tick_params(axis='both', width=1, pad = 2)
              for tick in ax.xaxis.get_major_ticks():
                  tick.set_pad(2)
  
              plt.savefig(newFilename,format=FigFormat,dpi=300)
              ax.cla()
              ax2.cla()
  
      if animate:
          # animate the pngs using ffmpeg
          system_call = "ffmpeg -framerate 3 -pattern_type glob -i '"+MLE_directory+"MLE_profiles*.png' -y -vcodec libx264 -s 1230x566 -pix_fmt yuv420p "+MLE_directory+"MLE_profiles.mp4"
          print (system_call)
          subprocess.call(system_call, shell=True)
          # delete the pngs if you want
          if not keep_pngs:
              system_call = "rm "+MLE_directory+"MLE_profiles*.png"
              subprocess.call(system_call, shell=True)
      plt.close(fig)
  
      
      
  def MakeChiPlotsChi(DataDirectory, fname_prefix, basin_list=[0], start_movern=0.2, d_movern=0.1, n_movern=7,
                      size_format='ESURF', FigFormat='png', animate=False, keep_pngs=False, parallel=False):
      """
      This function makes chi-elevation plots for each basin and each value of m/n
      where the channels are coloured by the chi value compared to the main stem.
      The main stem is plotted in black.
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_list: a list of the basins to make the plots for. If an empty list is passed then
          all the basins will be analysed. Default = basin 0.
          start_movern (float): the starting m/n value. Default is 0.2
          d_movern (float): the increment between the m/n values. Default is 0.1
          n_movern (float): the number of m/n values analysed. Default is 7.
          size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
          FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
          animate (bool): If this is true then it creates a movie of the chi-elevation plots coloured by MLE.
          keep_pngs (bool): If this is false and the animation flag is true, then the pngs are deleted and just the video is kept.
  
      Returns:
          Plot of each m/n value for each basin.
  
      Author: SMM, from FJC code
      """
      from matplotlib.ticker import FormatStrFormatter
  
      # check if a directory exists for the chi plots. If not then make it.
      MLE_directory = DataDirectory+'chi_plots/'
      if not os.path.isdir(MLE_directory):
          os.makedirs(MLE_directory)
  
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # make a figure
      fig = makefigure(size_format)
  
      gs = plt.GridSpec(100,100,bottom=0.15,left=0.1,right=1.0,top=1.0)
      ax = fig.add_subplot(gs[10:95,5:90])
      
      #colorbar axis
      #ax2 = fig.add_subplot(gs[10:95,82:85])
  
      # read in the csv files
      if not parallel:
          ProfileDF = Helper.ReadChiProfileCSV(DataDirectory, fname_prefix)
          DisorderDF = Helper.ReadDisorderUncertCSV(DataDirectory, fname_prefix)
      else:
          ProfileDF = Helper.AppendMovernCSV(DataDirectory, fname_prefix)
          DisorderDF = Helper.AppendDisorderCSV(DataDirectory, fname_prefix)
      best_fit_moverns = DisorderDF['median'].tolist()
      
      # get the number of basins
      basin_keys = list(DisorderDF['basin_key'])
      basin_keys = [int(x) for x in basin_keys]
      
      MOverNDict = dict(zip(basin_keys,best_fit_moverns))
  
      # get the list of basins
      if basin_list == []:
          print("You didn't give me a list of basins, so I'll just run the analysis on all of them!")
          basin_list = basin_keys
  
      # loop through each m over n value
      end_movern = float(start_movern)+float(d_movern)*(float(n_movern)-1)
      m_over_n_values = np.linspace(start_movern,end_movern,n_movern)
  
      for m_over_n in m_over_n_values:
          # read in the full stats file
  
          #Stupid floating point representation issues
          movern_str = "%.2f" % round(m_over_n,2)
          if movern_str.endswith('0'):
              movern_str = movern_str[:-1]
  
          print("This concavity is: "+movern_str)
  
          # loop through all the basins in the basin list
          for basin_key in basin_list:
              print("This basin key is: "+str(basin_key))
              #print("The best fit concavity is:")
              #print(MOverNDict[basin_key])
              
              # Format the best fit concavity string
              bf_movernstr = "%.2f" % round(MOverNDict[basin_key],2)
              if bf_movernstr.endswith('0'):
                  bf_movernstr = bf_movernstr[:-1]
                  
              #print("This concavity is:")
              #print(movern_str)
              #print("Best fit concavity in string format is:")
              #print(bf_movernstr)
              
              this_is_bf_concavity = False
              if (movern_str == bf_movernstr):
                  #print("==================================")
                  #print("This is the best fitting concavity")
                  #print("==================================")
                  this_is_bf_concavity = True
                  
              # mask the data frames for this basin
              ProfileDF_basin = ProfileDF[ProfileDF['basin_key'] == basin_key]
  
              # get the chi and elevation data for the main stem
              movern_key = 'm_over_n = %s' % movern_str
              X = list(ProfileDF_basin[movern_key])
              Elevation = list(ProfileDF_basin['elevation'])
  
              # get the colourmap to colour channels by the MLE value
              #NUM_COLORS = len(MLE)
              MLE_array = np.asarray(X)
              this_cmap = plt.cm.coolwarm
              cNorm  = colors.Normalize(vmin=np.min(X), vmax=np.max(X))
              plt.cm.ScalarMappable(norm=cNorm, cmap=this_cmap)
  
              # now plot the data with a colourmap
              if this_is_bf_concavity:
                  sc = ax.scatter(X,Elevation,c='r', s=2.5, edgecolors='none')
                  #sc = ax.scatter(X,Elevation,c=X,cmap=this_cmap, norm=cNorm, s=2.5, edgecolors='none')
              else:
                  sc = ax.scatter(X,Elevation,c='k', s=2.5, edgecolors='none')
                  #sc = ax.scatter(X,Elevation,c=X,cmap=this_cmap, norm=cNorm, s=2.5, edgecolors='none')
                  
              
  
              # some formatting of the figure
              ax.spines['top'].set_linewidth(1)
              ax.spines['left'].set_linewidth(1)
              ax.spines['right'].set_linewidth(1)
              ax.spines['bottom'].set_linewidth(1)
  
              # make the lables
              ax.set_xlabel("$\chi$ (m)")
              ax.set_ylabel("Elevation (m)")
  
              # the best fit m/n
              best_fit_movern = best_fit_moverns[basin_key]
              
              #print("The best fit concavity is: "+str(best_fit_movern))
  
              # label with the basin and m/n
              title_string = "Basin "+str(basin_key)+", "+ r"$\theta$ = "+movern_str
              if this_is_bf_concavity:
                  ax.text(0.05, 0.95, title_string,
                          verticalalignment='top', horizontalalignment='left',
                          transform=ax.transAxes,
                          color='red', fontsize=10)
                  ax.text(0.05, 0.85, "Best fit concavity",
                          verticalalignment='top', horizontalalignment='left',
                          transform=ax.transAxes,
                          color='red', fontsize=10)
              else:
                  ax.text(0.05, 0.95, title_string,
                          verticalalignment='top', horizontalalignment='left',
                          transform=ax.transAxes,
                          color='black', fontsize=10)
  
              # add the colorbar
              #colorbarlabel = "$\chi$ (m)"
              #cbar = plt.colorbar(sc,cmap=this_cmap,spacing='uniform', orientation='vertical',cax=ax2)
              #cbar.set_label(colorbarlabel, fontsize=10)
              #ax2.set_ylabel(colorbarlabel, fontname='Liberation Sans', fontsize=10)
              #ax2.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
  
              #save the plot
              newFilename = MLE_directory+"chi_profiles"+str(basin_key)+"_"+movern_str+"."+str(FigFormat)
  
              # This gets all the ticks, and pads them away from the axis so that the corners don't overlap
              ax.tick_params(axis='both', width=1, pad = 2)
              for tick in ax.xaxis.get_major_ticks():
                  tick.set_pad(2)
  
              plt.savefig(newFilename,format=FigFormat,dpi=300)
              ax.cla()
              #ax2.cla()
  
      if animate:
          # animate the pngs using ffmpeg
          system_call = "ffmpeg -framerate 3 -pattern_type glob -i '"+MLE_directory+"chi_profiles*.png' -y -vcodec libx264 -s 1230x566 -pix_fmt yuv420p "+MLE_directory+"chi_profiles.mp4"
          print (system_call)
          subprocess.call(system_call, shell=True)
          # delete the pngs if you want
          if not keep_pngs:
              system_call = "rm "+MLE_directory+"chi_profiles*.png"
              subprocess.call(system_call, shell=True)
      plt.close(fig)    
      
      
  def MakeChiPlotsColouredByK(DataDirectory, fname_prefix, basin_list=[0], start_movern=0.2, d_movern=0.1, n_movern=7,size_format='ESURF', FigFormat='png', animate=False, keep_pngs=False, parallel=False):
      """
      This function makes chi-elevation plots for each basin and each value of m/n
      where the channels are coloured by the K value (for model runs with spatially varying K).
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_list: a list of the basins to make the plots for. If an empty list is passed then
          all the basins will be analysed. Default = basin 0.
          start_movern (float): the starting m/n value. Default is 0.2
          d_movern (float): the increment between the m/n values. Default is 0.1
          n_movern (float): the number of m/n values analysed. Default is 7.
          size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
          FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
          animate (bool): If this is true then it creates a movie of the chi-elevation plots coloured by MLE.
          keep_pngs (bool): If this is false and the animation flag is true, then the pngs are deleted and just the video is kept.
  
      Returns:
          Plot of each m/n value for each basin.
  
      Author: FJC
      """
      from LSDPlottingTools import colours
  
      # check if a directory exists for the chi plots. If not then make it.
      K_directory = DataDirectory+'chi_plots_K/'
      if not os.path.isdir(K_directory):
          os.makedirs(K_directory)
  
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # make a figure
      fig = makefigure(size_format)
  
  
      gs = plt.GridSpec(100,100,bottom=0.15,left=0.1,right=0.95,top=1.0)
      ax = fig.add_subplot(gs[10:95,5:80])
      #colorbar axis
      ax2 = fig.add_subplot(gs[10:95,82:85])
  
      # read in the csv files
      if not parallel:
          ProfileDF = Helper.ReadChiProfileCSV(DataDirectory, fname_prefix)
          BasinStatsDF = Helper.ReadBasinStatsCSV(DataDirectory, fname_prefix)
      else:
          ProfileDF = Helper.AppendMovernCSV(DataDirectory,fname_prefix)
          BasinStatsDF = Helper.AppendBasinCSVs(DataDirectory,fname_prefix)
  
      # get the number of basins
      basin_keys = list(BasinStatsDF['basin_key'])
      basin_keys = [int(x) for x in basin_keys]
  
      # get the list of basins
      if basin_list == []:
          print("You didn't give me a list of basins, so I'll just run the analysis on all of them!")
          basin_list = basin_keys
  
      # loop through each m over n value
      end_movern = start_movern+d_movern*(n_movern-1)
      m_over_n_values = np.linspace(start_movern,end_movern,n_movern)
  
      # best fit moverns
      best_fit_moverns = SimpleMaxMLECheck(BasinStatsDF)
  
      for m_over_n in m_over_n_values:
          # read in the full stats file
          print("This m/n is: "+str(m_over_n))
          if not parallel:
              FullStatsDF = Helper.ReadFullStatsCSV(DataDirectory,fname_prefix,m_over_n)
          else:
              FullStatsDF = Helper.AppendFullStatsCSVs(DataDirectory,m_over_n,fname_prefix)
  
          # loop through all the basins in the basin list
          for basin_key in basin_list:
              print("This basin key is %s") %str(basin_key)
  
              # mask the data frames for this basin
              ProfileDF_basin = ProfileDF[ProfileDF['basin_key'] == basin_key]
              FullStatsDF_basin = FullStatsDF[FullStatsDF['basin_key'] == basin_key]
  
              #print FullStatsDF_basin
              print ("Getting the reference_source_key")
  
              print(FullStatsDF_basin.iloc[0]['reference_source_key'])
  
              # get the data frame for the main stem
              ProfileDF_MS = ProfileDF_basin[ProfileDF_basin['source_key'] == FullStatsDF_basin.iloc[0]['reference_source_key']]
  
              # get the data frame for the tributaries
              ProfileDF_basin = ProfileDF_basin[ProfileDF_basin['source_key'] != FullStatsDF_basin.iloc[0]['reference_source_key']]
              # merge with the full data to get the MLE for the tributaries
              ProfileDF_tribs = ProfileDF_basin.merge(FullStatsDF_basin, left_on = "source_key", right_on = "test_source_key")
  
              # get the chi and elevation data for the main stem
              movern_key = 'm_over_n = %s' %(str(m_over_n))
              MainStemX = list(ProfileDF_MS[movern_key])
              MainStemElevation = list(ProfileDF_MS['elevation'])
              MainStemK = list(ProfileDF_MS['K_value'])
  
              # get the chi, elevation, and MLE for the tributaries
              TributariesX = list(ProfileDF_tribs[movern_key])
              TributariesElevation = list(ProfileDF_tribs['elevation'])
              TributariesK = list(ProfileDF_tribs['K_value'])
  
              # get the colourmap to colour channels by the MLE value
              #NUM_COLORS = len(MLE)
              K_array = np.asarray(TributariesK)
              min_K = np.min(K_array)
              max_K = np.max(K_array)
              this_cmap = plt.cm.Spectral
              n_colours = 10
              this_cmap = colours.cmap_discretize(n_colours, this_cmap)
              cNorm  = colors.Normalize(vmin=min_K, vmax=max_K)
              plt.cm.ScalarMappable(norm=cNorm, cmap=this_cmap)
  
              # now plot the data with a colourmap
              sc = ax.scatter(TributariesX,TributariesElevation,c=TributariesK,cmap=this_cmap, norm=cNorm, s=2.5, edgecolors='none')
              sc = ax.scatter(MainStemX, MainStemElevation,c=MainStemK,cmap=this_cmap, norm=cNorm, s=2.5, edgecolors='none')
  
              # some formatting of the figure
              ax.spines['top'].set_linewidth(1)
              ax.spines['left'].set_linewidth(1)
              ax.spines['right'].set_linewidth(1)
              ax.spines['bottom'].set_linewidth(1)
  
              # make the labels
              ax.set_xlabel("$\chi$ (m)")
              ax.set_ylabel("Elevation (m)")
  
              # the best fit m/n
              best_fit_movern = best_fit_moverns[basin_key]
              print ("BEST FIT M/N IS: "+ str(best_fit_movern))
              print ("THIS M/N IS: "+str(m_over_n))
  
              # label with the basin and m/n
              title_string = "Basin "+str(basin_key)+", $m/n$ = "+str(m_over_n)
              if best_fit_movern == m_over_n:
                  ax.text(0.05, 0.95, title_string,
                          verticalalignment='top', horizontalalignment='left',
                          transform=ax.transAxes,
                          color='red', fontsize=10)
              else:
                  ax.text(0.05, 0.95, title_string,
                          verticalalignment='top', horizontalalignment='left',
                          transform=ax.transAxes,
                          color='black', fontsize=10)
  
              # add the colorbar
              colorbarlabel = "$K$"
              cbar = plt.colorbar(sc,cmap=this_cmap,spacing='uniform', orientation='vertical',cax=ax2)
              cbar.set_label(colorbarlabel, fontsize=10)
              ax2.set_ylabel(colorbarlabel, fontname='Liberation Sans', fontsize=10)
  
              #change labels to scientific notation
              colours.fix_colourbar_ticks(cbar,n_colours, cbar_type=float, min_value = min_K, max_value = max_K, cbar_label_rotation=0, cbar_orientation='vertical')
              # we need to get linear values between min and max K
              these_labels = np.linspace(min_K,max_K,n_colours)
              # now round these and convert to scientific notation
              these_labels = [str('{:.2e}'.format(float(x))) for x in these_labels]
              new_labels = []
              for label in these_labels:
                  a,b = label.split("e")
                  b = b.replace("0", "")
                  new_labels.append(a+' x 10$^{%s}$' % b)
  
              ax2.set_yticklabels(new_labels, fontsize=8)
  
              #save the plot
              newFilename = K_directory+"Chi_profiles_by_K_"+str(basin_key)+"_"+str(m_over_n)+"."+str(FigFormat)
  
              # This gets all the ticks, and pads them away from the axis so that the corners don't overlap
              ax.tick_params(axis='both', width=1, pad = 2)
              for tick in ax.xaxis.get_major_ticks():
                  tick.set_pad(2)
  
              plt.savefig(newFilename,format=FigFormat,dpi=300)
              ax.cla()
              ax2.cla()
  
      if animate:
          # animate the pngs using ffmpeg
          system_call = "ffmpeg -framerate 3 -pattern_type glob -i '"+K_directory+"Chi_profiles_by_K*.png' -y -vcodec libx264 -s 1230x566 -pix_fmt yuv420p "+K_directory+"Chi_profiles_by_K.mp4"
          print (system_call)
          subprocess.call(system_call, shell=True)
          # delete the pngs if you want
          if not keep_pngs:
              system_call = "rm "+K_directory+"Chi_profiles_by_K*.png"
              subprocess.call(system_call, shell=True)
      plt.close(fig)
  
  def MakeChiPlotsColouredByLith(DataDirectory, fname_prefix, basin_list=[0], start_movern=0.2, d_movern=0.1, n_movern=7,
                      size_format='ESURF', FigFormat='png', animate=False, keep_pngs=False, parallel=False):
      """
      This function makes chi-elevation plots for each basin and each value of m/n
      where the channels are coloured by the K value (for model runs with spatially varying K).
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_list: a list of the basins to make the plots for. If an empty list is passed then
          all the basins will be analysed. Default = basin 0.
          start_movern (float): the starting m/n value. Default is 0.2
          d_movern (float): the increment between the m/n values. Default is 0.1
          n_movern (float): the number of m/n values analysed. Default is 7.
          size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
          FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
          animate (bool): If this is true then it creates a movie of the chi-elevation plots coloured by MLE.
          keep_pngs (bool): If this is false and the animation flag is true, then the pngs are deleted and just the video is kept.
  
      Returns:
          Plot of each m/n value for each basin.
  
      Author: FJC
      """
      from LSDPlottingTools import colours
      print("WARNING DEPRECATED FUNCTION, USE THE MakeChiPlotsByLith FROM LSDMapLithoPlotting")
  
      # check if a directory exists for the chi plots. If not then make it.
      K_directory = DataDirectory+'chi_plots_Lith/'
      if not os.path.isdir(K_directory):
          os.makedirs(K_directory)
  
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # make a figure
      fig = makefigure(size_format)
  
  
      gs = plt.GridSpec(100,100,bottom=0.15,left=0.1,right=0.95,top=1.0)
      ax = fig.add_subplot(gs[10:95,5:80])
      #colorbar axis
      ax2 = fig.add_subplot(gs[10:95,82:85])
  
      # read in the csv files
      if not parallel:
          ProfileDF = Helper.ReadChiProfileCSV(DataDirectory, fname_prefix)
          BasinStatsDF = Helper.ReadBasinStatsCSV(DataDirectory, fname_prefix)
      else:
          ProfileDF = Helper.AppendMovernCSV(DataDirectory,fname_prefix)
          BasinStatsDF = Helper.AppendBasinCSVs(DataDirectory,fname_prefix)
  
      # get the number of basins
      basin_keys = list(BasinStatsDF['basin_key'])
      basin_keys = [int(x) for x in basin_keys]
  
      # get the list of basins
      if basin_list == []:
          print("You didn't give me a list of basins, so I'll just run the analysis on all of them!")
          basin_list = basin_keys
  
      # loop through each m over n value
      end_movern = start_movern+d_movern*(n_movern-1)
      m_over_n_values = np.linspace(start_movern,end_movern,n_movern)
  
      # best fit moverns
      best_fit_moverns = SimpleMaxMLECheck(BasinStatsDF)
  
      for m_over_n in m_over_n_values:
          # read in the full stats file
          print("This m/n is: "+str(m_over_n))
          if not parallel:
              FullStatsDF = Helper.ReadFullStatsCSV(DataDirectory,fname_prefix,m_over_n)
          else:
              FullStatsDF = Helper.AppendFullStatsCSVs(DataDirectory,m_over_n,fname_prefix)
  
          # loop through all the basins in the basin list
          for basin_key in basin_list:
              print("This basin key is %s") %str(basin_key)
  
              # mask the data frames for this basin
              ProfileDF_basin = ProfileDF[ProfileDF['basin_key'] == basin_key]
              FullStatsDF_basin = FullStatsDF[FullStatsDF['basin_key'] == basin_key]
  
              #print FullStatsDF_basin
              print ("Getting the reference_source_key")
  
              print(FullStatsDF_basin.iloc[0]['reference_source_key'])
  
              # get the data frame for the main stem
              ProfileDF_MS = ProfileDF_basin[ProfileDF_basin['source_key'] == FullStatsDF_basin.iloc[0]['reference_source_key']]
  
              # get the data frame for the tributaries
              ProfileDF_basin = ProfileDF_basin[ProfileDF_basin['source_key'] != FullStatsDF_basin.iloc[0]['reference_source_key']]
              # merge with the full data to get the MLE for the tributaries
              ProfileDF_tribs = ProfileDF_basin.merge(FullStatsDF_basin, left_on = "source_key", right_on = "test_source_key")
  
              # get the chi and elevation data for the main stem
              movern_key = 'm_over_n = %s' %(str(m_over_n))
              MainStemX = list(ProfileDF_MS[movern_key])
              MainStemElevation = list(ProfileDF_MS['elevation'])
              MainStemK = list(ProfileDF_MS[fname_prefix+"_geol"])
  
  
              # get the chi, elevation, and MLE for the tributaries
              TributariesX = list(ProfileDF_tribs[movern_key])
              TributariesElevation = list(ProfileDF_tribs['elevation'])
              TributariesK = list(ProfileDF_tribs[fname_prefix+"_geol"])
  
  
              # get the colourmap to colour channels by the MLE value
              #NUM_COLORS = len(MLE)
              K_array = np.asarray(TributariesK)
              min_K = np.min(K_array)
              max_K = np.max(K_array)
              this_cmap = plt.cm.Spectral
              n_colours = 10
              this_cmap = colours.cmap_discretize(n_colours, this_cmap)
              cNorm  = colors.Normalize(vmin=min_K, vmax=max_K)
              plt.cm.ScalarMappable(norm=cNorm, cmap=this_cmap)
  
              # now plot the data with a colourmap
              sc = ax.scatter(TributariesX,TributariesElevation,c=TributariesK,cmap=this_cmap, norm=cNorm, s=2.5, edgecolors='none')
              sc = ax.scatter(MainStemX, MainStemElevation,c=MainStemK,cmap=this_cmap, norm=cNorm, s=2.5, edgecolors='none')
  
              # some formatting of the figure
              ax.spines['top'].set_linewidth(1)
              ax.spines['left'].set_linewidth(1)
              ax.spines['right'].set_linewidth(1)
              ax.spines['bottom'].set_linewidth(1)
  
              # make the labels
              ax.set_xlabel("$\chi$ (m)")
              ax.set_ylabel("Elevation (m)")
  
              # the best fit m/n
              best_fit_movern = best_fit_moverns[basin_key]
              print ("BEST FIT M/N IS: "+ str(best_fit_movern))
              print ("THIS M/N IS: "+str(m_over_n))
  
              # label with the basin and m/n
              title_string = "Basin "+str(basin_key)+", $m/n$ = "+str(m_over_n)
              if best_fit_movern == m_over_n:
                  ax.text(0.05, 0.95, title_string,
                          verticalalignment='top', horizontalalignment='left',
                          transform=ax.transAxes,
                          color='red', fontsize=10)
              else:
                  ax.text(0.05, 0.95, title_string,
                          verticalalignment='top', horizontalalignment='left',
                          transform=ax.transAxes,
                          color='black', fontsize=10)
  
              # add the colorbar
              colorbarlabel = "$Lith$"
              cbar = plt.colorbar(sc,cmap=this_cmap,spacing='uniform', orientation='vertical',cax=ax2)
              cbar.set_label(colorbarlabel, fontsize=10)
              ax2.set_ylabel(colorbarlabel, fontname='Liberation Sans', fontsize=10)
  
              #change labels to scientific notation
              colours.fix_colourbar_ticks(cbar,n_colours, cbar_type=float, min_value = min_K, max_value = max_K, cbar_label_rotation=0, cbar_orientation='vertical')
              # we need to get linear values between min and max K
              these_labels = np.linspace(min_K,max_K,n_colours)
              # now round these and convert to scientific notation
              these_labels = [str('{:.2e}'.format(float(x))) for x in these_labels]
              new_labels = []
              for label in these_labels:
                  a,b = label.split("e")
                  b = b.replace("0", "")
                  new_labels.append(a+' x 10$^{%s}$' % b)
  
              ax2.set_yticklabels(new_labels, fontsize=8)
  
              #save the plot
              newFilename = K_directory+"Chi_profiles_by_Lith_"+str(basin_key)+"_"+str(m_over_n)+"."+str(FigFormat)
  
              # This gets all the ticks, and pads them away from the axis so that the corners don't overlap
              ax.tick_params(axis='both', width=1, pad = 2)
              for tick in ax.xaxis.get_major_ticks():
                  tick.set_pad(2)
  
              plt.savefig(newFilename,format=FigFormat,dpi=300)
              ax.cla()
              ax2.cla()
  
      if animate:
          # animate the pngs using ffmpeg
          system_call = "ffmpeg -framerate 3 -pattern_type glob -i '"+K_directory+"Chi_profiles_by_Lith*.png' -y -vcodec libx264 -s 1230x566 -pix_fmt yuv420p "+K_directory+"Chi_profiles_by_Lith.mp4"
          print (system_call)
          subprocess.call(system_call, shell=True)
          # delete the pngs if you want
          if not keep_pngs:
              system_call = "rm "+K_directory+"Chi_profiles_by_Lith*.png"
              subprocess.call(system_call, shell=True)
      plt.close(fig)
  
  
  def PlotProfilesRemovingOutliers(DataDirectory, fname_prefix, basin_list=[0], start_movern=0.2, d_movern=0.1, n_movern=7, size_format = "geomorphology", FigFormat="png", parallel=False):
      """
      This function is used to plot the chi profiles as they have outliers removed.
      It calls thefunction CheckMLEOutliers, which you should read to get details
      on how outliers are calculated and removed
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_list: a list of the basins to make the plots for. If an empty list is passed then
          all the basins will be analysed. Default = basin 0.
          start_movern (float): the starting m/n value. Default is 0.2
          d_movern (float): the increment between the m/n values. Default is 0.1
          n_movern (float): the number of m/n values analysed. Default is 7.
  
      Returns:
          Plots of chi profiles with basins removed
  
      Author: SMM
      """
  
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # make a figure
      fig = makefigure(size_format)
  
      gs = plt.GridSpec(100,100,bottom=0.15,left=0.1,right=1.0,top=1.0)
      ax = fig.add_subplot(gs[10:95,5:80])
      #colorbar axis
      ax2 = fig.add_subplot(gs[10:95,82:85])
  
      # we open the first file just so that we can get a counter list
      if not parallel:
          FirstDF = Helper.ReadFullStatsCSV(DataDirectory,fname_prefix,start_movern)
      else:
          FirstDF = Helper.AppendFullStatsCSVs(DataDirectory,start_movern,fname_prefix)
  
      # get the number of basins
      basin_keys = list(FirstDF['basin_key'])
      basin_keys = [float(x) for x in basin_keys]
  
      # get the list of basins
      if basin_list == []:
          print("You didn't give me a list of basins, so I'll just run the analysis on all of them!")
          basin_list = basin_keys
          basin_set = set(basin_list)
          basin_list = list(basin_set)
          basin_list = [int(i) for i in basin_list]
  
      # First we get all the information about outliers, m/n values and MLE
      # values from the CheckMLEOutliers function
      Outlier_counter, removed_sources_dict, best_fit_movern_dict, MLEs_dict = CheckMLEOutliers(DataDirectory, fname_prefix, basin_list, start_movern, d_movern, n_movern, parallel=parallel)
  
      # Now get the chi profiles of all the basins and channels
      # Load from file and put into a pandas data frame
      if not parallel:
          ProfileDF = Helper.ReadChiProfileCSV(DataDirectory,fname_prefix)
      else:
          ProfileDF = Helper.AppendMovernCSV(DataDirectory,fname_prefix)
  
      # now we need to get a plot for each basin, showing the incremental removal of outlying tribs
      for basin_number in basin_list:
  
          # Get the removed sources indices and the MLEs for this particular basin
          these_removed_sources = removed_sources_dict[basin_number]
          these_MLEs = MLEs_dict[basin_number]
  
          # loop through the best fit moverns
          # each value represents the MLE for a given number of removed outlying tributaries
          removed_sources_list = []
          for idx,best_fit_movern in enumerate(best_fit_movern_dict[basin_number]):
  
              print("The best fit m/n is: "+ str(best_fit_movern)+" and the index is "+str(idx))
  
              # mask the data frames for this basin
              ProfileDF_basin = ProfileDF[ProfileDF['basin_key'] == basin_number]
  
              # We also need to get the fullstats MLE file. This will be used
              # to colour tribs as well as get the source numbers from
              # the outlier list
              if not parallel:
                  FullStatsDF = Helper.ReadFullStatsCSV(DataDirectory,fname_prefix,best_fit_movern)
              else:
                  FullStatsDF = Helper.AppendFullStatsCSVs(DataDirectory,best_fit_movern,fname_prefix)
  
              # mask the data so you only get the correct basin
              FullStatsDF_basin = FullStatsDF[FullStatsDF['basin_key'] == basin_number]
  
              # extract the relevant data
              #MLE_values = list(FullStatsDF_basin['MLE'])
              #RMSE_values = list(FullStatsDF_basin['RMSE'])
              trib_values = list(FullStatsDF_basin['test_source_key'])
              ref_values = list(FullStatsDF_basin['reference_source_key'])
  
              # Now get the excluded tributaries for this iteration
              #removed_MLE = these_MLEs[idx]
  
              # Note that index 0 is the basin with no removed tributaries.
              if idx != 0:
                  removed_sources_list.extend(these_removed_sources[idx-1])
  
              # now you need to get the actual source numbers by indexing into the source list
              the_removed_sources = []
              for source_index in removed_sources_list:
                  the_removed_sources.append( trib_values[source_index]  )
  
              print("The main stem is: ")
              print( ref_values[0])
  
              print("The removed tribs are: ")
              print(the_removed_sources)
  
  
  
              # get the data frame for the main stem
              # It does this because the main stem source is always the 0 element in the trib_values list
              ProfileDF_MS = ProfileDF_basin[ProfileDF_basin['source_key'] == ref_values[0]]
  
              # get the data frame for the tributaries
              ProfileDF_basin = ProfileDF_basin[ProfileDF_basin['source_key'] != ref_values[0]]
  
              # now split the tributaries into exluded and non excluded tribs
              #ProfileDF_outliers = ProfileDF_basin.filter(items=removed_sources_list)
              #ProfileDF_kept = ProfileDF_basin.filter(items=removed_sources_list)
              ProfileDF_outliers = ProfileDF_basin[ProfileDF_basin.source_key.isin(the_removed_sources)]
              ProfileDF_kept = ProfileDF_basin[~ProfileDF_basin.source_key.isin(the_removed_sources)]
  
              #ProfileDF_basin = ProfileDF_basin[ProfileDF_basin.source_key.isin(removed_sources_list)]
              #ProfileDF_basin = ProfileDF_basin[ProfileDF_basin.source_key.isin([2,3,4])]
              #print(ProfileDF_basin)
  
              # merge with the full data to get the MLE for the tributaries
              ProfileDF_trib_outliers = ProfileDF_outliers.merge(FullStatsDF_basin, left_on = "source_key", right_on = "test_source_key")
              ProfileDF_trib_kept = ProfileDF_kept.merge(FullStatsDF_basin, left_on = "source_key", right_on = "test_source_key")
  
              # get the chi and elevation data for the main stem
              movern_key = 'm_over_n = %s' %(str(best_fit_movern))
              MainStemX = list(ProfileDF_MS[movern_key])
              MainStemElevation = list(ProfileDF_MS['elevation'])
  
              # get the chi, elevation, and MLE for the tributaries
              TributariesX_outliers = list(ProfileDF_trib_outliers[movern_key])
              TributariesElevation_outliers = list(ProfileDF_trib_outliers['elevation'])
              #TributariesMLE_outliers = list(ProfileDF_trib_outliers['MLE'])
  
              TributariesX_kept = list(ProfileDF_trib_kept[movern_key])
              TributariesElevation_kept = list(ProfileDF_trib_kept['elevation'])
              TributariesMLE_kept = list(ProfileDF_trib_kept['MLE'])
  
              # now reset the
  
              # get the colourmap to colour channels by the MLE value
              #NUM_COLORS = len(MLE)
              MLE_array = np.asarray(TributariesMLE_kept)
              this_cmap = plt.cm.Reds
              cNorm  = colors.Normalize(vmin=np.min(MLE_array), vmax=np.max(MLE_array))
              plt.cm.ScalarMappable(norm=cNorm, cmap=this_cmap)
  
              # now plot the data with a colourmap
              sc = ax.scatter(TributariesX_kept,TributariesElevation_kept,c=TributariesMLE_kept,cmap=this_cmap, norm=cNorm, s=2.5, edgecolors='none')
  
              # Add the outliers if the basin has them
              if(len(removed_sources_list)>0):
                  ax.scatter(TributariesX_outliers,TributariesElevation_outliers,c="b", norm=cNorm, s=2.5, edgecolors='none', alpha = 0.3)
  
              ax.plot(MainStemX,MainStemElevation,lw=2, c='k')
  
              # some formatting of the figure
              ax.spines['top'].set_linewidth(1)
              ax.spines['left'].set_linewidth(1)
              ax.spines['right'].set_linewidth(1)
              ax.spines['bottom'].set_linewidth(1)
  
              # make the lables
              ax.set_xlabel("$\chi$ (m)")
              ax.set_ylabel("Elevation (m)")
  
              # label with the basin and m/n
              title_string = "Basin "+str(basin_number)+", best fit "+r'$\theta$' + " = str(best_fit_movern)"
              ax.text(0.05, 0.95, title_string,
                      verticalalignment='top', horizontalalignment='left',
                      transform=ax.transAxes,
                      color='black', fontsize=10)
  
              # add the colorbar
              colorbarlabel = "$MLE$"
              cbar = plt.colorbar(sc,cmap=this_cmap,spacing='uniform', orientation='vertical',cax=ax2)
              cbar.set_label(colorbarlabel, fontsize=10)
              ax2.set_ylabel(colorbarlabel, fontname='Liberation Sans', fontsize=10)
  
              #save the plot
              newFilename = DataDirectory+"MLE_profiles"+str(basin_number)+"_"+str(best_fit_movern)+"_removed_"+str(idx)+".png"
  
              # This gets all the ticks, and pads them away from the axis so that the corners don't overlap
              ax.tick_params(axis='both', width=1, pad = 2)
              for tick in ax.xaxis.get_major_ticks():
                  tick.set_pad(2)
  
              plt.savefig(newFilename,format=FigFormat,dpi=300)
              ax.cla()
              ax2.cla()
  
  def PlotMLEWithMOverN(DataDirectory, fname_prefix, basin_list = [0], size_format='ESURF', FigFormat='png', start_movern=0.2, d_movern = 0.1, n_movern = 7, parallel=False):
      """
      This function makes a plot of the MLE values for each m/n showing how the MLE values change
      as you remove the tributaries.
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_list: a list of the basins to make the plots for. If an empty list is passed then
          all the basins will be analysed. Default = basin 0.
          size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
          FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
          start_movern (float): the starting m/n value. Default is 0.2
          d_movern (float): the increment between the m/n values. Default is 0.1
          n_movern (float): the number of m/n values analysed. Default is 7.
  
      Returns:
          Plots of MLE values for each m/n
  
      Author: FJC
      """
      from cycler import cycler
      from matplotlib import lines
      import matplotlib.patches as patches
  
      # check if a directory exists for the MLE plots. If not then make it.
      MLE_directory = DataDirectory+'MLE_plots/'
      if not os.path.isdir(MLE_directory):
          os.makedirs(MLE_directory)
  
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # set colours for tributaries
      # plt.rc('axes', prop_cycle=(cycler('color', ['k', '0.5', '0.5', '0.5', '0.5']) +
      #                        cycler('linestyle', ['-', '--', ':', '-.', '--'])))
  
      # make a figure
      fig = makefigure(size_format)
  
  
      gs = plt.GridSpec(100,100,bottom=0.15,left=0.1,right=0.85,top=0.9)
      ax = fig.add_subplot(gs[5:100,10:95])
  
      # we open the first file just so that we can get a counter list
      if not parallel:
          FirstDF = Helper.ReadFullStatsCSV(DataDirectory,fname_prefix,start_movern)
      else:
          FirstDF = Helper.AppendFullStatsCSVs(DataDirectory,start_movern,fname_prefix)
  
      # get the list of m over n values
      end_movern = start_movern+d_movern*(n_movern-1)
      m_over_n_values = np.linspace(start_movern,end_movern,n_movern)
  
      # get the number of basins
      basin_keys = list(FirstDF['basin_key'])
      basin_keys = [float(x) for x in basin_keys]
  
      # get the list of basins
      if basin_list == []:
          print("You didn't give me a list of basins, so I'll just run the analysis on all of them!")
          basin_list = basin_keys
          basin_set = set(basin_list)
          basin_list = list(basin_set)
          basin_list = [int(i) for i in basin_list]
  
      # First we get all the information about outliers, m/n values and MLE
      # values from the CheckMLEOutliers function
      Outlier_counter, removed_sources_dict, best_fit_movern_dict, MLEs_dict = CheckMLEOutliers(DataDirectory, fname_prefix, basin_list, start_movern, d_movern, n_movern, parallel=parallel)
  
      # Now get the chi profiles of all the basins and channels
      # Load from file and put into a pandas data frame
      if not parallel:
          ProfileDF = Helper.ReadChiProfileCSV(DataDirectory,fname_prefix)
      else:
          ProfileDF = Helper.AppendMovernCSV(DataDirectory,fname_prefix)
  
      # get list of line styles for plotting. This is hacky but not sure how else to do this.
      ls = lines.lineStyles.keys()
      ls = ls[3:]
      ls = ls + ls
  
      # loop through each basin and each number of removed tributaries
      for basin_number in basin_list:
  
          print ("This basin is: " +str(basin_number))
  
          # Get the removed sources indices and the MLEs for this particular basin
          basin_removed_sources = removed_sources_dict[basin_number]
          n_removed_sources = len(basin_removed_sources)
          basin_MLEs = MLEs_dict[basin_number]
  
          # get the best fit m over ns for this basin
          best_fit_moverns = best_fit_movern_dict[basin_number]
          print (best_fit_moverns)
  
          # colours for each iteration
          #colours = np.linspace(0.9,0.2,len(n_removed_sources)+1)
  
  
          # loop through the number of removed tributaires and get the MLE for each m/n for each iteration
          for i in range(n_removed_sources+1):
  
              # get the MLE of the best fit m over n
              best_fit_movern = best_fit_moverns[i]
              # get the index in the MLE list
              idx = int(round((best_fit_movern - start_movern)/d_movern,0))
              print (idx)
              best_fit_MLE = basin_MLEs[idx][i]
              print ("The best fit MLE is: "+str(best_fit_MLE)+", where "+ r'$\theta$'  ' = '+str(best_fit_movern))
  
              # get the MLEs for this iteration
              these_MLEs = basin_MLEs[:,i]
  
              # get the ratio of these MLEs to the best fit
              ratio_MLEs = [x/best_fit_MLE for x in these_MLEs]
  
              # no removed tributaries
              if i == 0:
                  # plot the data
                  ax.scatter(m_over_n_values,ratio_MLEs, label = str(i), c='k', s=5, zorder=100)
                  ax.plot(m_over_n_values,ratio_MLEs, c='0.75', ls="--")
  
                  # get the limits for the arrow
                  max_MLE = max(ratio_MLEs)
                  min_MLE = min(ratio_MLEs)
                  dy = (max_MLE-min_MLE)/8
                  spacing = 1.5
                  # add arrow at best fit m/n
                  ax.add_patch(
                      patches.Arrow(
                          best_fit_movern, #x
                          max_MLE-(dy*spacing), #y
                          0, #dx
                          dy, #dy
                          width = 0.05,
                          facecolor = 'r',
                          edgecolor = 'r')
                      )
                  ax.text(best_fit_movern-0.09, max_MLE-1.5*(dy*spacing), "Best-fit " + r"$\theta$ = "+str(best_fit_movern),fontsize=8, color="r")
              #remove tribs
              #else:
                  # plot the data
                  #ax.scatter(m_over_n_values,ratio_MLEs, label = str(i), s=5, c=colours[i]) # different linestyle for each iteration?
  
          # set the axes labels
          ax.set_xlabel(r'$\theta$')
          ax.set_ylabel('MLE')
  
          # set the ylim
          ax.set_ylim(min_MLE,max_MLE+(dy*spacing))
  
          # add the legend
          #ax.legend(loc='right', bbox_to_anchor=(1.25,0.5), title = 'Iterations', frameon=False)
  
          # some formatting of the figure
          ax.spines['top'].set_linewidth(1)
          ax.spines['left'].set_linewidth(1)
          ax.spines['right'].set_linewidth(1)
          ax.spines['bottom'].set_linewidth(1)
  
          # label with the basin and m/n
          best_fit_movern = best_fit_movern_dict[basin_number][0]
          #title_string = "Basin "+str(basin_number)+"; Best fit $m/n$: "+str(best_fit_movern)
          #ax.set_title(title_string)
          # ax.text(0, 1.1, title_string,
          #     verticalalignment='top', horizontalalignment='left',
          #     transform=ax.transAxes,
          #     color='red', fontsize=10)
  
          #save the plot
          newFilename = MLE_directory+"MLE_fxn_movern_"+str(basin_number)+"."+FigFormat
  
          # This gets all the ticks, and pads them away from the axis so that the corners don't overlap
          ax.tick_params(axis='both', width=1, pad = 2)
          for tick in ax.xaxis.get_major_ticks():
              tick.set_pad(2)
          ax.xaxis.set_major_locator(ticker.MultipleLocator(base=0.1))
  
          plt.savefig(newFilename,format=FigFormat,dpi=300)
          ax.cla()
      plt.close(fig)
  
  def MakeMOverNSummaryPlot(DataDirectory, fname_prefix, basin_list=[], start_movern=0.2, d_movern=0.1, n_movern=7,
                            size_format='ESURF', FigFormat='png',
                            SA_channels=False, show_legend=True,parallel=False,Chi_disorder=False,
                            Chi_all=True,Chi_bootstrap=True, SA_raw=True, SA_segmented=True):
      """
      This function makes a summary plot of the best fit m/n from the different
      methods.
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_list: list of basin keys to analyse, default = [] (all basins)
          start_movern (float): the starting m/n value. Default is 0.2
          d_movern (float): the increment between the m/n values. Default is 0.1
          n_movern (float): the number of m/n values analysed. Default is 7.
          size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
          FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
          SA_channels (bool): If true, will include the SA data separated by channel
          show_legend (bool): If true, display the legend with the plot
          Chi_disorder (bool): if true, will include the data from the chi disorder method
  
      Returns:
          Makes a summary plot
  
      Author: FJC, edited SMM 25/05/2018
      """
      # check if a directory exists for the summary plots. If not then make it.
      summary_directory = DataDirectory+'summary_plots/'
      if not os.path.isdir(summary_directory):
          os.makedirs(summary_directory)
  
      from matplotlib.ticker import FuncFormatter, MaxNLocator
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # make a figure
      fig = makefigure(size_format)
  
  
      print("SHOW LEGEND", show_legend)
  
      if show_legend:
          gs = plt.GridSpec(100,100,bottom=0.15,left=0.05,right=0.75,top=0.9)
      else:
          gs = plt.GridSpec(100,100,bottom=0.15,left=0.05,right=0.95,top=0.95)
  
      ax = fig.add_subplot(gs[5:100,10:95])
  
      # read in the summary csv
      df = Helper.ReadMOverNSummaryCSV(summary_directory,fname_prefix)
      
      
  
      if basin_list != []:
          basin_keys = basin_list
      else:
          # get the basin keys
          basin_keys = df['basin_key'].tolist()
          print (basin_keys)
  
      df = df[df['basin_key'].isin(basin_keys)]
  
      # plot the full chi data
      full_chi_keys = df['basin_key'].values
      full_chi_keys = full_chi_keys.astype(float) - 0.2
      print("Full chi keys are: ")
      print(full_chi_keys)
      
      # Now check what is in the dataframe
      if 'Chi_MLE_full' in df:
          if Chi_all:
              ax.scatter(full_chi_keys, df['Chi_MLE_full'],marker='o', edgecolors='k', lw=0.5, facecolors='#e34a33', s=15, zorder=200, label='Chi all data')
  
      # Now see if there is points data
      if 'Chi_MLE_points' in df:
          # plot the points data
          median_movern = df['Chi_MLE_points'].values
          points_max_err = df['Chi_MLE_points_max'].values
          points_max_err = points_max_err.astype(float)-median_movern.astype(float)
          points_min_err = df['Chi_MLE_points_min'].values
          points_min_err = median_movern.astype(float)-points_min_err.astype(float)
          errors = np.array(list(zip(points_min_err, points_max_err))).T
      
          print("The errors or the point data are")
          print(errors)
  
          print(df['basin_key'].values)
          points_chi_keys = df['basin_key'].values
          points_chi_keys = points_chi_keys.astype(float) - 0.1
  
          if Chi_bootstrap:
              ax.errorbar(points_chi_keys, df['Chi_MLE_points'], s=15, marker='o', xerr=None, yerr=errors,
                      ecolor='#fdbb84', fmt='none', elinewidth=1,label='_nolegend_')
              ax.scatter(points_chi_keys, df['Chi_MLE_points'], s=15, c='#fdbb84', marker='o', edgecolors='k',
                     lw=0.5,facecolors='#fdbb84', label='Chi bootstrap',zorder=200)
  
      # plot the chi disorder data if you want it
      if 'Chi_disorder' in df:
          if Chi_disorder:
              median_movern = df['Chi_disorder'].values
              points_max_err = df['Chi_disorder_max'].values
              points_max_err = points_max_err.astype(float)-median_movern.astype(float)
              points_min_err = df['Chi_disorder_min'].values
              points_min_err = median_movern.astype(float)-points_min_err.astype(float)
              errors = np.array(list(zip(points_min_err, points_max_err))).T
  
              disorder_chi_keys = df['basin_key'].values
              disorder_chi_keys = disorder_chi_keys.astype(float)-0.3
              ax.errorbar(disorder_chi_keys, df['Chi_disorder'], s=15, marker='o', xerr=None, yerr=errors, ecolor='#F06292', fmt='none', elinewidth=1,label='_nolegend_')
              ax.scatter(disorder_chi_keys, df['Chi_disorder'],marker='o', edgecolors='k', lw=0.5, facecolors='#F06292', s=15, zorder=100, label='Chi disorder')
  
  
      # plot the SA data
      SA_keys = df['basin_key'].values
      SA_sterr = df['SA_raw_sterr'].values
  
      if SA_raw:
          ax.errorbar(SA_keys, df['SA_raw'], yerr=SA_sterr, c='#2b8cbe', elinewidth=1, fmt='none',label='_nolegend_')
          ax.scatter(SA_keys, df['SA_raw'], s=15, c='#2b8cbe', edgecolors='k',lw=0.5, label='S-A all data', zorder=100)
  
      if SA_channels:
          # plot the SA data by tribs
          median_movern = df['SA_tribs'].values
          points_max_err = df['SA_tribs_max'].values
          points_max_err = points_max_err.astype(float)-median_movern.astype(float)
          points_min_err = df['SA_tribs_min'].values
          points_min_err = median_movern.astype(float)-points_min_err.astype(float)
          errors = np.array(list(zip(points_min_err, points_max_err))).T
  
          SA_tribs_keys = df['basin_key'].values
          SA_tribs_keys = SA_tribs_keys.astype(float)+0.1
          ax.errorbar(SA_tribs_keys, df['SA_tribs'], s=15, marker='D', facecolors='white', xerr=None, yerr=errors, edgecolors='r', fmt='none', elinewidth=1, linestyle = ":", ecolor='r',label='_nolegend_')
          ax.scatter(SA_tribs_keys, df['SA_tribs'], s=15, marker='D', facecolors='white', edgecolors='r', label='S-A by channel',zorder=100)
  
      # plot the segmented SA data
      median_movern = df['SA_segments'].values
      points_max_err = df['SA_segments_max'].values
      points_max_err = points_max_err.astype(float)-median_movern.astype(float)
      points_min_err = df['SA_segments_min'].values
      points_min_err = median_movern.astype(float)-points_min_err.astype(float)
      errors = np.array(list(zip(points_min_err, points_max_err))).T
  
      SA_segment_keys = df['basin_key'].values
      SA_segment_keys = SA_segment_keys.astype(float)+0.2
      if SA_segmented:
          ax.errorbar(SA_segment_keys, df['SA_segments'], s=15, marker='o', facecolors='#a6bddb', xerr=None, yerr=errors,
                  edgecolors='#a6bddb', fmt='none', elinewidth=1, linestyle = ":", ecolor='#a6bddb',label='_nolegend_')
          ax.scatter(SA_segment_keys, df['SA_segments'], s=15, marker='o', facecolors='#a6bddb', edgecolors='k',
                 lw=0.5, label='Segmented S-A', zorder=100)
  
      # set the axis labels
      ax.set_xlabel('Basin key')
      ax.set_ylabel('Best fit '+r'$\theta$')
  
      if show_legend:
          print ("ADDING THE LEGEND")
          # sort both labels and handles by labels
          handles, labels = ax.get_legend_handles_labels()
          labels, handles = list(zip(*sorted(list(zip(labels, handles)), key=lambda t: t[0])))
          # add the legend
          ax.legend(handles, labels,fontsize=8, bbox_to_anchor=(1.0,0.7),bbox_transform=plt.gcf().transFigure)
  
      # This gets all the ticks, and pads them away from the axis so that the corners don't overlap
      ax.tick_params(axis='both', width=1, pad = 2)
      for tick in ax.xaxis.get_major_ticks():
          tick.set_pad(2)
  
      # change tick spacing
      ax.xaxis.set_major_locator(ticker.MultipleLocator(base=2))
      ax.set_xlim(-1,)
      #ax.yaxis.set_major_locator(ticker.MultipleLocator(base=d_movern))
  
      # remove first tick from the x axis
      #ax.xaxis.set_major_locator(MaxNLocator(prune='lower'))
  
      #set y axis lims
      #end_movern = end_movern = start_movern+d_movern*(n_movern-1)
      #ax.set_ylim(start_movern-d_movern,1.3)
  
      # change y axis to the moverns tested
      # end_movern = end_movern = start_movern+d_movern*(n_movern-1)
      # print end_movern
      # ax.yaxis.set_ticks(np.arange(start_movern, end_movern, d_movern))
  
  
      newFilename = summary_directory+fname_prefix+"_"
  
      if Chi_all:
          newFilename = newFilename+"A"
      if Chi_bootstrap:
          newFilename = newFilename+"B"
      if Chi_disorder:
          newFilename = newFilename+"D"
      if SA_raw:
          newFilename = newFilename+"S"
      if SA_channels:
          newFilename = newFilename+"C"
      if SA_segmented:
          newFilename = newFilename+"G"
  
      newFilename = newFilename+"_movern_summary."+FigFormat
  
  
      plt.savefig(newFilename,format=FigFormat,dpi=300)
      ax.cla()
      plt.close(fig)
  
  def MakeMOverNPlotOneMethod(DataDirectory, fname_prefix, basin_list=[], start_movern=0.2, d_movern=0.1, n_movern=7,
                              size_format='ESURF', FigFormat='png', movern_method='chi_points'):
          """
          This function makes a summary plot of the best fit m/n, you choose which method
          you want to plot.
  
          Args:
              DataDirectory (str): the data directory with the m/n csv files
              fname_prefix (str): The prefix for the m/n csv files
              basin_list: list of basin keys to analyse, default = [] (all basins)
              start_movern (float): the starting m/n value. Default is 0.2
              d_movern (float): the increment between the m/n values. Default is 0.1
              n_movern (float): the number of m/n values analysed. Default is 7.
              size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
              FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
              movern_method (str): the method you want to plot. Can be 'chi_all', 'chi_points', 'chi_disorder', SA_raw', or 'SA_segments'. Default 'chi_points'
          Returns:
              Makes a summary plot
  
          Author: FJC edit SMM 15/11/2019
          """
          # check if a directory exists for the summary plots. If not then make it.
          summary_directory = DataDirectory+'summary_plots/'
          if not os.path.isdir(summary_directory):
              os.makedirs(summary_directory)
  
          from matplotlib.ticker import FuncFormatter, MaxNLocator
          # Set up fonts for plots
          label_size = 10
          rcParams['font.family'] = 'sans-serif'
          rcParams['font.sans-serif'] = ['Liberation Sans']
          rcParams['font.size'] = label_size
  
          # make a figure
          fig = makefigure(size_format)
  
  
          gs = plt.GridSpec(100,100,bottom=0.15,left=0.05,right=0.95,top=0.95)
  
          ax = fig.add_subplot(gs[5:100,10:95])
  
          # read in the summary csv
          df = Helper.ReadMOverNSummaryCSV(summary_directory,fname_prefix)
          print (df)
  
          if basin_list != []:
              basin_keys = basin_list
          else:
              # get the basin keys
              basin_keys = df['basin_key'].tolist()
              print (basin_keys)
  
          df = df[df['basin_key'].isin(basin_keys)]
  
          # plot the full chi data
          if movern_method=='chi_all':
              full_chi_keys = df['basin_key'].values
              full_chi_keys = full_chi_keys.astype(float)-0.2
              ax.scatter(full_chi_keys, df['Chi_MLE_full'],marker='o', edgecolors='k', lw=0.5, facecolors='#e34a33', s=15, zorder=200, label='Chi all data')
  
          elif movern_method=='chi_points':
              # plot the points data
              median_movern = df['Chi_MLE_points'].values
              points_max_err = df['Chi_MLE_points_max'].values
              points_max_err = points_max_err-median_movern
              points_min_err = df['Chi_MLE_points_min'].values
              points_min_err = median_movern-points_min_err
              errors = np.array(zip(points_min_err, points_max_err)).T
  
              points_chi_keys = df['basin_key'].values # -0.1 I removed that to fix the movern plots - Boris
              ax.errorbar(points_chi_keys, df['Chi_MLE_points'], s=15, marker='o', xerr=None, yerr=errors, ecolor='#fdbb84', fmt='none', elinewidth=1,label='_nolegend_')
              ax.scatter(points_chi_keys, df['Chi_MLE_points'], s=15, c='#fdbb84', marker='o', edgecolors='k', lw=0.5,facecolors='#fdbb84', label='Chi bootstrap',zorder=200)
  
              # plot the SA data
          elif movern_method=='SA_raw':
                  SA_keys = df['basin_key'].values
                  SA_sterr = df['SA_raw_sterr'].values
                  ax.errorbar(SA_keys, df['SA_raw'], yerr=SA_sterr, c='#2b8cbe', elinewidth=1, fmt='none',label='_nolegend_')
                  ax.scatter(SA_keys, df['SA_raw'], s=15, c='#2b8cbe', edgecolors='k',lw=0.5, label='S-A all data', zorder=100)
  
                  
          # plot the chi disorder data if you want it
          elif movern_method=='chi_disorder':
              median_movern = df['Chi_disorder'].values
              points_max_err = df['Chi_disorder_max'].values
              points_max_err = points_max_err.astype(float)-median_movern.astype(float)
              points_min_err = df['Chi_disorder_min'].values
              points_min_err = median_movern.astype(float)-points_min_err.astype(float)
              errors = np.array(list(zip(points_min_err, points_max_err))).T
  
              disorder_chi_keys = df['basin_key'].values
              disorder_chi_keys = disorder_chi_keys.astype(float)-0.3
              ax.errorbar(disorder_chi_keys, df['Chi_disorder'], s=15, marker='o', xerr=None, yerr=errors, ecolor='#F06292', fmt='none', elinewidth=1,label='_nolegend_')
              ax.scatter(disorder_chi_keys, df['Chi_disorder'],marker='o', edgecolors='k', lw=0.5, facecolors='#F06292', s=15, zorder=100, label='Chi disorder')                
                  
          else:
              # plot the segmented SA data
              median_movern = df['SA_segments'].values
              points_max_err = df['SA_segments_max'].values
              points_max_err = points_max_err-median_movern
              points_min_err = df['SA_segments_min'].values
              points_min_err = median_movern-points_min_err
              errors = np.array(zip(points_min_err, points_max_err)).T
  
              SA_segment_keys = df['basin_key'].values
              SA_segment_keys=SA_segment_keys.astype(float)+0.2
              ax.errorbar(SA_segment_keys, df['SA_segments'], s=15, marker='o', facecolors='#a6bddb', xerr=None, yerr=errors, edgecolors='#a6bddb', fmt='none', elinewidth=1, linestyle = ":", ecolor='#a6bddb',label='_nolegend_')
              ax.scatter(SA_segment_keys, df['SA_segments'], s=15, marker='o', facecolors='#a6bddb', edgecolors='k', lw=0.5, label='Segmented S-A', zorder=100)
  
          # set the axis labels
          ax.set_xlabel('Basin key')
          ax.set_ylabel('Best fit '+r'$\theta$')
  
          # This gets all the ticks, and pads them away from the axis so that the corners don't overlap
          ax.tick_params(axis='both', width=1, pad = 2)
          for tick in ax.xaxis.get_major_ticks():
              tick.set_pad(2)
  
          # change tick spacing
          ax.xaxis.set_major_locator(ticker.MultipleLocator(base=1))
          ax.set_xlim(-0.5,df['basin_key'].max()+0.5)
  
          newFilename = summary_directory+fname_prefix+"_movern_"+movern_method+"."+FigFormat
          plt.savefig(newFilename,format=FigFormat,dpi=300)
          ax.cla()
          plt.close(fig)
  
          
          
  def MakeMOverNDisorderDistancePlot(DataDirectory, fname_prefix, basin_list_list=[], start_movern=0.2, d_movern=0.1, n_movern=7,
                              size_format='ESURF', FigFormat='png', show_legend=True,parallel=False, group_names=[]): 
      """
      This function makes a summary plot of the best fit m/n, you choose which method
      you want to plot.
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_list: list of basin keys to analyse, default = [] (all basins)
          start_movern (float): the starting m/n value. Default is 0.2
          d_movern (float): the increment between the m/n values. Default is 0.1
          n_movern (float): the number of m/n values analysed. Default is 7.
          size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
          FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
          movern_method (str): the method you want to plot. Can be 'chi_all', 'chi_points', 'chi_disorder', SA_raw', or 'SA_segments'. Default 'chi_points'
      Returns:
          Makes a summary plot
  
      Author: SMM 15/11/2019
      """
      # check if a directory exists for the summary plots. If not then make it.
      summary_directory = DataDirectory+'summary_plots/'
      if not os.path.isdir(summary_directory):
          os.makedirs(summary_directory)
  
      from matplotlib.ticker import FuncFormatter, MaxNLocator
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # make a figure
      fig = makefigure(size_format)
  
  
      print("SHOW LEGEND", show_legend)
  
      if show_legend:
          gs = plt.GridSpec(100,100,bottom=0.15,left=0.05,right=0.75,top=0.9)
      else:
          gs = plt.GridSpec(100,100,bottom=0.15,left=0.05,right=0.95,top=0.95)
  
      ax = fig.add_subplot(gs[5:100,10:95])
  
      # read in the summary csv
      df = Helper.ReadMOverNSummaryCSV(summary_directory,fname_prefix)
      
      # you also need the basin info
      df_basin_info = Helper.ReadBasinInfoCSV(DataDirectory,fname_prefix)
      
      #print("The data frame is")
      #print(df)
      
      #print("The basin info df is:")
      #print(df_basin_info)
      
      # Now merge the two dataframes
      df_new = pd.merge(df,df_basin_info,on='basin_key')
      df = df_new
      
      #print("The new dataframe is")
      #print(df)
  
      if basin_list_list != []:
          basin_keys_list = basin_list_list
      else:
          # get the basin keys
          basin_keys_list = []
          basin_keys = df['basin_key'].tolist()
          basin_keys_list.append(basin_keys)
          print("You didn't select basin keys. I am picking all of them:")
          print (basin_keys_list)
  
      # This just makes a colour list for plotting
      colour_list = []
      colour_list.append('#e41a1c')
      colour_list.append('#377eb8')
      colour_list.append('#4daf4a')
      colour_list.append('#984ea3')
      colour_list.append('#ff7f00')
      colour_list.append('#ffff33')
      colour_list.append('#a65628')
      colour_list.append('#f781bf')
      colour_list.append('#999999')
      
      tab20_cm = cm.get_cmap('tab20')
      
          
      # Loop through the basin key list
      colour_index = 0;
      for basin_keys in basin_list_list:
          
          
          print("\n\n\nI am selecting the following basins")
          print(basin_keys)
          print("The colour indexs is: "+str(colour_index))
          
          this_df = df[df['basin_key'].isin(basin_keys)]
  
          # plot the chi disorder data if you want it. This will fail if the 
          
          if 'Chi_disorder' in this_df:
              colour_index_mod = colour_index % 9
              if group_names == []:
                  label_name = "Group "+str(colour_index)
              else:
                  label_name = group_names[colour_index]
              median_movern = this_df['Chi_disorder'].values
              points_max_err = this_df['Chi_disorder_max'].values
              points_max_err = points_max_err.astype(float)-median_movern.astype(float)
              points_min_err = this_df['Chi_disorder_min'].values
              points_min_err = median_movern.astype(float)-points_min_err.astype(float)
              errors = np.array(list(zip(points_min_err, points_max_err))).T
              lat_vals =  this_df['latitude'].values
              long_vals =  this_df['longitude'].values
  
              disorder_chi_keys = this_df['basin_key'].values
              disorder_chi_keys = disorder_chi_keys.astype(float)-0.3
              ax.errorbar(long_vals, this_df['Chi_disorder'], s=15, marker='o', xerr=None, yerr=errors, ecolor=colour_list[colour_index_mod], fmt='none', elinewidth=1,label='_nolegend_')
              ax.scatter(long_vals, this_df['Chi_disorder'],marker='o', edgecolors='k', lw=0.5, facecolors=colour_list[colour_index_mod], s=15, zorder=100, label=label_name)
          colour_index = colour_index+1
  
  
  
      # set the axis labels
      ax.set_xlabel('Longitude (decimal degrees)')
      ax.set_ylabel('Best fit '+r'$\theta$')
  
      if show_legend:
          print ("ADDING THE LEGEND")
          # sort both labels and handles by labels
          handles, labels = ax.get_legend_handles_labels()
          labels, handles = list(zip(*sorted(list(zip(labels, handles)), key=lambda t: t[0])))
          # add the legend
          ax.legend(handles, labels,fontsize=8, bbox_to_anchor=(1.0,0.7),bbox_transform=plt.gcf().transFigure)
  
      # This gets all the ticks, and pads them away from the axis so that the corners don't overlap
      ax.tick_params(axis='both', width=1, pad = 2)
      for tick in ax.xaxis.get_major_ticks():
          tick.set_pad(2)
  
      # change tick spacing
      #ax.xaxis.set_major_locator(ticker.MultipleLocator(base=2))
      #ax.set_xlim(-1,)
      #ax.yaxis.set_major_locator(ticker.MultipleLocator(base=d_movern))
  
      # remove first tick from the x axis
      #ax.xaxis.set_major_locator(MaxNLocator(prune='lower'))
  
      #set y axis lims
      #end_movern = end_movern = start_movern+d_movern*(n_movern-1)
      #ax.set_ylim(start_movern-d_movern,1.3)
  
      # change y axis to the moverns tested
      # end_movern = end_movern = start_movern+d_movern*(n_movern-1)
      # print end_movern
      # ax.yaxis.set_ticks(np.arange(start_movern, end_movern, d_movern))
  
  
      newFilename = summary_directory+fname_prefix+"_concavity_dist."+FigFormat
      plt.savefig(newFilename,format=FigFormat,dpi=300)
      ax.cla()
      plt.close(fig)    
          
          
  def MakeMOverNSummaryHistogram(DataDirectory, fname_prefix, basin_list=[], size_format='ESURF', FigFormat='png', start_movern=0.1, n_movern=7, d_movern=0.1, mn_method = "Chi", show_legend=True, parallel=False, Chi_disorder=False):
      """
      This function makes a stacked histogram showing the distribution of concavity (~m/n) values for each method
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
          FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
          start_movern (float): the starting m/n value. Default is 0.2
          d_movern (float): the increment between the m/n values. Default is 0.1
          n_movern (float): the number of m/n values analysed. Default is 7.
          mn_method (str): the method of calculating m/n, either "Chi", or "SA". Default "Chi"
          show_legend (bool): if true, add a legend to the plot (default=True)
          Chi_disorder(bool): if true, add the chi disorder analysis
  
      Returns:
          Histogram of the m/n values
  
      Author: FJC
      """
  
      # check if a directory exists for the summary plots. If not then make it.
      summary_directory = DataDirectory+'summary_plots/'
      if not os.path.isdir(summary_directory):
          os.makedirs(summary_directory)
  
      from matplotlib.ticker import FuncFormatter, MaxNLocator
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # make a figure
      fig = makefigure(size_format)
      figsize = fig.get_size_inches()
  
      # if show_legend:
      #     gs = plt.GridSpec(100,100,bottom=0.15,left=0.05,right=0.75,top=0.9)
      # else:
      #     gs = plt.GridSpec(100,100,bottom=0.1,left=0.05,right=0.95,top=0.95)
      #
      # ax = fig.add_subplot(gs[5:100,10:95])
  
      # read in the summary csv
      if not parallel:
        df = Helper.ReadMOverNSummaryCSV(summary_directory,fname_prefix)
      else:
        print ("MakeMOverNSummaryHistogram not parallelised yet")
        return
      
      Chi_disorder_only = False
      if Chi_disorder: 
          if 'Chi_MLE_points' not in df:
              print("I didn't find a bootstrap data line. Switching to disorder only")
              Chi_disorder_only= True
  
      # get the basin keys
      basin_keys = df['basin_key'].tolist()
      print (basin_keys)
      if Chi_disorder_only:
          columns = ['Chi_disorder', 'SA_raw', 'SA_segments']
          these_labels = ['Chi disorder', 'S-A all data', 'Segmented S-A']
          colours = ['#F06292', '#2b8cbe', '#a6bddb']
      else:
          if Chi_disorder:
              columns = ['Chi_MLE_full', 'Chi_MLE_points', 'Chi_disorder', 'SA_raw', 'SA_segments']
              these_labels = ['Chi all data', 'Chi bootstrap', 'Chi disorder', 'S-A all data', 'Segmented S-A']
              colours = ['#e34a33', '#fdbb84', '#F06292', '#2b8cbe', '#a6bddb']
  
          else:
              columns = ['Chi_MLE_full', 'Chi_MLE_points', 'SA_raw', 'SA_segments']
              these_labels = ['Chi all data', 'Chi bootstrap', 'S-A all data', 'Segmented S-A']
              colours = ['#e34a33', '#fdbb84', '#2b8cbe', '#a6bddb']
              
      x_spacing = 0.05
      fig, ax = joyplot.joyplot(df, figsize=figsize, column=columns, label_strings=these_labels, x_range=[0,1],grid="x",color=colours,x_title='Best fit '+r'$\theta$' +' distribution',x_spacing=x_spacing)
      #plt.xlabel('Best fit $m/n$')
      
      # change tick spacing
      x_loc = np.arange(0,1.01,0.05)
      print("x_loc is:")
      print(x_loc)
      labels = ["0","","","","0.2","","","","0.4","","","","0.6","","","","0.8","","","","1"]
      plt.xticks(x_loc, labels, rotation='vertical')
      
  
      newFilename = summary_directory+fname_prefix+"_movern_hist."+FigFormat
  
      plt.savefig(newFilename,format=FigFormat,dpi=300)
      plt.close(fig)
  
  def PlotMOverNByBasin(DataDirectory, fname_prefix, basin_list = [], size_format='ESURF', FigFormat='png', colour_points=False, column_header="Lithology"):
      """
      This function makes a plot of the best-fit m/n value calculated using the chi Monte Carlo points analysis
      for each basin. The points are coloured by another value which is specified by the user (e.g. lithology).
      The column header must be specified by the user.
  
      This is a work in progress. I haven't tested it yet and wrote it on the plane.
      I need to think about what is the best way to read in the lithology - will we have a string
      for each lithology? Then we want to map each string to a different colour.
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          basin_list: a list of the basins to analyse
          size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
          FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
          colour_points (bool): boolean to decide whether to colour the points by another value, e.g. lithology. Default = False
          column_header (str): if colour_points=True, then use this to tell us what column header in the csv you want to use
  
      Returns:
          m/n plot by basin_list
  
      Author: FJC
      """
  
      # check if a directory exists for the summary plots. If not then make it.
      summary_directory = DataDirectory+'summary_plots/'
      if not os.path.isdir(summary_directory):
          os.makedirs(summary_directory)
  
      # read in the summary csv file
      df = Helper.ReadMOverNSummaryCSV(DataDirectory,fname_prefix)
  
      # get the basin keys
      if basin_list == []:
          basin_list = df['basin_keys'].tolist()
  
      # mask the dataframe for the required basin keys
      df = df[df['basin_key'].isin(basin_list)]
  
      columns = df.columns()
  
      # set up the figure
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # make a figure
      fig = makefigure(size_format)
  
  
      gs = plt.GridSpec(100,100,bottom=0.1,left=0.05,right=0.95,top=0.95)
      ax = fig.add_subplot(gs[5:100,10:95])
  
      # check if the user wants to colour the points by a specific columns
      if colour_points:
          if not column_header.isin(columns):
              print ("The column header name does not exist! Please check your dataframe")
          else:
              values_to_colour = df[column_header].tolist()
              this_cmap = plt.cm.Set2
              ax.errorbar(basin_list, df['Chi_MLE_points'], s=20, marker='o', xerr=None, yerr=errors, ecolor=values_to_colour, cmap=this_cmap, fmt='none', elinewidth=1,label='_nolegend_')
              ax.scatter(basin_list, df['Chi_MLE_points'], s=20, c=values_to_colour, cmap=this_cmap, marker='o', edgecolors=None, lw=0.5, label='Chi bootstrap',zorder=200)
      else:
          ax.errorbar(basin_list, df['Chi_MLE_points'], s=20, marker='o', xerr=None, yerr=errors, ecolor='#fdbb84', fmt='none', elinewidth=1,label='_nolegend_')
          ax.scatter(basin_list, df['Chi_MLE_points'], s=20, c='#fdbb84', marker='o', edgecolors='k', lw=0.5,facecolors='#fdbb84', label='Chi bootstrap',zorder=200)
  
      # set the axis labels
      ax.set_xlabel('Basin key')
      ax.set_ylabel('Best fit '+r'$\theta$')
  
      # This gets all the ticks, and pads them away from the axis so that the corners don't overlap
      ax.tick_params(axis='both', width=1, pad = 2)
      for tick in ax.xaxis.get_major_ticks():
          tick.set_pad(2)
  
      # change tick spacing
      ax.xaxis.set_major_locator(ticker.MultipleLocator(base=1))
  
      newFilename = summary_directory+fname_prefix+"_movern_summary."+FigFormat
      if SA_channels:
          newFilename = summary_directory+fname_prefix+"_movern_summary_with_SA_tribs."+FigFormat
      plt.savefig(newFilename,format=FigFormat,dpi=300)
      ax.cla()
      plt.close(fig)
  
  
  
  # def MakeBasinJoyplot(DataDirectory, fname_prefix, basin_list=[], size_format='ESURF', FigFormat='png'):
  #     """
  #     This function makes a joyplot showing the m/n value for a list of basins. m/n value
  #     is calculated using the chi points method.
  #
  #     Args:
  #         DataDirectory (str): the data directory with the m/n csv files
  #         fname_prefix (str): The prefix for the m/n csv files
  #         basin_list: a list of the basins to analyse
  #         size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
  #         FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
  #
  #     Returns:
  #         Joyplot of the m/n values for each basin
  #
  #     Author: FJC
  #     """
  #
  #     # check if a directory exists for the summary plots. If not then make it.
  #     summary_directory = DataDirectory+'summary_plots/'
  #     if not os.path.isdir(summary_directory):
  #         os.makedirs(summary_directory)
  #
  #     from matplotlib.ticker import FuncFormatter, MaxNLocator
  #     # Set up fonts for plots
  #     label_size = 10
  #     rcParams['font.family'] = 'sans-serif'
  #     rcParams['font.sans-serif'] = ['Liberation Sans']
  #     rcParams['font.size'] = label_size
  #
  #     # make a figure
  #     if size_format == "geomorphology":
  #         #fig = plt.figure(1, facecolor='white',figsize=(6.25,3.5))
  #         figsize=(6.25,3.5)
  #         #l_pad = -40
  #     elif size_format == "big":
  #         #fig = plt.figure(1, facecolor='white',figsize=(16,9))
  #         figsize=(16,9)
  #         #l_pad = -50
  #     else:
  #         #fig = plt.figure(1, facecolor='white',figsize=(4.92126,3.2))
  #         figsize=(4.92126,3.2)
  #         #l_pad = -35
  #
  #     # if show_legend:
  #     #     gs = plt.GridSpec(100,100,bottom=0.15,left=0.05,right=0.75,top=0.9)
  #     # else:
  #     #     gs = plt.GridSpec(100,100,bottom=0.1,left=0.05,right=0.95,top=0.95)
  #     #
  #     # ax = fig.add_subplot(gs[5:100,10:95])
  #
  #     # read in the summary csv
  #     df = Helper.ReadMOverNSummaryCSV(summary_directory,fname_prefix)
  #     print df
  #
  #     # get the basin keys
  #     if basin_list == []:
  #         print "You didn't give me a list of basins so I'll analyse all of them"
  #         basin_list = df['basin_key'].tolist()
  #
  #     df = df[df['basin_key'].isin(basin_list)]
  #
  #     columns = ['basin_key']
  #     #these_labels = ['Chi all data', 'Chi Monte Carlo', 'S-A all data', 'Segmented S-A']
  #     #colours = ['#e34a33', '#fdbb84', '#2b8cbe', '#a6bddb']
  #     x_spacing = 0.2
  #     fig, ax = joyplot.joyplot(df, figsize=figsize, column=columns, x_range=[0,1],grid="x",x_title='Best fit $m/n$ distribution',x_spacing=x_spacing)
  #     #plt.xlabel('Best fit $m/n$')
  #
  #     newFilename = summary_directory+fname_prefix+"_basin_joyplots."+FigFormat
  #
  #     plt.savefig(newFilename,format=FigFormat,dpi=300)
  #     plt.close(fig)
  
  #=============================================================================
  # RASTER PLOTTING FUNCTIONS
  # Functions that interface with LSDMapFigure to plot the m/n analysis with
  # spatial data.
  #=============================================================================
  def MakeRasterPlotsBasins(DataDirectory, fname_prefix, size_format='ESURF', FigFormat='png', parallel=False):
      """
      This function makes a shaded relief plot of the DEM with the basins coloured
      by the basin ID.
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
          FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
  
      Returns:
          Shaded relief plot with the basins coloured by basin ID
  
      Author: FJC
      """
      # check if a directory exists for the chi plots. If not then make it.
      raster_directory = DataDirectory+'raster_plots/'
      if not os.path.isdir(raster_directory):
          os.makedirs(raster_directory)
  
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # set figure sizes based on format
      # make a figure
      fig = makefigure(size_format)
      fig_width_inches = fig.get_size_inches()[0]
  
      # get the basin IDs to make a discrete colourmap for each ID
      if not parallel:
          BasinInfoDF = Helper.ReadBasinInfoCSV(DataDirectory, fname_prefix)
      else:
          BasinInfoDF = Helper.AppendBasinInfoCSVs(DataDirectory,fname_prefix)
  
      print (BasinInfoDF)
  
      basin_keys = list(BasinInfoDF['basin_key'])
      basin_keys = [int(x) for x in basin_keys]
  
      basin_junctions = list(BasinInfoDF['outlet_junction'])
      basin_junctions = [int(x) for x in basin_junctions]
  
      print ('Basin keys are: ')
      print (basin_keys)
  
      # get a discrete colormap
      cmap = plt.cm.jet
  
      # going to make the basin plots - need to have bil extensions.
      print("I'm going to make the basin plots. Your topographic data must be in ENVI bil format or I'll break!!")
  
      # get the rasters
      raster_ext = '.bil'
      BackgroundRasterName = fname_prefix+raster_ext
      HillshadeName = fname_prefix+'_hs'+raster_ext
      BasinsName = fname_prefix+'_AllBasins'+raster_ext
  
      # create the map figure
      MF = MapFigure(HillshadeName, DataDirectory,coord_type="UTM_km", colourbar_location='None')
  
      # add the basins drape
      MF.add_basin_plot(BasinsName,fname_prefix,DataDirectory,
                        use_keys_not_junctions = True, show_colourbar = True,
                        discrete_cmap=True, n_colours=len(basin_keys), colorbarlabel = "Basin ID",
                        colourmap = cmap, adjust_text = False, parallel=parallel)
  
      # add the basin outlines ### need to parallelise
      if not parallel:
        Basins = LSDP.GetBasinOutlines(DataDirectory, BasinsName)
      else:
        Basins = LSDP.GetMultipleBasinOutlines(DataDirectory)
  
      MF.plot_polygon_outlines(Basins, linewidth=0.8)
  
      # add the channel network
      if not parallel:
          ChannelDF = Helper.ReadChiDataMapCSV(DataDirectory,fname_prefix)
      else:
          ChannelDF = Helper.AppendChiDataMapCSVs(DataDirectory,fname_prefix)
      ChannelPoints = LSDP.LSDMap_PointData(ChannelDF, data_type = "pandas", PANDEX = True)
      MF.add_point_data(ChannelPoints,show_colourbar="False", scale_points=True, column_for_scaling='drainage_area',alpha=0.5,zorder=100)
  
      # add the basin labelling
      label_dict = dict(zip(basin_junctions,basin_keys))
  
      if not parallel:
        Points = LSDP.GetPointWithinBasins(DataDirectory, BasinsName)
      else:
        Points = LSDP.GetPointsWithinMultipleBasins(DataDirectory, BasinsName)
  
      MF.add_text_annotation_from_shapely_points(Points, text_colour='k', label_dict=label_dict, zorder=200)
  
      # Save the figure
      ImageName = raster_directory+fname_prefix+'_basin_keys.'+FigFormat
      MF.save_fig(fig_width_inches = fig_width_inches, FigFileName = ImageName, FigFormat=FigFormat, Fig_dpi = 300)
  
  def MakeRasterPlotsMOverN(DataDirectory, fname_prefix, start_movern=0.2, n_movern=7, d_movern=0.1, movern_method='Chi_full', size_format='ESURF', FigFormat='png',lith = False, parallel=False):
      """
      This function makes a shaded relief plot of the DEM with the basins coloured
      by the best fit m/n
  
      Args:
          DataDirectory (str): the data directory with the m/n csv files
          fname_prefix (str): The prefix for the m/n csv files
          start_movern (int): The starting m/n, default = 0.2
          n_movern (int): The number of m/n values tested, default = 7.
          movern_method: the method of estimating m over n. Options are full chi "Chi_full", points "Chi_points", or slope-area "SA". Default is "Chi_full".
          size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
          FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
          lith (bool): the lithologic information if available. You need to follow the lithologic documentation to generate the right files. (This last does not exists yet).
  
      Returns:
          Shaded relief plot with the basins coloured by best fit m/n
  
      Author: FJC
      """
      # check if a directory exists for the chi plots. If not then make it.
      raster_directory = DataDirectory+'raster_plots/'
      if not os.path.isdir(raster_directory):
          os.makedirs(raster_directory)
  
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # make a figure
      fig = makefigure(size_format)
      fig_width_inches = fig.get_size_inches()[0]
  
      # get the basin IDs to make a discrete colourmap for each ID
      if not parallel:
          BasinDF = Helper.ReadBasinStatsCSV(DataDirectory,fname_prefix)
          # get the basin IDs to make a discrete colourmap for each ID
          BasinInfoDF = Helper.ReadBasinInfoCSV(DataDirectory, fname_prefix)
      else:
          BasinDF = Helper.AppendBasinCSVs(DataDirectory,fname_prefix)
          BasinInfoDF = Helper.AppendBasinInfoCSVs(DataDirectory,fname_prefix)
  
      basin_keys = list(BasinInfoDF['basin_key'])
      basin_keys = [int(x) for x in basin_keys]
  
      basin_junctions = list(BasinInfoDF['outlet_junction'])
      basin_junctions = [int(x) for x in basin_junctions]
  
      print ('Basin keys are: ')
      print (basin_keys)
  
      labeldict = {}
      # get the best fit m/n for each basin
      if movern_method == "Chi_full":
          print ('Making a plot for the full chi method')
          Outlier_counter, removed_sources_dict, best_fit_movern_dict, MLEs_dict = CheckMLEOutliers(DataDirectory, fname_prefix, basin_list=basin_keys, start_movern=start_movern, d_movern=d_movern, n_movern=n_movern, parallel=parallel)
          #MOverNDict = SimpleMaxMLECheck(BasinDF)
          m_over_ns = [round(i[0],2) for i in best_fit_movern_dict.values()]
          print ('You tested these m_over_ns')
          print (m_over_ns)
          MOverNDict = dict(zip(basin_keys,m_over_ns))
          label_list = [str(i)+": "+str(j) for i,j in zip(basin_keys,m_over_ns)]
          labeldict = dict(zip(basin_junctions,label_list))
          ImageName = raster_directory+fname_prefix+'_basins_movern_chi_full.'+FigFormat
          title = "Chi analysis (all data)"
      elif movern_method == "Chi_points":
          print ('Making a plot for the chi points method')
          if not parallel:
              PointsChiBasinDF = Helper.ReadMCPointsCSV(DataDirectory,fname_prefix)
          else:
              PointsChiBasinDF = Helper.AppendBasinPointCSVs(DataDirectory,fname_prefix)
  
          PointsDF = GetMOverNRangeMCPoints(PointsChiBasinDF,start_movern,d_movern,n_movern)
          moverns = PointsDF['Median_MOverNs'].tolist()
          MOverNDict = dict(zip(basin_keys,moverns))
          # labelling the basins
          moverns = [round(i,2) for i in moverns]
          label_list = [str(i)+": "+str(j) for i,j in zip(basin_keys,moverns)]
          labeldict = dict(zip(basin_junctions,label_list))
          ImageName = raster_directory+fname_prefix+'_basins_movern_chi_points.'+FigFormat
          title = "Chi bootstrap analysis"
      elif movern_method == "SA":
          SlopeAreaDF = SA.LinearRegressionRawData(DataDirectory,fname_prefix,parallel=parallel)
          moverns = SlopeAreaDF['regression_slope'].tolist()
          MOverNDict = dict(zip(basin_keys,moverns))
          # labelling the basins
          moverns = [round(i,2) for i in moverns]
          label_list = [str(i)+": "+str(j) for i,j in zip(basin_keys,moverns)]
          labeldict = dict(zip(basin_junctions,label_list))
          ImageName = raster_directory+fname_prefix+'_basins_movern_SA.'+FigFormat
          title = "Slope-area analysis"
      elif movern_method == "Chi_disorder":
          #DisorderDF = Helper.ReadDisorderCSV(DataDirectory, fname_prefix)
          if not parallel:
              DisorderDF = Helper.ReadDisorderUncertCSV(DataDirectory, fname_prefix)
          else:
              DisorderDF = Helper.AppendDisorderCSV(DataDirectory, fname_prefix)
          m_over_ns = DisorderDF['median'].tolist()
          MOverNDict = dict(zip(basin_keys,m_over_ns))
          moverns = [round(i,2) for i in m_over_ns]
          label_list = [str(i)+": "+str(j) for i,j in zip(basin_keys,moverns)]
          labeldict = dict(zip(basin_junctions,label_list))
          ImageName = raster_directory+fname_prefix+'_basins_movern_chi_disorder.'+FigFormat
          title = "Chi disorder analysis"
  
      else:
          print ("You didn't select an appropriate movern method. Please choose either 'Chi_full', 'Chi_points, or 'SA'.")
          sys.exit()
  
      # get moverns for cbar plotting. We always want a spacing of 0.1.
      min_max_str = ['min', 'max']
      end_movern = float(start_movern)+float(d_movern)*(float(n_movern)-1)
      all_moverns = np.linspace(start_movern,end_movern,n_movern)
      print ("END MOVERN:")
      print (end_movern)
      min_max_moverns = [start_movern, end_movern]
      cbar_dict = dict(zip(min_max_str,min_max_moverns))
  
      # work out how many colours we need for the colormap
      n_colours = len(all_moverns)
      print ("N colours is: "+str(n_colours))
  
      # get a discrete colormap
      mn_cmap = plt.cm.Reds
  
      # going to make the basin plots - need to have bil extensions.
      print("I'm going to make the basin m/n plots. Your topographic data must be in ENVI bil format or I'll break!!")
  
      # get the rasters
      raster_ext = '.bil'
      BackgroundRasterName = fname_prefix+raster_ext
      HillshadeName = fname_prefix+'_hs'+raster_ext
      BasinsName = fname_prefix+'_AllBasins'+raster_ext
  
      print ("THE TITLE IS:"+title)
      # create the map figure
      MF = MapFigure(HillshadeName, DataDirectory,coord_type="UTM_km", plot_title=title)
      # add the basins drape
      if lith:
  
          LithoMap = fname_prefix+"_LITHRAST"+raster_ext
          df_litho = pd.read_csv(DataDirectory+fname_prefix+"_lithokey.csv")
          MF.add_drape_image(LithoMap,DataDirectory,colourmap = plt.cm.jet,
                                  alpha=0.4,
                                  show_colourbar = False,
                                  colorbarlabel = "Colourbar", discrete_cmap=True, n_colours=df_litho.shape[0],
                                  norm = "None",
                                  colour_min_max = [],
                                  modify_raster_values=False,
                                  old_values=[], new_values=[], cbar_type=int,
                                  NFF_opti = False, custom_min_max = [])
  
      else:
          MF.add_basin_plot(BasinsName,fname_prefix,DataDirectory, value_dict = MOverNDict,
                            use_keys_not_junctions = True, show_colourbar = False,
                            discrete_cmap=True, n_colours=n_colours, colorbarlabel = "$m/n$",
                            colourmap = mn_cmap, adjust_text = False, cbar_dict=cbar_dict, parallel=parallel)
  
      # add the basin outlines
      if not parallel:
        Basins = LSDP.GetBasinOutlines(DataDirectory, BasinsName)
      else:
        Basins = LSDP.GetMultipleBasinOutlines(DataDirectory)
  
      MF.plot_polygon_outlines(Basins, linewidth=0.8)
  
      # add the channel network
      if not parallel:
          ChannelDF = Helper.ReadChiDataMapCSV(DataDirectory,fname_prefix)
      else:
          ChannelDF = Helper.AppendChiDataMapCSVs(DataDirectory,fname_prefix)
      ChannelPoints = LSDP.LSDMap_PointData(ChannelDF, data_type = "pandas", PANDEX = True)
      MF.add_point_data(ChannelPoints,show_colourbar="False", scale_points=True, column_for_scaling='drainage_area',alpha=0.1,zorder=100)
  
      # add the basin labelling
      if not parallel:
        Points = LSDP.GetPointWithinBasins(DataDirectory, BasinsName)
      else:
        Points = LSDP.GetPointsWithinMultipleBasins(DataDirectory, BasinsName)
  
      print ("Adding labels, the label dict is:", labeldict)
      MF.add_text_annotation_from_shapely_points_v2(Points, text_colour='k', label_dict=labeldict, zorder=200)
  
      if(lith):
          ImageName = ImageName[0:-(len(FigFormat)+1)]+"_lith."+FigFormat # adding a particle to the file
      MF.save_fig(fig_width_inches = fig_width_inches, FigFileName = ImageName, FigFormat=FigFormat, Fig_dpi = 300) # Save the figure
  
  
  #=============================================================================
  # UNCERTAINTY
  # Functions that make plots of uncertainties on the m/n analysis and comparison
  # between the different methods
  #=============================================================================
  
  def PlotMOverNDicts(DataDirectory,fname_prefix,SA_based_dict,Chi_based_dict, FigFormat = "png", size_format = "ESURF"):
      """
      This plots the estimates of m over n for basins contained in the SA and Chi based dicts
  
      Args:
          SA_based_dict (dict): A dictionary where the key is the basin and the values are lists of the best fit m/n for S-A plots
          Chi_based_dict (dict):  A dictionary where the key is the basin and the values are lists of the best fit m/n after removing outliers
  
      Author: SMM
      """
  
      # We need to deal with the colours. Each colour will indicate a different
      # Data point
  
  
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # make a figure
      fig = makefigure(size_format)
  
      gs = plt.GridSpec(100,100,bottom=0.15,left=0.1,right=0.85,top=0.9)
      ax = fig.add_subplot(gs[5:100,10:95])
      ax.grid(zorder=-100)
  
      # Loop through the basins adding data
      for basin in SA_based_dict:
          this_SA_list = SA_based_dict[basin]
          this_chi_list = Chi_based_dict[basin]
  
          this_SA_list = np.multiply(this_SA_list,-1)
          full_list = np.concatenate((this_SA_list,this_chi_list))
          basin_list = [basin]*len(full_list)
  
          colours = []
          this_num = -2
          for element in this_SA_list:
              colours.append(this_num)
              this_num = this_num+0.25
  
          this_num = 1
          for element in this_chi_list:
              colours.append(this_num)
              this_num = this_num+0.25
  
  
          #print("Hey there pal, let me print out the two lists for you")
          #print(this_SA_list)
          #print(this_chi_list)
  
          #print("These are the colour and basin lists")
          #print(full_list)
          #print(basin_list)
  
  
          #colours = range(len(full_list))
  
  
          ax.scatter(basin_list,full_list,c=colours,edgecolors='k',s=15,cmap = plt.cm.PuOr,zorder = 100)
  
  
  
      # set the axes labels
      ax.set_xlabel('basin')
      ax.set_ylabel('$m/n$')
      # Save the figure
      ImageName = DataDirectory+fname_prefix+'_MOverNbasins.'+FigFormat
      plt.savefig(ImageName, format=FigFormat, dpi=300)
      fig.clf()
  
  def plot_MCMC_analysis(DataDirectory,fname_prefix,basin_list=[],FigFormat='png',size_format='ESURF', parallel=False):
      """
      This function makes a plot of the MCMC uncertainty analysis for the MLE collinearity
      NOTE: haven't set this up for basins in parallel yet (FJC 19/10/17)
  
      Args:
          DataDirectory (str): the data directory
          fname_prefix (str): the DEM name without extension
          basin_list: list of basins to analyse. If none are passed then all basins are plotted.
          FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
          size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
          parallel (bool): If true the data exists in multiple files for each basin and must be recombined (Default is False)
  
      Author: FJC
      """
      # check if a directory exists for the chi plots. If not then make it.
      MCMC_directory = DataDirectory+'MCMC_plots/'
      if not os.path.isdir(MCMC_directory):
          os.makedirs(MCMC_directory)
  
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # make a figure
      fig = makefigure(size_format)
  
      gs = plt.GridSpec(100,100,bottom=0.15,left=0.1,right=0.85,top=0.9)
      ax = fig.add_subplot(gs[5:100,10:95])
  
      # get the basin info
      if not parallel:
        basin_df = Helper.ReadBasinInfoCSV(DataDirectory,fname_prefix)
        basin_stats_df = Helper.ReadBasinStatsCSV(DataDirectory,fname_prefix)
      else:
        basin_df = Helper.AppendBasinInfoCSVs(DataDirectory,fname_prefix)
        basin_stats_df = Helper.AppendBasinStatsCSVs(DataDirectory,fname_prefix)
  
      basin_keys = basin_df['basin_key']
  
      MOverNDict = SimpleMaxMLECheck(basin_stats_df)
  
      if basin_list == []:
          print ("You didn't give me a basin list so I'm going to plot all of them!")
          basin_list = basin_keys
  
      for basin in basin_list:
          # read in the chain csv file
          chain_df = Helper.ReadChainCSV(DataDirectory,fname_prefix,int(basin))
  
          best_fit_movern = MOverNDict[int(basin)]
          print (("The best fit m/n is:", best_fit_movern))
          n_iterations = len(chain_df['i'])
          #print n_iterations
          best_fit_line = [best_fit_movern] * len(chain_df['i'])
  
          # find the accepted values
          chain_df['AcceptedMask'] = chain_df['NAccepted'].diff()
          masked_df = chain_df.where(chain_df['AcceptedMask'] == 1.0)
          #print chain_df
  
          # plot the parameter with number of iterations
          ax.plot(masked_df['i'], masked_df['movern_New'], c='0.5', lw=0.5, zorder=1)
          ax.plot(chain_df['i'], best_fit_line, 'k--', zorder=100)
  
          #set plot labels
          ax.set_xlabel('N iterations')
          ax.set_ylabel('$m/n$')
          ax.set_title('Basin '+str(basin))
  
          # save the figure
          ImageName = MCMC_directory+fname_prefix+'_MCMC_basin' +str(basin)+'.'+FigFormat
          plt.savefig(ImageName, format=FigFormat, dpi=300)
          ax.cla()
  
  def PlotMCPointsUncertainty(DataDirectory,fname_prefix, basin_list=[0], FigFormat='png',size_format='ESURF', start_movern=0.2,d_movern=0.1,n_movern=7,parallel=False):
      """
      This function makes a plot showing the range in m/n calculated
      using the MC points analysis. Makes a separate plot for each basin.
  
      Args:
          DataDirectory (str): the data directory
          fname_prefix (str): the DEM name without extension
          basin_list: list of basins to be analysed. If empty then will analyse all of them.
          FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
          size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
  
      Returns:
          plot of uncertainty in m/n
  
      Author:
          FJC
      """
      import matplotlib.patches as patches
  
      # check if a directory exists for the chi plots. If not then make it.
      points_directory = DataDirectory+'MC_points_plots/'
      if not os.path.isdir(points_directory):
          os.makedirs(points_directory)
  
      # get the basin info
      if not parallel:
          basin_df = Helper.ReadBasinInfoCSV(DataDirectory,fname_prefix)
      else:
          basin_df = Helper.AppendBasinInfoCSVs(DataDirectory,fname_prefix)
  
      basin_keys = basin_df['basin_key']
  
      if basin_list == []:
          print ("You didn't give me a basin list so I'm going to plot all of them!")
          basin_list = basin_keys
  
      # get the uncertainty df
      if not parallel:
          PointsChiBasinDF = Helper.ReadMCPointsCSV(DataDirectory,fname_prefix)
      else:
          PointsChiBasinDF = Helper.AppendBasinPointCSVs(DataDirectory,fname_prefix)
  
      PointsChiBasinDF = PointsChiBasinDF[PointsChiBasinDF['basin_key'].isin(basin_list)]
  
      UncertaintyDF = GetMOverNRangeMCPoints(PointsChiBasinDF,start_movern,d_movern,n_movern)
  
      # set up the plot
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # make a figure
      fig = makefigure(size_format)
  
      gs = plt.GridSpec(100,100,bottom=0.15,left=0.1,right=0.85,top=0.9)
      ax = fig.add_subplot(gs[5:100,10:95])
  
      if not parallel:
          BasinStatsDF = Helper.ReadBasinStatsCSV(DataDirectory,fname_prefix)
      else:
          BasinStatsDF = Helper.AppendBasinCSVs(DataDirectory,fname_prefix)
      # best fit moverns
      best_fit_moverns = SimpleMaxMLECheck(BasinStatsDF)
  
      for basin_key in basin_keys:
          ThisBasinDF = PointsChiBasinDF[PointsChiBasinDF['basin_key'] == basin_key]
  
          # get the m/ns tested
          end_movern = start_movern+d_movern*(n_movern-1)
          all_moverns = np.linspace(start_movern,end_movern,n_movern)
  
          # get the median m/ns for this basin
          Medians = (ThisBasinDF.filter(regex='median')).values.tolist()[-1]
          FirstQs = (ThisBasinDF.filter(regex='FQ')).values.tolist()[-1]
          ThirdQs = (ThisBasinDF.filter(regex='TQ')).values.tolist()[-1]
  
          # now get the threshold value
          ThisUncertaintyDF = UncertaintyDF[UncertaintyDF['basin_key'] == basin_key]
          print (ThisUncertaintyDF)
  
          #plot the median and quartiles
          ax.plot(all_moverns,Medians,c="k",lw=1,zorder=2)
          ax.plot(all_moverns,FirstQs,c="k", lw=0.5, ls="--",zorder=2)
          ax.plot(all_moverns,ThirdQs,c="k", lw=0.5, ls="--",zorder=2)
          ax.fill_between(all_moverns, FirstQs, ThirdQs, color="0.8",alpha=0.75,zorder=1)
  
          #add a line for the threshold
          threshold = ThisUncertaintyDF.iloc[0]['FirstQ_threshold']
          min_movern = ThisUncertaintyDF.iloc[0]['Min_MOverNs']
          max_movern = ThisUncertaintyDF.iloc[0]['Max_MOverNs']
          threshold_moverns = np.linspace(min_movern,max_movern, 4)
          threshold_MLEs = np.full((n_movern,),threshold)
          ax.plot(all_moverns,threshold_MLEs,c='r',zorder=3, ls="--",lw=1)
  
          # add shaded background over range of m/n values
          ax.axvspan(min_movern,max_movern, alpha=0.1, color='red',zorder=0.2)
  
          # add arrow at best fit m/n
          # get the limits for the arrow
          max_MLE = max(Medians)
          min_MLE = min(FirstQs)
          dy = (max_MLE-min_MLE)/8
          spacing = 1.5
          # add arrow at best fit m/n
          ax.add_patch(
              patches.Arrow(
                  best_fit_moverns[basin_key], #x
                  max_MLE-(dy*spacing), #y
                  0, #dx
                  dy, #dy
                  width = 0.05,
                  facecolor = 'k',
                  edgecolor = 'k')
              )
          ax.text(best_fit_moverns[basin_key]-0.075, max_MLE-1.5*(dy*spacing), "Best-fit "+r"$\theta$",fontsize=8)
  
          # set the axes labels
          ax.set_xlabel(r'$\theta$')
          ax.set_ylabel('MLE')
          #ax.set_title('Basin '+str(basin_key))
          ax.set_xlim(start_movern,end_movern)
          ax.set_ylim(min(FirstQs),max(ThirdQs)+0.0005)
  
          # save the figure
          ImageName = points_directory+fname_prefix+'_MC_points' +str(basin_key)+'.'+FigFormat
          plt.savefig(ImageName, format=FigFormat, dpi=300)
          ax.cla()
  
      fig.clf()
      plt.close(fig)
  
  #=============================================================================
  # SENSITIVITY FUNCTIONS
  # Functions that make plots of sensitivity tests on the m/n analysis
  #=============================================================================
  def PlotSensitivityResultsSigma(DataDirectory,fname_prefix, FigFormat = "png", size_format = "ESURF", movern_method = "all"):
      """
      This function makes a plot of the results of a sensitivity analysis on sigma in the MLE method
      of calculating m/n.
      You need to specify the base directory - will look in every folder in this base directory
      with the prefix "Chi_analysis_sigma_" and get the m/n analysis in each sub-directory.
      NOTE - not set up for parallel (FJC 19/10/17)
  
      Args:
          DataDirectory (str): the root data directory
          fname_prefix (str): the DEM name without extension
          FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
          size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
          movern_method (str): the movern method you want to test. Can be "all" or "points"
  
      Author: FJC
      """
      import os
  
      # Try doing this as a dataframe and masking?!?!
      columns = ["sigma", "basin_key", "m_over_n"]
      combined_DF = pd.DataFrame(columns=columns)
  
      # loop through each sub-directory with the sensitivity results
      MLE_str = "sigma_"
      for subdir, dirs, files in os.walk(DataDirectory):
          for dir in dirs:
              if MLE_str in dir:
                  this_dir = DataDirectory+"/"+dir+'/'
                  print (this_dir)
  
                  this_df = pd.DataFrame(columns=columns)
  
                  # get this value of sigma
                  this_sigma = (dir.split("_"))[-1]
  
                  # get the best fit m/n dataframe
                  if movern_method == 'all':
                      BasinDF = Helper.ReadBasinStatsCSV(this_dir,fname_prefix)
                  elif movern_method == 'points':
                      BasinDF = Helper.ReadBasinStatsPointCSV(this_dir,fname_prefix)
                  else:
                      print ("You didn't specify a correct m/n method, you need to specify either 'all' or 'points'")
                  MOverNDict = SimpleMaxMLECheck(BasinDF)
                  this_df['sigma'] = [int(this_sigma)] * len(MOverNDict)
                  this_df['basin_key'] = MOverNDict.keys()
                  this_df['m_over_n'] = MOverNDict.values()
                  combined_DF = combined_DF.append(this_df, ignore_index=True)
  
      # Set up fonts for plots
      label_size = 10
      rcParams['font.family'] = 'sans-serif'
      rcParams['font.sans-serif'] = ['Liberation Sans']
      rcParams['font.size'] = label_size
  
      # make a figure
      fig = makefigure(size_format)
  
      gs = plt.GridSpec(100,100,bottom=0.15,left=0.1,right=0.85,top=0.9)
      ax = fig.add_subplot(gs[5:100,10:95])
      #ax.grid(zorder=-100)
  
      # now get the basin keys
      basin_keys = combined_DF['basin_key'].unique()
      basin_keys = [int(x) for x in basin_keys]
  
      keys = []
      sigmas = []
      # loop through the basin keys and plot sigma and m/n for each basin
      for key in basin_keys:
          this_df = combined_DF.loc[combined_DF['basin_key'] == key]
          #sort the data for plotting
          this_df.sort_values('sigma',inplace=True)
          this_df = this_df.loc[combined_DF['m_over_n'] > 0.2]
          if not this_df.empty:
              keys.append(int(key))
              sigmas.append(this_df['sigma'].iloc[0])
  
      ax.scatter(keys, sigmas, c='0.75', edgecolor='k')
      #
      # set the axes labels
      ax.set_xlabel('Basin key')
      ax.set_ylabel('$\sigma$ where $m/n$ is invariant')
      ax.set_xticks(basin_keys)
  
      # Save the figure
      ImageName = DataDirectory+fname_prefix+'_sensitivity.'+FigFormat
      plt.savefig(ImageName, format=FigFormat, dpi=300)
      fig.clf()
  
  # def PlotSensitivityResultsMLEWithSigma(DataDirectory,fname_prefix, FigFormat = "png", size_format = "ESURF"):
  #     """
  #     This function makes a plot of the results of a sensitivity analysis on sigma in the MLE method
  #     of calculating m/n.
  #     You need to specify the base directory - will look in every folder in this base directory
  #     with the prefix "sigma_" and get the m/n analysis in each sub-directory.
  #     Makes a plot of MLE with increasing sigma to show where it becomes invariant
  #
  #     Args:
  #         DataDirectory (str): the root data directory
  #         fname_prefix (str): the DEM name without extension
  #         FigFormat (str): The format of the figure. Usually 'png' or 'pdf'. If "show" then it calls the matplotlib show() command.
  #         size_format (str): Can be "big" (16 inches wide), "geomorphology" (6.25 inches wide), or "ESURF" (4.92 inches wide) (defualt esurf).
  #
  #     Author: FJC
  #     """
  #     import os
  #
  #     # Try doing this as a dataframe and masking?!?!
  #     columns = ["sigma", "basin_key", "m_over_n"]
  #     combined_DF = pd.DataFrame(columns=columns)
  #
  #     # loop through each sub-directory with the sensitivity results
  #     MLE_str = "sigma_"
  #     for subdir, dirs, files in os.walk(DataDirectory):
  #         for dir in dirs:
  #             if MLE_str in dir:
  #                 this_dir = DataDirectory+"/"+dir+'/'
  #                 print this_dir
  #
  #                 this_df = pd.DataFrame(columns=columns)
  #
  #                 # get this value of sigma
  #                 this_sigma = (dir.split("_"))[-1]
  #
  #                 # get the best fit m/n dataframe
  #                 BasinDF = Helper.ReadBasinStatsCSV(this_dir,fname_prefix)
  #                 MOverNDict = SimpleMaxMLECheck(BasinDF)
  #                 this_df['sigma'] = [int(this_sigma)] * len(MOverNDict)
  #                 this_df['basin_key'] = MOverNDict.keys()
  #                 this_df['m_over_n'] = MOverNDict.values()
  #                 combined_DF = combined_DF.append(this_df, ignore_index=True)
  #
  #     # Set up fonts for plots
  #     label_size = 10
  #     rcParams['font.family'] = 'sans-serif'
  #     rcParams['font.sans-serif'] = ['Liberation Sans']
  #     rcParams['font.size'] = label_size
  #
  #     # make a figure
  #     if size_format == "geomorphology":
  #         fig = plt.figure(1, facecolor='white',figsize=(6.25,3.5))
  #         #l_pad = -40
  #     elif size_format == "big":
  #         fig = plt.figure(1, facecolor='white',figsize=(16,9))
  #         #l_pad = -50
  #     else:
  #         fig = plt.figure(1, facecolor='white',figsize=(4.92126,3.2))
  #         #l_pad = -35
  #
  #     gs = plt.GridSpec(100,100,bottom=0.15,left=0.1,right=0.85,top=0.9)
  #     ax = fig.add_subplot(gs[5:100,10:95])
  #     #ax.grid(zorder=-100)
  #
  #     # now get the basin keys
  #     basin_keys = combined_DF['basin_key'].unique()
  #     basin_keys = [int(x) for x in basin_keys]
  #
  #     keys = []
  #     sigmas = []
  #     MLEs = []
  #     # loop through the basin keys and plot sigma and m/n for each basin
  #     for key in basin_keys:
  #         this_df = combined_DF.loc[combined_DF['basin_key'] == key]
  #         #sort the data for plotting
  #         this_df.sort_values('sigma',inplace=True)
  #         this_df = this_df.loc[combined_DF['m_over_n'] > 0.2]
  #         if not this_df.empty:
  #             keys.append(int(key))
  #             sigmas.append(this_df['sigma'].iloc[0])
  #
  #     ax.scatter(keys, sigmas, c='0.75', edgecolor='k')
  #     #
  #     # set the axes labels
  #     ax.set_xlabel('Basin key')
  #     ax.set_ylabel('$\sigma$ where $m/n$ is invariant')
  #     ax.set_xticks(basin_keys)
  #
  #     # Save the figure
  #     ImageName = DataDirectory+fname_prefix+'_sensitivity.'+FigFormat
  #     plt.savefig(ImageName, format=FigFormat, dpi=300)
  #     fig.clf()
  