diff --git a/deerlab/deerload.py b/deerlab/deerload.py
index b3dc432e4..2e7bd5a67 100644
--- a/deerlab/deerload.py
+++ b/deerlab/deerload.py
@@ -63,7 +63,6 @@ def deerload(fullbasename, plot=False, full_output=False, *args,**kwargs):
     parameters = read_description_file(filename_dsc)
     parDESC = parameters["DESC"]
     parSPL = parameters["SPL"]
-    
     # XPTS, YPTS, ZPTS specify the number of data points along x, y and z.
     if 'XPTS' in parDESC:
         nx = int(parDESC['XPTS'])
@@ -167,12 +166,9 @@ def deerload(fullbasename, plot=False, full_output=False, *args,**kwargs):
             abscissa[:Dimensions[index],index] = np.linspace(minimum,minimum+width,npts)
         if axistype == 'NTUP':
             raise ValueError('Cannot read data with NTUP axes.')
-    
-    # In case of column filled with NaN, erase the column in the array
-    abscissa = abscissa[:,~np.isnan(abscissa).all(axis=0)]
+
     dt_data = dt_spc
     dt_spc = dt_spc.newbyteorder(byteorder)
-
     # Read data matrix and separate complex case from real case.
     data = np.full((nx,ny,nz),np.nan)
     # reorganize the data in a "complex" way as the real part and the imaginary part are separated
@@ -182,11 +178,10 @@ def deerload(fullbasename, plot=False, full_output=False, *args,**kwargs):
         if parDESC['IKKF'] == 'REAL':
             data = np.full((nx,ny,nz),np.nan) 
             with open(filename_dta,'rb') as fp:
-                 data = np.frombuffer(fp.read(),dtype=dt_spc).reshape(nx,ny,nz)
+                 data = np.frombuffer(fp.read(),dtype=dt_spc)
             data = np.copy(data)
         elif parDESC['IKKF'] == 'CPLX':
             dt_new = np.dtype('complex')
-            data = np.full((2*nx*ny*nz,1),np.nan)
             with open(filename_dta,'rb') as fp:
                 data = np.frombuffer(fp.read(),dtype=dt_spc)
                 # Check if there is multiple harmonics (High field ESR quadrature detection)
@@ -199,25 +194,40 @@ def deerload(fullbasename, plot=False, full_output=False, *args,**kwargs):
                 n_harmonics = sum(harmonics)[0]
                 if n_harmonics != 0:
                     ny = int(len(data)/nx/n_harmonics)
-            
+
             # copy the data to a writable numpy array
-            data = np.copy(data.astype(dtype=dt_data).view(dtype=dt_new).reshape(nx,ny,nz))
+            data = np.copy(data.astype(dtype=dt_data).view(dtype=dt_new))
         else:
             raise ValueError("Unknown value for keyword IKKF in .DSC file!")
     else:
         warn("Keyword IKKF not found in .DSC file! Assuming IKKF=REAL.")
     
-    if nz == 1:
-        data = data.reshape(nx,ny)
-            
-     # ns -> us converesion
-    abscissa /= 1e3
+    # Split 1D-array into 3D-array according to XPTS/YPTS/ZPTS 
+    data = np.array_split(data,nz)
+    data = np.array(data).T
+    data = np.array_split(data,ny)
+    data = np.array(data).T
 
     # Ensue proper numpy formatting
-    abscissa,data = np.atleast_1d(abscissa,data)
-    abscissa = np.squeeze(abscissa)
+    data = np.atleast_1d(data)
     data = np.squeeze(data)
 
+    # Abscissa formatting
+    abscissa = np.atleast_1d(abscissa)
+    abscissa = np.squeeze(abscissa)
+    abscissas = []
+    # Convert to list of abscissas 
+    for absc in abscissa.T: 
+        # Do not include abcissas full of NaNs
+        if not all(np.isnan(absc)):
+            # ns -> us converesion
+            absc /= 1e3
+            # Remove nan values to ensure proper length of abscissa
+            abscissas.append(absc[~np.isnan(absc)])
+    # If 1D-dataset, return array instead of single-element list
+    if len(abscissas)==1:
+        abscissas = abscissas[0]
+
     if plot:
         plt.plot(abscissa,np.real(data),abscissa,np.imag(data))
         plt.xlabel("Time (μs)")
@@ -226,9 +236,9 @@ def deerload(fullbasename, plot=False, full_output=False, *args,**kwargs):
         plt.show()
     
     if full_output:
-        return abscissa, data, parameters
+        return abscissas, data, parameters
     else:
-        return abscissa, data
+        return abscissas, data
 
 
 def read_description_file(DSCFileName):