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()