<jupyter_start><jupyter_text>AOIWe first start with a spatial AOI, which should be a single GeoJSON Feature with a geometry type of Point, LineString, Polygon, MultiPoint, MultiLineString, or MultiPolygon. A fast and simple tool to create such as an AOI is http://geojson.io/. Save a GeoJSON Feature (not a FeatureCollection!) in a file accessible by this notebook.<jupyter_code>aoi_filename = "./lakota.geojson"

# read in AOI as a GeoDataFrame
import geopandas as gpd
aoi = gpd.read_file(aoi_filename)

# get the geometry of the AOI as a dictionary for use with PySTAC Client
from shapely.geometry import mapping
geom = mapping(aoi.to_dict()['geometry'][0])
geom<jupyter_output><empty_output><jupyter_text>Searching a STAC APIWe can search any STAC API by providing the root URL to the Client.open() function.We can now use the AOI geometry for the intersects parameter with **pystac-client**. Edit/add parameters below.<jupyter_code># STAC API - Landsat Collection 2
url = "https://landsatlook.usgs.gov/stac-server"

# Search parameters
params = {
    "collections": ["landsat-c2l2-sr"],
    "intersects": geom,
    "datetime": "2015-01-01/2020-12-31",
    "limit": 100,
    "query": ["platform=LANDSAT_8", "eo:cloud_cover<10"]
}

from pystac_client import Client
cat = Client.open(url)
search = cat.search(**params)

matched = search.matched()
print(f"{search.matched()} scenes found")
#! conda install -y pystac-client
%%time
from pystac import ItemCollection

# get all items found in search
items_dict = []
for item in search.get_all_items_as_dict()['features']:
    for a in item['assets']:
        if 'alternate' in item['assets'][a] and 's3' in item['assets'][a]['alternate']:
            item['assets'][a]['href'] = item['assets'][a]['alternate']['s3']['href']
        item['assets'][a]['href'] = item['assets'][a]['href'].replace('usgs-landsat-ard', 'usgs-landsat')
    items_dict.append(item)

# Create GeoDataFrame from resulting Items
items_gdf = items_to_geodataframe(items_dict)
item_collection = ItemCollection(items_dict)<jupyter_output><empty_output><jupyter_text>Using Pandas/GeoPandasWith the STAC Items in a GeoDataFrame, we can view the geometries of the Items GeoDataFrame and the AOI GeoDataFrame together.<jupyter_code>plot_polygons(items_gdf) * aoi.hvplot.paths(geo=True, line_width=3, line_color='red')<jupyter_output><empty_output><jupyter_text>STAC AssetsWe haven't looked at any data yet, as the actual data in STAC Items is stored in assets so we are limited to exploring the metadata about the data. To use the assets a user should first understand what assets are available.  Here we use Pandas to take a look at the Asset dictionary of the first Item (note that not all Items may have the same assets).<jupyter_code># use the first Item to see what assets are available

assets = pd.DataFrame.from_dict(items_dict[0]['assets'], orient='index')

for f in ['alternate', 'file:checksum', 'proj:transform', 'rel']:
    if f in assets:
        del assets[f]

assets<jupyter_output><empty_output><jupyter_text>OpenDataCube<jupyter_code># Convert the STAC item(s) to ODC datasets
import yaml

from odc import stac
from pyproj import CRS
from pystac.extensions.projection import ProjectionExtension

def open_odc(items, crs=None, resolution=None):
    configuration_str = """---
        landsat-c2l2-sr:
          measurements:
            'red':
              dtype: float32
              nodata: 0
              units: 'm'
        """
    configuration = yaml.load(configuration_str, Loader=yaml.CSafeLoader)
    datasets = list(stac.stac2ds(items, configuration))
    
    proj = ProjectionExtension.ext(items[0])
    if crs is None:
        crs = CRS.from_epsg(proj.epsg)
    if resolution is None:
        resolution = (proj.transform[4], proj.transform[0])

    #data = stac.dc_load(datasets, measurements=['red','green','blue', 'nir08', 'qa_pixel'], chunks={"x": 1024, "y": 1024}, output_crs=crs, resolution=resolution)
    data = stac.dc_load(datasets, measurements=['red','green','blue', 'nir08', ], 
      chunks={"x": 1024, "y": 1024}, output_crs=crs, resolution=resolution)
    return data
help(stac.stac2ds)
# hello
# open found items as an OpenDataCube

_datacube = open_odc(item_collection) #, 'epsg:4326', resolution='0.002')
_datacube.to_array(dim='bands')
import rioxarray

datacube = _datacube.rio.clip([geom], crs='epsg:4326')
datacube.to_array(dim='bands')
# create RGB image for display
from odc.algo import to_rgba

rgba = to_rgba(datacube, clamp=(1000,20000), bands=['red', 'green', 'blue'])
rgba
from dask.distributed import Client, progress
client = Client(processes=False, threads_per_worker=8,
                n_workers=2, memory_limit='8GB', dashboard_address=':8080')
client
rgba.compute()
rgba
rgba.hvplot.rgb('x', 'y', bands='band', groupby='time', aspect=1, frame_width=400)
datacube['red'][1].plot()
datacube['green'][10].plot(cmap='BuPu')
# %%time
# # process with dask

# from dask.distributed import wait

# rgba = cluster.persist(rgba)
# _ = wait(rgba)
# rgba
#help (Client.persist)
datacube<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1267937880872813, 'fraction_numerical': 0.018675054059367015, 'mean_word_length': 4.93006993006993, 'pattern_counts': {'":': 9, '<': 22, '<?xml version=': 0, '>': 21, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '99699', 'license': 'MIT', 'n_tokens_mistral': 1786, 'n_tokens_neox': 1691, 'n_words': 470, 'path': '00_ODC/000_so_this_works_wonderfully_save_for ugly_plots_using_ODC_for_kirill_and_damien_perhaps_missing_common_names_is_the_issue-Copy1.ipynb', 'repo_name': 'tonybutzer/i-experiments', 'size': '5087'}