In [1]:

from pystac_client import Client

from odc.stac import stac_load, configure_rio


config = {
    "s2_l2a": {
        "assets": {
            "*": {
                "data_type": "uint16",
                "nodata": 0,
                "unit": "1",
            },
            "SCL": {
                "data_type": "uint8",
                "nodata": 0,
                "unit": "1",
            },
        },
        "aliases": {
            "costal_aerosol": "B01",
            "blue": "B02",
            "green": "B03",
            "red": "B04",
            "red_edge_1": "B05",
            "red_edge_2": "B06",
            "red_edge_3": "B07",
            "nir": "B08",
            "nir_narrow": "B08A",
            "water_vapour": "B09",
            "swir_1": "B11",
            "swir_2": "B12",
            "mask": "SCL",
            "aerosol_optical_thickness": "AOT",
            "scene_average_water_vapour": "WVP",
        },
    }
}



configure_rio(
    cloud_defaults=True,
    aws={"aws_unsigned": True},
    AWS_S3_ENDPOINT="s3.af-south-1.amazonaws.com",
)


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

In [2]:

for collection in catalog.get_all_collections():
    print(collection)

<CollectionClient id=Arrivals>
<CollectionClient id=alos_palsar_mosaic>
<CollectionClient id=cci_landcover>
<CollectionClient id=cgls_landcover>
<CollectionClient id=crop_mask_eastern>
<CollectionClient id=crop_mask_northern>
<CollectionClient id=crop_mask_sahel>
<CollectionClient id=crop_mask_southern>
<CollectionClient id=crop_mask_western>
<CollectionClient id=dem_cop_30>
<CollectionClient id=dem_cop_90>
<CollectionClient id=dem_srtm>
<CollectionClient id=dem_srtm_deriv>
<CollectionClient id=esa_worldcover>
<CollectionClient id=fc_ls>
<CollectionClient id=fc_ls_summary_annual>
<CollectionClient id=gm_ls5_ls7_annual>
<CollectionClient id=gm_ls5_ls7_annual_lowres>
<CollectionClient id=gm_ls8_annual>
<CollectionClient id=gm_ls8_annual_lowres>
<CollectionClient id=gm_s2_annual>
<CollectionClient id=gm_s2_annual_lowres>
<CollectionClient id=gm_s2_semiannual>
<CollectionClient id=gm_s2_semiannual_lowres>
<CollectionClient id=gmw>
<CollectionClient id=io_lulc>
<CollectionClient id=jers_sar_mosaic>
<CollectionClient id=ls5_sr>
<CollectionClient id=ls5_st>
<CollectionClient id=ls7_sr>
<CollectionClient id=ls7_st>
<CollectionClient id=ls8_sr>
<CollectionClient id=ls8_st>
<CollectionClient id=ls9_sr>
<CollectionClient id=ls9_st>
<CollectionClient id=nasadem>
<CollectionClient id=ndvi_climatology_ls>
<CollectionClient id=pc_s2_annual>
<CollectionClient id=rainfall_chirps_daily>
<CollectionClient id=rainfall_chirps_monthly>
<CollectionClient id=s1_rtc>
<CollectionClient id=s2_l2a>
<CollectionClient id=wofs_ls>
<CollectionClient id=wofs_ls_summary_alltime>
<CollectionClient id=wofs_ls_summary_annual>

In [3]:

# %% [markdown]
# ## Find STAC Items to Load
#
# ### Define query parameters

# %%
# Set a bounding box
# [xmin, ymin, xmax, ymax] in latitude and longitude
bbox = [37.76, 12.49, 37.77, 12.50]

# Set a start and end date
start_date = "2020-09-01"
end_date = "2020-12-01"

# Set the STAC collections
collections = ["s1_rtc"]


# %% [markdown]
# ### Construct query and get items from catalog

# %%
# Build a query with the set parameters
query = catalog.search(
    bbox=bbox, collections=collections, datetime=f"{start_date}/{end_date}"
)

In [4]:

# Search the STAC catalog for all items matching the query
items = list(query.get_items())
print(f"Found: {len(items):d} datasets")

Found: 8 datasets

In [5]:

# %% [markdown]
# ## Load the Data
#
# In this step, we specify the desired coordinate system, resolution (here 20m), and bands to load. We also pass the bounding box to the `stac_load` function to only load the requested data. Since the band aliases are contained in the `config` dictionary, bands can be loaded using these aliaes (e.g. `"red"` instead of `"B04"` below).
#
# The data will be lazy-loaded with dask, meaning that is won't be loaded into memory until necessary, such as when it is displayed.

# %%
crs = "EPSG:6933"
resolution = 20

ds = stac_load(
    items,
    bands=("vh", "vv"),
    crs=crs,
    resolution=resolution,
    chunks={},
    groupby="solar_day",
    stac_cfg=config,
    bbox=bbox,
)

