diff --git a/observations.py b/observations.py
index ea45957db8a9351d07ad5a60f5ea8a4e9cd6cd03..f56413d53f1c3886539b1518300817c8b08c8f5b 100644
--- a/observations.py
+++ b/observations.py
@@ -14,7 +14,7 @@ def decode_observations(observation_files):
     # https://stackoverflow.com/questions/20906474/import-multiple-csv-files-into-pandas-and-concatenate-into-one-dataframe    
 
     # Create a list of dataframes, one per file
-    all_files = glob.glob(observation_files)
+    all_files = sorted(glob.glob(observation_files))
     li = []
     for filename in all_files:
         df = pd.read_csv(filename, index_col=None, header=0)
@@ -39,6 +39,64 @@ def decode_observations(observation_files):
 
     return observations,ocoords
 
+def random_sort_observations(observations,ocoords):
+
+    # Check that input arrays have consistent dimensions
+    assert (observations.shape == ocoords.shape),\
+        f'observation and coordinate arrays must have the same size'
+
+    # Preallocate arrays
+    nassim,nobs = observations.shape
+    shuffled_obs = np.zeros(observations.shape)+np.nan
+    shuffled_coords = np.zeros(observations.shape)+np.nan
+
+    # Randomize observation order, differently at each time
+    indices = np.arange(nobs)
+    for i in range(nassim):
+        np.random.shuffle(indices)
+        shuffled_obs[i,:] = observations[i,indices]
+        shuffled_coords[i,:] = ocoords[i,indices]
+        
+    return shuffled_obs,shuffled_coords
+
 if __name__ == '__main__':
 
+    printout = False
+
+    # Read the observations
     observations,ocoords = decode_observations('./LES_data/Theta/*csv')
+
+    if printout:
+        ntim = 0
+        for i in range(observations.shape[1]):
+            print(ocoords[ntim,i],observations[ntim,i])
+
+    # Sort the observations randomly
+    observations,ocoords = random_sort_observations(observations,ocoords)
+    if printout:
+        for i in range(observations.shape[1]):
+            print(ocoords[ntim,i],observations[ntim,i])
+
+    makeplots = False
+    # Make plots
+    if makeplots:
+        fig, ax1 = p.subplots(1,1,constrained_layout=True)
+        fig.set_size_inches(4,3)
+
+        tint = 300
+        t = np.linspace(0,tint*(ocoords.shape[0]-1),ocoords.shape[0])[:,None]+np.zeros(ocoords.shape)
+        t = t/3600.
+        c = ax1.pcolormesh(t,ocoords,observations,vmin=290,vmax=296)
+        ax1.contour(t,ocoords,observations,
+                np.linspace(290,290+0.003*2000,13),
+                colors='white',
+                linestyles='--',
+                linewidths=0.75)
+        ax1.set_ylim([0,2000])
+        ax1.set_ylabel(r'Height (m)')
+        ax1.set_xlabel(r'Time (s)')
+        ax1.set_title(r'$\overline{\theta}$ (K) obs in t=(%u-%u) h'%(t.min(),t.max()))
+        ax1.legend(frameon=False)
+        p.colorbar(c,orientation='horizontal')
+        fig.savefig('LES_obs.png',format='png',dpi=300)
+        p.close(fig)