diff --git a/skais/ais/ais_trajectory.py b/skais/ais/ais_trajectory.py
index f8364722878ef74364cb05f363da1ed002425b83..e123aff7898a03ad1eff6d03c94f70b3734d1e71 100644
--- a/skais/ais/ais_trajectory.py
+++ b/skais/ais/ais_trajectory.py
@@ -5,8 +5,9 @@ import numpy as np
 from numba import jit
 from scipy.interpolate import interp1d
 
-from skais.utils.geography import great_circle, position_from_distance
+from skais.utils.geography import great_circle, position_from_distance, get_coord
 from skais.ais.ais_points import AISPoints
+from skais.utils.geometry import bresenham
 
 
 @jit(nopython=True)
@@ -41,7 +42,7 @@ def apply_func_on_window(dat, func, radius, on_edge='copy'):
         return result
     elif on_edge == 'ignore':
         for i in range(0, dat.shape[0]):
-            lower_bound = max(0, i-radius)
+            lower_bound = max(0, i - radius)
             upper_bound = min(dat.shape[0], i + radius + 1)
             data = dat[lower_bound:upper_bound]
             result[i] = func(data)
@@ -222,4 +223,39 @@ class AISTrajectory(AISPoints):
             if current_label != row[label_column]:
                 current_label = row[label_column]
                 result.append((row['ts_sec'], current_label))
-        return result
\ No newline at end of file
+        return result
+
+    def generate_array_from_positions(self, height=256, width=256, link=True, bounding_box='fit', features=None,
+                                      node_size=0):
+        nb_channels = 1
+        if features is not None:
+            nb_channels = len(features)
+        data = np.zeros((height, width, nb_channels), dtype=np.uint8)
+        if bounding_box != 'fit':
+            raise ValueError("feature not implemented")
+        positions = self.df[['longitude', 'latitude']].to_numpy()
+        min_lon, max_lon = (min(positions[:, 0]), max(positions[:, 0]))
+        min_lat, max_lat = (min(positions[:, 1]), max(positions[:, 1]))
+        for longitude, latitude in positions:
+            x_coord, y_coord = get_coord(latitude, longitude, height, width, min_lat, max_lat, min_lon, max_lon)
+
+            x_lower_bound = max(0, x_coord - node_size)
+            x_upper_bound = min(height - 1, x_coord + node_size)
+
+            y_lower_bound = max(0, y_coord - node_size)
+            y_upper_bound = min(width - 1, y_coord + node_size)
+
+            for x in range(x_lower_bound, x_upper_bound + 1):
+                for y in range(y_lower_bound, y_upper_bound + 1):
+                    data[x, y] = [1]
+
+        if link:
+            lon, lat = positions[0, 0], positions[0, 1]
+            for longitude, latitude in positions[1:]:
+                x_prv, y_prev = get_coord(lat, lon, height, width, min_lat, max_lat, min_lon, max_lon)
+                x_nxt, y_nxt = get_coord(latitude, longitude, height, width, min_lat, max_lat, min_lon, max_lon)
+
+                lon, lat = longitude, latitude
+                for x, y in bresenham(x_prv, y_prev, x_nxt, y_nxt):
+                    data[x, y] = [1]
+        return data
diff --git a/skais/tests/ais/test_ais_trajectory.py b/skais/tests/ais/test_ais_trajectory.py
index bea4d8843011eaabc1ae0b81bf1dc023b3e62706..fe8efc4c087776701515b8036ae9ecf4939e3b27 100644
--- a/skais/tests/ais/test_ais_trajectory.py
+++ b/skais/tests/ais/test_ais_trajectory.py
@@ -40,7 +40,7 @@ class TestAISTrajectory(unittest.TestCase):
 
         for r, e in zip(result, expected):
             np.testing.assert_array_equal(r, e)
