import numpy as npimport geopandas as gpdimport rioxarray as rioxrimport matplotlib.pyplot as pltfrom shapely.geometry import Polygon# Used to access STAC catalogsfrom pystac_client import Client# Used to sign items from the MPC STAC catalogimport planetary_computer# other libraries for nice outputsfrom IPython.display import Image
mpc_env
We use the Client function from the pystac_client package to acess the catalog
The modifier parameter is needed to access the data in the MPC catalog.
Exploration
Let’s check out some of the catalog’s metadata:
# metadata from the catalogprint('Title:', catalog.title)print('Description:',catalog.description)
Title: Microsoft Planetary Computer STAC API
Description: Searchable spatiotemporal metadata describing Earth science datasets hosted by the Microsoft Planetary Computer
We can access the catalog collection by using the get_collections() method
catalog.get_collections() # Nothing happens# Generator objects are lazy objects you need to ask it to print it out
<generator object Client.get_collections at 0x7fd1394efbc0>
Notice the output of get_collections() is a generator.
This is a special kind of lazy object in Python over which you can loop over like a list. Unlike a list, the items in a generator do not exist in memory until you explicitely iterate over them or convert them to a list. Let’s try getting the collections from the catalog again:
# get collections and print their namescollections =list(catalog.get_collections())print('Number of Collections:', len(collections))print('Collection IDs:')for collection in collections:print('-',collection.id)
description: The [National Agriculture Imagery Program](https://www.fsa.usda.gov/programs-and-services/aerial-photography/imagery-programs/naip-imagery/) (NAIP) provides U.S.-wide, high-resolution aerial imagery, with four spectral bands (R, G, B, IR). NAIP is administered by the [Aerial Field Photography Office](https://www.fsa.usda.gov/programs-and-services/aerial-photography/) (AFPO) within the [US Department of Agriculture](https://www.usda.gov/) (USDA). Data are captured at least once every three years for each state. This dataset represents NAIP data from 2010-present, in [cloud-optimized GeoTIFF](https://www.cogeo.org/) format.
We can narrow the search withing the catalog by specifying a time range, an area of interest, and the collection name. The simplest ways to define the area of interest to look for data in the catalog are:
a GeoJSON-type dictionary with the coordinates of the bounding box,
as a list [xmin, ymin, xmax, ymax] with the coordinate values defining the four corners of the bounding box.
In this lesson we will look for the NAIP scenes over Santa Barbara from 2018 to 2023. We’ll use the GeoJSON method to define the area of interest:
# temporal range of interesttime_range ="2018-01-01/2023-01-01"#NCEAS bounding box# NCEAS bounding box (as a GeoJSON)# Dictionary and the second key has a list of tupplesbbox = {"type": "Polygon","coordinates":[ [ [-119.70608227128903, 34.426300194372274], [-119.70608227128903, 34.42041139020533], [-119.6967885126002, 34.42041139020533], [-119.6967885126002, 34.426300194372274], [-119.70608227128903, 34.426300194372274] ] ],}# catalog searchsearch = catalog.search( collections=['naip'], intersects=bbox, datetime=time_range)search
<pystac_client.item_search.ItemSearch at 0x7fd13313b650>
To get the items found in the search (or check if there were any matches in the search) we use the item_collection() method:
items = search.item_collection()len(items)
2
This output tells us there were two items in the catalog that matched our search!
items# Assets in the items you see the url for getting the data
Notice each asset has an href, which is a link to the asset object (i.e. the data). For example, we can use the URL for the rendered preview asset to plot it:
The raster data in our current item is in the image asset. Again, we access this data via its URL. This time, we open it using rioxr.open_rasterio() directly:
Notice this raster has four bands. So we cannot use the .plot.imshow() method directly (as this function only works when we have three bands). Thus we need select the bands we want to plot (RGB) before plotting:
# plot raster with correct ratiosize =6# height in in of plot heightaspect = sb.rio.width / sb.rio.height # select R,G,B bands and plotsb.sel(band=[1,2,3]).plot.imshow(size=size, aspect=aspect)
<matplotlib.image.AxesImage at 0x7fd129f1c250>
Exercise
The ‘cop-dem-glo-90’ collection contains the Copernicus DEM at 90m resolution (the data we previously used for the Grand Canyon).
Reuse the bbox for Santa Barbara to look for items in this collection. Get the first item in the search and check its assets. Check the item’s rendered preview asset by clicking on it’s URL. Open the item’s data using rioxarray
Use the bbox for Santa Barbara to look for items in this collection/
Get the first item in the search and check its assets
