diff --git a/skais/ais/ais_trajectory.py b/skais/ais/ais_trajectory.py
index f4ec956b6930cca43121640fc4e2f93a698ce915..336f6da6d695870ec2f2af9dddd203042fde6d1c 100644
--- a/skais/ais/ais_trajectory.py
+++ b/skais/ais/ais_trajectory.py
@@ -60,13 +60,125 @@ def apply_time_sequence(dat, time, func):
     return result
 
 
+@jit(nopython=True)
 def __get_image_value__(features, bounds):
-    value = []
+    value = np.zeros((len(features)))
+    i = 0
     for f, b in zip(features, bounds):
-        value.append(1 - (b[1] - f - b[0]) / (b[1] - b[0]))
+        value[i] = (1 - (b[1] - f - b[0]) / (b[1] - b[0]))
+        i = i + 1
     return value
 
 
+def __get_bounding_box__(bounding_box, positions, ref_index):
+    if bounding_box == 'fit':
+        lower_lon, upper_lon = (min(positions[:, 0]), max(positions[:, 0]))
+        lower_lat, upper_lat = (min(positions[:, 1]), max(positions[:, 1]))
+    elif bounding_box == 'centered':
+        center_lon, center_lat = positions[ref_index]
+        min_lon, max_lon = (min(positions[:, 0]), max(positions[:, 0]))
+        min_lat, max_lat = (min(positions[:, 1]), max(positions[:, 1]))
+
+        distance_to_center = max(center_lon - min_lon, max_lon - center_lon, center_lat - min_lat,
+                                 max_lat - center_lat)
+
+        upper_lat = center_lat + distance_to_center
+        lower_lat = center_lat - distance_to_center
+        upper_lon = center_lon + distance_to_center
+        lower_lon = center_lon - distance_to_center
+    elif type(bounding_box) is list:
+        if type(bounding_box[0]) is not numbers.Number:
+            upper_lon = bounding_box[1][0]
+            lower_lon = bounding_box[0][0]
+            upper_lat = bounding_box[1][1]
+            lower_lat = bounding_box[0][1]
+        else:
+            center_lon, center_lat = positions[ref_index]
+            distance_to_center_lon = bounding_box[0]
+            distance_to_center_lat = bounding_box[1]
+            upper_lat = center_lat + distance_to_center_lat
+            lower_lat = center_lat - distance_to_center_lat
+            upper_lon = center_lon + distance_to_center_lon
+            lower_lon = center_lon - distance_to_center_lon
+    else:
+        raise ValueError(f"Option not supported: {bounding_box}")
+
+    if lower_lat == upper_lat:
+        lower_lat -= 1
+        upper_lat += 1
+    if lower_lon == upper_lon:
+        lower_lon -= 1
+        upper_lon += 1
+
+    return lower_lon, upper_lon, lower_lat, upper_lat
+
+
+@jit(nopython=True)
+def generate_points(data, positions, height, width, node_size, lower_lat, upper_lat, lower_lon, upper_lon):
+    for longitude, latitude in positions:
+        x_coord, y_coord = get_coord(latitude, longitude, height, width, lower_lat, upper_lat, lower_lon,
+                                     upper_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]
+
+
+@jit(nopython=True)
+def generate_links(data, positions, height, width, lower_lat, upper_lat, lower_lon, upper_lon):
+    lon, lat = positions[0, 0], positions[0, 1]
+    for longitude, latitude in positions[1:]:
+        x_prv, y_prev = get_coord(lat, lon, height, width, lower_lat, upper_lat, lower_lon, upper_lon)
+        x_nxt, y_nxt = get_coord(latitude, longitude, height, width, lower_lat, upper_lat, lower_lon,
+                                 upper_lon)
+
+        lon, lat = longitude, latitude
+        for x, y in bresenham(x_prv, y_prev, x_nxt, y_nxt):
+            data[x, y] = [1]
+
+
+@jit(nopython=True)
+def generate_points_with_features(data, positions, features_vectors, bounds, node_size, height, width,
+                                  lower_lat, upper_lat, lower_lon, upper_lon, link):
+    for pos, f in zip(positions, features_vectors):
+        latitude = pos[1]
+        longitude = pos[0]
+        x_coord, y_coord = get_coord(latitude, longitude, height, width, lower_lat, upper_lat, lower_lon,
+                                     upper_lon)
+        value = __get_image_value__(f, bounds)
+        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):
+                for i, v in enumerate(value):
+                    data[x, y, i] = v
+    if link:
+        lon, lat = positions[0, 0], positions[0, 1]
+
+        value = __get_image_value__(features_vectors[0], bounds)
+        for pos, f in zip(positions[1:], features_vectors[1:]):
+            latitude = pos[1]
+            longitude = pos[0]
+            x_prv, y_prev = get_coord(lat, lon, height, width, lower_lat, upper_lat, lower_lon, upper_lon)
+            x_nxt, y_nxt = get_coord(latitude, longitude, height, width, lower_lat, upper_lat, lower_lon,
+                                     upper_lon)
+            lon, lat = longitude, latitude
+            for x, y in bresenham(x_prv, y_prev, x_nxt, y_nxt):
+                for i, v in enumerate(value):
+                    data[x, y, i] = v
+            value = __get_image_value__(f, bounds)
+
+
 class AISTrajectory(AISPoints):
     def __init__(self, df, mmsi=0, interpolation_time=None):
         df = df.drop_duplicates(subset=['ts_sec'])
