NAIP search by D2S projects and clip¶

This guide will walk you through the steps to search a STAC Catalog using the STAC API to find NAIP items within a D2S project and clip an item by your D2S project boundary.

To get started, you'll need access to a D2S instance with existing projects. In addition to d2spy, this guide will use Python packages leafmap, numpy, and pystac_client.

In [ ]:

# Uncomment and run the following line if working out of Google Colab
# !pip install numpy
# !pip install pystac_client
# !pip install d2spy

# Uncomment and run the following line if working out of Google Colab # !pip install numpy # !pip install pystac_client # !pip install d2spy

In [ ]:

import leafmap
import numpy as np
from pystac_client import Client

from d2spy.extras.utils import clip_by_mask
from d2spy.workspace import Workspace

import leafmap import numpy as np from pystac_client import Client from d2spy.extras.utils import clip_by_mask from d2spy.workspace import Workspace

You must connect to your D2S workspace before you can request any data. The Workspace module's connect method can be used to login to a D2S instance and connect to your workspace in one go.

Note: This tutorial uses a D2S instance hosted at https://ps2.d2s.org. You will need to have an account and access to data on this instance to use it. Change the URL if you are self-hosting an instance or using an instance hosted elsewhere.

In [ ]:

# Connect to D2S workspace
workspace = Workspace.connect("https://ps2.d2s.org", "yourD2Semail@example.com")

# Connect to D2S workspace workspace = Workspace.connect("https://ps2.d2s.org", "yourD2Semail@example.com")

Locate the project you're interested in for NAIP data.

In [ ]:

# Change the search term in `.filter_by_title` to match your project
project = workspace.get_projects().filter_by_title("INDOT")[0]
project_boundary = project.get_project_boundary()

# Change the search term in `.filter_by_title` to match your project project = workspace.get_projects().filter_by_title("INDOT")[0] project_boundary = project.get_project_boundary()

In an upcoming step, you will perform a spatial query using STAC API to find NAIP items located within your project. This query will require providing STAC API with a bounding box, [xmin, ymin, xmax, ymax] for the project.

In [ ]:

# Load project boundary coordinates as numpy array
boundary_arr = np.array(project_boundary["geometry"]["coordinates"][0])
# Create a bounding box for the project
bounding_box = [
    boundary_arr[:, 0].min(),
    boundary_arr[:, 1].min(),
    boundary_arr[:, 0].max(),
    boundary_arr[:, 1].max()
]

# Load project boundary coordinates as numpy array boundary_arr = np.array(project_boundary["geometry"]["coordinates"][0]) # Create a bounding box for the project bounding_box = [ boundary_arr[:, 0].min(), boundary_arr[:, 1].min(), boundary_arr[:, 0].max(), boundary_arr[:, 1].max() ]

D2S provides a STAC API that hosts a NAIP collection for you to search. The API is accessible at https://stac-api.d2s.org. For a more user-friendly interface to browse the data exposed by the API, visit https://stac.d2s.org. In the following cells, you'll connect to the STAC API using pystac_client and search the NAIP collection using your bounding box.

In [ ]:

# Connect to STAC API
client = Client.open("https://stac-api.d2s.org")

# Connect to STAC API client = Client.open("https://stac-api.d2s.org")

In [ ]:

# Search for items from 2020 in the NAIP collection
search = client.search(
    max_items=10,
    collections=["naip"],
    bbox=bounding_box,
    datetime=['2020-01-01', '2020-12-31'],
)

# Search for items from 2020 in the NAIP collection search = client.search( max_items=10, collections=["naip"], bbox=bounding_box, datetime=['2020-01-01', '2020-12-31'], )

In [ ]:

# Print STAC Item ID and STAC Browser URL for search results
stac_browser_base_item_url = "https://stac.d2s.org/collections/naip/items"
items = []
for item in search.items():
    print(f"ID: {item.id}, URL: {stac_browser_base_item_url}/{item.id}")
    items.append(item)

# A single item is returned by this example query
item = items[0]

# Print STAC Item ID and STAC Browser URL for search results stac_browser_base_item_url = "https://stac.d2s.org/collections/naip/items" items = [] for item in search.items(): print(f"ID: {item.id}, URL: {stac_browser_base_item_url}/{item.id}") items.append(item) # A single item is returned by this example query item = items[0]

Clip NAIP Raster and Save Locally¶

First you need to find the asset within the STAC Item that represents the NAIP image/raster.

In [ ]:

# Print name of assets available in the STAC Item
for asset in item.assets:
    print(asset)

# Print name of assets available in the STAC Item for asset in item.assets: print(asset)

In [ ]:

# The URL for the NAIP raster is in the "image" asset
naip_url = item.assets["image"].href
print(naip_url)

# The URL for the NAIP raster is in the "image" asset naip_url = item.assets["image"].href print(naip_url)

The NAIP raster is a Cloud Optimized GeoTIFF, allowing us to stream the dataset directly, rather than downloading it in full before clipping it.

In [ ]:

# Desired location and name for clipped raster
out_filename = "/tmp/clipped_naip.tif"

# Clip the raster
clip_by_mask(in_raster=naip_url, geojson=project_boundary, out_raster=out_filename)

# Desired location and name for clipped raster out_filename = "/tmp/clipped_naip.tif" # Clip the raster clip_by_mask(in_raster=naip_url, geojson=project_boundary, out_raster=out_filename)

The results can be visualized with leafmap.

In [ ]:

m = leafmap.Map()
m.add_raster(out_filename)
m.add_geojson(project_boundary, style={"color": "black", "weight": 3, "fill": False})
m

m = leafmap.Map() m.add_raster(out_filename) m.add_geojson(project_boundary, style={"color": "black", "weight": 3, "fill": False}) m

Once finished viewing your data, you can revoke your authorization session by logging out.

In [ ]:

workspace.logout()

workspace.logout()