STAC across Languages

Some background

Over a month ago I was invited as a lecturer at the OpenGeoHub 2023 Summer School in Poznan, Poland. It was a great experience, meeting brilliant people, learning from other lecturers and enjoying the city. Many thanks to the organizers!

My lecture was on Tools and packages to query and process Sentinel-1 and Sentinel-2 data with R and Python. Following on the lecture and inspired by the geocompx blogpost on geographic data analysis in R and Python, I started the quest to do a similar post with the basic content of my lecture.

Why do I call it a quest? Well, it took me a while. Before this post, my website was based on distill, so the migration process from a distill-based website to a quarto-based one was the first task. Why did I have to do the migration? Because Quarto allows the combination of two or more programming languages in a single file, crucial for data science across languages.

I managed to migrate thanks to JaKub Nowosad’s talk at OpenGeoHub23 (Nowosad 2023), the Quarto documentation and to posts from Danielle Navarro (Navarro 2022). Big thanks!

But in any case, I can now inaugurate my new website with a post on how to query STAC data and create on-demand data cubes with R and Python! So let’s get to it.

Libraries

Let’s first load the libraries needed. Here you can get familiar with the tab setting of this post. Clicking on the R or Python tab will switch between languages. Clicking on either tab will switch the panel in the whole document.

Python

import dask.array as da # handle dask arrays
import geopandas as gpd # handle geospatial data frames
from IPython.display import Image # visualize URLs
import pystac_client # connecting to the STAC API
from rasterio.enums import Resampling # perform re-sampling operations
import shapely # create vector objects
import stackstac # build an on-demand STAC data cube

R

library(dplyr) # data wrangling
library(gdalcubes) # on-demand data cubes
library(knitr) # visualize URLs
library(rstac) # connecting to the STAC API
library(sf) # handle geospatial data frames

Querying STAC data

STAC stands for SpatioTemporal Asset Catalog and it is “a common language to describe geospatial information, so it can more easily be worked with, indexed, and discovered”.

There are several APIs available to query data, you can browse them all in the STAC catalog index. Some of these APIs will require authentication to access the data. We will use the Earth Search catalog for this post, which allows querying data on Amazon Web Services (AWS). The data we will fetch (Sentinel-2) does not require authentication.

In Python we can use the pystac-client library (Duckworth 2021), while in R we can use the rstac library (Simoes et al. 2021). Refer to their documentation pages for more information.

To start fetching data, we will open a client. We can then see the collections available for this API. In Python we call a for loop on the client collections, while in R we use the get_request() function. Both methods make an HTTP GET request to STAC web services.

Python

api_url = 'https://earth-search.aws.element84.com/v1'
client = pystac_client.Client.open(api_url)
for collection in client.get_collections():
    print(collection)

<CollectionClient id=sentinel-2-c1-l2a-preview>
<CollectionClient id=cop-dem-glo-30>
<CollectionClient id=naip>
<CollectionClient id=sentinel-2-l2a>
<CollectionClient id=sentinel-2-l1c>
<CollectionClient id=cop-dem-glo-90>
<CollectionClient id=sentinel-1-grd>
<CollectionClient id=landsat-c2-l2>

R

api_url = "https://earth-search.aws.element84.com/v1"
client = stac(api_url)
collections(client) |> get_request()

###STACCollectionList
- collections (8 item(s)):
  - sentinel-2-c1-l2a-preview
  - cop-dem-glo-30
  - naip
  - sentinel-2-l2a
  - sentinel-2-l1c
  - cop-dem-glo-90
  - sentinel-1-grd
  - landsat-c2-l2
- field(s): collections, links, context

We will choose here Sentinel-2 Level 2a data, which provides atmospherically corrected data representing surface reflectance.

We can also set the coordinates for our query, in this case a POINT geometry, which will help to retrieve scenes that intersect with it. At last we can also set a time extent, where we define the start and end of our query.

Python

# collection ID
collection = 'sentinel-2-l2a'
# coordinates
lon = 16.9
lat = 52.4
point = shapely.Point(lon, lat)
# date range
time_extent = '2022-05-01/2022-10-01'

R

# collection ID
collection = 'sentinel-2-l2a'
# coordinates
lon = 16.9
lat = 52.4
point = st_point(c(lon, lat))
# date range
time_extent = c("2022-05-01", "2022-10-01")

Now, we can pass the arguments to our search. Note for the R search we had to slightly format the datetime parameter. This is because of a bug when fetching the data from STAC.

You will see we can also query based on scene metadata. In this example we look for scenes with cloud cover lower than 10% for the whole scene. In Python we do the filtering within the search, while in R we perform it after the HTTP POST request.

Finally, we check how many items we fetched with our query, also to compare that R and Python give us the same results.

Python

search = client.search(
    collections=[collection],
    intersects=point,
    datetime=time_extent,
    query=["eo:cloud_cover<10"],
)
items = search.item_collection()
len(items)

8

R