@@ -233,125 +345,44 @@ class AISTrajectory(AISPoints):
                 result.append((row['ts_sec'], current_label))
         return result
 
-    def generate_array_from_positions(self, height=256, width=256, link=True, bounding_box='fit', ref_index=-1, features=None,
-                                      node_size=0):
-        nb_channels = 1
+    def generate_array_from_positions(self, height=256, width=256, link=True, bounding_box='fit', ref_index=-1,
+                                      features=None, node_size=0):
 
         positions = self.df[['longitude', 'latitude']].to_numpy()
-        if bounding_box == 'fit':
-            lower_lon, upper_lon = (min(positions[:, 0]), max(positions[:, 0]))
-            lower_lat, upper_lat = (min(positions[:, 1]), max(positions[:, 1]))
-        elif bounding_box == 'centered':
-            center_lon, center_lat = positions[ref_index]
-            min_lon, max_lon = (min(positions[:, 0]), max(positions[:, 0]))
-            min_lat, max_lat = (min(positions[:, 1]), max(positions[:, 1]))
-
-            distance_to_center = max(center_lon - min_lon, max_lon - center_lon, center_lat - min_lat,
-                                     max_lat - center_lat)
-
-            upper_lat = center_lat + distance_to_center
-            lower_lat = center_lat - distance_to_center
-            upper_lon = center_lon + distance_to_center
-            lower_lon = center_lon - distance_to_center
-        elif type(bounding_box) is list:
-            if type(bounding_box[0]) is not numbers.Number:
-                upper_lon = bounding_box[1][0]
-                lower_lon = bounding_box[0][0]
-                upper_lat = bounding_box[1][1]
-                lower_lat = bounding_box[0][1]
-            else:
-                center_lon, center_lat = positions[ref_index]
-                distance_to_center_lon = bounding_box[0]
-                distance_to_center_lat = bounding_box[1]
-                upper_lat = center_lat + distance_to_center_lat
-                lower_lat = center_lat - distance_to_center_lat
-                upper_lon = center_lon + distance_to_center_lon
-                lower_lon = center_lon - distance_to_center_lon
-        else:
-            raise ValueError(f"Option not supported: {bounding_box}")
 
-        if lower_lat == upper_lat:
-            lower_lat -= 1
-            upper_lat += 1
-        if lower_lon == upper_lon:
-            lower_lon -= 1
-            upper_lon += 1
+        lower_lon, upper_lon, lower_lat, upper_lat = __get_bounding_box__(bounding_box, positions, ref_index)
 
+        bounds = []
         if features is None:
-            data = np.zeros((height, width, nb_channels), dtype=np.uint8)
-            for longitude, latitude in positions:
-                x_coord, y_coord = get_coord(latitude, longitude, height, width, lower_lat, upper_lat, lower_lon,
-                                             upper_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, lower_lat, upper_lat, lower_lon, upper_lon)
-                    x_nxt, y_nxt = get_coord(latitude, longitude, height, width, lower_lat, upper_lat, lower_lon,
-                                             upper_lon)
-
-                    lon, lat = longitude, latitude
-                    for x, y in bresenham(x_prv, y_prev, x_nxt, y_nxt):
-                        data[x, y] = [1]
-
+            features_vectors = None
+        elif type(features) is list:
+            features_vectors = self.df[features].to_numpy()
+            for c in features_vectors.T:
+                bounds.append((0, max(c)))
+        elif type(features) is str:
+            features_vectors = self.df[[features]].to_numpy()
+            for c in features_vectors.T:
+                bounds.append((0, max(c)))
+        elif type(features) is dict:
+            bounds = list(features.values())
+            features_vectors = self.df[features.keys()].to_numpy()
         else:
