Get src_crs from last index position of gcps tuple

Need to reproject from a proper source coordinate reference system (src_crs) instead of arbitrarily.
gjoseph92 · Sep 23, 2022 · 4f656fe · 4f656fe
1 parent 66f3781
commit 4f656fe
Showing 1 changed file with 16 additions and 10 deletions.
diff --git a/stackstac/rio_reader.py b/stackstac/rio_reader.py
@@ -344,23 +344,29 @@ def _open(self) -> ThreadsafeRioDataset:
                 )
 
             # Only make a VRT if the dataset doesn't match the spatial spec we want
-            if (
-                hasattr(ds.crs, "to_epsg")
-                and self.spec.vrt_params
-                != {
-                    "crs": ds.crs.to_epsg(),
-                    "transform": ds.transform,
-                    "height": ds.height,
-                    "width": ds.width,
-                }
-            ) or ds.gcps is not None:
+            if hasattr(ds.crs, "to_epsg") and self.spec.vrt_params != {
+                "crs": ds.crs.to_epsg(),
+                "transform": ds.transform,
+                "height": ds.height,
+                "width": ds.width,
+            }:
                 with self.gdal_env.open_vrt:
                     vrt = WarpedVRT(
                         ds,
                         sharing=False,
                         resampling=self.resampling,
                         **self.spec.vrt_params,
                     )
+            elif ds.gcps is not None:
+                with self.gdal_env.open_vrt:
+                    src_crs = ds.gcps[-1]
+                    vrt = WarpedVRT(
+                        ds,
+                        src_crs=src_crs,
+                        sharing=False,
+                        resampling=self.resampling,
+                        **self.spec.vrt_params,
+                    )
             else:
                 logger.info(f"Skipping VRT for {self.url!r}")
                 vrt = None