Index: trunk/SDToolBox/data_processing.py
===================================================================
diff -u -r64 -r66
--- trunk/SDToolBox/data_processing.py	(.../data_processing.py)	(revision 64)
+++ trunk/SDToolBox/data_processing.py	(.../data_processing.py)	(revision 66)
@@ -13,6 +13,7 @@
 from SDToolBox.output_data import OutputData
 from sklearn.decomposition import PCA as skPCA
 from sklearn.preprocessing import StandardScaler
+import math
 
 import numpy as np
 import xarray as xr
@@ -143,58 +144,56 @@
         times_size = len(times)
 
         derivative_points = \
-            np.zeros((times_size, latitude_size, longitude_size), dtype=float)
-
+            np.nan*np.ones((times_size, latitude_size, longitude_size), dtype=float) # iniitialize with NaNs
         gradient_slp = []
 
         for time in range(1, len(times)):
-            for longitude_idx in range(2, len(longitude)-1): # x
-                for latitude_idx in range(2, len(latitude)-1): # y
+            for longitude_idx in range(2, longitude_size-1):
+                for latitude_idx in range(2, latitude_size-1):
 
                     phi = np.pi *  \
                         np.abs(meshed_latitudes[latitude_idx, longitude_idx])/180
                     derivative_longitude_first = \
-                        (slp[:, latitude_idx, longitude_idx] - slp[:, latitude_idx, longitude_idx - 1]) / np.cos(phi)
+                        (slp[time, latitude_idx, longitude_idx] - slp[time, latitude_idx, longitude_idx - 1]) / np.cos(phi)
                     derivative_longitude_second = \
-                        (slp[:, latitude_idx, longitude_idx + 1] - slp[:, latitude_idx, longitude_idx]) / np.cos(phi)
+                        (slp[time, latitude_idx, longitude_idx + 1] - slp[time, latitude_idx, longitude_idx]) / np.cos(phi)
                     derivative_latitude_first = \
-                        (slp[:, latitude_idx, longitude_idx] - slp[:, latitude_idx - 1, longitude_idx]) / np.cos(phi)
+                        (slp[time, latitude_idx, longitude_idx] - slp[time, latitude_idx - 1, longitude_idx]) / np.cos(phi)
                     derivative_latitude_second = \
-                        (slp[:, latitude_idx + 1, longitude_idx] - slp[:, latitude_idx, longitude_idx]) / np.cos(phi)
+                        (slp[time, latitude_idx + 1, longitude_idx] - slp[time, latitude_idx, longitude_idx]) / np.cos(phi)
 
-                    derivative_points[:, latitude_idx, longitude_idx] = np.square(derivative_longitude_first) + np.square(derivative_longitude_second)/2 + np.square(derivative_latitude_first) + np.square(derivative_latitude_second)/2
-        gradient_slp.append(derivative_points)
+                    derivative_points[time - 1, latitude_idx - 2, longitude_idx - 2] =  \
+                        (np.square(derivative_longitude_first) +
+                            np.square(derivative_longitude_second))/2 \
+                        + (np.square(derivative_latitude_first) +
+                            np.square(derivative_latitude_second))/2
+        slp.attrs['gradient'] = derivative_points
+
+        return slp
         # probably here we need to use .assign_coords function to add grad2slp to the existing array
 
-        data_array = xr.DataArray(
-            gradient_slp,
-            coords=[time, latitude, longitude],
-            dims=['times', 'longitude', 'latitude'])
-
-        return data_array
-
     def compute_PCA(self, slp: xr.DataArray):
 
-        grad2slp = slp.coords['grad2slp']
+        grad2slp = slp.attrs['grad2slp']
 
         if grad2slp is None:
             raise Exception(om.error_no_gradient_was_calculated)
 
         scaler = StandardScaler(with_mean=True, with_std=True)
         scaler.fit(grad2slp)
-        transformed_grad2slp = scaler.transform(grad2slp)
+        transformed_grad2slp = scaler.transform(grad2slp) #normalizing the features
         exvar = 99.88/100
         pca = skPCA(n_components=exvar)
         pca.fit(transformed_grad2slp)
         principal_components = pca.components_
         temporal_indices = \
             np.dot(transformed_grad2slp, principal_components.transpose())
         standard_temporal_indices = temporal_indices.std(axis=0)
-        pca.mean = scaler.mean_()
-        pca.std = np.std(scaler.var_())
+        pca.mean = scaler.mean_ # probably the mean is on the transformed 
+        pca.std = np.std(scaler.var_)
         pca.eof = principal_components
         pca.pc = temporal_indices
         pca.variance = \
             standard_temporal_indices**2/np.sum(standard_temporal_indices**2)*100
-        pca.variance = \
-            pca.explained_variance_ratio_()
+
+        return pca.mean, pca.std, pca.variance