diff --git a/__pycache__/heatmap_visualisation.cpython-39.pyc b/__pycache__/heatmap_visualisation.cpython-39.pyc
index 8161ab61c997fbecdf290aa568589925a223abfd..dee5032c377d62ecfafacddfce7a8db612154c25 100644
Binary files a/__pycache__/heatmap_visualisation.cpython-39.pyc and b/__pycache__/heatmap_visualisation.cpython-39.pyc differ
diff --git a/behaviours.py b/behaviours.py
new file mode 100644
index 0000000000000000000000000000000000000000..5cf8c1b867b8364e4a79a2199cf280bae7561305
--- /dev/null
+++ b/behaviours.py
@@ -0,0 +1,124 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Feb  8 16:47:19 2023
+
+@author: adria036
+"""
+
+
+import os
+os.chdir(r"C:\Users\adria036\OneDrive - Wageningen University & Research\iAdriaens_doc\Projects\cKamphuis\nlas\scripts\uwb") 
+
+
+#%% import modules
+
+
+from datetime import date, datetime
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib
+import seaborn as sns
+#%matplotlib qt
+
+
+#%% import data, seth paths and constants
+
+
+# path_out 
+svpath = r"C:\Users\adria036\OneDrive - Wageningen University & Research\iAdriaens_doc\Projects\cKamphuis\nlas\results\exploration"
+
+# path to per barn areas
+path = os.path.join("W:","\ASG","WLR_Dataopslag","DairyCampus","3406_Nlas","uwb_processed")
+path_zones = os.path.join("W:","\ASG","WLR_Dataopslag","DairyCampus","3406_Nlas","raw","dc_sewio")
+fna = "\\Barn_areas.xlsx"
+area_zones = pd.read_excel(path_zones+fna, sheet_name = "areas")    
+# barn_edges = pd.read_excel(path_zones+fna, sheet_name = "edges")
+del fna, path_zones
+
+
+# settings
+settings = {'barn' : [60,61], # [60,61,62,70,71,72,73],
+            'startdate' : date(2022,10,1),
+            'enddate' :  date(2022,10,10),
+            'cows' : [2658,2155], # or specific cow number
+            }
+
+# files that comply with settings
+fn = []
+for b in range(0,len(settings["barn"])):
+    print("barn = " + str(settings["barn"][b]))
+    if settings["cows"] == 0:
+        fbarn = [f for f in os.listdir(path + "/barn" + str(settings["barn"][b])) \
+                if os.path.isfile(os.path.join(path,"barn"+str(settings["barn"][b]),f)) \
+                    and (datetime.strptime(f[5:13], '%Y%m%d').date() >= settings["startdate"]) \
+                    and (datetime.strptime(f[5:13], '%Y%m%d').date() <= settings["enddate"])]
+        fbarn.sort()
+    else:
+        fbarn = [f for f in os.listdir(path + "/barn" + str(settings["barn"][b])) \
+                if os.path.isfile(os.path.join(path,"barn"+str(settings["barn"][b]),f)) \
+                    and (int(f[26:-4]) in settings["cows"]) \
+                    and (datetime.strptime(f[5:13], '%Y%m%d').date() >= settings["startdate"]) \
+                    and (datetime.strptime(f[5:13], '%Y%m%d').date() <= settings["enddate"])]
+        fbarn.sort()
+    fn.extend(fbarn)
+    fn.sort()
+
+# find unique cows
+cows = list(set([int(f[26:-4]) for f in fn]))
+
+# read data
+data = pd.DataFrame([])
+for f in fn:
+    barn = f[19:21]
+    sub = pd.read_csv(path + "/barn" + barn + "/" + f, 
+                      usecols = ["cowid","barn","date","t","xnew","ynew","area","zone","X","y"],
+                      dtype = {"cowid" : "int64","barn" : "int64","date" : "object",
+                               "t" : "int64", "xnew":"float64","ynew":"float64",
+                               "area":"object","zone":"float64","X":"float64","y" : "float64"})
+    sub["date"] = pd.to_datetime(sub["date"], format = "%Y-%m-%d") 
+    data = pd.concat([data,sub])
+data = data.sort_values(by = ["cowid","date","t"])
+data = data.reset_index(drop=1)
+
+# correct the areas
+cows = data["cowid"].drop_duplicates().reset_index(drop=1)
+for cow in cows:
+    print(cow)
+    days = data.loc[(data["cowid"]==cow),["date"]].drop_duplicates() 
+    days = days.sort_values(by = ["date"]).reset_index(drop=1)
+    for dd in days["date"]:
+        print(dd)
+        barnnr = int(data.loc[(data["cowid"] == cow) &
+                          (data["date"] == dd),"barn"].mean()) 
+        myarea = area_zones.loc[(area_zones["barn"] == barnnr) | \
+                            (area_zones["barn"].isna() == True) \
+                            ,:].reset_index(drop=1)
+        for j in range(0,len(myarea)):
+            x1,y1,x2,y2,x3,y3,x4,y4 = myarea[["x1","y1","x2","y2","x3","y3","x4","y4"]].values[j]
+            data.loc[((data["cowid"] == cow) & (data["date"] == dd) & \
+                      (data["xnew"] >= x1) & (data["xnew"] <= x2) & \
+                      (data["xnew"] >= x3) & (data["xnew"] <= x4) & \
+                      (data["ynew"] >= y1) & (data["ynew"] >= y2) & \
+                      (data["ynew"] <= y3) & (data["ynew"] <= y4)),"area"] = myarea["area"][j]
+            data.loc[((data["cowid"] == cow) & (data["date"] == dd) & \
+                      (data["xnew"] >= x1) & (data["xnew"] <= x2) & \
+                      (data["xnew"] >= x3) & (data["xnew"] <= x4) & \
+                      (data["ynew"] >= y1) & (data["ynew"] >= y2) & \
+                      (data["ynew"] <= y3) & (data["ynew"] <= y4)),"zone"] = myarea["zone"][j]
+        del x1,y1,x2,y2,x3,y3,x4,y4
+del cows, dd, days
+data.loc[data["zone"].isna(),"zone"] = 8
+data.loc[data["zone"] == 8 ,"area"] = "unknown"
+
+data["date"] = data["date"].dt.date
+
+# check correctly read ids
+print(data[["cowid","barn"]].drop_duplicates())
+
+
+#%% RESTING BEHAVIOUR
+
+cows = data["cowid"].drop_duplicates().reset_index(drop=1)
+for cow in cows:
+    subset = 
diff --git a/exploration.py b/exploration.py
index d7f5e529ecd64e327bbe4cd8ce235496d4cefb0b..9e4988f9cd80845a8d31de776987677ebc1b9fee 100644
--- a/exploration.py
+++ b/exploration.py
@@ -16,8 +16,9 @@ import pandas as pd
 import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib
+from matplotlib.patches import Rectangle
 import seaborn as sns
-from heatmap_visualisation import heatmap_barn
+from heatmap_visualisation import heatmap_barn,heatmap_barn2
 #%matplotlib qt
 
 #%% set paths and constants and load data
@@ -39,7 +40,7 @@ del fna, path_zones
 # settings
 settings = {'barn' : [60,61,62,70,71,72,73], # [60,61,62,70,71,72,73],
             'startdate' : date(2022,10,1),
-            'enddate' :  date(2022,10,15),
+            'enddate' :  date(2022,10,2),
             'cows' : [907,89,98,517,1748,495,3107], # or specific cow number
             }
 
@@ -260,8 +261,6 @@ for cow in cows:
 # ------------------------- individual cow heatmaps ---------------------------
 cows = data["cowid"].drop_duplicates().reset_index(drop=1)
 
-
-
 #barn limits
 barn = {"60":[64.8,75.6],
         "61":[54,64.8],
@@ -339,6 +338,11 @@ for cow in cows:
         print(dd)
         subset = data.loc[(data["cowid"]==cow) & \
                           (data["date"]==dd),:].copy()
+        # set limits (based on information barn)
+        barnnr = int(data.loc[(data["cowid"] == cow) &
+                          (data["date"] == dd),"barn"].mean())
+        x_lim = barn[str(barnnr)]
+        y_lim = [-18,0]
         
         fig, ax = plt.subplots(1, 2, gridspec_kw={'width_ratios': [3, 1.5]}, figsize = (22,7))
         # plot scatter + line
@@ -358,15 +362,10 @@ for cow in cows:
         ax[0].set_xticks([0,14400,28800,43200,57600,72000,86400])
         ax[0].set_xticklabels(['00:00','04:00','08:00','12:00','16:00','20:00','24:00'])
         ax[0].set_xlabel("time [h]")
-        ax[0].set_title("cow = " + str(cow) + ', date = ' + str(dd))
+        ax[0].set_title("barn = " +str(barnnr) + ', cow = ' + str(cow) + \
+                     ", date = " + str(dd))
         ax[0].set_xlim([-1000,87400])
         
-        # set limits (based on information barn)
-        barnnr = int(data.loc[(data["cowid"] == cow) &
-                          (data["date"] == dd),"barn"].mean())
-        x_lim = barn[str(barnnr)]
-        y_lim = [-18,0]
-        
         # plot heatmap
         ax[1] = heatmap_barn(subset["xnew"], subset["ynew"], 0.1, 0.1, x_lim, y_lim)
         ax[1].set_xlabel('y coordinate')
@@ -376,6 +375,37 @@ for cow in cows:
         ax[1].xaxis.tick_top()
         ax[1].yaxis.tick_right()
         ax[1].yaxis.set_label_position("right")
+        # ax[1].invert_xaxis()
         
+        # plot rectangles ~ barn areas
+        myarea = area_zones.loc[(area_zones["barn"] == barnnr) | \
+                            (area_zones["barn"].isna() == True) \
+                            ,:].reset_index(drop=1)
+        myarea = myarea.drop((myarea.loc[myarea["zone"]==7]).index)
+        myarea = myarea.drop((myarea.loc[myarea["zone"]==0]).index)
+        myarea = myarea.reset_index(drop=1)
+        yax = ax[1].get_ylim()
+        ycor = barn[str(barnnr)]
+        colors = []
+        for j in range(0,len(myarea)):
+            x1,y1,x2,y2,x3,y3,x4,y4 = myarea[["x1","y1","x2","y2","x3","y3","x4","y4"]].values[j]
+            color = myarea["area"][j]
+            X = y3*(-10)
+            W = (y2*(-10)) - X
+            Y = round((x1-ycor[0]),2)*10
+            H = round((x4-x3),2)*10
+            ax[1].add_patch(Rectangle((X,Y),W,H, color = palette[color],
+                                    fc = palette[color],
+                                    alpha = 0.15,
+                                    linestyle = 'dotted',
+                                    linewidth= 2
+                                    ))
+            ax[1].add_patch(Rectangle((X,Y),W,H, color = palette[color],
+                                    fc = 'none',
+                                    linestyle = 'dotted',
+                                    linewidth= 2
+                                    ))
+        del x2,y2,x3,y3,x4,y4,X,Y,H,W
+
         plt.savefig(svpath+"\\spatial_cow_" + str(cow) + "_" + str(dd).replace("-","") + ".png")
         plt.close()
\ No newline at end of file
diff --git a/heatmap_visualisation.py b/heatmap_visualisation.py
index d49d6d232408acbd8f814523b027bcaef25b0292..e9c7247a018745d2436457fae795f537479fdc40 100644
--- a/heatmap_visualisation.py
+++ b/heatmap_visualisation.py
@@ -115,7 +115,9 @@ def heatmap_barn(data_x,data_y,xstep,ystep,x_lim,y_lim):
                      yticklabels = round(1/xstep),
                      xticklabels = round(1/ystep*2),
                      robust = True,
-                     cbar = False
+                     cbar = False,
+                     cmap = "Greys_r"
+                     # alpha = 0.5
                      )
 
     return ax
@@ -179,3 +181,142 @@ for cow in cows["cowid"]:
         plt.close()
 """
 
