diff --git a/openeo_driver/save_result.py b/openeo_driver/save_result.py
index cfab273b..bd5a1719 100644
--- a/openeo_driver/save_result.py
+++ b/openeo_driver/save_result.py
@@ -314,8 +314,7 @@ def __init__(
         """
         super().__init__(data=timeseries)
         if not isinstance(regions, (GeometryCollection, DriverVectorCube)):
-            # TODO: raise exception instead of warning?
-            _log.warning(
+            raise ValueError(
                 f"AggregatePolygonResult: GeometryCollection or DriverVectorCube expected but got {type(regions)}"
             )
         self._regions = regions
@@ -445,17 +444,21 @@ def _create_point_timeseries_xarray(self, feature_ids, timestamps, lats, lons, a
         return the_array
 
     def to_netcdf(self, destination: Optional[str] = None) -> str:
-        def features_ids_from_index(geometries):
-            return ["feature_%d" % i for i in range(len(geometries.geoms))]
-
         if isinstance(self._regions, GeometryCollection):
             points = [r.representative_point() for r in self._regions.geoms]
-            feature_ids = features_ids_from_index(self._regions)
-        else:
+            feature_ids = ["feature_%d" % i for i in range(len(self._regions.geoms))]
+        elif isinstance(self._regions, DriverVectorCube):
             points = [r.representative_point() for r in self._regions.get_geometries()]
             feature_ids = self._regions.get_ids()
-            feature_ids = (list(feature_ids) if feature_ids is not None
-                           else features_ids_from_index(self._regions.get_geometries()))
+            feature_ids = (
+                list(feature_ids)
+                if feature_ids is not None
+                else ["feature_%d" % i for i in range(self._regions.geometry_count())]
+            )
+        else:
+            raise ValueError(
+                f"AggregatePolygonResult: GeometryCollection or DriverVectorCube expected but got {type(self._regions)}"
+            )
 
         lats = [p.y for p in points]
         lons = [p.x for p in points]
diff --git a/tests/test_save_result_netcdf.py b/tests/test_save_result_netcdf.py
index 9c19847f..8bc7a73b 100644
--- a/tests/test_save_result_netcdf.py
+++ b/tests/test_save_result_netcdf.py
@@ -1,9 +1,11 @@
 from pathlib import Path
 
+import geopandas as gpd
 import numpy as np
 import pytest
 import xarray as xr
 from numpy.testing import assert_array_equal, assert_allclose
+from openeo_driver.datacube import DriverVectorCube
 from openeo_driver.save_result import AggregatePolygonResult, AggregatePolygonResultCSV
 from shapely.geometry import GeometryCollection, Polygon
 
@@ -150,3 +152,27 @@ def test_aggregate_polygon_result_CSV(tmp_path):
     timeseries_ds.red.sel( t='2017-09-05')
 
     assert_allclose( timeseries_ds.red.sel(feature=2).sel( t='2017-09-06').data,4645.719597)
+
+
+def test_aggregate_polygon_result_vector_cube_no_ids_to_netcdf(tmp_path):
+    """Test for https://github.com/Open-EO/openeo-python-driver/issues/519"""
+    timeseries = {
+        "2020-06-01T00:00:00Z": [[77.2], [77.15]],
+        "2020-06-02T00:00:00Z": [[77.21], [77.16]],
+    }
+    # GeoDataFrame without an 'id' column so that get_ids() returns None
+    gdf = gpd.GeoDataFrame(
+        geometry=[
+            Polygon([(77.11, 28.56), (77.29, 28.56), (77.29, 28.69), (77.11, 28.69)]),
+            Polygon([(77.11, 28.56), (77.20, 28.56), (77.20, 28.62), (77.11, 28.62)]),
+        ]
+    )
+    regions = DriverVectorCube(gdf)
+
+    result = AggregatePolygonResult(timeseries, regions=regions)
+    filename = result.to_netcdf(tmp_path / "timeseries_vector_cube_no_ids.nc")
+
+    timeseries_ds = xr.open_dataset(filename)
+    # Feature IDs should be auto-generated as "feature_0", "feature_1", ...
+    assert list(timeseries_ds.coords["feature_names"].data) == ["feature_0", "feature_1"]
+    assert_array_equal(timeseries_ds.band_0.sel(feature=0).data, [77.2, 77.21])