Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Using DigitalEarthAfrica data #132

Closed
vitusbenson opened this issue Feb 11, 2022 · 6 comments
Closed

Using DigitalEarthAfrica data #132

vitusbenson opened this issue Feb 11, 2022 · 6 comments

Comments

@vitusbenson
Copy link

vitusbenson commented Feb 11, 2022

Hi @gjoseph92

(solved some of my issues myself, thats why I edit again)

I would like to use Sentinel 1 RTC from Digital Earth Africa on AWS.

I try the following:

import stackstac
import pystac_client
import rasterio

URL = "https://explorer.digitalearth.africa/stac/"
catalog = pystac_client.Client.open(URL)

bbox = (-3, 15.25, -2.9, 15.35)

items = catalog.search(
    bbox = bbox,
    collections=["s1_rtc"],
    datetime="2021-11-01/2021-12-01"
).get_all_items()

with rasterio.Env(aws_unsigned = True,AWS_S3_ENDPOINT= 's3.af-south-1.amazonaws.com'):
     stack = stackstac.stack(items)
     stack.compute()

When I do this I get the error

[NotImplementedError: Cannot automatically compute the resolution, since asset 'vh' on item 0 '620fe0bb-2e29-5ac3-86c1-7bcc9698da29' has a non-rectilinear geotrans (its data is is not axis-aligned, or "north-up"): [-4.0, 0.0002, 0.0, 16.0, 0.0, -0.0002, 0.0, 0.0, 1.0]. We should be able to handle this but just don't want to deal with it right now.

Please specify the `resolution=` argument.]()

However when opening such Sentinel 1 tif manually:

test = xr.open_rasterio("https://deafrica-sentinel-1.s3.af-south-1.amazonaws.com/s1_rtc/N14W009/2018/01/06/022245/s1_rtc_022245_N14W009_2018_01_06_VH.tif")

I get test.transform == (0.0002, 0.0, -9.0, 0.0, -0.0002, 15.0), which to me seems like a legit Geotransform..

Do you have any idea?

Thanks(:

@Kirill888
Copy link

@vitusbenson
Transform in the file is fine, but transform as recorded in the STAC document uses GDAL geotransform order, which is not correct proj extension order which instead expects affine matrix order. So I'm guessing that's why you see this issue.

https://explorer.digitalearth.africa/stac/collections/s1_rtc/items/46bbaf19-62fd-58c0-9538-9a48c2aec7cc

Digital Earth Africa team has confirmed an issue with metadata, problem is upstream from them but hopefully will be addressed some time soon.

@alexgleith
Copy link

Thanks for reporting @vitusbenson, I'm working on resolving this by updating the DE Africa metadata documents.

@gjoseph92
Copy link
Owner

Thanks @Kirill888, was literally writing out the same thing. You beat me by a minute!

Yeah, the proj extension spec could be a little clearer here that the transform field should not be in GDAL order (cc @matthewhanson). It's clear from the matrix-math listed, but the fact that "GDAL(GetGeoTransform)" is referenced is a little misleading.

So @vitusbenson, just to be clear: this isn't a bug with stackstac, it's just the that the STAC metadata is wrong. Because stackstac is careful to never open the underlying GeoTIFFs directly (docs needed to explain this clearly, xref #112), it doesn't matter that xr.open_rasterio gives you the correct geotrans, because stackstac only gets to look at the metadata.

@Kirill888
Copy link

@gjoseph92 agree about the wording in the stac extension spec: stac-extensions/projection#9

@vitusbenson
Copy link
Author

So @vitusbenson, just to be clear: this isn't a bug with stackstac, it's just the that the STAC metadata is wrong. Because stackstac is careful to never open the underlying GeoTIFFs directly (docs needed to explain this clearly, xref #112), it doesn't matter that xr.open_rasterio gives you the correct geotrans, because stackstac only gets to look at the metadata.

Thanks guys, I got a dirty bugfix working by manually changing the geotransform from GDAL to Proj format in the pystac catalog that I then give to stacstack.

@alexgleith
Copy link

I believe this is fixed in our Sentinel-1 data, here's an updated scene: https://explorer.digitalearth.africa/stac/collections/s1_rtc/items/0a461211-1d56-5bb0-8baf-bb9917d8f77f

Transform is now:

"proj:transform": [
0.0002,
0,
50,
0,
-0.0002,
12,
0,
0,
1
],

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants