CDIP banner
CDIP banner recent historic documents
sub-menu
Documentation
 
FAQs & Summaries
Glossary
Publications
 
Introduction
  History and Funding
  Program Goals
Wave Measurement
  Wave Generation
  Wave Dynamics
  Irregular Waves
  Spectral Analysis
  Gauging Waves
  Hurricane Events
  Tsunami Events
Instrumentation
  Underwater Sensors
  Surface Buoys
  Meteorological
Data Acquisition
  System Organization
  Hardware
  Software
Data Processing
  System Organization
  Software
  Quality Control
Data Management
  Stations and Sets
  Files and Storage
CDIP Products
  Data Formats
  Web Products
  COOS Integration
  WIS
  QARTOD
  Wave Eval Tool
  Metadata
  Custom Products
  NDBC XML/NWS Format
  NDBC Dial-A-Buoy
  Access Instructions
 
Related Links
 
c-- EDITOR ---------------------------------------------------------------------
c
c Contains our time series quality control checks:
c
c   value_test, delta_test, flat_episodes_test, mean_shift_test,
c   equal_peaks_test, acceleration_test, mean_crossing_test,
c   period_distribution_test and spike_edit.
c
c The chi-square test is not included here yet.
c 
c Used by: .wv/displacement/source/disp_editor.f
c          .wv/temperature/source/temp_editor.f
c          .wv/barometer/source/barm_editor.f
c
c-------------------------------------------------------------------------------

      module editor

      use elemental_stats

      contains
 

c-- VALUE_TEST -----------------------------------------------------------------
c
c   Tests whether time series values are within min_val and max_val. 
c   If a time series value is less than min_val, returns 1.
c   If a time series value is greater than max_val, returns 2.
c   If there are time series values exceeding both limits, returns 3.
c   Also, min_val and max_val are set to the minimum and maximum values
c   in the time series respectively for error diagnostics.
c
c-------------------------------------------------------------------------------
      integer function value_test(ts_vals, ts_start, ts_end, ts_size, 
     *                                min_val, max_val)

        integer ts_start, ts_end, ts_size
        real ts_vals(ts_size)
        real max_val, min_val, min_limit, max_limit

        value_test = 0

        min_limit = min_val
        max_limit = max_val

        min_val = MINVAL(ts_vals(ts_start:ts_end))
        max_val = MAXVAL(ts_vals(ts_start:ts_end))

        if ( min_val .lt. min_limit .and. max_val .gt. max_limit ) then
          value_test = 3
        else if ( max_val .gt. max_limit ) then
          value_test = 2
        else if ( min_val .lt. min_limit ) then
          value_test = 1
        end if

      end function


c-- DELTA_TEST -----------------------------------------------------------------
c
c   Checks the difference between a point with index i, and its successor. 
c   If greater than max difference, returns i, else returns 0.
c   We assume there are no missing values within the ts window.
c
c   Used by: .wv/temperature/
c
c-------------------------------------------------------------------------------
      integer function delta_test(ts_vals, ts_start, ts_end, ts_size, 
     *                                max_difference)

        integer ts_start, ts_end, ts_size
        real ts_vals(ts_size)
        real max_difference

        delta_test = 0   
        do i=ts_start, ts_end-1
            if ( abs(ts_vals(i+1) - ts_vals(i)) .gt. max_difference ) then
              delta_test = i
              return
            end if
        end do

      end function