-            bounds = []
-            if type(features) is list:
-                features_vectors = self.df[features].to_numpy()
-                for c in features_vectors.T:
-                    bounds.append((0, max(c)))
-            elif type(features) is str:
-                features_vectors = self.df[[features]].to_numpy()
-                for c in features_vectors.T:
-                    bounds.append((0, max(c)))
-            elif type(features) is dict:
-                bounds = list(features.values())
-                features_vectors = self.df[features.keys()].to_numpy()
-            else:
-                raise TypeError("Type not supported")
+            raise TypeError("Type not supported")
 
+        if features_vectors is not None:
             nb_channels = len(features_vectors.T)
-            data = np.zeros((height, width, nb_channels), dtype=np.float)
-
+        else:
+            nb_channels = 1
+        data = np.zeros((height, width, nb_channels), dtype=np.float)
 
-            for pos, f in zip(positions, features_vectors):
-                latitude = pos[1]
-                longitude = pos[0]
-                x_coord, y_coord = get_coord(latitude, longitude, height, width, lower_lat, upper_lat, lower_lon,
-                                             upper_lon)
-                value = __get_image_value__(f, bounds)
-                x_lower_bound = max(0, x_coord - node_size)
-                x_upper_bound = min(height - 1, x_coord + node_size)
+        if features_vectors is None:
 
-                y_lower_bound = max(0, y_coord - node_size)
-                y_upper_bound = min(width - 1, y_coord + node_size)
+            generate_points(data, positions, height, width, node_size, lower_lat, upper_lat, lower_lon, upper_lon)
 
-                for x in range(x_lower_bound, x_upper_bound + 1):
-                    for y in range(y_lower_bound, y_upper_bound + 1):
-                        for i, v in enumerate(value):
-                            data[x, y, i] = v
             if link:
-                lon, lat = positions[0, 0], positions[0, 1]
-
-                value = __get_image_value__(features_vectors[0], bounds)
-                for pos, f in zip(positions[1:], features_vectors[1:]):
-                    latitude = pos[1]
-                    longitude = pos[0]
-                    x_prv, y_prev = get_coord(lat, lon, height, width, lower_lat, upper_lat, lower_lon, upper_lon)
-                    x_nxt, y_nxt = get_coord(latitude, longitude, height, width, lower_lat, upper_lat, lower_lon,
-                                             upper_lon)
-                    lon, lat = longitude, latitude
-                    for x, y in bresenham(x_prv, y_prev, x_nxt, y_nxt):
-                        for i, v in enumerate(value):
-                            data[x, y, i] = v
-                    value = __get_image_value__(f, bounds)
+                generate_links(data, positions, height, width, lower_lat, upper_lat, lower_lon, upper_lon)
+
+        else:
+            generate_points_with_features(data, positions, features_vectors, bounds, node_size, height, width,
+                                          lower_lat, upper_lat, lower_lon, upper_lon, link)
         return data
diff --git a/skais/utils/geography.py b/skais/utils/geography.py
index de21a9aed417ba8e959aa0bb5a5780f9b879634f..ae539e1ddac0640c78c1e2b266601db67c2dd251 100644
--- a/skais/utils/geography.py
+++ b/skais/utils/geography.py
@@ -23,6 +23,7 @@ def great_circle(lat1, lat2, long1, long2):
     return d
 
 
+@jit(nopython=True)
 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)
diff --git a/skais/utils/geometry.py b/skais/utils/geometry.py
index 0919af4a5b6dd1501a9b5fcb0f3bf746a336cec3..7a111d0f32f291f673eff5dca58607c80e8fb446 100644
--- a/skais/utils/geometry.py
+++ b/skais/utils/geometry.py
@@ -1,3 +1,7 @@
+from numba import jit
+
+
+@jit(nopython=True)
 def bresenham(x1, y1, x2, y2):
     dx = int(x2 - x1)
     dy = int(y2 - y1)