Enable opening datasets with gcps even if no valid crs is present #182

weiji14 · 2022-09-23T18:09:51Z

Some datasets like Sentinel-1 GRD files don't have a coordinate reference system (crs), but they may have ground control points (gcps) and can still be opened by rasterio.

Test this branch using pip install git+https://github.com/weiji14/stackstac.git@open_vrt_with_gcps

Fixes #181.

Some datasets like Sentinel-1 GRD files don't have a coordinate reference system (crs), but they may have ground control points (gcps) and can still be opened by rasterio.

weiji14 · 2022-09-23T18:20:25Z

Probably should write a regression test. @scottyhq, any good sample Sentinel-1 GRD files with a stable link that you know of?

Need to reproject from a proper source coordinate reference system (src_crs) instead of arbitrarily.

scottyhq · 2022-09-23T19:05:21Z

Probably should write a regression test. @scottyhq, any good sample Sentinel-1 GRD files with a stable link that you know of?

GRDs are really big. xarray-sentinel has a minimized one here for testing https://github.com/bopen/xarray-sentinel/tree/main/tests/data/S1B_IW_GRDH_1SDV_20210401T052623_20210401T052648_026269_032297_ECC8.SAFE.

Or rioxarray just constructs a synthetic dataset with GCPs for a test: https://github.com/corteva/rioxarray/blob/bfd616594bb82886133b5f4b9356ed16d39014fc/test/integration/test_integration_rioxarray.py#L1053

Not very well written, but it's a start I guess. Synthetic GeoTIFF generation code modified from https://github.com/corteva/rioxarray/blob/bfd616594bb82886133b5f4b9356ed16d39014fc/test/integration/test_integration_rioxarray.py#L1053.

gjoseph92

Thanks a bunch for working on this @weiji14, and thanks for the test too!

stackstac/rio_reader.py

Put the src_crs in the vrt_params dict to avoid an elif statement.

weiji14 · 2022-09-23T22:36:03Z

stackstac/tests/test_rio_reader.py

+    np.testing.assert_allclose(
+        actual=array, desired=np.array([[0.0, 0.0], [0.0, 0.0]], dtype=np.float32)
+    )


Thanks a bunch for working on this @weiji14, and thanks for the test too!

Don't thank me yet, this is a really bad test 😅 It only checks that things can run (and that a GeoTIFF with gcps works without an error message), but I don't like this assert statement... Any pointers?

Hm, maybe we could read it back using rasterio with (what we think are) the right parameters, and check that the results match?

Only other thing I can think of would be to write some data into the file (probably like np.linspace(...).reshape(...)). Then we could make a more meaningful assertion about the values we get out.

gjoseph92 · 2022-10-24T17:10:49Z

@weiji14 did you get any further with this? Is this actually working as is now?

weiji14 · 2022-10-24T17:23:23Z

@weiji14 did you get any further with this? Is this actually working as is now?

Technically, I think the implementation is sound. But the unit test could be much improved. Ok if you want to merge this first and improve the unit test at a later date.

It's not a good idea in dask to mutate inputs. This object could be shared with other tasks. Also, I think we just need to check `ds.crs is None`—if `ds.crs` is a CRS, it will always have a `to_epsg` method. (Though that could fail, so I've added handling for that case too.)

gjoseph92

@weiji14 I pushed a change here to avoid mutating spec.vrt_params, since that object could possibly be shared with other dask tasks.

To be clear, in your example in #181, you're still getting the purple/green/yellow output? That makes me think this might not be working correctly yet?

gjoseph92 · 2022-10-24T19:07:33Z

stackstac/tests/test_rio_reader.py

+    np.testing.assert_allclose(
+        actual=array, desired=np.array([[0.0, 0.0], [0.0, 0.0]], dtype=np.float32)
+    )


Hm, maybe we could read it back using rasterio with (what we think are) the right parameters, and check that the results match?

Only other thing I can think of would be to write some data into the file (probably like np.linspace(...).reshape(...)). Then we could make a more meaningful assertion about the values we get out.

The test is still failing for me though.

weiji14 · 2022-10-24T20:29:10Z

@weiji14 I pushed a change here to avoid mutating spec.vrt_params, since that object could possibly be shared with other dask tasks.

I'm not seeing any new commits in this PR. Could you try again?

To be clear, in your example in #181, you're still getting the purple/green/yellow output? That makes me think this might not be working correctly yet?

If I recall correctly, there's actually 2 things needed to fix the purple/green/yellow output.

