diff --git a/stackstac/rio_reader.py b/stackstac/rio_reader.py
index 226c455..8771ebf 100644
--- a/stackstac/rio_reader.py
+++ b/stackstac/rio_reader.py
@@ -7,6 +7,7 @@
 
 import numpy as np
 import rasterio as rio
+import rasterio.errors
 from rasterio.vrt import WarpedVRT
 
 from .rio_env import LayeredEnv
@@ -323,9 +324,7 @@ def _open(self) -> ThreadsafeRioDataset:
         with self.gdal_env.open:
             with time(f"Initial read for {self.url!r} on {_curthread()}: {{t}}"):
                 try:
-                    ds = SelfCleaningDatasetReader(
-                        self.url, sharing=False
-                    )
+                    ds = SelfCleaningDatasetReader(self.url, sharing=False)
                 except Exception as e:
                     msg = f"Error opening {self.url!r}: {e!r}"
                     if exception_matches(e, self.errors_as_nodata):
@@ -343,18 +342,30 @@ def _open(self) -> ThreadsafeRioDataset:
                     "a separate STAC asset), so you'll need to exclude this asset from your analysis."
                 )
 
+            ds_epsg: int | None
+            try:
+                ds_epsg = ds.crs.to_epsg() if ds.crs is not None else None
+            except rasterio.errors.CRSError:
+                ds_epsg = None
+
             # Only make a VRT if the dataset doesn't match the spatial spec we want
-            if self.spec.vrt_params != {
-                "crs": ds.crs.to_epsg(),
+            if ds_epsg is None or self.spec.vrt_params != {
+                "crs": ds_epsg,
                 "transform": ds.transform,
                 "height": ds.height,
                 "width": ds.width,
             }:
+                # Set source crs from ground control points (gcps) if present
+                src_crs = (
+                    ds.gcps[-1] if ds.crs is None and ds.gcps is not None else None
+                )
+
                 with self.gdal_env.open_vrt:
                     vrt = WarpedVRT(
                         ds,
                         sharing=False,
                         resampling=self.resampling,
+                        src_crs=src_crs,
                         **self.spec.vrt_params,
                     )
             else:
diff --git a/stackstac/tests/test_rio_reader.py b/stackstac/tests/test_rio_reader.py
new file mode 100644
index 0000000..1735515
--- /dev/null
+++ b/stackstac/tests/test_rio_reader.py
@@ -0,0 +1,53 @@
+import tempfile
+
+import numpy as np
+import rasterio
+import rasterio.enums
+from rasterio.control import GroundControlPoint
+from rasterio.crs import CRS
+from rasterio.windows import Window
+from stackstac.raster_spec import RasterSpec
+from stackstac.rio_reader import AutoParallelRioReader
+
+
+def test_dataset_read_with_gcps():
+    """
+    Ensure that GeoTIFFs with ground control points (gcps) can be read using
+    AutoParallelRioReader.
+
+    Regression test for https://github.com/gjoseph92/stackstac/issues/181.
+    """
+    with tempfile.NamedTemporaryFile(suffix=".tif") as tmpfile:
+        src_gcps = [
+            GroundControlPoint(row=0, col=0, x=156113, y=2818720, z=0),
+            GroundControlPoint(row=0, col=800, x=338353, y=2785790, z=0),
+            GroundControlPoint(row=800, col=800, x=297939, y=2618518, z=0),
+            GroundControlPoint(row=800, col=0, x=115698, y=2651448, z=0),
+        ]
+        crs = CRS.from_epsg(32618)
+        with rasterio.open(
+            tmpfile.name,
+            mode="w",
+            height=800,
+            width=800,
+            count=1,
+            dtype=np.uint8,
+            driver="GTiff",
+        ) as source:
+            source.gcps = (src_gcps, crs)
+
+        reader = AutoParallelRioReader(
+            url=tmpfile.name,
+            spec=RasterSpec(
+                epsg=4326, bounds=(90, -10, 110, 10), resolutions_xy=(10, 10)
+            ),
+            resampling=rasterio.enums.Resampling.bilinear,
+            dtype=np.dtype("float32"),
+            fill_value=np.nan,
+            rescale=True,
+        )
+        array = reader.read(window=Window(col_off=0, row_off=0, width=10, height=10))
+
+    np.testing.assert_allclose(
+        actual=array, desired=np.array([[0.0, 0.0], [0.0, 0.0]], dtype=np.float32)
+    )