diff --git a/pyflwdir/gis_utils.py b/pyflwdir/gis_utils.py
index e669648..4effaaf 100644
--- a/pyflwdir/gis_utils.py
+++ b/pyflwdir/gis_utils.py
@@ -487,7 +487,7 @@ def distance(idx0, idx1, ncol, latlon=False, transform=IDENTITY):
 
 
 ## VECTORIZE
-def features(flowpaths, xs, ys, **kwargs):
+def features(flowpaths, xs=None, ys=None, transform=None, shape=None, **kwargs):
     """Returns a LineString feature for each stream
 
     Parameters
@@ -496,8 +496,10 @@ def features(flowpaths, xs, ys, **kwargs):
         linear indices of flowpaths
     xs, ys : 1D-array of float
         x, y coordinates
-    idxs_ds: list of intp
-        linear index of first cell on next downstream flow path
+    transform : Affine
+        Coefficients mapping pixel coordinates to coordinate reference system.
+    shape : tuple of int
+        The height, width  of the raster.
     kwargs : extra sample maps key-word arguments
         optional maps to sample from
         e.g.: strord=flw.stream_order()
@@ -507,8 +509,17 @@ def features(flowpaths, xs, ys, **kwargs):
     feats : list of dict
         Geofeatures, to be parsed by e.g. geopandas.GeoDataFrame.from_features
     """
+    if xs is None or ys is None:
+        if transform is None or shape is None:
+            raise ValueError(
+                "transform and shape should be provided if xs and ys are None"
+            )
+        _size = shape[0] * shape[1]
+    else:
+        _size = xs.size
+
     for key in kwargs:
-        if not isinstance(kwargs[key], np.ndarray) or kwargs[key].size != xs.size:
+        if not isinstance(kwargs[key], np.ndarray) or kwargs[key].size != _size:
             raise ValueError(
                 f'Kwargs map "{key}" should be ndarrays of same size as coordinates'
             )
@@ -520,12 +531,17 @@ def features(flowpaths, xs, ys, **kwargs):
         idx0 = idxs[0]
         pit = idxs[-1] == idxs[-2]
         props = {key: kwargs[key].flat[idx0] for key in kwargs}
+        if xs is None or ys is None:
+            xi, yi = idxs_to_coords(idxs, transform, shape)
+            coordinates = list(zip(xi, yi))
+        else:
+            coordinates = [(xs[i], ys[i]) for i in idxs]
         feats.append(
             {
                 "type": "Feature",
                 "geometry": {
                     "type": "LineString",
-                    "coordinates": [(xs[i], ys[i]) for i in idxs],
+                    "coordinates": coordinates,
                 },
                 "properties": {"idx": idx0, "idx_ds": idxs[-1], "pit": pit, **props},
             }
diff --git a/pyflwdir/pyflwdir.py b/pyflwdir/pyflwdir.py
index 12cd6eb..c19eef4 100644
--- a/pyflwdir/pyflwdir.py
+++ b/pyflwdir/pyflwdir.py
@@ -974,13 +974,12 @@ def geofeatures(self, flowpaths, xs=None, ys=None, **kwargs):
             Geofeatures, to be parsed by e.g. geopandas.GeoDataFrame.from_features
         """
         # get geoms and return features
-        if xs is None or ys is None:
-            idxs0 = np.arange(self.size, dtype=np.intp)
-            xs, ys = gis.idxs_to_coords(idxs0, self.transform, self.shape)
         feats = gis.features(
             flowpaths=flowpaths,
-            xs=self._check_data(xs, "xs"),
-            ys=self._check_data(ys, "ys"),
+            xs=self._check_data(xs, "xs", optional=True),
+            ys=self._check_data(ys, "ys", optional=True),
+            transform=self.transform,
+            shape=self.shape,
             **kwargs,
         )
         return feats