diff --git a/source/campo/op_experimental/randommove.py b/source/campo/op_experimental/randommove.py
new file mode 100644
index 0000000..0bef1d4
--- /dev/null
+++ b/source/campo/op_experimental/randommove.py
@@ -0,0 +1,49 @@
+import numpy as np
+import random 
+from pathlib import Path 
+import sys
+from field_agent.field_agent_interactions import raster_values_to_feature, zonal_values_to_feature
+import campo
+import math
+import matplotlib.pyplot as plt
+
+def randommove_to_boolean (boolean_fieldprop, field_pset, point_prop):
+    '''
+    - boolean_fieldprop = a property that is a field and contains true values for places where an agent may move to   
+    - field_pset = a property set describing the domain of the field of the study area (Type: propertyset)
+    - point_propset = the property set that has the move (= point agents)
+    - nragents = the number of point agents 
+    return lists of coordinates with the same length as the number of agents in a point propertyset (sueful wehn you want to make them move there) 
+    '''
+    # need to flip again because when calculating with this map points have different orientation than field (see rerasterize)
+    if isinstance(boolean_fieldprop, np.ndarray):
+        map_flipped = np.flip (boolean_fieldprop, axis=0)
+    elif isinstance (boolean_fieldprop, campo.property.Property):
+        map_flipped = np.flip (boolean_fieldprop.values()[0], axis=0) 
+    else:
+        raise TypeError ('boolean_fieldprop needs to be of type campo.property or of a numpy array with same dimensions as field property values')
+    
+    nragents = len (point_prop.values().values.values())
+    for fidx, area in enumerate (field_pset.space_domain): 
+        nr_cols = int(area[5])
+        xmin = area [0]
+        ymin = area [1] 
+        resolution = math.fabs (area[2] - xmin)/nr_cols
+    # finding the indices of the places where the fieldcondition is true 
+    coords_idx = np.argwhere (map_flipped) #coordinates of the spawning grounds in [y,x]
+    # collecting the coordinate combination in a tuple so as to prevent them from being 'disconnected' from eachother 
+    coords_list = [tuple(row) for row in coords_idx] 
+    random_newindex = random.sample (coords_list, nragents) # nr of agents is the subsetsize
+    # seperating the tuples in the list in two seperate list
+    yindex, xindex = zip(*random_newindex) #assuming that tuple list is reversed, first gives y then x as in column = x and row = y 
+    xindex = np.array (xindex)
+    yindex = np.array (yindex)
+    xcoords = np.zeros(nragents)
+    ycoords = np.zeros(nragents)
+
+    # make from x index a x coordinate by using resolution and bounding box information
+    for i, xvalue in enumerate(xindex):
+        xcoords [i] = (xvalue*resolution + xmin) 
+    for j, yvalue in enumerate(yindex): 
+        ycoords[j] = (yvalue *resolution + ymin)
+    return xcoords, ycoords 
\ No newline at end of file
diff --git a/source/campo/op_experimental/zonal.py b/source/campo/op_experimental/zonal.py
new file mode 100644
index 0000000..8531420
--- /dev/null
+++ b/source/campo/op_experimental/zonal.py
@@ -0,0 +1,58 @@
+from osgeo import gdal
+from osgeo import osr
+from osgeo import ogr
+
+gdal.UseExceptions()
+
+import math
+import numpy as np
+
+from campo.points import Points
+from campo.areas import Areas
+from campo.property import Property
+
+def zonal_values_to_feature(point_pset, field_pset, field_prop_class, field_prop_var, operation):
+  ''' Queries aggregative raster values of a zone of a field property based on the location
+    of point agents within a classification map, given a certain aggregative operation. 
+    Writes this as a new point agent property. Only works for one fieldagent existing at the location. 
+    Fieldagents with overlapping domains cannot generate an output
+    E.g.: According to both the agent's location and a soil map (field_prop_class), the agent is positioned in 
+        soil '2' , which is clay. With the operation 'mean', the mean rainfall (from field_prop_var) 
+        is calculated and attributed to the agent
+  Parameters:
+    point_pset: property set of the point agents to be attributed the location 
+    field_pset: property set of a single field 
+    field_prop_class: property describing classes or groups of cells of which the zonal extent shall be the windowsize of the aggregration
+    field_prop_var: property describing the variable which needs to be aggregated,
+      in case of a boolean map, 'True' values are 1
+    operation: operation: aggregative numpy operation as a string ('sum', 'mean', 'min', 'max', 'etc')
+  Returns: 
+    point_prop: property of point agents with values of aggregated field values'''
+  
+  # Generate operator, first checking if operation is available in numpy
+  if not hasattr (np, operation):
+    raise ValueError (f'Unsupported numpy operation, {operation}')
+  operator = getattr(np, operation)
+
+  # Identifying the zone the agent is in
+  agents_zoneIDs = raster_values_to_feature(point_pset,field_pset, field_prop_class)
+  # Generate empty point property given point property space definitions
+  point_prop = Property('emptycreatename', point_pset.uuid, point_pset.space_domain, point_pset.shapes)
+  # make as many field properties as there are agents:  
+  # Loop over space attributes of the different points in the point agents propertyset
+  for pidx, ID in enumerate(agents_zoneIDs.values()):
+    # Making a boolean map concerning the extent of the zone for each agent
+    aggr_zone_var = np.zeros(1)
+    for fidx,area in enumerate(field_pset.space_domain):
+        zone_extent = np.where (field_prop_class.values()[fidx] == ID, 1, 0)
+        variable_zone_array = np.multiply (zone_extent, field_prop_var.values()[fidx])
+        # we don't need to flip this time, since the raster_values_to_feature already gave 
+        # the right topological relationship between the field and the agent: 
+        # the zone_extent array describes the zone in which the agent is positioned. 
+        # This array might be flipped, but this won't lead to any different outcomes of aggregative operations
+        aggr_zone_var[fidx] = operator (variable_zone_array)
+    #field_prop_array [pidx] = field_prop# array as long as the number of agents filled with a field prop for each agent
+    # Write the value to the point property for each point agent
+    point_prop.values()[pidx] = aggr_zone_var.item() 
+    
+  return point_prop
\ No newline at end of file