# View the Xarray Dataset
ds

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Input In [5], in <module>
      9 crs = "EPSG:6933"
     10 resolution = 20
---> 12 ds = stac_load(
     13     items,
     14     bands=("vh", "vv"),
     15     crs=crs,
     16     resolution=resolution,
     17     chunks={},
     18     groupby="solar_day",
     19     stac_cfg=config,
     20     bbox=bbox,
     21 )
     23 # View the Xarray Dataset
     24 ds

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/site-packages/odc/stac/_load.py:454, in load(items, bands, groupby, resampling, chunks, crs, resolution, geobox, bbox, lon, lat, x, y, align, like, geopolygon, stac_cfg, patch_url, product_cache, skip_broken_datasets, progress_cbk, fuse_func, **kw)
    450     if "resolution" not in geo:
    451         # Need to pick resolution
    452         geo["resolution"] = pick_best_resolution(dss, bands)
--> 454 dss = list(stac2ds(items, stac_cfg, product_cache=product_cache))
    455 auto_fill_geo(geo, dss, bands)
    457 if patch_url is not None:

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/site-packages/odc/stac/_eo3.py:735, in stac2ds(items, cfg, product_cache)
    733 # Have not seen this collection yet, figure it out
    734 if product is None:
--> 735     product = infer_dc_product(item, cfg)
    736     products[collection_id] = product
    738 yield item_to_ds(item, product)

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/functools.py:888, in singledispatch.<locals>.wrapper(*args, **kw)
    884 if not args:
    885     raise TypeError(f'{funcname} requires at least '
    886                     '1 positional argument')
--> 888 return dispatch(args[0].__class__)(*args, **kw)

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/site-packages/odc/stac/_eo3.py:543, in infer_dc_product_from_item(item, cfg)
    539 # We assume that grouping of data bands into grids is consistent across
    540 # entire collection, so we compute it once and keep it on a product object
    541 # at least for now.
    542 if has_proj:
--> 543     _, band2grid = compute_eo3_grids(data_bands)
    544 else:
    545     band2grid = _band2grid_from_gsd(data_bands)

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/site-packages/odc/stac/_eo3.py:275, in compute_eo3_grids(assets)
    272     gsd = geobox_gsd(geobox)
    273     return f"g{gsd:g}"
--> 275 geoboxes = dicttoolz.valmap(asset_geobox, assets)
    277 # GeoBox to list of bands that share same footprint
    278 grids: Dict[GeoBox, List[str]] = {}

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/site-packages/toolz/dicttoolz.py:83, in valmap(func, d, factory)
     72 """ Apply function to values of dictionary
     73 
     74 >>> bills = {"Alice": [20, 15, 30], "Bob": [10, 35]}
   (...)
     80     itemmap
     81 """
     82 rv = factory()
---> 83 rv.update(zip(d.keys(), map(func, d.values())))
     84 return rv

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/site-packages/odc/stac/_eo3.py:242, in asset_geobox(asset)
    239     raise ValueError(f"Asset transform is not affine: {_proj.transform}")
    241 affine = Affine(*_proj.transform[:6])
--> 242 return GeoBox(w, h, affine, _proj.crs_string)

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/site-packages/datacube/utils/geometry/_base.py:989, in GeoBox.__init__(self, width, height, affine, crs)
    988 def __init__(self, width: int, height: int, affine: Affine, crs: MaybeCRS):
--> 989     assert is_affine_st(affine), "Only axis-aligned geoboxes are currently supported"
    990     self.width = width
    991     self.height = height

AssertionError: Only axis-aligned geoboxes are currently supported

In [6]:

items_dict = query.get_all_items_as_dict()['features']

In [7]:

len(items_dict)

Out[7]:

In [8]:

items_dict[0]

Out[8]:

