Clipping rasters with d2spy¶
This guide will walk you through the steps for clipping a Cloud Optimized GeoTIFF hosted on D2S with d2spy. We will visualize the results with leafmap.
To get started, you will need to have access to a D2S instance where you have created projects and the open-source library leafmap added to your Python environment.
# Uncomment and run the following line if working out of Google Colab
# !pip install d2spy
# !pip install leafmap
import os
from datetime import date
import leafmap
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.
# Connect to D2S workspace
workspace = Workspace.connect("https://ps2.d2s.org", "yourD2Semail@example.com")
First we will need to find the project that contains the Cloud Optimized GeoTIFF (COG) we will be clipping.
# Change the search term in `.filter_by_title` to match your project
project = workspace.get_projects().filter_by_title("clip demo")[0]
print(project)
Next, find the flight that the COG was uploaded to within the project. If you have multiple flights with the same acquisition date you may need to change the index at the end of the first line of code.
# Change the date range in `filter_by_date` to match the acquistion date of the flight in your project
flight = project.get_flights().filter_by_date(date(2022, 6, 23), date(2022, 6, 23))[0]
print(flight)
Now that we have the correct flight selected, we can find the data product. If you have multiple data products with the same data type within this flight you may need to change the index at the end of the first line of code.
# Change the search term in `.filter_by_data_type` to match your COG's data type
data_product = flight.get_data_products().filter_by_data_type("dsm")[0]
print(data_product)
The polygon feature we will be using to clip the COG needs to be in a Python dictionary that matches the GeoJSON polygon feature format. In this example, we will fetch the polygon feature from a map layer previously uploaded to this project.
# This project only has one map layer so we can safely access it from the first index position
map_layer = project.get_map_layers()[0]
# The polygon feature is currently inside a GeoJSON Feature Collection
# Extract the GeoJSON Polygon Feature from the Feature Collection
clip_feature = map_layer["features"][0]
print(clip_feature)
If the COG is protected we will need to set our D2S API key in the environment. If it is public, this step can be skipped.
You can create an API key from the Profile page of your D2S instance (e.g., https://ps2.d2s.org/auth/profile). Be careful with this API key as it can be used to access any of your data products. You can revoke your current API key at any time from the Profile page.
if workspace.api_key:
os.environ["D2S_API_KEY"] = workspace.api_key
else:
# Note, you will need to log in again to fetch the key from workspace.api_key or you can manually enter the key above
print("Please create an API key from your D2S Profile page")
We can use leafmap to quickly visualize our COG and clip boundary.
os.environ["TITILER_ENDPOINT"] = "https://titiler.d2s.org"
m = leafmap.Map()
m.add_cog_layer(data_product.url, colormap_name="rainbow", name="DSM")
m.add_geojson(map_layer, style={"color": "black", "weight": 3, "fill": False})
m
out_filename = "/tmp/clipped_cog.tif"
data_product.clip(geojson_feature=clip_feature, out_raster=out_filename)
We will use leafmap once again to visualize the results.
m = leafmap.Map()
m.add_raster(out_raster, colormap="rainbow")
m.add_geojson(map_layer, style={"color": "black", "weight": 3, "fill": False})
m
Once finished viewing your data, you can revoke your authorization session by logging out.
workspace.logout()