diff --git a/src/covariance.py b/src/covariance.py
index 5a2791f..8bc34a3 100644
--- a/src/covariance.py
+++ b/src/covariance.py
@@ -547,11 +547,11 @@ def _calcCovarianceMatrix(self, method='focal', nthreads=None):
                             self.config.adaptive_subsampling)\
                 .reshape(nleaves, nleaves)
 
-            if self.quadtree.leaf_mean_px_var is not None:
+            if self.quadtree.leaf_mat_mean_px_var is not None:
                 self._log.debug(
                     'Adding variance from scene.displacement_px_var')
                 cov_matrix[num.diag_indices_from(cov_matrix)] +=\
-                    self.quadtree.leaf_mean_px_var
+                    self.quadtree.leaf_mat_mean_px_var
 
         else:
             raise TypeError('Covariance calculation %s method not defined!'
diff --git a/src/quadtree.py b/src/quadtree.py
index e0a82f0..39458c7 100644
--- a/src/quadtree.py
+++ b/src/quadtree.py
@@ -84,13 +84,32 @@ def var(self):
 
     @property_cached
     def mean_px_var(self):
-        """ Variance of displacement
+        """ Mean pixel variance of displacement in leaf
         :type: float
         """
         if self.displacement_px_var is not None:
             return float(num.nanmean(self.displacement_px_var))
+            
+        return None
+    
+    @property_cached
+    def mat_mean_px_var(self):
+        """ Mean matrix pixel variance in leaf
+        
+        Mean value of Matrix with pixel variance in
+        :type: float
+        """
+        if self.displacement_px_var is not None:
+            nvalid_leaf = num.sum(num.isfinite(self.displacement_px_var))
+            nmat = nvalid_leaf * nvalid_leaf
+            if nmat>0.:
+                return float(num.nansum(self.displacement_px_var)/nmat)
+            else:
+                return float(0.)
+            
         return None
 
+
     @property_cached
     def corr_median(self):
         """ Standard deviation of node's displacement corrected for median
@@ -760,6 +779,17 @@ def leaf_mean_px_var(self):
         if self.scene.displacement_px_var is not None:
             return num.array([lf.mean_px_var for lf in self.leaves])
         return None
+    
+    @property
+    def leaf_mat_mean_px_var(self):
+        """
+        :getter: Mean pixel variance in each quadtree,
+            if :attr:`kite.Scene.displacement_px_var` is set.
+        :type: :class:`numpy.ndarray`, size ``N``.
+        """
+        if self.scene.displacement_px_var is not None:
+            return num.array([lf.mat_mean_px_var for lf in self.leaves])
+        return None
 
     @property_cached
     def leaf_means(self):
diff --git a/src/scene.py b/src/scene.py
index 8a90642..8e32bab 100644
--- a/src/scene.py
+++ b/src/scene.py
@@ -920,7 +920,7 @@ def save(self, filename=None):
         _file, ext = op.splitext(filename)
         filename = _file if ext in ['.yml', '.npz'] else filename
 
-        components = ['_displacement', 'theta', 'phi']
+        components = ['_displacement', 'theta', 'phi', 'displacement_px_var']
         self._log.debug('Saving scene data to %s.npz' % filename)
 
         num.savez('%s.npz' % (filename),
@@ -956,6 +956,9 @@ def load(cls, filename):
             displacement = data['arr_0']
             theta = data['arr_1']
             phi = data['arr_2']
+            if num.shape(data)[0]>3:
+                displacement_px_var = data['arr_3']
+
         except IOError:
             raise UserIOWarning('Could not load data from %s.npz' % basename)
 
@@ -965,11 +968,20 @@ def load(cls, filename):
         except IOError:
             raise UserIOWarning('Could not load %s.yml' % basename)
 
-        scene = cls(
-            displacement=displacement,
-            theta=theta,
-            phi=phi,
-            config=config)
+        if num.shape(data)[0]>3:
+            scene = cls(
+                displacement=displacement,
+                theta=theta,
+                phi=phi,
+                displacement_px_var = displacement_px_var,
+                config=config)
+        else:
+            scene = cls(
+                displacement=displacement,
+                theta=theta,
+                phi=phi,
+                config=config)
+                
         scene._log.debug('Loading from %s[.npz,.yml]', basename)
 
         scene.meta.filename = filename