Getting Started

Objective: Find a place, search for imagery, and rasterize it.

This tutorial demonstrates a few of the basic APIs in the Descartes Labs Platform. We’ll start by figuring out where we want to look. From there, we’ll briefly cover how to search for products using the Metadata API and what a GeoContext is. Next, we will demonstrate how to search for products and filter imagery using the Scenes API. Finally, we’ll rasterize the available imagery into a numpy array and display it using their built in methods.

In [1]: import descarteslabs as dl

In [2]: import matplotlib.pyplot as plt

In [3]: from descarteslabs.vectors import FeatureCollection, properties as p

In [4]: metadata_client = dl.Metadata()

So you now have access to a giant archive of imagery. First question is – where do you want to look? In general you’ll have to define a geometry of the area that you’re interested in. There are some ways to help you define that area. One is the website Another one is the vectors service which contains various products with geometries, such as states, cities, counties, and more. Some examples:

  • States of the USA
    • name: us_states_area
    • id: 2c3a5d20b8204023ab07218ae37ee31
  • Cities of the USA
    • name: us_cities_area
    • id: d1349cc2d8854d998aa6da92dc2bd24
  • Counties of the USA
    • name: us_county_area
    • id: 1f360e60a2ae4e3687e7d04eea2faaf

As an example, let’s try to find Taos, New Mexico (a favorite place for Cartesians to go hiking, biking, camping, and skiing).

# Find potential matches
In [5]: fc = FeatureCollection('1f360e60a2ae4e3687e7d04eea2faaf')

In [6]: matches = [f for f in fc.filter(properties=(p.county_name == 'Taos')).features()]

