Ctrl+K

Quickstart with PySTAC-Client

Quickstart with PySTAC-Client#

In this guide, we will show you how to use the pystac-client to:

  • List the available data collections in the catalog

  • Search the catalog using your AOI and TOI

  • Download data to your local machine

  • Access a subset of data from a cloud optimized geotiff

Setting up Dependencies and Auth#

First we need to install the pystac-client and rasterio, which we will use later for accessing data directly from a cloud optimized geotiff.

%pip install pystac-client
%pip install rasterio

import json
import pystac
from pystac_client import Client

Next we will set up basic authentication and open a connection to the STAC API.

import base64
# Replace with your username and password
userpass = 'username:password'
b64 = base64.b64encode(userpass.encode()).decode()
headers = {'Authorization':'Basic ' + b64}

catalog = Client.open('https://fusion-stac.hydrosat.com/', headers)

List the Available Collections#

Now let’s see what collections are available through the API.

collections = catalog.get_children()
for collection in collections:
    print(f"{collection.id} - {collection.title}")

Search a Collection with an AOI and Time Range#

You can find fusion imagery in [collection name]. To search within this collection, first set up your query parameters. For your AOI, you can either provide a valid geometry (e.g. point, polygon) or define a bounding box.

aoi = {
        "type": "Point",
        "coordinates": [
          -120.21102905273436,
          36.535019209789
        ]
      }
      
bbaoi = [-120.5141,36.4316,-120.4541,36.4827]

For your time of interest (TOI), you can either specify a specific datetime or a range. Date and time expressions adhere to RFC 3339.

toi = ["2021-09-14T00:00:00Z","2021-09-16T00:00:00Z"]

These can then be used in a simple search.

search = catalog.search(
    collections = "starfm_predictions_modis_landsat",
    intersects = aoi,
    datetime = toi,
    max_items = 50
)
items = search.get_all_items()
len(items)

or

search = catalog.search(
    collections = "starfm_predictions_modis_landsat",
    bbox = bbaoi,
    datetime = toi,
    max_items = 5
)
items = search.get_all_items()
len(items)

Now let’s have a look at the items that were returned in the search.

import IPython.display
itemjson = items.to_dict()
IPython.display.JSON(itemjson)

Download Data to your Local Machine#

We’ll now loop through the returned items and download them to our local machine. You must be authenticated to see data assets. For unauthenticated users, the list of assets for any item will be empty.

import requests
basefp = <your-local-path>

for item in itemjson["features"]:

    filepath = item["assets"]["lst"]["href"]
    filename = filepath[filepath.rfind("/")+1:]
    sep = "?"
    truncname = filename.split(sep, 1)[0]
    outfp = basefp + truncname
   
    with requests.get(filepath, stream=True) as result:
        result.raise_for_status()
        with open(outfp, 'wb') as f:
            for chunk in result.iter_content(chunk_size=10000000):
                f.write(chunk)

Access a Subset of Data from the Cloud Optimized Geotiff#

A key advantage of hosting the data as cloud optimized geotiffs (COG) is that you can access only the subset of a file that you need. Here we will use rasterio to open a fusion product and display the result.

import matplotlib.pyplot as plt

filepath = itemjson["features"][1]["assets"]["lst"]["href"]

# Open the file with Rasterio
import rasterio
from rasterio.windows import Window
from matplotlib import pyplot as plt
Session = rasterio.Env()

with Session:
    with rasterio.open(filepath) as src:
        meta = src.meta
        tags = src.tags()
        print(meta)
        w = meta.get("width")
        h = meta.get("height")
        win = Window(w/6, h/6, 2000, 2000)
        data = src.read(1, window=win)
        data[data == meta['nodata']] = 0
        fig, ax = plt.subplots(1, 1, figsize=(12,12))
        plt.imshow(data, cmap='magma')
        plt.colorbar(shrink=0.5)

Creating a Time Series Animation#

%pip install pystac-client
%pip install rasterio
%pip install opencv-python

import json
import pystac
from pystac_client import Client

import base64
userpass = '<username>:<password>'
b64 = base64.b64encode(userpass.encode()).decode()
headers = {'Authorization':'Basic ' + b64}

catalog = Client.open("https://fusion-stac.hydrosat.com/", headers)

geom = {'type': 'Polygon', 'coordinates': [[[-119.71183776855467, 36.34776181462616], [-119.68643188476562, 36.34776181462616], [-119.68643188476562, 36.363245618141725], [-119.71183776855467, 36.363245618141725], [-119.71183776855467, 36.34776181462616]]]}

search = catalog.search(
    collections = "starfm_predictions_modis_landsat",
    intersects = geom,
    datetime = ["2021-08-17T00:00:00Z","2021-09-10T00:00:00Z"],
    max_items = 20
)

items = search.get_all_items()
len(items)

import IPython.display
#Sort the results in time
itemjson = items.to_dict()
features = itemjson['features']
IPython.display.JSON(features)
assetlist = {}
for f in features:
    assetlist.update( {f["properties"]["datetime"] : f["assets"]["lst"]["href"]})  

%matplotlib notebook
import numpy as np
import rasterio
from rasterio.windows import Window
import matplotlib.animation as animation
import matplotlib.pyplot as plt

fig = plt.figure()

# Open the file with Rasterio
Session = rasterio.Env(
            GDAL_DISABLE_READDIR_ON_OPEN='YES',
            CPL_VSIL_CURL_USE_HEAD='NO',
            )
ims=[]
for key in sorted(assetlist.keys()):   
    with Session:
        filepath = assetlist[key]
        timestamp = key
        with rasterio.open(filepath) as src:
            meta = src.meta
            w = meta.get("width")
            h = meta.get("height")
            win = Window(0, 0, w, h)
            img1 = src.read(1, window=win)
            img1[img1 == meta['nodata']] = 0
            im = plt.imshow(img1, animated=True)
            ims.append([im])
            
ani = animation.ArtistAnimation(fig, ims, interval=50, blit=True, repeat_delay=500)

ani.save('anim.gif', writer="imagemagick")
plt.show()
Contents