Getting the image plotted at (roughly) the correct location, and this requires reading the coordinate reference system (CRS) and ground control points (GCP) which this Pull Request has fixed (at least I think it's implemented correctly).
Resampling/Interpolation - this is something this PR won't be handling, but which the user would still need to set to get a proper Sentinel-1 GRD image plotted (that is not purple/green/yellow).

I'm having trouble with the second part (resampling), which would require someone a bit smarter on the SAR data structure to get right (e.g. @scottyhq). I wish I had a proper end-to-end validated example where someone has taken a Sentinel-1 GRD image (from Microsoft Planetary Computer), read it properly with rasterio (with some set of parameters) and plotted it, but I'm missing that picture right now.

this also gets the test passing

gjoseph92 · 2022-10-24T21:06:00Z

Oops maybe I forgot to push. I see them now in https://github.com/gjoseph92/stackstac/pull/182/files/1b3c47c6bd0d3ec9dd17b25dc8cb1c8496a131a7..d7b08203ab754295bf5644fa0a661df9fa94e8b9.

I'd be a little surprised to find that it's just a resampling issue. As long as the data types are right and we're not truncating floats to ints or something, I don't quite see why a different resampling method would change meaningful data into something completely meaningless-looking.

Furthermore, if you look at the example scene from #181 in PC explorer, you can see that what we're seeing from stackstac.show also has a different spatial extent, as well as looking meaningless. I suppose that could be a STAC metadata issue like #152 though (I haven't looked at the metadata at all).

I guess where I'm leaning right now is that if we can't get this to seem to work properly and return meaningful data, I'd rather raise a NotImplementedError if ds.crs is None than try to use GCPs but possibly return bad values.

weiji14 · 2023-01-06T19:10:58Z

stackstac/rio_reader.py

-            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,


MIght need to set src_transform in addition to transform as mentioned by @vincentsarago in #181 (comment)?

So looking at the code it might be a bit more complex.

What I do in rio-tiler, is automatically create a WarpedVRT when we have gcps (which in stackstac case should replace the ds) and then I create another WarpedVRT on top of it if I need warping

Enable opening datasets with gcps even if no valid crs is present

66f3781

Some datasets like Sentinel-1 GRD files don't have a coordinate reference system (crs), but they may have ground control points (gcps) and can still be opened by rasterio.

weiji14 marked this pull request as ready for review September 23, 2022 18:17

gjoseph92 mentioned this pull request Sep 23, 2022

AttributeError: 'NoneType' object has no attribute 'to_epsg' #181

Open

weiji14 marked this pull request as draft September 23, 2022 18:29

Get src_crs from last index position of gcps tuple

4f656fe

Need to reproject from a proper source coordinate reference system (src_crs) instead of arbitrarily.

weiji14 added 2 commits September 23, 2022 15:44

Add an inline comment about trying to make VRT when gcps are present

cc0cd1f

gjoseph92 reviewed Sep 23, 2022

View reviewed changes

stackstac/rio_reader.py Show resolved Hide resolved

Remove duplicated WarpedVRT instantiation

1b3c47c

Put the src_crs in the vrt_params dict to avoid an elif statement.

weiji14 commented Sep 23, 2022

View reviewed changes

gjoseph92 reviewed Oct 24, 2022

View reviewed changes

fix up test imports

fca7961

The test is still failing for me though.

fix conditional

d7b0820

this also gets the test passing

weiji14 commented Jan 6, 2023

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Enable opening datasets with gcps even if no valid crs is present #182

Enable opening datasets with gcps even if no valid crs is present #182

weiji14 commented Sep 23, 2022 •

edited

weiji14 commented Sep 23, 2022

scottyhq commented Sep 23, 2022

gjoseph92 left a comment

weiji14 Sep 23, 2022 •

edited

gjoseph92 Oct 24, 2022

gjoseph92 commented Oct 24, 2022

weiji14 commented Oct 24, 2022 •

edited

gjoseph92 left a comment

gjoseph92 Oct 24, 2022

weiji14 commented Oct 24, 2022

gjoseph92 commented Oct 24, 2022

weiji14 Jan 6, 2023

vincentsarago Jan 13, 2023

Enable opening datasets with gcps even if no valid crs is present #182

Are you sure you want to change the base?

Enable opening datasets with gcps even if no valid crs is present #182

Conversation

weiji14 commented Sep 23, 2022 • edited

weiji14 commented Sep 23, 2022

scottyhq commented Sep 23, 2022

gjoseph92 left a comment

Choose a reason for hiding this comment

weiji14 Sep 23, 2022 • edited

Choose a reason for hiding this comment

gjoseph92 Oct 24, 2022

Choose a reason for hiding this comment

gjoseph92 commented Oct 24, 2022

weiji14 commented Oct 24, 2022 • edited

gjoseph92 left a comment

Choose a reason for hiding this comment

gjoseph92 Oct 24, 2022

Choose a reason for hiding this comment

weiji14 commented Oct 24, 2022

gjoseph92 commented Oct 24, 2022

weiji14 Jan 6, 2023

Choose a reason for hiding this comment

vincentsarago Jan 13, 2023

Choose a reason for hiding this comment

weiji14 commented Sep 23, 2022 •

edited

weiji14 Sep 23, 2022 •

edited

weiji14 commented Oct 24, 2022 •

edited