+def heatmap_barn2(
+    data_x, data_y, xstep, ystep, x_lim, y_lim, image_data=False, **kwargs
+):
+    """
+    Plot the frequency of positions registered, given data_x and data_y
+    xstep and ystep determine how fine these positions are taken
+    x_lim and y_lim set boundaries (e.g. incl feeding lane)
+
+    Parameters
+    ----------
+    data_x : TYPE = Pandas series
+        Data of the x-coordinate as registered by the sewio system
+    data_y : TYPE = Pandas series
+        Data of the y-coordinate as registered by the sewio system
+    xstep : TYPE = float
+        How fine you want to determine location in x-direction (e.g. 1 is one
+            meter, 0.1 = 10 centimeters)
+    ystep : TYPE = float
+        How fine you want to determine location in x-direction (e.g. 1 is one
+            meter, 0.1 = 10 centimeters)
+    x_lim : TYPE = list of floats
+        Physical boundaries of plot you want to make (e.g. cut off feeding lane
+            then indicate lower boundary in y direction = -4)
+    y_lim : TYPE = list of floats
+        Physical boundaries of plot you want to make (e.g. cut off feeding lane
+            then indicate lower boundary in y direction = -4)
+    image_data : TYPE = bool
+        Do X,Y data originate from an image?
+
+    Returns
+    -------
+    ax : TYPE = sns heatmap
+        figure axes with heatmap
+
+    """
+    # transform data_x and data_y for count
+    x_transf = (data_x * (1 / xstep)).round() * xstep
+    y_transf = (data_y * (1 / ystep)).round() * ystep
+
+    # put x_transf and y_transf in DataFrame
+    frame = {"xT": x_transf.round(decimals=4), "yT": y_transf.round(decimals=4)}
+    df = pd.DataFrame(frame)
+
+    # count and groupby
+    heatdata = df.groupby(["xT", "yT"]).size()
+    heatdata = heatdata.to_frame(name="counts")
+    heatdata = heatdata.reset_index()
+
+    # set limits to agree with steps
+    x_low = round(round(x_lim[0] * (1 / xstep)) * xstep, ndigits=4)
+    x_high = round(round(x_lim[1] * (1 / xstep)) * xstep, ndigits=4)
+    y_low = round(round(y_lim[0] * (1 / ystep)) * ystep, ndigits=4)
+    y_high = round(round(y_lim[1] * (1 / ystep)) * ystep, ndigits=4)
+
+    # delete the rows for which values are outside x_lim or y_lim
+    heatdata = heatdata.loc[
+        (heatdata["xT"] >= x_lim[0])
+        & (heatdata["xT"] <= x_lim[1])
+        & (heatdata["yT"] >= y_lim[0])
+        & (heatdata["yT"] <= y_lim[1])
+    ]
+    heatdata.reset_index(drop=True, inplace=True)
+
+    # expected axes limits based on lim and step
+    x = np.arange(x_low, x_high + xstep, xstep, dtype=float)
+    y = np.arange(y_low, y_high + ystep, ystep, dtype=float)
+
+    # Create rectangular dataframe (x_high, y_high)
+    if image_data:
+        temp_holder = pd.DataFrame(
+            data={
+                "xT": pd.Series(
+                    [
+                        int(i)
+                        for i in np.arange(x_low, x_high + 1, xstep, dtype=float)
+                        for _ in np.arange(y_low, y_high + 1, ystep, dtype=int)
+                    ]
+                ),
+                "yT": np.tile(
+                    np.arange(y_low, y_high + 1, ystep, dtype=int),
+                    len(np.arange(x_low, x_high + 1, xstep, dtype=float)),
+                ),
+            }
+        )
+    else:
+        temp_holder = pd.DataFrame(
+            data={
+                "xT": pd.Series(
+                    [
+                        i
+                        for i in np.arange(x_low, x_high, xstep, dtype=float)
+                        for _ in np.arange(y_low, y_high, ystep, dtype=float)
+                    ]
+                ),
+                "yT": np.tile(
+                    np.arange(y_low, y_high, ystep, dtype=float),
+                    len(np.arange(x_low, x_high, xstep, dtype=float)),
+                ),
+            }
+        )
+        
+    # 
+    temp_holder["yT"] = round(temp_holder["yT"],1)
+    temp_holder["xT"] = round(temp_holder["xT"],1)
+    heatdata["yT"] = round(heatdata["yT"],1)
+    heatdata["xT"] = round(heatdata["xT"],1)
+    
+
+    # merge into frame
+    tdf = pd.merge(temp_holder, heatdata, how="left", on=["xT", "yT"])
+    tdf["counts"] = tdf["counts"].fillna(0)
+    tdf = tdf.pivot(index="yT", columns="xT", values="counts").T
+    
+    # convert to df with right colnames and indices
+    xcols = np.array(x.round(decimals=3),dtype = str)
+    ycols = np.flip(np.array(y.round(decimals=3),dtype = str))
+    # df = pd.DataFrame(tdf, columns = ycols, index = xcols)
+    
+
+    defaultkwargs = {'yticklabels' : round(1 / ystep),
+                     'xticklabels' : round(1 / xstep * 2),
+                     'cmap' : plt.cm.YlOrRd_r,}
+
+    kwargs = {**defaultkwargs, **kwargs}
+
+    # make heatmap
+    ax = sns.heatmap(
+        tdf,
+        yticklabels = round(1/xstep),
+        xticklabels = round(1/ystep*2),
+        robust = True,
+        cbar = False,
+        cmap = "Greys"
+    )
+
+
+
+
+    return ax