Modification algorithme corTsgLinear.m :

1 - Ajustement lineaire permis pour 3 points et plus
2 - si seulement 2 points calcul de l'ecart moyen

tsg_util/corTsgLinear.m

@@ -16,6 +16,11 @@ function [error] = corTsgLinear(hMainFig, PARA, dateMin, dateMax)
 % less than 4 samples
 %
 % The error minimum cannot be lower than 0.01
+%
+% 2009/03/23 Y. Gouriou
+% 1) Compute linear fit witha minimum of 3 points.
+% 2) If only 2 points: no linear fit, but compute the average difference
+%    with the 2 points.
 
 % Get application data
 % --------------------
@@ -68,11 +73,11 @@ if isempty( tsg.([PARA '_ADJUSTED']) )
 end
 
 if dateMax > dateMin
- 
+
   % Find samples within TIME_WINDOWS with Good, probably Good, QC
   % -------------------------------------------------------------
   ind = find( tsg.DAYD_EXT    >= dateMin &  tsg.DAYD_EXT    <= dateMax &...
-              tsg.([SAMPLE '_EXT_QC']) <= PROBABLY_GOOD);
+    tsg.([SAMPLE '_EXT_QC']) <= PROBABLY_GOOD);
 
   if ~isempty(ind)
 
@@ -82,51 +87,74 @@ if dateMax > dateMin
 
     % Compute linear fit of the TSG/SAMPLE difference
     % -----------------------------------------------
-    if ~isempty(ind2) && length(ind2) > 1
-      if ~isempty(tsg.EXT_DIF(ind(ind2)))
-
-        % Linear fit applied to the difference tsg-sample
-        % -----------------------------------------------
-        X = tsg.DAYD_EXT(ind(ind2));
-        Y = tsg.EXT_DIF(ind(ind2));
-
-        [p, S, mu] = polyfit( X, Y, 1);
-
-        % The correction is applied to the TSG between dateMin and dateMax using
-        % a linear interpolation only on measurements with keptCode Quality
-        % Codes
-        % ----------------------------------------------------------------------
-        for icode = 1 : length( keptCode )
-          dtTsg = find( tsg.DAYD    >= dateMin  & tsg.DAYD <= dateMax &...
-            tsg.([PARA '_QC']) == keptCode( icode ));
-
-          if ~isempty( dtTsg )
-            [tsg.([PARA '_ADJUSTED'])(dtTsg),...
-                                    tsg.([PARA '_ADJUSTED_ERROR'])(dtTsg)] =...
-                                    polyval( p, tsg.DAYD(dtTsg), S, mu);
-            tsg.([PARA '_ADJUSTED'])(dtTsg) = ...
+    if ~isempty(ind2) && length(ind2) > 2
+      %       if ~isempty(tsg.EXT_DIF(ind(ind2)))    Is this line useful?
+
+      % Linear fit applied to the difference tsg-sample
+      % -----------------------------------------------
+      X = tsg.DAYD_EXT(ind(ind2));
+      Y = tsg.EXT_DIF(ind(ind2));
+
+      [p, S, mu] = polyfit( X, Y, 1);
+
+      % The correction is applied to the TSG between dateMin and dateMax using
+      % a linear interpolation only on measurements with keptCode Quality Codes
+      % ------------------------------------------------------------------------
+      for icode = 1 : length( keptCode )
+        dtTsg = find( tsg.DAYD    >= dateMin  & tsg.DAYD <= dateMax &...
+          tsg.([PARA '_QC']) == keptCode( icode ));
+
+        if ~isempty( dtTsg )
+
+          % Compute the correction + the error
+          % ----------------------------------
+          [tsg.([PARA '_ADJUSTED'])(dtTsg),...
+            tsg.([PARA '_ADJUSTED_ERROR'])(dtTsg)] =...
+                                            polyval( p, tsg.DAYD(dtTsg), S, mu);
+
+          % Compute the corrected value : orignal value + correction
+          % --------------------------------------------------------
+          tsg.([PARA '_ADJUSTED'])(dtTsg) =...
                             tsg.(PARA)(dtTsg) + tsg.([PARA '_ADJUSTED'])(dtTsg);
-            tsg.([PARA '_ADJUSTED_QC'])(dtTsg) = tsg.([PARA '_QC'])(dtTsg);
 
-            % The error is maximum, equal to 1, if the median is computed with
-            % less than 4 samples
-            % -------------------------------------------------------------
-            if length( ind2 ) < 4
-              tsg.([PARA '_ADJUSTED_ERROR'])(dtTsg) = 1;
-            end
-          end
+          % Transfer the QC
+          % ---------------
+          tsg.([PARA '_ADJUSTED_QC'])(dtTsg) = tsg.([PARA '_QC'])(dtTsg);
+        end
+      end
+
+      % If there is only 2 points, compute the mean difference. Apply an error max.
+      % --------------------------------------------------------------------------
+    elseif ~isempty(ind2) && length(ind2) == 2
+
+      meanDif = mean( tsg.EXT_DIF(ind(ind2)) );
 
+      % The correction is applied to the TSG between dateMin and dateMax using
+      % a linear interpolation only on measurements with keptCode Quality Codes
+      % ------------------------------------------------------------------------
+      for icode = 1 : length( keptCode )
+        dtTsg = find( tsg.DAYD    >= dateMin  & tsg.DAYD <= dateMax &...
+                      tsg.([PARA '_QC']) == keptCode( icode ));
+
+        if ~isempty( dtTsg )
+          tsg.([PARA '_ADJUSTED'])(dtTsg)       = tsg.(PARA)(dtTsg) + meanDif;
+          tsg.([PARA '_ADJUSTED_ERROR'])(dtTsg) = meanDif / 2;
+          tsg.([PARA '_ADJUSTED_QC'])(dtTsg)    = tsg.([PARA '_QC'])(dtTsg);
         end
+      end
 
-        % The error minimum cannot be lower than 0.01 or equal to Inf
-        % -----------------------------------------------------------
-        tsg.([PARA '_ADJUSTED_ERROR'])...
-                            ( tsg.([PARA '_ADJUSTED_ERROR']) < 0.01 ) = 0.01;
+    end
 
-        tsg.([PARA '_ADJUSTED_ERROR'])...
+    % The error minimum cannot be lower than 0.01 or equal to Inf
+    % -----------------------------------------------------------
+    if ~isempty(ind2) && length(ind2) >= 2
+      tsg.([PARA '_ADJUSTED_ERROR'])...
+                               ( tsg.([PARA '_ADJUSTED_ERROR']) < 0.01 ) = 0.01;
+
+      tsg.([PARA '_ADJUSTED_ERROR'])...
                             (isinf(tsg.([PARA '_ADJUSTED_ERROR'])) == 1) = 0.01;
-      end
     end
+    %       end
   end
 
   % Update tsg application data

tsgqc.m

@@ -1927,8 +1927,8 @@ hrbInterpCancel = uicontrol( ...
           % --------------------------
           if ~isempty(tsg.DAYD_EXT)
 
-            ind = find(tsg.DAYD_EXT > p1(1,1) &  tsg.DAYD_EXT < p2(1,1) & ...
-              tsg.EXT_DIF  > p1(1,2) &  tsg.EXT_DIF  < p2(1,2));
+            ind = find(tsg.DAYD_EXT > p1(1) &  tsg.DAYD_EXT < p2(1) & ...
+              tsg.EXT_DIF  > p1(2) &  tsg.EXT_DIF  < p2(2));
 
             % Keep the information on the indices of the selected zone
             % --------------------------------------------------------