-            
+
     def test_sliding_window_too_short(self):
         trajectory = AISTrajectory(
             pd.DataFrame(
@@ -319,7 +319,7 @@ class TestAISTrajectory(unittest.TestCase):
         self.assertEqual(0, compute_trajectory.py_func(times, 800))
 
     def test_apply_func_on_window(self):
-        self.assertRaises(ValueError, apply_func_on_window,np.arange(10), 0, 0, 'not valid string')
+        self.assertRaises(ValueError, apply_func_on_window, np.arange(10), 0, 0, 'not valid string')
 
     def test_apply_func_on_window_ignore(self):
         result = apply_func_on_window(np.arange(10), np.mean, 2, 'ignore')
@@ -369,7 +369,7 @@ class TestAISTrajectory(unittest.TestCase):
         trajectory = AISTrajectory(
             pd.DataFrame(
                 {
-                    "label": [1 for _ in range(11)] + [2 for _ in range(10)]+ [1 for _ in range(10)],
+                    "label": [1 for _ in range(11)] + [2 for _ in range(10)] + [1 for _ in range(10)],
                     "ts_sec": [i for i in range(0, 18001, 600)]
                 }
             )
@@ -378,4 +378,79 @@ class TestAISTrajectory(unittest.TestCase):
         result = trajectory.get_time_per_label_shift()
         expected = [(0, 1), (6600, 2), (12600, 1)]
 
-        self.assertListEqual(result, expected)
\ No newline at end of file
+        self.assertListEqual(result, expected)
+
+    def test_generate_array_from_positions(self):
+        trajectory = AISTrajectory(
+            pd.DataFrame(
+                {
+                    "latitude": [0, 10, 0, -10],
+                    "longitude": [0, 10, 10, -10],
+                    "ts_sec": [i for i in range(4)]
+                }
+            )
+        )
+
+        result = trajectory.generate_array_from_positions(height=9, width=9, link=False, bounding_box='fit',
+                                                          features=None, node_size=0).reshape((9, 9))
+        expected = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 1],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0],
+                             [0, 0, 0, 0, 1, 0, 0, 0, 1],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0],
+                             [1, 0, 0, 0, 0, 0, 0, 0, 0]])
+
+        np.testing.assert_array_equal(result, expected)
+
+    def test_generate_array_from_positions_node_size(self):
+        trajectory = AISTrajectory(
+            pd.DataFrame(
+                {
+                    "latitude": [0, 10, 0, -10],
+                    "longitude": [0, 10, 10, -10],
+                    "ts_sec": [i for i in range(4)]
+                }
+            )
+        )
+
+        result = trajectory.generate_array_from_positions(height=9, width=9, link=False, bounding_box='fit',
+                                                          features=None, node_size=1).reshape((9, 9))
+        expected = np.array([[0, 0, 0, 0, 0, 0, 0, 1, 1],
+                             [0, 0, 0, 0, 0, 0, 0, 1, 1],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0],
+                             [0, 0, 0, 1, 1, 1, 0, 1, 1],
+                             [0, 0, 0, 1, 1, 1, 0, 1, 1],
+                             [0, 0, 0, 1, 1, 1, 0, 1, 1],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0],
+                             [1, 1, 0, 0, 0, 0, 0, 0, 0],
+                             [1, 1, 0, 0, 0, 0, 0, 0, 0]])
+
+        np.testing.assert_array_equal(result, expected)
+
+    def test_generate_array_from_positions_with_line(self):
+        trajectory = AISTrajectory(
+            pd.DataFrame(
+                {
+                    "latitude": [0, 10, 0, 20],
+                    "longitude": [0, 10, 20, 20],
+                    "ts_sec": [i for i in range(4)]
+                }
+            )
+        )
+
+        result = trajectory.generate_array_from_positions(height=9, width=18, link=True, bounding_box='fit',
+                                                          features=None, node_size=0).reshape((9, 18))
+        expected = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+                             [0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1],
+                             [0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1],
+                             [0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1],
+                             [0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1],
+                             [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1]])
+
+        np.testing.assert_array_equal(result, expected)
\ No newline at end of file
diff --git a/skais/tests/utils/test_geometry.py b/skais/tests/utils/test_geometry.py
new file mode 100644
index 0000000000000000000000000000000000000000..7a0b76ddb703f43a136c315e5abcf1833460a03d
--- /dev/null
+++ b/skais/tests/utils/test_geometry.py
@@ -0,0 +1,42 @@
+import unittest
+
+from skais.utils.geometry import bresenham
+
+
+class TestGeometry(unittest.TestCase):
+    def test_bresenham(self):
+        result = bresenham(3, 4, 16, 9)
+        expected = [(3, 4), (4, 4), (5, 5), (6, 5), (7, 6), (8, 6), (9, 6), (10, 7), (11, 7), (12, 7), (13, 8), (14, 8),
+                    (15, 9), (16, 9)]
+
+        self.assertListEqual(result, expected)
+
+    def test_bresenham_inverted(self):
+        result = bresenham(16, 9, 3, 4)
+        expected = [(3, 4), (4, 4), (5, 5), (6, 5), (7, 6), (8, 6), (9, 6), (10, 7), (11, 7), (12, 7), (13, 8),
+                    (14, 8), (15, 9), (16, 9)]
+
+        self.assertListEqual(result, expected)
+
+
+    def test_bresenham_inverted_2(self):
+        result = bresenham(16, 4, 3, 9)
+        expected = [(3, 9), (4, 9), (5, 8), (6, 8), (7, 7), (8, 7), (9, 7), (10, 6), (11, 6), (12, 6), (13, 5), (14, 5),
+                    (15, 4), (16, 4)]
+        self.assertListEqual(result, expected)
+
+    def test_bresenham_same_line(self):
+        result = bresenham(3, 4, 10, 4)
+        expected = [(3, 4), (4, 4), (5, 4), (6, 4), (7, 4), (8, 4), (9, 4), (10, 4)]
+
+        self.assertListEqual(result, expected)
+
+    def test_bresenham_same_column(self):
+        result = bresenham(3, 4, 3, 10)
+        expected = [(3, 4), (3, 5), (3, 6), (3, 7), (3, 8), (3, 9), (3, 10)]
+
+        self.assertListEqual(result, expected)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/skais/utils/geography.py b/skais/utils/geography.py
index 79b110194c7da6bd332b310c1e46a70beb93a2dd..de21a9aed417ba8e959aa0bb5a5780f9b879634f 100644
--- a/skais/utils/geography.py
+++ b/skais/utils/geography.py
@@ -23,6 +23,13 @@ def great_circle(lat1, lat2, long1, long2):
     return d
 
 
