Datacube Indexing RuntimeError: mismatch of attribute 'resolution'

[suhendra0812]

Hi, all.

I encountered an issue when executing the datacube indexing process. It raised RuntimeError: mismatch of attribute 'resolution'. I used the dataset which is processed using geocode function with EPSG:4326. Here the detail of the geocode script:

from pathlib import Path
import time

from pyroSAR.snap.util import geocode

paths = sorted(Path("sentinel/SAR/GRD").glob("*/*/*/*.zip"))

for path in paths:
    print(path)
    infile = str(path)
    outdir = str(path.parent)

    start_time = time.time()

    result = geocode(
        infile,
        outdir,
        t_srs=4326,
        tr=10,
        nodataValueAtSea=False,
        test=False
    )

    end_time = time.time()
    elapse_time = end_time - start_time
    print(f"Elapse time: {elapse_time/60} minutes")

After finished the geocode process, I try to execute datacube indexing process with the following script:

from pathlib import Path

from pyroSAR.datacube_util import Product, Dataset
from pyroSAR.ancillary import find_datasets, groupby

# find pyroSAR files by metadata attributes
archive_dir = Path("sentinel")
scenes = find_datasets(str(archive_dir), recursive=True)

# group the found files by their file basenames
# files with the same basename are considered to belong to the same dataset
grouped = groupby(scenes, "outname_base")

# define the polarization units describing the data sets
units = {"VV": "backscatter VV", "VH": "backscatter VH"}

# create a new product
product_name = "s1_gamma0_rtc"
product_type = "gamma0"
with Product(
    name=product_name,
    product_type=product_type,
    description="Gamma Naught RTC backscatter",
) as prod:
    for dataset in grouped:
        with Dataset(dataset, units=units) as ds:

            # add the dataset to the product
            prod.add(ds)

            # parse datacube indexing YMLs from product and data set metadata
            yml_index_oudir = Path(dataset[0]).parent
            yml_index = sorted(yml_index_oudir.glob("*.yml"))
            if len(yml_index) == 0:
                prod.export_indexing_yml(ds, str(yml_index_oudir))

    # write the product YML
    yml_product = archive_dir.joinpath(f"{product_name}.yml")
    prod.write(str(yml_product))

    # print the product metadata which is written to the product YML
    print(prod)

But I encountered this error:

Traceback (most recent call last):
  File "E:\Data\eodatasets\s1_indexing_pyrosar.py", line 29, in <module>
    prod.add(ds)
  File "C:\Users\suhen\miniconda3\envs\odc\lib\site-packages\pyroSAR\datacube_util.py", line 552, in add
    self.check_integrity(dataset, allow_new_measurements=True)
  File "C:\Users\suhen\miniconda3\envs\odc\lib\site-packages\pyroSAR\datacube_util.py", line 592, in check_integrity
    raise RuntimeError("mismatch of attribute '{0}': {1}, {2}".format(attr, val_ds, val_prod))
RuntimeError: mismatch of attribute 'resolution': {'longitude': 8.983152841195213e-05, 'latitude': 8.983152841195213e-05}, {'longitude': 8.983152841195205e-05, 'latitude': 8.983152841195205e-05}

Is there any mistake in my script or anything else?

These are datasets that I use:

>>> for dataset in grouped:
...     print(dataset)
... 
['E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2017\\07\\20\\S1A__IW___A_20170720T104110_VH_NR_Orb_Cal_TF_TC_dB.tif', 'E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2017\\07\\20\\S1A__IW___A_20170720T104110_VV_NR_Orb_Cal_TF_TC_dB.tif']
['E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2018\\08\\20\\S1A__IW___A_20180820T104118_VH_NR_Orb_Cal_TF_TC_dB.tif', 'E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2018\\08\\20\\S1A__IW___A_20180820T104118_VV_NR_Orb_Cal_TF_TC_dB.tif']
['E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2019\\07\\10\\S1A__IW___A_20190710T104122_VH_NR_Orb_Cal_TF_TC_dB.tif', 'E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2019\\07\\10\\S1A__IW___A_20190710T104122_VV_NR_Orb_Cal_TF_TC_dB.tif']
['E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2020\\02\\23\\S1A__IW___A_20200223T104123_VH_NR_Asm_Orb_Cal_TF_TC_dB.tif', 'E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2020\\02\\23\\S1A__IW___A_20200223T104123_VV_NR_Asm_Orb_Cal_TF_TC_dB.tif']
['E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2020\\07\\28\\S1A__IW___A_20200728T104130_VH_NR_Orb_Cal_TF_TC_dB.tif', 'E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2020\\07\\28\\S1A__IW___A_20200728T104130_VV_NR_Orb_Cal_TF_TC_dB.tif']
['E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2015\\03\\24\\S1A__IW___D_20150324T215246_VH_NR_Orb_Cal_TF_TC_dB.tif', 'E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2015\\03\\24\\S1A__IW___D_20150324T215246_VV_NR_Orb_Cal_TF_TC_dB.tif']
['E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2015\\03\\29\\S1A__IW___D_20150329T220041_VH_NR_Orb_Cal_TF_TC_dB.tif', 'E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2015\\03\\29\\S1A__IW___D_20150329T220041_VV_NR_Orb_Cal_TF_TC_dB.tif']
['E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2015\\07\\22\\S1A__IW___D_20150722T215252_VH_NR_Orb_Cal_TF_TC_dB.tif', 'E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2015\\07\\22\\S1A__IW___D_20150722T215252_VV_NR_Orb_Cal_TF_TC_dB.tif']
['E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2016\\08\\09\\S1A__IW___D_20160809T215257_VH_NR_Orb_Cal_TF_TC_dB.tif', 'E:\\Data\\eodatasets\\sentinel\\SAR\\GRD\\2016\\08\\09\\S1A__IW___D_20160809T215257_VV_NR_Orb_Cal_TF_TC_dB.tif']

Please assist me if there is something wrong with my script.
Thank you in advance.

Best Regards,
Suhendra

[johntruckenbrodt]

Hi @suhendra0812. Just as a word of warning... pyroSAR's Open Data Cube functionality has not seen an update since 2019 so it is very likely that it does not longer work as it is.
However, your error seems to be coming from the fact that I always worked with data in UTM projection (and thus resolution in meters) and you work with resolution in degrees. There is a slight mismatch between resolutions:

{'longitude': 8.983152841195213e-05, 'latitude': 8.983152841195213e-05}
{'longitude': 8.983152841195205e-05, 'latitude': 8.983152841195205e-05}

This test in check_integrity needs to be adjusted to account for this. Unfortunately I currently do not have the time to do this. You are very welcome to contribute..

[suhendra0812]

Ok, I will change the CRS to UTM.
Thanks @johntruckenbrodt for the clarification.

[johntruckenbrodt]

Sure, changing to UTM would also work.