c-- FLAT_EPISODES_TEST ---------------------------------------------------------
c
c   Tests for a slowly changing value over a specified number of points. If this
c   happens more than 5 times, we reject the data by returning 1, otherwise
c   we return 0.
c
c-------------------------------------------------------------------------------
      integer function flat_episodes_test(ts_vals, ts_start, ts_end, ts_size,
     *                                    paros, hs)

        integer ts_start, ts_end, ts_size
        real ts_vals(ts_size)
        real max_diff, hs
        integer num_points  !* from ref_pnt check num_points preceding.
        integer ref_pnt, flat_cnt, episode_cnt, predecessor
        logical paros

        flat_episodes_test = 0

        if( paros ) then       
          num_points = 10          
          max_diff = 1.0    !* cm reduced from 1.5cm 10/23/2001
        else
          num_points = 4
          max_diff = 0.5    !* cm
        end if

        episode_cnt = 0
        do ref_pnt = num_points+ts_start, ts_end
           flat_cnt = 0
           do i = 1, num_points       
              predecessor = ref_pnt - i
              if(abs(ts_vals(ref_pnt)-ts_vals(predecessor)) .lt. max_diff) then
                flat_cnt = flat_cnt + 1
              end if
           end do
           if (flat_cnt .ge. num_points) episode_cnt = episode_cnt + 1
           if(episode_cnt .gt. 5) then
             flat_episodes_test = 1
             return
           end if 
        end do

      end function


c-- MEAN_SHIFT_TEST ------------------------------------------------------------
c
c   Breaks the time series up into sections of 2048 points, and sub-sections
c   "intervals" of 256 points. The series interval means are compared to 
c   neigboring intervals. If the difference of any two adjacent means
c   exceeds 10% of the significant wave height of the 2048 section, 
c   the data is rejected, mean_shift_test returns 1. 
c
c   For non-standard series lengths such as 1024 points, the section
c   length is reduced to 1024, but the interval length stays the same.
c   If the time series is greater than 2048 points, the above is done for each
c   section separately.
c
c   Modifications:
c     o 14May2007 - set sect_limit to max(0.1*hs,7) rather than max(0.1*hs,15)
c       In the first week of april 2007 at 073 the mean was wandering, but
c       the hs was low, so the default of 15 was being used which did not
c       catch the problem. 
c
c-------------------------------------------------------------------------------
      integer function mean_shift_test(ts_vals, ts_start, ts_end, ts_size, 
     *                                   surge_filter, energy_basin, shift_diff)

        integer ts_start, ts_end, ts_size
        real ts_vals(ts_size)
        logical surge_filter, energy_basin
        integer :: num_intervals=8, prv_interval, series_length
        integer sctn_length, sctn_last, sctn_start, int_start, int_end, interval
        integer interval_span
        real sctn_limit, int_mean(8), shift_diff

        mean_shift_test = 0

        series_length = ts_end - ts_start + 1  !* number of points in the series
        sctn_length = 2048  !* initialize to get into while loop
        interval_span = 256  !* number of points in each interval

c       sctn_last = 0
        sctn_last = ts_start - 1

        do while ( sctn_length .eq. 2048 )

          sctn_length = min(ts_end-sctn_last,2048)

          if ( sctn_length .gt. interval_span ) then

            sctn_start = sctn_last + 1
            sctn_last = sctn_start + sctn_length - 1

            if ( surge_filter .or. energy_basin ) then
              sctn_limit = 60
            else 
              sctn_limit = max(
     *          0.10*4.0*stdev_lh(ts_vals, sctn_start, sctn_last, ts_size),
     *          7.0)
            end if

            num_intervals = FLOOR(sctn_length/real(interval_span))
            int_start = sctn_start
            do interval = 1, num_intervals
              int_end = int_start + interval_span - 1
              int_mean(interval) = mean_lh(ts_vals,int_start,int_end,ts_size)
              int_start = int_end + 1
            end do

            do interval = 2, num_intervals
               prv_interval = interval - 1
               shift_diff = abs( int_mean(interval) - int_mean(prv_interval) )
               if( shift_diff .gt. sctn_limit ) then
                  mean_shift_test = 1
                  return
               endif
            end do

          end if

        end do

      end function