{'type': 'Feature',
 'stac_version': '1.0.0',
 'id': '46bbaf19-62fd-58c0-9538-9a48c2aec7cc',
 'properties': {'title': 'N12E037_2020_09_02_03F841',
  'proj:epsg': 4326,
  'platform': 'sentinel-1a',
  'odc:product': 's1_rtc',
  'instruments': ['c-sar'],
  'odc:file_format': 'GeoTIFF',
  'odc:region_code': 'N12E037',
  'sat:orbit_state': 'descending',
  'end_datetime': '2020-09-02T03:17:10.538021Z',
  'constellation': 'sentinel-1',
  'sar:product_type': 'RTC',
  'sar:polarizations': ['VV', 'VH'],
  'start_datetime': '2020-09-02T03:16:45.539433Z',
  'sar:frequency_band': 'C',
  'sat:absolute_orbit': 34174,
  'sat:relative_orbit': 152,
  'sar:instrument_mode': 'IW',
  'sar:center_frequency': 5.405,
  'created': '2021-12-25T16:27:00.110433Z',
  'sar:observation_direction': 'right',
  'proj:shape': [5000, 5000],
  'proj:transform': [37.0, 0.0002, 0.0, 13.0, 0.0, -0.0002, 0.0, 0.0, 1.0],
  'datetime': '2020-09-02T03:16:45.539433Z',
  'cubedash:region_code': 'N12E037'},
 'geometry': {'type': 'Polygon',
  'coordinates': [[[37.017147804993954, 12.0],
    [38.0, 12.0],
    [38.0, 13.0],
    [37.22209236793799, 13.0],
    [37.14467726492754, 12.640145779886035],
    [37.017147804993954, 12.0]]]},
 'links': [{'rel': 'self',
   'href': 'https://explorer.digitalearth.africa/stac/collections/s1_rtc/items/46bbaf19-62fd-58c0-9538-9a48c2aec7cc',
   'type': 'application/json'},
  {'rel': 'odc_yaml',
   'href': 'https://explorer.digitalearth.africa/dataset/46bbaf19-62fd-58c0-9538-9a48c2aec7cc.odc-metadata.yaml',
   'type': 'text/yaml',
   'title': 'ODC Dataset YAML'},
  {'rel': 'collection',
   'href': 'https://explorer.digitalearth.africa/stac/collections/s1_rtc'},
  {'rel': 'product_overview',
   'href': 'https://explorer.digitalearth.africa/product/s1_rtc',
   'type': 'text/html',
   'title': 'ODC Product Overview'},
  {'rel': 'alternative',
   'href': 'https://explorer.digitalearth.africa/dataset/46bbaf19-62fd-58c0-9538-9a48c2aec7cc',
   'type': 'text/html',
   'title': 'ODC Dataset Overview'}],
 'assets': {'vh': {'href': 's3://deafrica-sentinel-1/s1_rtc/N12E037/2020/09/02/03F841/s1_rtc_03F841_N12E037_2020_09_02_VH.tif',
   'type': 'image/tiff; application=geotiff; profile=cloud-optimized',
   'title': 'vh',
   'eo:bands': [{'name': 'vh'}],
   'proj:epsg': 4326,
   'proj:shape': [5000, 5000],
   'proj:transform': [37.0, 0.0002, 0.0, 13.0, 0.0, -0.0002, 0.0, 0.0, 1.0],
   'roles': ['data']},
  'vv': {'href': 's3://deafrica-sentinel-1/s1_rtc/N12E037/2020/09/02/03F841/s1_rtc_03F841_N12E037_2020_09_02_VV.tif',
   'type': 'image/tiff; application=geotiff; profile=cloud-optimized',
   'title': 'vv',
   'eo:bands': [{'name': 'vv'}],
   'proj:epsg': 4326,
   'proj:shape': [5000, 5000],
   'proj:transform': [37.0, 0.0002, 0.0, 13.0, 0.0, -0.0002, 0.0, 0.0, 1.0],
   'roles': ['data']},
  'area': {'href': 's3://deafrica-sentinel-1/s1_rtc/N12E037/2020/09/02/03F841/s1_rtc_03F841_N12E037_2020_09_02_AREA.tif',
   'type': 'image/tiff; application=geotiff; profile=cloud-optimized',
   'title': 'area',
   'eo:bands': [{'name': 'area'}],
   'proj:epsg': 4326,
   'proj:shape': [5000, 5000],
   'proj:transform': [37.0, 0.0002, 0.0, 13.0, 0.0, -0.0002, 0.0, 0.0, 1.0],
   'roles': ['data']},
  'mask': {'href': 's3://deafrica-sentinel-1/s1_rtc/N12E037/2020/09/02/03F841/s1_rtc_03F841_N12E037_2020_09_02_MASK.tif',
   'type': 'image/tiff; application=geotiff; profile=cloud-optimized',
   'title': 'mask',
   'eo:bands': [{'name': 'mask'}],
   'proj:epsg': 4326,
   'proj:shape': [5000, 5000],
   'proj:transform': [37.0, 0.0002, 0.0, 13.0, 0.0, -0.0002, 0.0, 0.0, 1.0],
   'roles': ['data']},
  'angle': {'href': 's3://deafrica-sentinel-1/s1_rtc/N12E037/2020/09/02/03F841/s1_rtc_03F841_N12E037_2020_09_02_ANGLE.tif',
   'type': 'image/tiff; application=geotiff; profile=cloud-optimized',
   'title': 'angle',
   'eo:bands': [{'name': 'angle'}],
   'proj:epsg': 4326,
   'proj:shape': [5000, 5000],
   'proj:transform': [37.0, 0.0002, 0.0, 13.0, 0.0, -0.0002, 0.0, 0.0, 1.0],
   'roles': ['data']},
  'metadata': {'href': 's3://deafrica-sentinel-1/s1_rtc/N12E037/2020/09/02/03F841/s1_rtc_03F841_N12E037_2020_09_02_metadata.xml',
   'type': 'application/xml',
   'roles': ['metadata']}},
 'bbox': [37.017147804993954, 12.0, 38.0, 13.0],
 'stac_extensions': ['https://stac-extensions.github.io/eo/v1.0.0/schema.json',
  'https://stac-extensions.github.io/projection/v1.0.0/schema.json'],
 'collection': 's1_rtc'}

