Selecting satellites from landsat-c2-l2 collection

[cullen-molitor]

Hello,

I am trying to use the landsat-c2-l2 collection. I am looking for a strategy that would allow me to select a particular satellite from the collection, or at the very least exclude Landsat-7. My main issue is with the Scan Line Correction failure and the resulting stripes. I will be using a single image at a time so filling missing values with multiple images is not an option. I am open to suggestions for gap filling, but thought it would be easier to just use the sensors without the issue.

I have tried setting my query to query={"platform": {"in": ["landsat-5", "landsat-8"]}} or query={"platform": {"neq": 'landsat-7'}} but this did not seem to work. I think I have messed up the syntax but cannot figure out how.

I see that the data catalog lists the spectral bands from the different sensors, but I am unsure how to search with those names (TM_B1) instead of the STAC Key (blue).

Here is what I have tried so far, but with no luck:

import stackstac
import pystac_client
import pyproj
import planetary_computer as pc
from pystac.extensions.eo import EOExtension as eo
# from pystac.extensions.sat import SatExtension as sat

catalog = pystac_client.Client.open('https://planetarycomputer.microsoft.com/api/stac/v1')

search = catalog.search(
    collections=['landsat-c2-l2'],  
    intersects=dict(type="Point", coordinates=[20, 5]),
    datetime="2008-02-01/2008-03-01",
    query={"eo:cloud_cover": {"lt": 20}},  
    # query={"platform": {"in": ["landsat-5", "landsat-8"]}},  # Tried a few different ways to query
    # query={"platform": {"neq": 'landsat-7'}},
    # query={"eo:cloud_cover": {"lt": 10}, "platform": {"neq": 'landsat-7'}},
    limit=1,
)
items = pc.sign(search); items[0].properties

Output

{'gsd': 30,
 'created': '2022-05-06T16:39:43.054109Z',
 'sci:doi': '10.5066/P9C7I13B',
 'datetime': '2010-01-03T08:46:41.053573Z',
 'platform': 'landsat-7',
 'proj:epsg': 32634,
 'proj:shape': [6971, 7981],
 'description': 'Landsat Collection 2 Level-2',
 'instruments': ['etm+'],
 'eo:cloud_cover': 0.0,
 'proj:transform': [30.0, 0.0, 275985.0, 0.0, -30.0, 744615.0],
 'view:off_nadir': 0,
 'landsat:wrs_row': '056',
 'landsat:scene_id': 'LE71800562010003ASN00',
 'landsat:wrs_path': '180',
 'landsat:wrs_type': '2',
 'view:sun_azimuth': 135.9089404,
 'landsat:correction': 'L2SP',
 'view:sun_elevation': 49.47199538,
 'landsat:cloud_cover_land': 0.0,
 'landsat:collection_number': '02',
 'landsat:collection_category': 'T1'}

stack = stackstac.stack(items, assets=["red", "green", "blue",])
x_utm, y_utm = pyproj.Proj(stack.crs)(lon, lat)
buffer = 3000  # meters
aoi = stack.loc[..., y_utm+buffer:y_utm-buffer, x_utm-buffer:x_utm+buffer]
data = aoi.compute()
data.plot.imshow(row="time", rgb="band", robust=True, size=6)

Any suggestions or tips would be greatly appreciated!

Side note: I think a tutorial notebook to explain common ways of dealing with the SLC failure and how to implement them on MPC would be a great addition to the collection.

[TomAugspurger]

You can use the in operator to filter to specific platforms using the STAC query extension. Our example notebook now includes a section demonstrating this with pystac-client.

Using pystac-client that would be

import pystac_client

catalog = pystac_client.Client.open('https://planetarycomputer.microsoft.com/api/stac/v1')
search = catalog.search(
    collections=['landsat-c2-l2'],
    datetime="2022-02-01/2022-03-01",
    query={"platform": {"in": ["landsat-8", "landsat-9"]}},
)

next(search.get_items())

All the items returned in the search iterable will be from landsat-8 or landsat-9.

(orignal post below)

---

Agreed that a tutorial would be helpful here. I'm waiting for a couple of things to settle down before writing it up, and a couple of limitations to be removed.

Right now, I don't think that pgstac supports the {"platform": {"in": [values]}} style syntax: (edit: stac-utils/pgstac#124 is fixing this).

For now, as long as the result set isn't too large, I think you're best off filtering client-side:

search = catalog.search(
    collections=['landsat-c2-l2'],  
    intersects=dict(type="Point", coordinates=[20, 5]),
    datetime="2008-02-01/2008-03-01",
    query={"eo:cloud_cover": {"lt": 20}},  
)
items = [item for item in search.get_all_items() if item.platform != "landsat-7"]