c-- EQUAL_PEAKS_TEST -----------------------------------------------------------
c
c   Identifies the peaks and troughs of the time series ("peaks")
c   and compares those peak values. If the value of the current peak is
c   the same as one or more of the next 10 peaks, and this is
c   repeated for two other peaks, the data is rejected, and this
c   routine returns 1, otherwise it returns 0.
c
c-------------------------------------------------------------------------------
      integer function equal_peaks_test(ts_vals, ts_start, ts_end, ts_size, hs)

        integer ts_start, ts_end, ts_size
        real ts_vals(ts_size), hs
        integer prv_interval, num_episodes, peak
        integer num_peaks, num_equal, pnt, prev, next, test_interval
        integer int_start, int_end
        real peak_vals(8000)

        equal_peaks_test = 0   !* initialize as no error

        if ( hs .lt. 30.0 ) return !* 30cm

c -- identify peaks (see note above).

        num_peaks = 0
        do pnt = ts_start+1, ts_end-1
           prev = pnt - 1
           next = pnt + 1
           if( (ts_vals(pnt) .gt. ts_vals(prev) .and.
     *          ts_vals(pnt) .gt. ts_vals(next)) .or.
     *         (ts_vals(pnt) .lt. ts_vals(prev) .and.
     *          ts_vals(pnt) .lt. ts_vals(next)) ) then
               num_peaks = num_peaks + 1
               peak_vals(num_peaks) = ts_vals(pnt)
           end if
        end do

        test_interval = 10     !* number of points 
        num_episodes = 0       !* peak episode counter

        do peak = 1, num_peaks

           int_start = peak + 1
           int_end = min(peak+test_interval, num_peaks)
           num_equal = 0

           do i = int_start, int_end
              if(peak_vals(peak) .eq. peak_vals(i)) num_equal = num_equal + 1
           end do

           if( num_equal .ge. 2 ) num_episodes = num_episodes + 1
           if( num_episodes .gt. 2) then
             equal_peaks_test = 1   !* too many equal peaks
             return
           end if

        end do

      end function


c-- ACCELERATION_TEST ----------------------------------------------------------
c
c   The acceleration of a particle floating on the surface of the ocean
c   should never exceed g=980cm/s^2 in the downward direction, but due
c   to the symmetry of the swell it should also not exceed that limit
c   in the upward direction (correct logic?). In this test, we approximate
c   the second derivative for each point of the time series and test it
c   against this natural limit. If the acc is greater than g more than
c   three times, we reject and return 1, otherwise return 0.
c
c   2002/11/20 : Limit changed from g to (1/3)g
c-------------------------------------------------------------------------------
      integer function acceleration_test(ts_vals, ts_start, ts_end, ts_size, 
     *                                    sample_rate)

        integer ts_start, ts_end, ts_size
        real ts_vals(ts_size), sample_rate
        real prv_dv1, dv1, dv2       !* first and second derivatives
        integer times_exceeded, next1, pnt

        acceleration_test = 0
        do pnt= ts_start, ts_end-2
          next1 = pnt + 1
          dv1 = (ts_vals(next1) - ts_vals(pnt))*sample_rate
          if ( pnt .gt. ts_start ) then
            dv2 = (dv1 - prv_dv1)*sample_rate
            if ( abs(dv2) .gt. (980.0/3.0) ) then
              times_exceeded = times_exceeded + 1
              if ( times_exceeded .gt. 3 ) then
                acceleration_test = 1
                return
              end if
            end if
          end if
          prv_dv1 = dv1
        end do

      end function


c-- MEAN_CROSSING_TEST ---------------------------------------------------------
c
c   For a displacement time series, we expect many mean crossings. 
c   We reject if there are no mean crossings within max_points
c   points. We divide the series into sections of 1024 points, and for each
c   section calculate a mean for the test which is used to determine
c   whether or not we cross the mean.
c
c-------------------------------------------------------------------------------
      integer function mean_crossing_test(ts_vals, ts_start, ts_end, ts_size, 
     *                                   hs, max_points)

        integer ts_start, ts_end, ts_size
        real ts_vals(ts_size)
        real hs              !* 4*(standard deviation)
        integer sctn_length, prv_interval, counter, next1, pnt
        integer sctn_last, sctn_start, int_start, int_end, interval
        integer max_points, series_length
        real sctn_limit, sctn_mean, diff

        mean_crossing_test = 0