items = client  |> 
  stac_search(
      collections = collection,
      intersects = point, 
      datetime = paste0(time_extent,"T00:00:00Z", collapse = "/"),
      limit = 100
    ) |> 
  post_request() |> 
  items_filter(
    filter_fn = function(x) {x$properties$`eo:cloud_cover` < 10}
  )
items_length(items)

[1] 8

To explore the items as data frames we can either convert to a geopandas table in Python or to an sf table in R. We added an extra fid column to allow index matching to the original item list.

Python

items_df = gpd.GeoDataFrame.from_features(items.to_dict(), crs="epsg:4326")
items_df

                                            geometry  ...                   updated
0  POLYGON ((16.49848 53.24020, 17.39166 53.22562...  ...  2022-11-06T13:14:57.270Z
1  POLYGON ((16.49848 53.24020, 17.39346 53.22558...  ...  2022-11-06T12:53:48.090Z
2  POLYGON ((16.49848 53.24020, 18.14175 53.20819...  ...  2022-11-06T13:12:58.202Z
3  POLYGON ((16.49848 53.24020, 18.14175 53.20819...  ...  2022-11-06T13:14:37.621Z
4  POLYGON ((16.49848 53.24020, 18.14175 53.20819...  ...  2022-11-06T12:54:05.312Z
5  POLYGON ((16.49848 53.24020, 17.39346 53.22558...  ...  2022-11-06T12:55:07.931Z
6  POLYGON ((16.49848 53.24020, 18.14175 53.20819...  ...  2022-11-06T12:52:50.337Z
7  POLYGON ((16.49848 53.24020, 18.14175 53.20819...  ...  2022-11-06T13:16:38.195Z

[8 rows x 42 columns]

R

# we add an extra column to be able to subset the original item list
items_df = items_as_sf(items) |> 
  mutate(fid = row_number())
# wrap function around st_as_sf and as_tibble just for nicer printing
st_as_sf(as_tibble(items_df))

Simple feature collection with 8 features and 42 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 16.465 ymin: 52.22258 xmax: 18.14175 ymax: 53.2402
Geodetic CRS:  WGS 84
# A tibble: 8 × 43
  `s2:generation_time`        s2:medium_proba_clouds_percenta…¹ `s2:product_uri`
  <chr>                                                   <dbl> <chr>           
1 2022-08-26T16:20:59.000000Z                           0.00172 S2A_MSIL2A_2022…
2 2022-08-16T16:25:57.000000Z                           2.27    S2A_MSIL2A_2022…
3 2022-08-03T13:51:55.000000Z                           5.69    S2A_MSIL2A_2022…
4 2022-07-19T11:39:43.000000Z                           0.00007 S2B_MSIL2A_2022…
5 2022-07-14T17:50:57.000000Z                           1.29    S2A_MSIL2A_2022…
6 2022-06-27T16:28:10.000000Z                           0.732   S2A_MSIL2A_2022…
7 2022-06-24T14:39:14.000000Z                           0.00130 S2A_MSIL2A_2022…
8 2022-06-19T11:43:29.000000Z                           0.0514  S2B_MSIL2A_2022…
# ℹ abbreviated name: ¹​`s2:medium_proba_clouds_percentage`
# ℹ 40 more variables: `s2:processing_baseline` <chr>, `grid:code` <chr>,
#   `s2:product_type` <chr>, `s2:snow_ice_percentage` <dbl>,
#   `s2:vegetation_percentage` <dbl>, `processing:software` <chr>,
#   `s2:high_proba_clouds_percentage` <dbl>,
#   `s2:unclassified_percentage` <dbl>, `s2:water_percentage` <dbl>,
#   `s2:not_vegetated_percentage` <dbl>, `s2:cloud_shadow_percentage` <dbl>, …

When we have the items as data frame, we can further filter the data based on the table columns. Below, we filter for data with a cloud cover lower than 2% and with a no-data pixel percentage below 10%.

To illustrate how to fetch item properties, we can select the first item in our list and get the datetime of this scene.

Python

ids = items_df.loc[
  (items_df['eo:cloud_cover'] <= 2) &
  (items_df['s2:nodata_pixel_percentage'] <= 10)
]
item = items[ids.index[0]]
item.datetime

datetime.datetime(2022, 7, 19, 10, 6, 8, 645000, tzinfo=tzutc())

R

ids = items_df |> 
  filter(
    `eo:cloud_cover` <= 2, `s2:nodata_pixel_percentage` <= 10
  ) |> 
  pull(fid)
item = items$features[[ids[1]]]
item$properties$datetime

[1] "2022-07-19T10:06:08.645000Z"

We can also take a look at the selected scene by fetching the item thumbnail. Here you will see that both R and Python return the same scene, as they should! 😉.

Python

thumbnail = item.assets["thumbnail"].href
Image(url = thumbnail)

R

item |> 
  assets_url(asset_names = "thumbnail") |> 
  include_graphics()

Creating a STAC data cube

Fetching individual scenes is interesting, but the real power of working with the STAC API is the possibility to work with big EO data without the need to download a single pixel to your own computer!

For this, we can create an on-demand data cube that can help us organize and analyze our data, doing the data wrangling needed and only fetching the data we need as a final step.

To create on-demand data cubes, we will use R’s gdalcubes (Appel, Pebesma, and Mohr 2021) and Python’s stackstac (Joseph 2023).

We will create our data cube on a smaller subset. Here we load a GeoJSON file of Poznan, where the summer school took place. We read it in and extract its bounding box.

Python

poznan = "https://raw.githubusercontent.com/loreabad6/ogh23/main/data/poznan.geojson"
footprint = gpd.read_file(poznan)
bbox = footprint.total_bounds

R

poznan = "https://raw.githubusercontent.com/loreabad6/ogh23/main/data/poznan.geojson"
footprint = read_sf(poznan)
bbox = footprint |> st_transform(crs = 32633) |> st_bbox()

Now we can create our cube. The following code chunk has several steps depending if you are in the R or Python tab. I kept it as a single block since this is how we get an on-demand data cube on each end. In essence what we do is:

1. Define the assets we want.

I selected RGB+NIR bands to illustrate this process, but of course we can set any other band. I included the SCL band, which we will use to mask clouds.

2. Set-up the cube.

• In Python we do this with the stackstac.stack function, where we can set the resolution, assets, bounding box, and re-sampling method for the data. Separately we define a cloud mask to get remove invalid pixels.
• In R, this works in a lengthier way, using more intermediate functions. We first define the collection with the assets with the stac_image_collection function. Then, we define an extent in the form of a list. This extent will have the start and end date as well as the bounding box of our cube. Now we can define a view, which defines the CRS, spatial and temporal resolution, and the aggregation and re-sampling method. We finally define the cloud mask with image_mask and can then build a raster_cube.

Python

assets = ["red","green","blue","nir","scl"]
cube_all = stackstac.stack(
    items,
    assets=assets,
    resolution=200,
    bounds_latlon=bbox,
    resampling=Resampling.bilinear
)
scl = cube_all.sel(band=["scl"])
s2_mask = da.isin(scl, [3,8,9])
cube = cube_all.where(~s2_mask)
cube

<xarray.DataArray 'stackstac-488683ba46ca080d7bae3993b8b89142' (time: 8,
                                                                band: 5,
                                                                y: 140, x: 143)>
dask.array<where, shape=(8, 5, 140, 143), dtype=float64, chunksize=(1, 1, 140, 143), chunktype=numpy.ndarray>
Coordinates: (12/52)
  * time                                     (time) datetime64[ns] 2022-06-19...
    id                                       (time) <U24 'S2B_33UXU_20220619_...
  * band                                     (band) <U5 'red' 'green' ... 'scl'
  * x                                        (x) float64 6.19e+05 ... 6.474e+05
  * y                                        (y) float64 5.823e+06 ... 5.795e+06
    s2:generation_time                       (time) <U27 '2022-06-19T11:43:29...
    ...                                       ...
    raster:bands                             (band) object {'nodata': 0, 'dat...
    title                                    (band) <U30 'Red (band 4) - 10m'...
    common_name                              (band) object 'red' ... None
    center_wavelength                        (band) object 0.665 0.56 ... None
    full_width_half_max                      (band) object 0.038 0.045 ... None
    epsg                                     int32 32633
Attributes:
    spec:        RasterSpec(epsg=32633, bounds=(619000, 5795000, 647600, 5823...
    crs:         epsg:32633
    transform:   | 200.00, 0.00, 619000.00|\n| 0.00,-200.00, 5823000.00|\n| 0...
    resolution:  200

xarray.DataArray

'stackstac-488683ba46ca080d7bae3993b8b89142'

• time: 8
• band: 5
• y: 140
• x: 143

• dask.array<chunksize=(1, 1, 140, 143), meta=np.ndarray>

  Array Chunk Bytes 6.11 MiB 156.41 kiB Shape (8, 5, 140, 143) (1, 1, 140, 143) Dask graph 40 chunks in 10 graph layers Data type float64 numpy.ndarray

• Coordinates: (52)
  • time
    (time)
    datetime64[ns]
    2022-06-19T10:06:08.087000 ... 2...

    array(['2022-06-19T10:06:08.087000000', '2022-06-24T10:06:16.177000000',
           '2022-06-27T10:16:12.822000000', '2022-07-14T10:06:16.331000000',
           '2022-07-19T10:06:08.645000000', '2022-08-03T10:06:14.409000000',
           '2022-08-16T10:16:12.609000000', '2022-08-26T10:16:12.680000000'],
          dtype='datetime64[ns]')

  • id
    (time)
    <U24
    'S2B_33UXU_20220619_0_L2A' ... '...

    array(['S2B_33UXU_20220619_0_L2A', 'S2A_33UXU_20220624_0_L2A',
           'S2A_33UXU_20220627_0_L2A', 'S2A_33UXU_20220714_0_L2A',
           'S2B_33UXU_20220719_0_L2A', 'S2A_33UXU_20220803_0_L2A',
           'S2A_33UXU_20220816_0_L2A', 'S2A_33UXU_20220826_0_L2A'],
          dtype='<U24')

  • band
    (band)
    <U5
    'red' 'green' 'blue' 'nir' 'scl'

    array(['red', 'green', 'blue', 'nir', 'scl'], dtype='<U5')

  • x
    (x)
    float64
    6.19e+05 6.192e+05 ... 6.474e+05

    array([619000., 619200., 619400., 619600., 619800., 620000., 620200., 620400.,
           620600., 620800., 621000., 621200., 621400., 621600., 621800., 622000.,
           622200., 622400., 622600., 622800., 623000., 623200., 623400., 623600.,
           623800., 624000., 624200., 624400., 624600., 624800., 625000., 625200.,
           625400., 625600., 625800., 626000., 626200., 626400., 626600., 626800.,
           627000., 627200., 627400., 627600., 627800., 628000., 628200., 628400.,
           628600., 628800., 629000., 629200., 629400., 629600., 629800., 630000.,
           630200., 630400., 630600., 630800., 631000., 631200., 631400., 631600.,
           631800., 632000., 632200., 632400., 632600., 632800., 633000., 633200.,
           633400., 633600., 633800., 634000., 634200., 634400., 634600., 634800.,
           635000., 635200., 635400., 635600., 635800., 636000., 636200., 636400.,
           636600., 636800., 637000., 637200., 637400., 637600., 637800., 638000.,
           638200., 638400., 638600., 638800., 639000., 639200., 639400., 639600.,
           639800., 640000., 640200., 640400., 640600., 640800., 641000., 641200.,
           641400., 641600., 641800., 642000., 642200., 642400., 642600., 642800.,
           643000., 643200., 643400., 643600., 643800., 644000., 644200., 644400.,
           644600., 644800., 645000., 645200., 645400., 645600., 645800., 646000.,
           646200., 646400., 646600., 646800., 647000., 647200., 647400.])

  • y
    (y)
    float64
    5.823e+06 5.823e+06 ... 5.795e+06

    array([5823000., 5822800., 5822600., 5822400., 5822200., 5822000., 5821800.,
           5821600., 5821400., 5821200., 5821000., 5820800., 5820600., 5820400.,
           5820200., 5820000., 5819800., 5819600., 5819400., 5819200., 5819000.,
           5818800., 5818600., 5818400., 5818200., 5818000., 5817800., 5817600.,
           5817400., 5817200., 5817000., 5816800., 5816600., 5816400., 5816200.,
           5816000., 5815800., 5815600., 5815400., 5815200., 5815000., 5814800.,
           5814600., 5814400., 5814200., 5814000., 5813800., 5813600., 5813400.,
           5813200., 5813000., 5812800., 5812600., 5812400., 5812200., 5812000.,
           5811800., 5811600., 5811400., 5811200., 5811000., 5810800., 5810600.,
           5810400., 5810200., 5810000., 5809800., 5809600., 5809400., 5809200.,
           5809000., 5808800., 5808600., 5808400., 5808200., 5808000., 5807800.,
           5807600., 5807400., 5807200., 5807000., 5806800., 5806600., 5806400.,
           5806200., 5806000., 5805800., 5805600., 5805400., 5805200., 5805000.,
           5804800., 5804600., 5804400., 5804200., 5804000., 5803800., 5803600.,
           5803400., 5803200., 5803000., 5802800., 5802600., 5802400., 5802200.,
           5802000., 5801800., 5801600., 5801400., 5801200., 5801000., 5800800.,
           5800600., 5800400., 5800200., 5800000., 5799800., 5799600., 5799400.,
           5799200., 5799000., 5798800., 5798600., 5798400., 5798200., 5798000.,
           5797800., 5797600., 5797400., 5797200., 5797000., 5796800., 5796600.,
           5796400., 5796200., 5796000., 5795800., 5795600., 5795400., 5795200.])

  • s2:generation_time
    (time)
    <U27
    '2022-06-19T11:43:29.000000Z' .....

    array(['2022-06-19T11:43:29.000000Z', '2022-06-24T14:39:14.000000Z',
           '2022-06-27T16:28:10.000000Z', '2022-07-14T17:50:57.000000Z',
           '2022-07-19T11:39:43.000000Z', '2022-08-03T13:51:55.000000Z',
           '2022-08-16T16:25:57.000000Z', '2022-08-26T16:20:59.000000Z'],
          dtype='<U27')

  • instruments
    ()
    <U3
    'msi'

    array('msi', dtype='<U3')

  • s2:snow_ice_percentage
    ()
    int32
    0

    array(0)

  • s2:thin_cirrus_percentage
    (time)
    float64
    7.433 0.000713 ... 0.4116 0.001074

    array([7.43253e+00, 7.13000e-04, 5.00000e-05, 5.48640e-02, 5.43800e-03,
           1.80420e-02, 4.11558e-01, 1.07400e-03])

  • s2:datatake_type
    ()
    <U8
    'INS-NOBS'

    array('INS-NOBS', dtype='<U8')

  • earthsearch:s3_path
    (time)
    <U79
    's3://sentinel-cogs/sentinel-s2-...

    array(['s3://sentinel-cogs/sentinel-s2-l2a-cogs/33/U/XU/2022/6/S2B_33UXU_20220619_0_L2A',
           's3://sentinel-cogs/sentinel-s2-l2a-cogs/33/U/XU/2022/6/S2A_33UXU_20220624_0_L2A',
           's3://sentinel-cogs/sentinel-s2-l2a-cogs/33/U/XU/2022/6/S2A_33UXU_20220627_0_L2A',
           's3://sentinel-cogs/sentinel-s2-l2a-cogs/33/U/XU/2022/7/S2A_33UXU_20220714_0_L2A',
           's3://sentinel-cogs/sentinel-s2-l2a-cogs/33/U/XU/2022/7/S2B_33UXU_20220719_0_L2A',
           's3://sentinel-cogs/sentinel-s2-l2a-cogs/33/U/XU/2022/8/S2A_33UXU_20220803_0_L2A',
           's3://sentinel-cogs/sentinel-s2-l2a-cogs/33/U/XU/2022/8/S2A_33UXU_20220816_0_L2A',
           's3://sentinel-cogs/sentinel-s2-l2a-cogs/33/U/XU/2022/8/S2A_33UXU_20220826_0_L2A'],
          dtype='<U79')

  • s2:cloud_shadow_percentage
    (time)
    object
    0 0.000989 1.168417 ... 1.743145 0

    array([0, 0.000989, 1.168417, 2.347726, 0, 4.482971, 1.743145, 0],
          dtype=object)

  • s2:water_percentage
    (time)
    float64
    1.067 1.266 ... 0.7679 0.7824

    array([1.066825, 1.266303, 0.75182 , 1.253815, 1.263951, 1.116363,
           0.767899, 0.782386])

  • grid:code
    ()
    <U10
    'MGRS-33UXU'

    array('MGRS-33UXU', dtype='<U10')

  • earthsearch:payload_id
    (time)
    <U74
    'roda-sentinel2/workflow-sentine...

    array(['roda-sentinel2/workflow-sentinel2-to-stac/58f82edc8c40c1d0d56d8a26269d02eb',
           'roda-sentinel2/workflow-sentinel2-to-stac/9b8a11349c8bf4f2437e596421530090',
           'roda-sentinel2/workflow-sentinel2-to-stac/fd889d2975690939ea2990c8b784fc31',
           'roda-sentinel2/workflow-sentinel2-to-stac/1cea63f7fa1542d8a1c7bcd0b774584d',
           'roda-sentinel2/workflow-sentinel2-to-stac/4f39107c8fe71ec8f29d674a811d6c21',
           'roda-sentinel2/workflow-sentinel2-to-stac/6245fb70f8c638a7681c672a7bf7a006',
           'roda-sentinel2/workflow-sentinel2-to-stac/b5b785b4fac74933d25c20aec8ce66bd',
           'roda-sentinel2/workflow-sentinel2-to-stac/bff15f451a9bb2107d01668bb557f526'],
          dtype='<U74')

  • mgrs:latitude_band
    ()
    <U1
    'U'

    array('U', dtype='<U1')

  • s2:nodata_pixel_percentage
    (time)
    object
    0 0 ... 60.335821 60.458738

    array([0, 0, 60.372794, 0, 3e-06, 3e-06, 60.335821, 60.458738],
          dtype=object)

  • mgrs:utm_zone
    ()
    int32
    33

    array(33)

  • mgrs:grid_square
    ()
    <U2
    'XU'

    array('XU', dtype='<U2')

  • earthsearch:boa_offset_applied
    ()
    bool
    True

    array(True)

  • s2:sequence
    ()
    <U1
    '0'

    array('0', dtype='<U1')

  • s2:processing_baseline
    ()
    <U5
    '04.00'

    array('04.00', dtype='<U5')

  • view:sun_azimuth
    (time)
    float64
    158.9 158.5 162.7 ... 165.2 166.8

    array([158.91324506, 158.50630661, 162.6663807 , 157.84964052,
           157.99361406, 159.43576449, 165.15376613, 166.81057205])

  • platform
    (time)
    <U11
    'sentinel-2b' ... 'sentinel-2a'

    array(['sentinel-2b', 'sentinel-2a', 'sentinel-2a', 'sentinel-2a',
           'sentinel-2b', 'sentinel-2a', 'sentinel-2a', 'sentinel-2a'],
          dtype='<U11')

  • s2:product_uri
    (time)
    <U65
    'S2B_MSIL2A_20220619T100029_N040...

    array(['S2B_MSIL2A_20220619T100029_N0400_R122_T33UXU_20220619T114329.SAFE',
           'S2A_MSIL2A_20220624T100041_N0400_R122_T33UXU_20220624T143914.SAFE',
           'S2A_MSIL2A_20220627T100611_N0400_R022_T33UXU_20220627T162810.SAFE',
           'S2A_MSIL2A_20220714T100041_N0400_R122_T33UXU_20220714T175057.SAFE',
           'S2B_MSIL2A_20220719T095559_N0400_R122_T33UXU_20220719T113943.SAFE',
           'S2A_MSIL2A_20220803T100041_N0400_R122_T33UXU_20220803T135155.SAFE',
           'S2A_MSIL2A_20220816T100611_N0400_R022_T33UXU_20220816T162557.SAFE',
           'S2A_MSIL2A_20220826T100611_N0400_R022_T33UXU_20220826T162059.SAFE'],
          dtype='<U65')

  • created
    (time)
    <U24
    '2022-11-06T13:16:38.195Z' ... '...

    array(['2022-11-06T13:16:38.195Z', '2022-11-06T12:52:50.337Z',
           '2022-11-06T12:55:07.931Z', '2022-11-06T12:54:05.312Z',
           '2022-11-06T13:14:37.621Z', '2022-11-06T13:12:58.202Z',
           '2022-11-06T12:53:48.090Z', '2022-11-06T13:14:57.270Z'],
          dtype='<U24')

  • processing:software
    ()
    object
    {'sentinel2-to-stac': '0.1.0'}

    array({'sentinel2-to-stac': '0.1.0'}, dtype=object)

  • s2:vegetation_percentage
    (time)
    float64
    80.88 87.55 84.59 ... 55.84 64.06

    array([80.8779  , 87.545991, 84.593809, 64.089584, 64.467937, 49.599385,
           55.840641, 64.060515])

  • s2:medium_proba_clouds_percentage
    (time)
    float64
    0.05143 0.001297 ... 2.271 0.00172

    array([5.143000e-02, 1.297000e-03, 7.322780e-01, 1.287753e+00,
           7.000000e-05, 5.689905e+00, 2.270966e+00, 1.720000e-03])

  • proj:epsg
    ()
    int32
    32633

    array(32633)

  • s2:product_type
    ()
    <U7
    'S2MSI2A'

    array('S2MSI2A', dtype='<U7')

  • s2:high_proba_clouds_percentage
    (time)
    float64
    0.007488 0.000484 ... 0.000537

    array([7.488000e-03, 4.840000e-04, 3.348470e-01, 4.827590e-01,
           3.000000e-06, 3.067140e+00, 2.117714e+00, 5.370000e-04])

  • constellation
    ()
    <U10
    'sentinel-2'

    array('sentinel-2', dtype='<U10')

  • s2:saturated_defective_pixel_percentage
    ()
    int32
    0

    array(0)

  • s2:dark_features_percentage
    (time)
    float64
    0.002057 0.004107 ... 0.01486

    array([0.002057, 0.004107, 0.006966, 0.00361 , 0.004429, 0.003646,
           0.010289, 0.01486 ])

  • s2:datatake_id
    (time)
    <U34
    'GS2B_20220619T100029_027604_N04...

    array(['GS2B_20220619T100029_027604_N04.00',
           'GS2A_20220624T100041_036584_N04.00',
           'GS2A_20220627T100611_036627_N04.00',
           'GS2A_20220714T100041_036870_N04.00',
           'GS2B_20220719T095559_028033_N04.00',
           'GS2A_20220803T100041_037156_N04.00',
           'GS2A_20220816T100611_037342_N04.00',
           'GS2A_20220826T100611_037485_N04.00'], dtype='<U34')

  • s2:reflectance_conversion_factor
    (time)
    float64
    0.9691 0.9683 ... 0.9741 0.9779

    array([0.96912679, 0.96833888, 0.96797547, 0.96750684, 0.96788495,
           0.97039583, 0.97414715, 0.9779224 ])

  • view:sun_elevation
    (time)
    float64
    59.41 59.34 59.71 ... 50.17 46.97

    array([59.4067386 , 59.33648111, 59.7142557 , 57.43426893, 56.58894221,
           53.33325433, 50.17160714, 46.96899755])

  • s2:datastrip_id
    (time)
    <U64
    'S2B_OPER_MSI_L2A_DS_2BPS_202206...

    array(['S2B_OPER_MSI_L2A_DS_2BPS_20220619T114329_S20220619T100031_N04.00',
           'S2A_OPER_MSI_L2A_DS_ATOS_20220624T143914_S20220624T100039_N04.00',
           'S2A_OPER_MSI_L2A_DS_ATOS_20220627T162810_S20220627T100933_N04.00',
           'S2A_OPER_MSI_L2A_DS_ATOS_20220714T175057_S20220714T100046_N04.00',
           'S2B_OPER_MSI_L2A_DS_2BPS_20220719T113943_S20220719T095747_N04.00',
           'S2A_OPER_MSI_L2A_DS_ATOS_20220803T135155_S20220803T100037_N04.00',
           'S2A_OPER_MSI_L2A_DS_ATOS_20220816T162557_S20220816T100922_N04.00',
           'S2A_OPER_MSI_L2A_DS_ATOS_20220826T162059_S20220826T100933_N04.00'],
          dtype='<U64')

  • updated
    (time)
    <U24
    '2022-11-06T13:16:38.195Z' ... '...

    array(['2022-11-06T13:16:38.195Z', '2022-11-06T12:52:50.337Z',
           '2022-11-06T12:55:07.931Z', '2022-11-06T12:54:05.312Z',
           '2022-11-06T13:14:37.621Z', '2022-11-06T13:12:58.202Z',
           '2022-11-06T12:53:48.090Z', '2022-11-06T13:14:57.270Z'],
          dtype='<U24')

  • s2:granule_id
    (time)
    <U62
    'S2B_OPER_MSI_L2A_TL_2BPS_202206...

    array(['S2B_OPER_MSI_L2A_TL_2BPS_20220619T114329_A027604_T33UXU_N04.00',
           'S2A_OPER_MSI_L2A_TL_ATOS_20220624T143914_A036584_T33UXU_N04.00',
           'S2A_OPER_MSI_L2A_TL_ATOS_20220627T162810_A036627_T33UXU_N04.00',
           'S2A_OPER_MSI_L2A_TL_ATOS_20220714T175057_A036870_T33UXU_N04.00',
           'S2B_OPER_MSI_L2A_TL_2BPS_20220719T113943_A028033_T33UXU_N04.00',
           'S2A_OPER_MSI_L2A_TL_ATOS_20220803T135155_A037156_T33UXU_N04.00',
           'S2A_OPER_MSI_L2A_TL_ATOS_20220816T162557_A037342_T33UXU_N04.00',
           'S2A_OPER_MSI_L2A_TL_ATOS_20220826T162059_A037485_T33UXU_N04.00'],
          dtype='<U62')

  • s2:not_vegetated_percentage
    (time)
    float64
    10.43 11.0 12.07 ... 35.87 34.91

    array([10.429534, 11.000518, 12.070826, 30.188122, 34.062394, 34.022531,
           35.867879, 34.906355])

  • s2:unclassified_percentage
    (time)
    float64
    0.1322 0.1796 ... 0.9699 0.2326

    array([0.132232, 0.179595, 0.340984, 0.291771, 0.195772, 2.000013,
           0.96991 , 0.232552])

  • s2:degraded_msi_data_percentage
    (time)
    object
    0 0 0.0045 0 0 0 0.0045 0.0044

    array([0, 0, 0.0045, 0, 0, 0, 0.0045, 0.0044], dtype=object)

  • eo:cloud_cover
    (time)
    float64
    7.491 0.002495 ... 4.8 0.003331

    array([7.491449e+00, 2.495000e-03, 1.067175e+00, 1.825375e+00,
           5.511000e-03, 8.775087e+00, 4.800238e+00, 3.331000e-03])

  • gsd
    (band)
    object
    10 10 10 10 None

    array([10, 10, 10, 10, None], dtype=object)

  • raster:bands
    (band)
    object
    {'nodata': 0, 'data_type': 'uint...

    array([{'nodata': 0, 'data_type': 'uint16', 'bits_per_sample': 15, 'spatial_resolution': 10, 'scale': 0.0001, 'offset': -0.1},
           {'nodata': 0, 'data_type': 'uint16', 'bits_per_sample': 15, 'spatial_resolution': 10, 'scale': 0.0001, 'offset': -0.1},
           {'nodata': 0, 'data_type': 'uint16', 'bits_per_sample': 15, 'spatial_resolution': 10, 'scale': 0.0001, 'offset': -0.1},
           {'nodata': 0, 'data_type': 'uint16', 'bits_per_sample': 15, 'spatial_resolution': 10, 'scale': 0.0001, 'offset': -0.1},
           {'nodata': 0, 'data_type': 'uint8', 'spatial_resolution': 20}],
          dtype=object)

  • title
    (band)
    <U30
    'Red (band 4) - 10m' ... 'Scene ...

    array(['Red (band 4) - 10m', 'Green (band 3) - 10m',
           'Blue (band 2) - 10m', 'NIR 1 (band 8) - 10m',
           'Scene classification map (SCL)'], dtype='<U30')

  • common_name
    (band)
    object
    'red' 'green' 'blue' 'nir' None

    array(['red', 'green', 'blue', 'nir', None], dtype=object)

  • center_wavelength
    (band)
    object
    0.665 0.56 0.49 0.842 None

    array([0.665, 0.56, 0.49, 0.842, None], dtype=object)

  • full_width_half_max
    (band)
    object
    0.038 0.045 0.098 0.145 None

    array([0.038, 0.045, 0.098, 0.145, None], dtype=object)

  • epsg
    ()
    int32
    32633

    array(32633)

• Indexes: (4)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2022-06-19 10:06:08.087000', '2022-06-24 10:06:16.177000',
                   '2022-06-27 10:16:12.822000', '2022-07-14 10:06:16.331000',
                   '2022-07-19 10:06:08.645000', '2022-08-03 10:06:14.409000',
                   '2022-08-16 10:16:12.609000', '2022-08-26 10:16:12.680000'],
                  dtype='datetime64[ns]', name='time', freq=None))

  • band
    PandasIndex

    PandasIndex(Index(['red', 'green', 'blue', 'nir', 'scl'], dtype='object', name='band'))

  • x
    PandasIndex

    PandasIndex(Index([619000.0, 619200.0, 619400.0, 619600.0, 619800.0, 620000.0, 620200.0,
           620400.0, 620600.0, 620800.0,
           ...
           645600.0, 645800.0, 646000.0, 646200.0, 646400.0, 646600.0, 646800.0,
           647000.0, 647200.0, 647400.0],
          dtype='float64', name='x', length=143))

  • y
    PandasIndex

    PandasIndex(Index([5823000.0, 5822800.0, 5822600.0, 5822400.0, 5822200.0, 5822000.0,
           5821800.0, 5821600.0, 5821400.0, 5821200.0,
           ...
           5797000.0, 5796800.0, 5796600.0, 5796400.0, 5796200.0, 5796000.0,
           5795800.0, 5795600.0, 5795400.0, 5795200.0],
          dtype='float64', name='y', length=140))

• Attributes: (4)

  spec :

  RasterSpec(epsg=32633, bounds=(619000, 5795000, 647600, 5823000), resolutions_xy=(200, 200))

  crs :

  epsg:32633

  transform :

  | 200.00, 0.00, 619000.00| | 0.00,-200.00, 5823000.00| | 0.00, 0.00, 1.00|

  resolution :

  200

R

assets = c("red","green","blue","nir","scl")
col = stac_image_collection(items$features, asset_names = assets)
extent = list(
  t0 = time_extent[1], t1 = time_extent[2],
  left = bbox["xmin"], bottom = bbox["ymin"], 
  right = bbox["xmax"], top = bbox["ymax"]
)
view = cube_view(
  srs = "epsg:32633",  
  extent = extent,
  dx = 200, dy = 200, dt = "P1D",
  aggregation = "median",
  resampling = "bilinear"
)
s2_mask = image_mask("scl", values = c(3,8,9))
cube = raster_cube(col, view, mask = s2_mask)
cube

A data cube proxy object

Dimensions:
               low             high count pixel_size chunk_size
t       2022-05-01       2022-10-01   154        P1D          1
y 5795085.08512264 5822885.08512264   139        200        128
x 619173.072575812 647573.072575812   142        200        128

Bands:
   name offset scale nodata unit
1  blue      0     1    NaN     
2 green      0     1    NaN     
3   nir      0     1    NaN     
4   red      0     1    NaN     
5   scl      0     1    NaN     

And now we have a data cube! Let’s so a simple reduction over time of our data to view an RGB composite.

Python

rgb = cube.sel(band=["red", "green", "blue"]).median("time", keep_attrs=True)
rgb = rgb.compute()
rgb.plot.imshow(
    rgb="band",
    robust=True
)

R

gdalcubes_options(parallel = 4, cache = TRUE) 
rgb = cube |> 
    select_bands(c("blue","green","red")) |> 
    reduce_time(c("median(blue)", "median(green)", "median(red)")) 
rgb |> 
    plot(rgb = 3:1) 

And there we have it, a side by side how-to create on-demand STAC data cubes. If you want to get more functionalities and examples on what can be done with the packages please refer to the materials of my lecture: Tools and packages to query and process Sentinel-1 and Sentinel-2 data with R and Python.

Happy STACing!

References

Appel, Marius, Edzer Pebesma, and Matthias Mohr. 2021. “Cloud-Based Processing of Satellite Image Collections in r Using STAC, COGs, and on-Demand Data Cubes.” 2021. https://r-spatial.org/r/2021/04/23/cloud-based-cubes.html.

Duckworth, Jon. 2021. “STAC Python Client.” 2021. https://pystac-client.readthedocs.io/en/stable/index.html.

Joseph, Gabe. 2023. “Stackstac: Easier Cloud-Native Geoprocessing.” 2023. https://stackstac.readthedocs.io/en/latest/.

Navarro, Danielle. 2022. “Porting a Distill Blog to Quarto.” April 20, 2022. https://blog.djnavarro.net/posts/2022-04-20_porting-to-quarto.

Nowosad, Jakub. 2023. “Sharing Your Geospatial Knowledge in the Open.” 2023. https://www.youtube.com/watch?v=AUt-GOZJQaY&ab_channel=TomislavHengl%28OpenGeoHubFoundation%29.

Simoes, Rolf, Felipe Souza, Matheus Zaglia, Gilberto Ribeiro Queiroz, Rafael Santos, and Karine Ferreira. 2021. “Rstac: An r Package to Access Spatiotemporal Asset Catalog Satellite Imagery.” In 2021 IEEE International Geoscience and Remote Sensing Symposium IGARSS, 7674–77. https://doi.org/10.1109/IGARSS47720.2021.9553518.