Planetary Computer
Reading Data from the STAC API.
from pystac_client import Client
import json
import planetary_computer as pc
catalog = Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")
collections = catalog.get_children()
for collection in collections:
print(f"{collection.id} - {collection.title}")
landsat = catalog.get_child("landsat-8-c2-l2")
for band in landsat.extra_fields["summaries"]["eo:bands"]:
name = band["name"]
description = band["description"]
common_name = "" if "common_name" not in band else f"({band['common_name']})"
ground_sample_distance = band["gsd"]
print(f"{name} {common_name}: {description} ({ground_sample_distance}m resolution)")
area_of_interest = {
"type": "Polygon",
"coordinates": [
[
[-122.27508544921875, 47.54687159892238],
[-121.96128845214844, 47.54687159892238],
[-121.96128845214844, 47.745787772920934],
[-122.27508544921875, 47.745787772920934],
[-122.27508544921875, 47.54687159892238],
]
],
}
time_range = "2020-12-01/2020-12-31"
search = catalog.search(
collections=["landsat-8-c2-l2"], intersects=area_of_interest, datetime=time_range
)
items = list(search.get_items())
for item in items:
print(f"{item.id}: {item.datetime}")
selected_item = sorted(items, key=lambda item: item.properties["eo:cloud_cover"])[0]
for asset_key, asset in selected_item.assets.items():
print(f"{asset_key:<25} - {asset.title}")
thumbnail_asset = selected_item.assets["thumbnail"]
print(json.dumps(thumbnail_asset.to_dict(), indent=2))
import planetary_computer as pc
signed_href = pc.sign(thumbnail_asset.href)
from PIL import Image
from urllib.request import urlopen
Image.open(urlopen(signed_href))
area_of_interest = {
"type": "Polygon",
"coordinates": [
[
[-111.9839859008789, 40.5389819819361],
[-111.90502166748045, 40.5389819819361],
[-111.90502166748045, 40.57015381856105],
[-111.9839859008789, 40.57015381856105],
[-111.9839859008789, 40.5389819819361],
]
],
}
range_old = "2010-01-01/2013-01-01"
range_new = "2018-01-01/2020-01-01"
catalog = Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")
search_old = catalog.search(
collections=["naip"], intersects=area_of_interest, datetime=range_old
)
search_new = catalog.search(
collections=["naip"], intersects=area_of_interest, datetime=range_new
)
items_old = list(search_old.get_items())
items_new = list(search_new.get_items())
print(f"{len(items_old)} Items found in the 'old' range")
print(f"{len(items_new)} Items found in the 'new' range")
from shapely.geometry import shape
area_shape = shape(area_of_interest)
target_area = area_shape.area
def area_of_overlap(item):
overlap_area = shape(item.geometry).intersection(shape(area_of_interest)).area
return overlap_area / target_area
item_old = sorted(items_old, key=area_of_overlap, reverse=True)[0]
item_new = sorted(items_new, key=area_of_overlap, reverse=True)[0]
import rasterio
from rasterio import windows
from rasterio import features
from rasterio import warp
import numpy as np
from PIL import Image
def create_image(item):
print(item.datetime)
href = pc.sign(item.assets["image"].href)
with rasterio.open(href) as ds:
aoi_bounds = features.bounds(area_of_interest)
warped_aoi_bounds = warp.transform_bounds("epsg:4326", ds.crs, *aoi_bounds)
aoi_window = windows.from_bounds(transform=ds.transform, *warped_aoi_bounds)
band_data = ds.read(indexes=[1, 2, 3], window=aoi_window)
img = Image.fromarray(np.transpose(band_data, axes=[1, 2, 0]))
w = img.size[0]
h = img.size[1]
aspect = w / h
# Downscale a bit for plotting
target_w = 800
target_h = (int)(target_w / aspect)
return img.resize((target_w, target_h), Image.BILINEAR)
create_image(item_old)
create_image(item_new)