c -- return if hs smaller than tidal range

        if ( hs .lt. 10.0 ) return  !* 10cm

c -- mean crossing test


        series_length = ts_end - ts_start + 1  !* number of points in the series
        max_sctn_length = 1024  !* max num points of a section
        sctn_length = min(series_length, max_sctn_length)   !* num of pnts in 1st section
        max_points = 0.15 * sctn_length !* test limit: 15% of the section.
        sctn_last = ts_start - 1

        do while ( sctn_length .gt. 0 )

            sctn_start = sctn_last + 1
            sctn_last = sctn_start + sctn_length - 1

            sctn_mean = mean_lh(ts_vals, sctn_start, sctn_last, ts_size)

c -- if mean crossing, reset counter

            counter = 0
            do pnt=sctn_start, sctn_last-1
              counter = counter + 1
              next1 = pnt + 1
              if( (ts_vals(pnt)-sctn_mean)*(ts_vals(next1)-sctn_mean) .lt. 0.0) then
                counter = 0
              end if
              if ( counter .ge. max_points ) then
                mean_crossing_test = 1
                return
              end if
            end do

            sctn_length = min(series_length-sctn_last, max_sctn_length)

        end do

      end function


c-- PERIOD_DISTRIBUTION_TEST ---------------------------------------------------
c
c   The distribution of the period (length of time between two successive 
c   upward mean crossings) is approximated by counting the number of periods
c   that fall within specified intervals (bins). If any bin with lower interval
c   greater than 22 seconds contains more than max_limit fraction of periods, 
c   we reject the data and return 1, else we return 0. 
c
c-------------------------------------------------------------------------------
      integer function period_distribution_test(
     *     ts_vals, ts_start, ts_end, ts_size, hs, sample_rate)

        integer ts_start, ts_end, ts_size
        real ts_vals(ts_size), sample_rate, fraction_in_bin
        real hs                     !* 4*(standard deviation)
        real :: max_limit=.2        !* more than 20% in one band
        real band(18), series_mean, adj_value, prv_adj_value, period
        integer cut, pnt, num_up_crossings, pnts_in_period
        integer bin(17)

        data band /  0.0,  4.0,  6.0,  8.0, 10.0, 12.0, 14.0,
     *              16.0, 18.0, 20.0, 22.0, 24.0, 26.0, 28.0,
     *              30.0, 32.0, 34.0, 1000.0/

        period_distribution_test = 0

c -- return if hs too small

        if ( hs .lt. 30.0 ) return  !* 30 cm

        series_mean = mean_lh(ts_vals,ts_start,ts_end,ts_size)

c -- bin period distribution

        num_up_crossings = 0
        pnts_in_period = 0
        do cut=1,17
          bin(cut) = 0
        end do
        prv_adj_value = ts_vals(ts_start) - series_mean

        do pnt = ts_start+1, ts_end

          pnts_in_period = pnts_in_period + 1
          adj_value = ts_vals(pnt) - series_mean

          if(prv_adj_value .lt. 0.0 .and. (prv_adj_value*adj_value) .le. 0.0) then
            num_up_crossings = num_up_crossings + 1
            period = real(pnts_in_period)/sample_rate
            do cut=1,17
              if ( period .ge. band(cut) .and. period .lt. band(cut+1) ) then
                bin(cut) = bin(cut) + 1
              end if
            end do
            pnts_in_period = 0
          end if

          prv_adj_value = adj_value

        end do

        if ( num_up_crossings .eq. 0 ) then
          period_distribution_test = 1
          return
        end if

c -- check if the fraction exceeds max_limit

        do cut=11,17   !* check only periods greater than 22 seconds.
          
          fraction_in_bin = bin(cut)/real(num_up_crossings)
          if ( fraction_in_bin .gt. max_limit ) then
            period_distribution_test = 1
            return
          end if

        end do

      end function