+def get_coord(lat, lon, height, width, min_lat, max_lat, min_lon, max_lon):
+    x_coord = max(min(height - int(height * (lat - min_lat) / (max_lat - min_lat)) - 1, height - 1), 0)
+    y_coord = max(min(int((width - 1) * (lon - min_lon) / (max_lon - min_lon)), width - 1), 0)
+
+    return x_coord, y_coord
+
+
 def position_from_distance(position, distances):
     lat = np.arcsin(
         np.sin(np.radians(position[0])) * np.cos(distances[0] / R) + np.cos(np.radians(position[0])) * np.sin(
diff --git a/skais/utils/geometry.py b/skais/utils/geometry.py
new file mode 100644
index 0000000000000000000000000000000000000000..adbd8a58b0d66373ebb2eecd6ca8f3bc8aaecefe
--- /dev/null
+++ b/skais/utils/geometry.py
@@ -0,0 +1,52 @@
+def bresenham(x1, y1, x2, y2):
+    dx = int(x2 - x1)
+    dy = int(y2 - y1)
+
+    sx = sy = 1
+    if dx < 0:
+        sx = -1
+    if dy < 0:
+        sy = -1
+
+    x = x1
+    y = y1
+
+    pixels = [(x1, y1)]
+    if abs(dx) > abs(dy): # slope < 1
+        if x1 > x2:
+            tmp = x2
+            x2 = x1
+            x1 = tmp
+
+            tmp = y2
+            y2 = y1
+            y1 = tmp
+        p = (2 * abs(dy)) - abs(dx)
+        
+        for x in range(x1 + 1, x2 + 1):
+            if p < 0:
+                p += 2 * abs(dy)
+            else:
+                y += sy
+                p += (2 * abs(dy)) - (2 * abs(dx))
+            pixels.append((x, y))
+    else: # slope >= 1
+        if y1 > y2:
+            tmp = x2
+            x2 = x1
+            x1 = tmp
+
+            tmp = y2
+            y2 = y1
+            y1 = tmp
+        p = (2 * abs(dx)) - abs(dy)
+        for y in range(y1 + 1, y2 + 1):
+            if p < 0:
+                p += 2 * abs(dx)
+            else:
+                x += sx
+                p += (2 * abs(dx)) - (2 * abs(dy))
+            pixels.append((x, y))
+
+    return pixels
+