In [7]: matches
   'geometry': {
     'coordinates': (
         (-106.058364, 36.296978),
         (-106.052349, 36.294994),
         (-106.00071, 36.27795),
         (-105.999389, 36.277728),
         (-105.992675, 36.275444),
         (-105.898493, 36.244354),
     'type': 'Polygon'
   'id': u'35055',
   'properties': {
     u'aff_geo_code': u'0500000US35055',
     u'area_land_meters': 5704119046,
     u'area_water_meters': 3383126,
     u'county_fips_code': u'055',
     u'county_gnis_code': u'00933056',
     u'county_name': u'Taos',
     u'geo_id': u'35055',
     u'legal_area_code': u'06',
     u'state_fips_code': u'35'
   'type': 'Feature'
# Let's make Taos county our area of interest
In [8]: aoi = matches[0]

# And convert it into a GeoJSON Feature object (which contains 'geometry')
In [9]: shape = aoi.geojson

Searching for available imagery

The Descartes Labs platform pulls together data from over 75 different sources. We call each of these distinct offerings products, and they all have their own unique IDs. Even derived or external data you save into the platform has a product ID. The first parameter when searching for data is a list of product IDs to specify what kind of imagery you want returned. To view products you have access to, you can use the catalog interface which allows you to search our products and displays all their metadata information including details about spectral bands, classes, lifespan, and more. The catalog is an interface for our Metadata API. You can get back all the same information programmatically using the Metadata API.

# Return all the products you have access to
In [10]: products = metadata_client.products()

In [11]: print("I currently have access to {} products. Some of these may not be publically available.".format(len(products)))
I currently have access to 200 products. Some of these may not be publically available.
# Get back all band names based on a product ID
In [12]: bands = metadata_client.bands('landsat:LC08:PRE:TOAR')

In [13]: [band['name'] for band in bands]

Lets find some Landsat 8 imagery over our AOI using the Scenes API

There are two ways to search for and obtain data: the Scenes API, or a combination of the Metadata API and the Raster API. The Scenes API handles the details of scaling, clipping, resolution, coordinate systems, and no data values through a paradigm we call a GeoContext, while the latter requires more knowledge of the imagery you are requesting. We will focus on the Scenes API for the remainder of this tutorial.

# Create a SceneCollection matching our search criteria and view its summary
In [14]: scenes, ctx =['geometry'],
   ....:                                products="landsat:LC08:01:RT:TOAR",
   ....:                                start_datetime="2017-11-01",
   ....:                                end_datetime="2018-07-01",
   ....:                                limit=500)

In [15]: scenes
SceneCollection of 60 scenes
  * Dates: Nov 08, 2017 to Jun 27, 2018
  * Products: landsat:LC08:01:RT:TOAR: 60

Inspect the GeoContext of the SceneCollection

When a SceneCollection is created, an associated GeoContext is instantiated which intelligently sets defaults for rasterizing the imagery. This includes: cutline or boundary, spatial resolution, and coordinate reference system.

A GeoContext is immutable, but you can overwrite the default values by creating a new context as demonstrated below.

 # View the context of our SceneCollection
In [16]: ctx
    bounds=(-106.058364, 36.013014, -105.200117, 36.996017),
# Create a copy of the original context, and set the resolution to 60m.
In [17]: lowres_context = ctx.assign(resolution=60)

The Scenes API uses Shapely under the hood. This means you have accesss to all the Shapely functionality such as plotting the geometry.

# Plot the shape of Taos inline using the geometry method.
In [18]: ctx.geometry
Out[18]: <shapely.geometry.polygon.Polygon at 0x7f1597f2df90>

Display a single image

The original Scene Collection has 41 Scenes, though no single Scene has data covering our entire AOI. Before we mosaic the data, let’s take a look at a single scene, accessing its data in the form of an ndarray.

# Pick a Scene that has more than a sliver of data within our AOI.
In [19]: scene = scenes[16]
# A Scene or SceneCollection only contains the Metadata for the imagery
# We have to use another method to get the data back, in this case, the ndarray method
# Note that we are passing in the low resolution GeoContext
arr = scene.ndarray("red green blue", lowres_context)
dl.scenes.display(arr, title="Partial Coverage")

Create a mosaic of the SceneCollection

Now that we understand what a single Scene is, let’s take a look at the SceneCollection. We know a single image provides insufficient data for our AOI. The image above is missing data in the upper right hand corner of the shape. The mosaic method will give us complete coverage by overlaying scenes into a single raster. Note that where the images overlap, the pixels from the last image in the list simply “cover up” the others. The method returns a NumPy masked array, where pixels outside our AOI or containing invalid data are masked out.

# Request a mosaic of data containing the red green and blue bands
mosaic = scenes.mosaic("red green blue", lowres_context)
dl.scenes.display(mosaic, title="Taos Mosaic")

Create a cloud free composite

The mosaic function simply takes the last valid data value for any given pixel location, sometimes resulting in cloudy or undesirable results as seen above. One way to avoid issues with data in individual scenes is to create a composite: a single image of data aggregated from many scenes. To remove clouds, a simple yet often-effective method is to take the median value of each pixel across many images.

We can do that using SceneCollection.stack and NumPy. The stack method returns a 4D masked ndarray, with each Scene in the SceneCollection as a separate “layer”. Then the compositing is a simple call to

# Request a stack of data containing the red green and blue bands.
In [20]: stack = scenes.stack("red green blue", lowres_context)

In [21]: print("The stack method returns a {} containing the data for the same {} Scenes returned by the search method above.".format(type(stack),len(stack)))
The stack method returns a <class ''> containing the data for the same 60 Scenes returned by the search method above.
# We import the NumPy library so we can call the median method on our masked array to produce the final dataset
import numpy as np
composite =, axis=0)
dl.scenes.display(composite, title="Taos Composite")

Of course, you can work with any number of bands available for analysis and visualization. To complete the tutorial, let’s use the same process to display a false color image. Areas that appear red indicate healthy vegetation.

false_color_stack = scenes.stack("nir swir1 swir2", lowres_context)

composite =, axis=0)
dl.scenes.display(composite, title="Taos Infrared Composite")