c -- SPIKE_EDIT ----------------------------------------------------------------
c
c Checks for spikes in a time series window. If finds a spike, sets
c ts value to the avg of the current and previous point. Returns 0 for
c success, 1 for too many spikes and 2 for too many loops. This check is
c identical to the original spike check in the .dsk2af/editors.
c
c-------------------------------------------------------------------------------
      integer function spike_edit(ts_vals, ts_start, ts_end, ts_size, 
     *                 std_mult, min_limit, percent_spikes_allowed, num_loops)

        integer loops, ts_start, ts_end, spike_cnt, ts_size
        real std_mult, min_limit, spike_limit, percent_spikes_allowed
        real spread, sigma, ts_mean, ts_vals(ts_size)

        spike_limit = percent_spikes_allowed*(ts_end-ts_start+1)/100.0

        do loops=1,num_loops

          sigma = stdev_lh(ts_vals,ts_start,ts_end,ts_size)
          spread = max(std_mult*sigma,min_limit)
          ts_mean = mean_lh(ts_vals, ts_start, ts_end, ts_size)
          spike_cnt = 0
          do i = ts_start, ts_end
            if( abs(ts_vals(i) - ts_mean) .gt. spread ) then
              if ( i .eq. ts_start ) then
                ts_vals(i) = ts_vals(i+1)
              else if ( i .eq. ts_end ) then
                ts_vals(i) = ts_vals(i-1)
              else
                ts_vals(i) = (ts_vals(i-1) + ts_vals(i+1))/2.0
              end if
              spike_cnt = spike_cnt + 1
              if ( spike_cnt .gt. spike_limit ) then
                spike_edit = 1  !* number of spikes exceed limit
                return
              end if
            end if
          end do

          if ( spike_cnt .eq. 0 ) then
            spike_edit = 0
            return
          end if

          if (loops .eq. num_loops ) then
            spike_edit = 2  !* still have spikes after num_loops
            return
          end if

        end do

      end function


c -- SURGE_SPIKE_EDIT ----------------------------------------------------------
c
c Checks for spikes in surge series. These series are long and
c have spikes due to phone transmission that need to be removed prior
c to any other editing.
c
c-------------------------------------------------------------------------------
      integer function surge_spike_edit(ts_vals, ts_start, ts_end, ts_size, spike_cnt)

        integer ts_start, ts_end, ts_size, spike_cnt, pnt
        real max_spikes, ts_vals(ts_size)

         surge_spike_edit = 0
         spike_cnt = 0       
         max_spikes = .01*ts_size !* 1% spikes allowed

         do pnt = ts_start+1, ts_end
            if(abs(ts_vals(pnt) - ts_vals(pnt-1)) .gt. 40.) then  !* 40 cm.
               spike_cnt = spike_cnt + 1
               if ( spike_cnt .gt. max_spikes) then
                 surge_spike_edit = 1
                 return
               end if
               ts_vals(pnt) = ts_vals(pnt-1)
            end if
         end do

      end function


c -- BUOY_MEAN_TEST ------------------------------------------------------------
c
c Tests whether the mean is within specs for a non-directional buoy.
c Returns 1 if mean is outside specs, 0 otherwise.
c
c NOTE: TEST DISABLED! This test only works with old-format df files,
c or rd files - i.e. time series that have not had the mean removed.
c It can not be applied to our current nondir buoy df files.
c-------------------------------------------------------------------------------
      integer function buoy_mean_test(ts_vals, ts_start, ts_end, ts_size,
     *            cal_factor, cal_offset, series_avg)

        integer ts_start, ts_end, ts_size
        real ts_vals(ts_size), cal_factor, cal_offset

        expected_avg = abs(-16500. * cal_factor + cal_offset)
c       expected_avg = abs(-16500. * cal_factor + cal_offset + 1005.)
        delta = 600. * cal_factor   !* acceptable deviation

        series_avg = mean_lh(ts_vals, ts_start, ts_end, ts_size)

        if ( abs(series_avg - expected_avg) .gt. delta ) then
          buoy_mean_test = 1	!* TEST DISABLED
        else
          buoy_mean_test = 0
        end if
        
      end function


      end
Official UCSD Web Page