---

(original post below)

STAC has an alternative, more powerful / verbose query language cql2. In the explorer, you can click the "Advanced" link to build query, and then "Explore results in the Hub" to copy the code. See https://planetarycomputer.microsoft.com/docs/overview/explorer/#custom-filters.

Unfortunately, due to stac-utils/pgstac#122, this will fail for more than one value. So you'll need to do two separate searches.

from pystac_client import Client
import planetary_computer as pc

# Search against the Planetary Computer STAC API
catalog = Client.open(
  "https://planetarycomputer.microsoft.com/api/stac/v1"
)

# Define your area of interest
aoi = {
  "type": "Polygon",
  "coordinates": [
    [
      [-108.92460847488842, 39.09935495125271],
      [-106.63067499990376, 39.09935495125271],
      [-106.63067499990376, 40.47212159588767],
      [-108.92460847488842, 40.47212159588767],
      [-108.92460847488842, 39.09935495125271]
    ]
  ]
}

# Define your temporal range
daterange = {"interval": ["1982-08-22T00:00:00Z", "2022-06-06T23:59:59Z"]}

# Define your search with CQL2 syntax
search = catalog.search(filter_lang="cql2-json", filter={
  "op": "and",
  "args": [
    {"op": "s_intersects", "args": [{"property": "geometry"}, aoi]},
    {"op": "anyinteracts", "args": [{"property": "datetime"}, daterange]},
    {"op": "=", "args": [{"property": "collection"}, "landsat-c2-l2"]},
    {"op": "=", "args": [{"property": "platform"}, "landsat-8"]}
  ]
})

# Grab the first item from the search results and sign the assets
first_item = next(search.get_items())
pc.sign_item(first_item).assets

Once that bug is fixed you'll be able to do select multiple and do it all in one search. And hopefully we can support the in shorthand at some point in the future.

[cullen-molitor]

Thanks, I didn't realize there was another option for structuring the search, or that the explorer would build the query for you! This is very helpful.

[TomAugspurger]

@cullen-molitor stac-utils/pgstac#124 has been merged and deployed now, so you can provided in queries now.

[gilbertocamara]

Dear @cullen-molitor @TomAugspurger, in the latest development version of the sits R package, we provide full support to use the Planetary Computer LANDSAT-C2-L2 collection, including platform selection. Even if you are not a R user, you may still run a few lines of R code to select the images you need. An exemple follows

# install sits version 1.1.0 current development package (soon on CRAN)
devtools::install_github("e-sensing/sits@dev")
# select a shapefile to serve as a region of interest
# ROI can also be defined as Lat/long boxes in WGS84
shp_file <- system.file("extdata/shapefiles/df_bsb/df_bsb.shp",
        package = "sits"
 )
# read the shapefile as an sf object
sf_bsb <- sf::read_sf(shp_file)
# select a subset of MSPC LANDSAT-C2-L2 collection
 landsat_cube <- sits::sits_cube(
                source = "MSPC",
                collection = "LANDSAT-C2-L2",
                roi = sf_bsb,
                platform = "LANDSAT-5", 
                bands = c("GREEN",  "NIR08", "SWIR22"),
                start_date = as.Date("2008-07-18"),
                end_date = as.Date("2008-10-23")
 )
# list the files you want
landsat_cube$file_info[[1]]

[TomAugspurger]

I suspect that would hit the same bug in the Planetary Computer's API server if you did platform = c("LANDSAT-5", "LANDSAT-6").

stac-utils/pgstac#124 is fixing this in pgstac, and then we'll get it deployed to staging and production.

[gilbertocamara]

Hi @TomAugspurger, indeed in sits one can select only one Landsat platform at a time. However, is this a real problem? In practice, we only had at maximum two Landsat satellites operating together at any given time, so being able to select one should solve 95% of requests.

[gilbertocamara]

The important bit (at least for me) is how concise and short should the API be. IMHO, Planetary Computer users would benefit from an API that is as clean (or cleaner) than the Google Earth Engine one.

[TomAugspurger]

Yeah, this post was specifically about filtering out Landsat-7 so I wanted to mention that the issue is server-side and will affect all STAC clients.

My earlier comment was out-of-date: pgstac now supports in-style queries, it's just hitting this same bug as the cql-2 queries. Once stac-utils/pgstac#124 is fixed and deployed, both sits (via rstac) and pystac-client will be able to do this in filter cleanly. No need for cql-2.