text: #!/usr/local/sci/bin/python
  #*****************************
  #
  # Detect outliers from station and neighbour
  #
  #
  #************************************************************************
  #                    SVN Info
  #$Rev:: 116                                           $:  Revision of last commit
  #$Author:: rdunn                                      $:  Author of last commit
  #$Date:: 2017-01-30 15:24:24 +0000 (Mon, 30 Jan 2017) $:  Date of last commit
  #************************************************************************
  
  
  import numpy as np
  import scipy as sp
  import math
  import os
  import gc
  
  # RJHD routines
  import qc_utils as utils
  import netcdf_procs as ncdfp
  
  
  N_NEIGHBOURS = 10
  
  UNFLAG_COL_DICT = {"spike":{"temperatures":27 ,"dewpoints":28, "slp":29},\
                         "climatological":{"temperatures":24 ,"dewpoints":25, "slp":26},\
                         "odd":{"temperatures":54 ,"dewpoints":55, "slp":57},\
                         "gap":{"slp":7},\
                         "dpd":{"dewpoints":31}}
  
  #*******************************************************
  def get_neighbour_parameters(station_info):
      '''
      Return two big symmetrical arrays of the station separations and bearings
      '''
  
      distances = np.zeros([len(station_info),len(station_info)])
      angles = np.zeros([len(station_info),len(station_info)])
      
      all_lats = np.array(station_info[:,1], dtype=float)
      all_lons = np.array(station_info[:,2], dtype=float)
  
  
      for st, stat in enumerate(station_info):
          distances[st,:], angles[st,:] = utils.get_dist_and_bearing([all_lats[st],all_lons[st]],[all_lats,all_lons])
          distances[st,st] = 0. # by hand set the diagonal to zero
          angles[st,st] = 0. # by hand set the diagonal to zero
  
      return distances, angles # get_neighbour_parameters
  
  #*******************************************************
  def get_neighbours(station_loc, st_elev, distances, bearings, elevations, sep_limit = 500., elev_limit = 300.):
      '''
      Get the neighbours for the station using distance, angles and elevations
      Will return 10 or less neighbours.
  
      :param int station_loc: sequence number of station
      :param float st_elev: station elevation
      :param array distances: distances of all other stations
      :param array bearings: bearings of all other stations
      :param array elevations: elevations of all other stations
      :param float sep_limit: separation limit (default 500km)
      :param float elev_limit: elevation limit (default 300m)
  
      :returns: final set of neighbours as numbers in station list
      '''
  
  
      # get bearings into quadrants to make searching easier
      quadrants = np.zeros(bearings.shape)
      for a,angle in enumerate([270,180,90,0]):
          locs = np.where((bearings >= angle) & (bearings < angle+90.))
          quadrants[locs] = a+1
  
      # start by sorting by distance
      ordering = np.argsort(distances)
  
      neighbour_locations = np.array([], dtype = int)
      neighbour_quadrants = np.array([])
      for index in ordering:
  
          if index == station_loc:
              # don't choose itself!
              continue
  
          # if within range of distance and elevation
          if distances[index] <= sep_limit:
              if np.abs(elevations[index] - st_elev) <= elev_limit:
  
                  neighbour_locations=np.append(neighbour_locations, index)
                  neighbour_quadrants = np.append(neighbour_quadrants, quadrants[index]) # combine into one data structure
                                 
                  # if have more than 10, see if angularly spread
                  if len(neighbour_locations) >= N_NEIGHBOURS:
                      
                      locs1 = neighbour_locations[neighbour_quadrants == 1]
                      locs2 = neighbour_locations[neighbour_quadrants == 2]
                      locs3 = neighbour_locations[neighbour_quadrants == 3]
                      locs4 = neighbour_locations[neighbour_quadrants == 4]
  
                      # if at least two per quadrant then stop looking
                      if len(locs1) >= 2 and len(locs2) >= 2 and len(locs3) >= 2 and len(locs4) >= 2:
                          break
  
  
      # if have more than 10
      if len(neighbour_locations) >= N_NEIGHBOURS:
  
          # now go through to ensure have only 10, and two per quadrant
          final_locs = np.concatenate((locs1[:2], locs2[:2], locs3[:2], locs4[:2]), axis = 0).reshape(-1)
  
          # and add the rest in order of distance
          for index in neighbour_locations:
              if index not in final_locs:
                  final_locs = np.append(final_locs, index)
  
              if len(final_locs) == N_NEIGHBOURS:
                  break
  
          return final_locs
      # else just return the ones found
      else:
          return neighbour_locations # get_neighbours
  
  #*******************************************************
  def get_all_neighbours(station_loc, st_elev, distances, bearings, elevations, sep_limit = 500., elev_limit = 300., max_neighbours = 20):
      '''
      Get the neighbours for the station using distance, angles and elevations
      Returns all neighbours within distance and elevation ranges 
  
      :param int station_loc: sequence number of station
      :param float st_elev: station elevation
      :param array distances: distances of all other stations
      :param array bearings: bearings of all other stations
      :param array elevations: elevations of all other stations
      :param float sep_limit: separation limit (default 500km)
      :param float elev_limit: elevation limit (default 300m)
      :param int max_neighbours: maximum number to return (default 50)
  
      :returns: final set of neighbours as numbers in station list
      '''
  
  
      # get bearings into quadrants to make searching easier
      quadrants = np.zeros(bearings.shape)
      for a,angle in enumerate([270,180,90,0]):
          locs = np.where((bearings >= angle) & (bearings < angle+90.))
          quadrants[locs] = a+1
  
      # start by sorting by distance
      ordering = np.argsort(distances)
  
      neighbour_locations = np.array([], dtype = int)
      neighbour_quadrants = np.array([])
      for index in ordering:
  
          if index == station_loc:
              # don't choose itself!
              continue
  
          # if within range of distance and elevation
          if distances[index] <= sep_limit:
              if np.abs(elevations[index] - st_elev) <= elev_limit:
  
                  neighbour_locations=np.append(neighbour_locations, index)
                  neighbour_quadrants = np.append(neighbour_quadrants, quadrants[index]) # combine into one data structure
                                 
  
      return neighbour_locations[:max_neighbours], neighbour_quadrants[:max_neighbours] # get_all_neighbours
  
  #*******************************************************
  def hourly_daily_anomalies(timeseries, obs_per_day = 6):
      '''
      Format the time series to get allow for sensible correlations
  
       - Process into 24h x N_days
       - Obtain daily average and hence hourly anomalies from daily average (removes annual cycle)
       - Obtain hourly average of anomalies to get double-anomalies (removes diurnal cycle)
  
      :param array timeseries: data to be processed in 1-D array
      :param int obs_per_day: number of observations per 24hr to get a daily average to process (default = 6)
  
      :returns: anomalies - timeseries with removed annual and diurnal cycles.
      '''
  
      timeseries = timeseries.reshape((-1,24))
  
      daily_mean = np.ma.mean(timeseries, axis = 1)
  
      not_mask_count = timeseries.count(axis = 1)
  
      daily_mean = np.ma.masked_where(not_mask_count < obs_per_day, daily_mean) # completeness check - only for correlations
  
      hourly_anomalies = timeseries - np.ma.repeat(daily_mean, 24).reshape((-1,24)) # removed annual cycle
  
      hourly_mean = np.ma.mean(hourly_anomalies, axis = 0)
  
      anomalies = hourly_anomalies - np.tile(hourly_mean, timeseries.shape[0]).reshape((-1,24)) # removed diurnal cycle
  
      return anomalies.ravel() # hourly_daily_anomalies
  
  #*******************************************************
  def select_neighbours(station, variable, neighbour_info, neighbours, neighbour_distances, neighbour_quadrants, data_locs, datastart, dataend, logfile, second = False, diagnostics = False, plots = False):
      '''
      From the list of nearby stations select the ones which will be good neighours for the test.
      Select on basis of correlation, overlap of data points and bearing (quadrants)
      
      :param object station: station object
      :param str variable: which variable to proces
      :param array neighbour_info: array of ID, lat, lon and elev
      :param array neighbours: which station sequence numbers are the nearby stations
      :param array neighbour_distances: distances to nearby stations
      :param array neighbour_quadrants: bearings to nearby stations (in 90deg bins)
      :param array data_locs: path to data files
      :param datetime datastart: start of data set
      :param datetime dataend: end of data set
      :param file logfile: logfile to store outputs
      :param boolean second: second run through
      :param boolean diagnostics: output diagnostic information
      :param boolean plots: make a plot
  
      :returns: final_locs - array of station sequence numbers to use.
      '''
  
      first = not second
  
      # set up storage arrays
      n_correlations = np.zeros(len(neighbours))
      n_distances = np.zeros(len(neighbours))
      n_quadrants = np.zeros(len(neighbours))
      n_overlaps = np.zeros(len(neighbours))
      combined_score = np.zeros(len(neighbours))
  
      # get station data
      st_var = getattr(station, variable)
      st_anomalies = hourly_daily_anomalies(st_var.data[:])
  
      # go through initial list and extract correlations and overlaps
      for nn, nn_loc in enumerate(neighbours):
  
          n_details = neighbour_info[nn]
          neigh = utils.Station(n_details[0], float(n_details[1]), float(n_details[2]), float(n_details[3]))
  
          if first:
              ncdfp.read(os.path.join(data_locs, neigh.id + "_internal.nc"), neigh, [variable], diagnostics = diagnostics, read_input_station_id = False)
          elif second:
              ncdfp.read(os.path.join(data_locs, neigh.id + "_internal2.nc"), neigh, [variable], diagnostics = diagnostics, read_input_station_id = False)
  
          dummy = utils.create_fulltimes(neigh, [variable], datastart, dataend, [], do_input_station_id = False)
  
          # get the correlations of data to this neighbour
          neigh_var = getattr(neigh, variable)
          neigh_anomalies = hourly_daily_anomalies(neigh_var.data[:])
          # correlation = np.ma.corrcoef(neigh_var.data, st_var.data)[1,0]
          correlation = np.ma.corrcoef(neigh_anomalies, st_anomalies)[1,0]
  
          overlap = len(np.where(np.logical_or(neigh_var.data.mask, st_var.data.mask) == False)[0])/float(len(st_var.data.compressed()))
  
          if not math.isnan(correlation):
              n_correlations[nn] = correlation
              n_overlaps[nn] = overlap
              combined_score[nn] = correlation + overlap
              n_distances[nn] = neighbour_distances[nn]
              n_quadrants[nn] = neighbour_quadrants[nn]
              
          # clear up to save memory
          del dummy
          del neigh_var
          del neigh_anomalies
          gc.collect()
      # sort in order of the combination of correlation and overlap
      sort_order = np.argsort(combined_score)[::-1]
  
      # and select the best 10
      # final_selection = neighbours[sort_order][:10]
  
      # sort out the quadrants
                      
      locs1 = neighbours[sort_order][n_quadrants[sort_order] == 1]
      locs2 = neighbours[sort_order][n_quadrants[sort_order] == 2]
      locs3 = neighbours[sort_order][n_quadrants[sort_order] == 3]
      locs4 = neighbours[sort_order][n_quadrants[sort_order] == 4]
  
      final_locs = np.concatenate((locs1[:2], locs2[:2], locs3[:2], locs4[:2]), axis = 0).reshape(-1)
  
      # and add the rest in order of combined score
      for index in neighbours[sort_order]:
          if index not in final_locs:
              final_locs = np.append(final_locs, index)
              
          if len(final_locs) == N_NEIGHBOURS:
              break
  
      # output table showing distances, correlations, overlaps, the combined score and which ones were selected
      if plots or diagnostics:
          print "{:14s} {:10s} {:10s} {:10s} {:10s} {:10s} {:10s} {:10s}".format("Neighbour","Distance","Elevation", "Correl'n", "Overlap", "Combined", "Quadrant","Selected")
      else:
          logfile.write("{:14s} {:10s} {:10s} {:10s} {:10s} {:10s} {:10s} {:10s}\n".format("Neighbour","Distance","Elevation", "Correl'n", "Overlap", "Combined", "Quadrant","Selected")
  )
  
      selected_correlations = []
      selected_overlaps = []
      for nn, nn_loc in enumerate(neighbours[sort_order]):
  
          selected = ""
          if nn_loc in final_locs: 
              selected = "Y"
              if plots:
                  selected_correlations += [n_correlations[sort_order[nn]]]
                  selected_overlaps += [n_overlaps[sort_order[nn]]]
  
  
          neigh_details = neighbour_info[sort_order][nn]
          if plots or diagnostics:
              print "{:14s} {:10.1f} {:10.1f} {:10.5f} {:10.3f} {:10.3f} {:10.0f} {:10s}".format(neigh_details[0], n_distances[sort_order][nn], float(neigh_details[3]), n_correlations[sort_order][nn], n_overlaps[sort_order][nn], combined_score[sort_order][nn], n_quadrants[sort_order][nn], selected)
          else:
              logfile.write("{:14s} {:10.1f} {:10.1f} {:10.5f} {:10.3f} {:10.3f} {:10.0f} {:10s}\n".format(neigh_details[0], n_distances[sort_order][nn], float(neigh_details[3]), n_correlations[sort_order][nn], n_overlaps[sort_order][nn], combined_score[sort_order][nn], n_quadrants[sort_order][nn], selected))
  
              
      # plot of correlations and overlaps, with selected stations highlighted
      if plots:
          import matplotlib.pyplot as plt
  
          plt.clf()
          plt.plot(n_correlations, n_overlaps, 'bo')
          plt.plot(selected_correlations, selected_overlaps, 'ro')
          plt.xlabel("correlations")
          plt.ylabel("data overlap")
          plt.title("{} - {}".format(station.id, variable))
          plt.show()
  
      return final_locs # select_neighbours
  
  #*******************************************************
  def plot_target_neigh_diffs_dist(differences, iqr):
      '''
      Plot the distribution of target-neighbour differences
      
      :param array differences: masked difference array
      :param float iqr: inter quartile range of differences
  
      :returns: 
      '''
      import matplotlib.pyplot as plt
      
      plt.clf()
      
      bins, bincenters = utils.create_bins(differences.compressed(), 1.0)
      
      hist, binEdges = np.histogram(differences.compressed(), bins=bins)
      plot_hist = np.array([float(x) if x != 0 else 1e-1 for x in hist])
      plt.step(bincenters, plot_hist, 'k-', label = 'observations', where='mid')
      
      fit = utils.fit_gaussian(bincenters, hist, max(hist), mu=np.mean(differences.compressed()), sig = np.std(differences.compressed()))
      plot_gaussian = utils.gaussian(bincenters, fit)
      plt.plot(bincenters, plot_gaussian, 'b-', label = 'Gaussian fit')
      
      plt.axvline(5.*iqr, c = 'r')
      plt.axvline(-5.*iqr, c = 'r')
      
      print "only shows lowest of monthly IQRs"
  
      plt.ylabel("Frequency")
      plt.gca().set_yscale('log')
      plt.ylim([0.1,2*max(hist)])
      
      plt.show()      
  
      return # plot_target_neigh_diffs_dist
  
  
  #*******************************************************
  def detect(station, neighbour, variable, flags, neighbour_count, start, end, distance = 0, diagnostics = False, plots = False):
      '''
      Detect which observations are outliers
  
      :param MetVar station: station object (target)
      :param MetVar neighbour: station object (neighbour)
      :param string variable: which variable to process
      :param array flags: array to store how many neighbours thing each obs is bad
      :param array neighbour_count: how many neighbours present at each obs
      :param datetime start: start of dataset
      :param datetime end: end of dataset
      :param int distance: separation of target and neighbour
      :param bool diagnostics: extra output
      :param bool plots: make figures
  
      :returns: None
      '''
  
      FILTERING_FLAG_COL = {"temperatures":[0,1,4,5,8,12,16,20,27,41,44,58],
                            "dewpoints":[0,2,4,6,8,9,13,17,21,28,30,31,32,42,45,59],
                            "slp":[0,3,4,7,11,15,19,23,29,43,46,60],
                            "windspeeds":[0,4,10,14,18,22,56,62,63,64]} # not used, but ready for it.
      
      
      st_var = getattr(station, variable)
      neigh_var = getattr(neighbour, variable)
  
      # filter by flags - not all (no Climatological [24,25], or Odd cluster [54,55,56,57]), T record check not in D, 
      total_flags = np.sum(station.qc_flags[:,FILTERING_FLAG_COL[variable]], axis = 1)
      st_filtered = np.ma.masked_where(total_flags == 1, st_var.data)
      neigh_filtered = np.ma.masked_where(total_flags == 1, neigh_var.data)
  
      # match the observation times
      match = np.where(np.logical_and((st_filtered.data != st_var.mdi), (neigh_filtered.data != neigh_var.mdi)))
  
      month_ranges = utils.month_starts_in_pairs(start, end).reshape(-1,12,2) # in year-long sets of pairs.        
  
      if len(match[0]) >= 100:
  
          neighbour_count[match] += 1 # number of neighbours with data present
  
          differences = np.ma.zeros(len(st_filtered))
          differences.fill(st_var.mdi)
          differences.mask = True
  
          differences[match] = st_filtered.data[match] - neigh_filtered.data[match]
          differences.mask[match] = False
  
          all_iqrs = np.zeros(len(differences))
          # get monthly IQR values
          for month in range(12):
    
              this_month, dummy1, dummy2 = utils.concatenate_months(month_ranges[:,month,:], differences, hours = False)
  
              if len(this_month.compressed()) > 4:
  
                  iqr = utils.IQR(this_month.compressed())
                  if iqr <= 2.: iqr = 2.
              else:
                  iqr = 2.
  
              # and copy back into the array
              for year in month_ranges[:,month,:]:
                  all_iqrs[year[0]:year[1]] = iqr
  
          if plots:
              plot_target_neigh_diffs_dist(differences, min(all_iqrs))
  
          dubious = np.ma.where(np.ma.abs(differences) > 5. * all_iqrs)
  
          if len(dubious[0]) >= 1.:
  
              if variable == "slp":
                  # check if they are storms
                  positive = np.ma.where(differences > 5. * iqr)
                  negative = np.ma.where(differences < -5. * iqr)
                  
                  # if majority negative (2/3) and separation > 100
  
                  if (distance > 100.) and (float(len(positive[0]))/len(dubious[0]) < 0.333):
  
                      if len(positive[0]) > 0:
                          flags[positive] += 1
                      if len(negative[0]) > 0:
                          neighbour_count[match] -= 1
                      
                  else:
                      flags[dubious] += 1
              else:
  
                  flags[dubious] += 1
  
      return # detect
  
  #*******************************************************
  def plot_outlier(station, variable, outlier_locs, all_data, datastart):
      '''
      Plot the outlier location (either to flag or unflag) with the target and all neighbours
  
      :param MetVar station: station object
      :param str variable: variable to process
      :param array outlier_locs: locations which are outliers to plot
      :param array all_data: all data from neighbours for plotting
      :param datetime datastart: start of dataset
  
      :returns: None
      '''
  
      import matplotlib.pyplot as plt
      import datetime as dt
      
      YLABELS = {"temperatures":"Temperature (C)", "dewpoints":"Dewpoints (C)", "slp":"SLP (hPa)", "windspeeds":"Wind Speed (m/s)"}
      extra = 48 # hours
  
      indata = getattr(station, variable).data
      
      for location in outlier_locs:
          
          plot_times = utils.times_hours_to_datetime(station.time.data[location-extra: location+extra], datastart)
          
          plt.clf()
          plt.plot(plot_times, indata[location-extra:location+extra], 'bo', ls='-')
          plt.plot(plot_times[extra], indata[location], 'ro', markersize=10)
          
          for nn in range(all_data.shape[0]):
              plt.plot(plot_times,all_data[nn,location-extra:location+extra] , c='0.5', ls='-')
              
          plt.ylabel(YLABELS[variable])
          plt.title("{:s} {:s}".format(station.id,dt.datetime.strftime(plot_times[extra], "%d/%m/%Y")))
          plt.show()
                      
      return # plot_outlier
  
  #*******************************************************
  def unflagging_locs(differences, flags, neigh_count, dpd_count = [], flag_value = 1):
      '''
      Return locations where flags to be set to zero
      Currently has deliberate bug, but this is to match IDL
  
      :param array differences: normalised differences of target - median of neighbours
      :param array flags: flag array to check where to unset
      :param array neigh_count: number of neighbour obs at each timestamp
      :param list dpd_count: number of DPD flags set at each timestamp
      :param int flag_value: can use to look for tentative flags
  
      :returns:
         unset_locs - list of locations to unset
      '''
  
      sufficient_neighbours, = np.where(neigh_count >= 3)
  
      if len(sufficient_neighbours) > 0:
  
          flag_locs, = np.where(flags[sufficient_neighbours] == flag_value)
  
          if len(flag_locs) > 0:
  
              if dpd_count != []:
                  # if at least 2/3 of neighbours also have DPD flag set, then unset these
                  dpd_proportion = dpd_count/neigh_count
                  
                  unset_locs = np.where(dpd_proportion[sufficient_neighbours][flag_locs] >= 2./3.)
                  
              else:
                  
                  unset_locs = np.where(differences[sufficient_neighbours][flag_locs] <= 4.5)
                  
              if len(unset_locs) > 0:
                  # have some values that fall within 4.5MAD of median or where DPD is
                  #  present in sufficiently many neighbours, remove these flags
  
                  return sufficient_neighbours[flag_locs][unset_locs]
              
  
      return [] # unflagging
  
  #*******************************************************
  def bn_median(masked_array, axis=None):
      """
      https://github.com/astropy/ccdproc/blob/122cdbd5713140174f057eaa8fdb6f9ce03312df/docs/ccdproc/bottleneck_example.rst
      Perform fast median on masked array
  
      Parameters
  
      masked_array : `numpy.ma.masked_array`
          Array of which to find the median.
  
      axis : int, optional
          Axis along which to perform the median. Default is to find the median of
          the flattened array.
      """
      import numpy as np
      import bottleneck as bn
      data = masked_array.filled(fill_value=np.NaN)
      med = bn.nanmedian(data, axis=axis)
      # construct a masked array result, setting the mask from any NaN entries
      return np.ma.array(med, mask=np.isnan(med))
  
  #*******************************************************
  def median_absolute_deviation(a, axis=None):
      """from astropy.stats, converted to use np.ma - 13-Oct-2014 RJHD
  
      Compute the median absolute deviation.
  
      Returns the median absolute deviation (MAD) of the array elements.
      The MAD is defined as ``median(abs(a - median(a)))``.
  
      Parameters
      ----------
      a : array_like
          Input array or object that can be converted to an array.
      axis : int, optional
          Axis along which the medians are computed. The default (axis=None)
          is to compute the median along a flattened version of the array.
  
      Returns
      -------
      median_absolute_deviation : ndarray
          A new array holding the result. If the input contains
          integers, or floats of smaller precision than 64, then the output
          data-type is float64.  Otherwise, the output data-type is the same
          as that of the input.
  
      Examples
      --------
  
      This will generate random variates from a Gaussian distribution and return
      the median absolute deviation for that distribution::
  
          >>> from astropy.stats import median_absolute_deviation
          >>> from numpy.random import randn
          >>> randvar = randn(10000)
          >>> mad = median_absolute_deviation(randvar)
  
      See Also
      --------
      numpy.median
  
      """
  
      a = np.ma.array(a, copy=False)
      a_median = bn_median(a, axis=axis)
  
      # re-broadcast the output median array to subtract it
      if axis is not None:
          a_median = np.ma.expand_dims(a_median, axis=axis)
  
      # calculated the median average deviation
      return bn_median(np.ma.abs(a - a_median), axis=axis) # median_absolute_deviation
  
  #*******************************************************
  def do_unflagging(station, variable, all_data, reporting_accuracies, neigh_count, dpd_flags, FLAG_COL_DICT, start, logfile, plots = False, diagnostics = False):
      '''
      Set up and run the unflagging process for the specified tests
  
      :param MetVar station: station object
      :param string variable: variable to process
      :param array all_data: array containing all neighbour obs for full time period
      :param array reporting accuracies: reporting accuracy for each neighbour
      :param array neigh_count: number of neighbours with data at each time stamp
      :param array dpd_flags: number of neighbours that have DPD set at each time stamp
      :param dict FLAG_COL_DICT: look up dictionary to 
      :param datetime start: start of dataset
      :param file logfile: logfile to store outputs
      :param bool plots: do plots
      '''
  
  
      # unflagging using neighbours
      '''This is slow - np.ma.median is known to be slow
      https://github.com/astropy/ccdproc/issues/74
      https://github.com/astropy/ccdproc/blob/122cdbd5713140174f057eaa8fdb6f9ce03312df/docs/ccdproc/bottleneck_example.rst'''
      mean_of_neighbours = bn_median(all_data, axis = 0)
      std_of_neighbours = median_absolute_deviation(all_data, axis = 0)
  
      # find where the spread of neighbour observations is less than 1/2
      #    of maximum reporting accuracy
      std_of_neighbours[std_of_neighbours < 0.5*max(reporting_accuracies)] = 0.5*max(reporting_accuracies)
  
      # create series of normalised differences of obs from neighbour mean
      st_var = getattr(station, variable)
      normalised_differences = np.ma.abs(st_var.data - mean_of_neighbours)/std_of_neighbours
              
      for qc_test in ["climatological","gap","odd","dpd"]:
          
          if qc_test == "dpd" and variable == "dewpoints":
              flags = station.qc_flags[:, UNFLAG_COL_DICT[qc_test][variable]] 
              unset_locs = unflagging_locs(normalised_differences, flags, neigh_count, dpd_count = dpd_flags)
          elif qc_test == "dpd":
              # only unflag DPD on dewpoints
              continue
  
          elif qc_test == "gap" and variable != "slp":
              # only unflag gap check on slp observations
              continue
  
          else:
              flags = station.qc_flags[:, UNFLAG_COL_DICT[qc_test][variable]] 
              if qc_test == "gap" or qc_test == "climatological":
                  # only tentative flags
                  unset_locs = unflagging_locs(normalised_differences, flags, neigh_count, flag_value = 2)
              else:
                  unset_locs = unflagging_locs(normalised_differences, flags, neigh_count)
  
          if len(unset_locs) > 0:
              station.qc_flags[unset_locs, UNFLAG_COL_DICT[qc_test][variable]] = 0
              
              # need to unflag attribute if and only if no other flags are set
              subset_flags = station.qc_flags[:, FLAG_COL_DICT[variable]]
              total_flags = np.sum(subset_flags[unset_locs, :], axis = 1)
              clean_locs = np.where(total_flags == 0)
              st_var.flags[unset_locs[clean_locs]] = 0
  
          # and print result
          if plots or diagnostics:
              utils.print_flagged_obs_number(logfile, "Unflagging "+qc_test, variable, len(unset_locs), noWrite = True)
          else:
              utils.print_flagged_obs_number(logfile, "Unflagging "+qc_test, variable, len(unset_locs))
              
  
          if plots:
              if len(unset_locs) > 0:
                  plot_outlier(station, variable, unset_locs, all_data, start)
  
      station = utils.append_history(station, "Unflagging - "+variable)  
      
      return # do_unflagging
  
  
  #************************************************************************
  if __name__ == "__main__":
      
      print "Utilities for neighbour check"
  