In [9]:

import rioxarray as rxr

In [10]:

url = 's3://deafrica-sentinel-1/s1_rtc/N12E037/2020/09/02/03F841/s1_rtc_03F841_N12E037_2020_09_02_VV.tif'

In [11]:

da = rxr.open_rasterio(url)

In [12]:

da

In [13]:

da.spatial_ref.crs_wkt

Out[13]:

'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'

In [14]:

crs = "EPSG:4326"
resolution = 0.002

ds = stac_load(
    items,
    bands=("vh", "vv"),
    crs=crs,
    resolution=resolution,
    chunks={},
    groupby="solar_day",
    stac_cfg=config,
    bbox=bbox
)

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Input In [14], in <module>
      1 crs = "EPSG:4326"
      2 resolution = 0.002
----> 4 ds = stac_load(
      5     items,
      6     bands=("vh", "vv"),
      7     crs=crs,
      8     resolution=resolution,
      9     chunks={},
     10     groupby="solar_day",
     11     stac_cfg=config,
     12     bbox=bbox
     13 )

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/site-packages/odc/stac/_load.py:454, in load(items, bands, groupby, resampling, chunks, crs, resolution, geobox, bbox, lon, lat, x, y, align, like, geopolygon, stac_cfg, patch_url, product_cache, skip_broken_datasets, progress_cbk, fuse_func, **kw)
    450     if "resolution" not in geo:
    451         # Need to pick resolution
    452         geo["resolution"] = pick_best_resolution(dss, bands)
--> 454 dss = list(stac2ds(items, stac_cfg, product_cache=product_cache))
    455 auto_fill_geo(geo, dss, bands)
    457 if patch_url is not None:

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/site-packages/odc/stac/_eo3.py:735, in stac2ds(items, cfg, product_cache)
    733 # Have not seen this collection yet, figure it out
    734 if product is None:
--> 735     product = infer_dc_product(item, cfg)
    736     products[collection_id] = product
    738 yield item_to_ds(item, product)

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/functools.py:888, in singledispatch.<locals>.wrapper(*args, **kw)
    884 if not args:
    885     raise TypeError(f'{funcname} requires at least '
    886                     '1 positional argument')
--> 888 return dispatch(args[0].__class__)(*args, **kw)

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/site-packages/odc/stac/_eo3.py:543, in infer_dc_product_from_item(item, cfg)
    539 # We assume that grouping of data bands into grids is consistent across
    540 # entire collection, so we compute it once and keep it on a product object
    541 # at least for now.
    542 if has_proj:
--> 543     _, band2grid = compute_eo3_grids(data_bands)
    544 else:
    545     band2grid = _band2grid_from_gsd(data_bands)

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/site-packages/odc/stac/_eo3.py:275, in compute_eo3_grids(assets)
    272     gsd = geobox_gsd(geobox)
    273     return f"g{gsd:g}"
--> 275 geoboxes = dicttoolz.valmap(asset_geobox, assets)
    277 # GeoBox to list of bands that share same footprint
    278 grids: Dict[GeoBox, List[str]] = {}

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/site-packages/toolz/dicttoolz.py:83, in valmap(func, d, factory)
     72 """ Apply function to values of dictionary
     73 
     74 >>> bills = {"Alice": [20, 15, 30], "Bob": [10, 35]}
   (...)
     80     itemmap
     81 """
     82 rv = factory()
---> 83 rv.update(zip(d.keys(), map(func, d.values())))
     84 return rv

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/site-packages/odc/stac/_eo3.py:242, in asset_geobox(asset)
    239     raise ValueError(f"Asset transform is not affine: {_proj.transform}")
    241 affine = Affine(*_proj.transform[:6])
--> 242 return GeoBox(w, h, affine, _proj.crs_string)

File /home/conda/store/26851070ae96f50ac395255946c7876bb86f109ca9ae22af854d21e68393d114-pangeo/lib/python3.9/site-packages/datacube/utils/geometry/_base.py:989, in GeoBox.__init__(self, width, height, affine, crs)
    988 def __init__(self, width: int, height: int, affine: Affine, crs: MaybeCRS):
--> 989     assert is_affine_st(affine), "Only axis-aligned geoboxes are currently supported"
    990     self.width = width
    991     self.height = height

AssertionError: Only axis-aligned geoboxes are currently supported

In [ ]: