From 6820c3c8753ab5a1db72535a2ffba0e923ebc0ff Mon Sep 17 00:00:00 2001 From: Marcin <106817115+MathewNWSH@users.noreply.github.com> Date: Tue, 11 Mar 2025 11:10:53 +0100 Subject: [PATCH 1/6] https://github.com/gjoseph92/stackstac/issues/181#issuecomment-1374026051 --- stackstac/rio_reader.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/stackstac/rio_reader.py b/stackstac/rio_reader.py index 13b31fb..a58bcba 100644 --- a/stackstac/rio_reader.py +++ b/stackstac/rio_reader.py @@ -7,6 +7,7 @@ import numpy as np import rasterio as rio +from rasterio.transform import from_gcps from rasterio.vrt import WarpedVRT from .rio_env import LayeredEnv @@ -343,10 +344,10 @@ def _open(self) -> ThreadsafeRioDataset: # 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(), - "transform": ds.transform, - "height": ds.height, - "width": ds.width, + "crs": ds.gcps[1], + "transform": from_gcps(ds.gcps[0]), + # "height": ds.height, + # "width": ds.width, }: with self.gdal_env.open_vrt: vrt = WarpedVRT( From 5d9af583306fa6dd0edebefb40d813ba9f671630 Mon Sep 17 00:00:00 2001 From: Marcin <106817115+MathewNWSH@users.noreply.github.com> Date: Tue, 11 Mar 2025 12:29:57 +0100 Subject: [PATCH 2/6] Update rio_reader.py --- stackstac/rio_reader.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/stackstac/rio_reader.py b/stackstac/rio_reader.py index a58bcba..35d173f 100644 --- a/stackstac/rio_reader.py +++ b/stackstac/rio_reader.py @@ -7,7 +7,7 @@ import numpy as np import rasterio as rio -from rasterio.transform import from_gcps +from rasterio.transform import transform, from_gcps from rasterio.vrt import WarpedVRT from .rio_env import LayeredEnv @@ -345,7 +345,10 @@ def _open(self) -> ThreadsafeRioDataset: # Only make a VRT if the dataset doesn't match the spatial spec we want if self.spec.vrt_params != { "crs": ds.gcps[1], - "transform": from_gcps(ds.gcps[0]), + "transform": transform.from_gcps(ds.gcps[0]), + "nodata": ds.nodata, + "add_alpha": False, + "src_nodata": ds.nodata, # "height": ds.height, # "width": ds.width, }: From 3cfe7479775d3317950238fa63a4d68484d3c9f2 Mon Sep 17 00:00:00 2001 From: Marcin <106817115+MathewNWSH@users.noreply.github.com> Date: Tue, 11 Mar 2025 12:32:09 +0100 Subject: [PATCH 3/6] Update rio_reader.py --- stackstac/rio_reader.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/stackstac/rio_reader.py b/stackstac/rio_reader.py index 35d173f..38abdf0 100644 --- a/stackstac/rio_reader.py +++ b/stackstac/rio_reader.py @@ -7,7 +7,8 @@ import numpy as np import rasterio as rio -from rasterio.transform import transform, from_gcps +from rasterio import transform +from rasterio.transform import from_gcps from rasterio.vrt import WarpedVRT from .rio_env import LayeredEnv From 3ebefba09accd0b1184f65963a8cc6f1477433ad Mon Sep 17 00:00:00 2001 From: Marcin <106817115+MathewNWSH@users.noreply.github.com> Date: Tue, 11 Mar 2025 13:52:01 +0100 Subject: [PATCH 4/6] Update raster_spec.py --- stackstac/raster_spec.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/stackstac/raster_spec.py b/stackstac/raster_spec.py index 2a0f190..83c6493 100644 --- a/stackstac/raster_spec.py +++ b/stackstac/raster_spec.py @@ -54,10 +54,10 @@ def shape(self) -> Tuple[int, int]: @cached_property def vrt_params(self) -> dict: - height, width = self.shape return { - "crs": self.epsg, - "transform": self.transform, - "height": height, - "width": width, + "crs": self.gcps[1], + "transform": transform.from_gcps(self.gcps[0]), + "nodata": self.nodata, + "add_alpha": False, + "src_nodata": self.nodata, } From 314248d10d3594cbe70e40f34c02f908d8f441da Mon Sep 17 00:00:00 2001 From: Marcin <106817115+MathewNWSH@users.noreply.github.com> Date: Tue, 11 Mar 2025 13:52:28 +0100 Subject: [PATCH 5/6] Update rio_reader.py --- stackstac/rio_reader.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/stackstac/rio_reader.py b/stackstac/rio_reader.py index 38abdf0..83ec6e1 100644 --- a/stackstac/rio_reader.py +++ b/stackstac/rio_reader.py @@ -345,20 +345,12 @@ def _open(self) -> ThreadsafeRioDataset: # Only make a VRT if the dataset doesn't match the spatial spec we want if self.spec.vrt_params != { - "crs": ds.gcps[1], - "transform": transform.from_gcps(ds.gcps[0]), - "nodata": ds.nodata, - "add_alpha": False, - "src_nodata": ds.nodata, - # "height": ds.height, - # "width": ds.width, }: with self.gdal_env.open_vrt: vrt = WarpedVRT( ds, sharing=False, resampling=self.resampling, - add_alpha=ds.nodata is None, **self.spec.vrt_params, ) else: From 05d148e5cfe845e2ea7d6195eb3defea99675e1f Mon Sep 17 00:00:00 2001 From: Marcin Niemyjski Date: Tue, 11 Mar 2025 16:09:18 +0100 Subject: [PATCH 6/6] handle both S1 and S2 --- stackstac/raster_spec.py | 10 +++++----- stackstac/rio_reader.py | 31 +++++++++++++++++++++++++++---- 2 files changed, 32 insertions(+), 9 deletions(-) diff --git a/stackstac/raster_spec.py b/stackstac/raster_spec.py index 83c6493..2a0f190 100644 --- a/stackstac/raster_spec.py +++ b/stackstac/raster_spec.py @@ -54,10 +54,10 @@ def shape(self) -> Tuple[int, int]: @cached_property def vrt_params(self) -> dict: + height, width = self.shape return { - "crs": self.gcps[1], - "transform": transform.from_gcps(self.gcps[0]), - "nodata": self.nodata, - "add_alpha": False, - "src_nodata": self.nodata, + "crs": self.epsg, + "transform": self.transform, + "height": height, + "width": width, } diff --git a/stackstac/rio_reader.py b/stackstac/rio_reader.py index 83ec6e1..e69482f 100644 --- a/stackstac/rio_reader.py +++ b/stackstac/rio_reader.py @@ -344,18 +344,41 @@ def _open(self) -> ThreadsafeRioDataset: ) # Only make a VRT if the dataset doesn't match the spatial spec we want - if self.spec.vrt_params != { - }: + # Determine current spatial parameters based on GCP presence + if ds.gcps[0]: # Check if there are GCPs + current_crs = ds.gcps[1].to_epsg() + current_transform = transform.from_gcps(ds.gcps[0]) + else: + current_crs = ds.crs.to_epsg() + current_transform = ds.transform + + current_params = { + "crs": current_crs, + "transform": current_transform, + "height": ds.height, + "width": ds.width, + } + + if self.spec.vrt_params != current_params: + # Prepare VRT parameters + vrt_kwargs = dict(self.spec.vrt_params) + if ds.gcps[0]: + # Add GCP-based source parameters + vrt_kwargs.update({ + "src_crs": ds.gcps[1], + "src_transform": current_transform, + "src_nodata": ds.nodata, + }) with self.gdal_env.open_vrt: vrt = WarpedVRT( ds, sharing=False, resampling=self.resampling, - **self.spec.vrt_params, + add_alpha=ds.nodata is None, + **vrt_kwargs, ) else: vrt = None - if ds.driver in MULTITHREADED_DRIVER_ALLOWLIST: return ThreadLocalRioDataset(self.gdal_env, ds, vrt=vrt) # ^ NOTE: this forces all threads to wait for the `open()` we just did before they can open their