diff --git a/skais/ais/ais_trajectory.py b/skais/ais/ais_trajectory.py
index db1cd8f7311fbdd1e17cb2e6afc863d1bdc997c6..e6ba885dcafca82f0a023eee0196003472ce28b4 100644
--- a/skais/ais/ais_trajectory.py
+++ b/skais/ais/ais_trajectory.py
@@ -59,6 +59,13 @@ def apply_time_sequence(dat, time, func):
     return result
 
 
+def __get_image_value__(features, bounds):
+    value = []
+    for f, b in zip(features, bounds):
+        value.append(1 - (b[1] - f) / b[1])
+    return value
+
+
 class AISTrajectory(AISPoints):
     def __init__(self, df, mmsi=0, interpolation_time=None):
         df = df.drop_duplicates(subset=['ts_sec'])
@@ -228,9 +235,7 @@ class AISTrajectory(AISPoints):
     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()
@@ -243,26 +248,70 @@ class AISTrajectory(AISPoints):
             min_lon -= 1
             max_lon += 1
 
-        for longitude, latitude in positions:
-            x_coord, y_coord = get_coord(latitude, longitude, height, width, min_lat, max_lat, min_lon, max_lon)
+        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, 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)
+                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)
+                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]
+                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)
 
-        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]
 
-                lon, lat = longitude, latitude
-                for x, y in bresenham(x_prv, y_prev, x_nxt, y_nxt):
-                    data[x, y] = [1]
+        else:
+            if type(features) is list:
+                nb_channels = len(features)
+            elif type(features) is str:
+                features = [features]
+            else:
+                raise TypeError("Type not supported")
+            data = np.zeros((height, width, nb_channels), dtype=np.float)
+            features_vectors = self.df[features].to_numpy()
+            bounds = []
+            for c in features_vectors.T:
+                bounds.append((min(c), max(c)))
+
+            for pos, f in zip(positions, features_vectors):
+                latitude = pos[1]
+                longitude = pos[0]
+                x_coord, y_coord = get_coord(latitude, longitude, height, width, min_lat, max_lat, min_lon, max_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, 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):
+                        for i, v in enumerate(value):
+                            data[x, y, i] = v
+                    value = __get_image_value__(f, bounds)
         return data
diff --git a/skais/tests/ais/test_ais_trajectory.py b/skais/tests/ais/test_ais_trajectory.py
index 9c92727e540588e34886f79fc46aa817a7a7c244..d5e64a0e87a67f15f19289c2f1d527f215b178db 100644
--- a/skais/tests/ais/test_ais_trajectory.py
+++ b/skais/tests/ais/test_ais_trajectory.py
@@ -545,7 +545,7 @@ class TestAISTrajectory(unittest.TestCase):
         )
 
         result = trajectory.generate_array_from_positions(height=9, width=18, link=True, bounding_box='fit',
-                                                          features=['sog', 'cog'], node_size=0).reshape((9, 18))
+                                                          features=['sog', 'cog'], node_size=0)
         expected = np.array([[[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,0], [0,0], [0.5,0.25]],
                              [[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,0], [0,0], [0.5,0.25]],
                              [[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,0], [0,0], [0.5,0.25]],